matlab大地测量高斯投影正反算程序设计实验

合集下载

概述基于matlab的坐标正反算.doc

测量程序设计实验报告实验名称：坐标正反算实验三坐标正反算一、实验目的编写坐标正反算程序，并对格式化文件数据进行计算，验证程序。

二、实验内容1、编写坐标正算程序1）建立以xy_direct命名的函数，函数输入输出格式为[x2,y2] = xy_direct(x1,y1,distance, azimuth)度转度分秒：>> function dms= degree2dms(jiaodu)>>degree = fix(jiaodu);>>mimute = fix((jiaodu-degree)*60);>>second = ((jiaodu-degree)*60-mimute)*60;>>dms = degree+mimute/100+second/10000;度分秒转度：>> function degree = dms2degree(jiaodu)>>degree = fix(jiaodu);>> mimute = fix((jiaodu-degree)*100);>>second = (jiaodu-degree-mimute/100)*10000;>>degree = degree+mimute/60+second/3600;弧度转度：>> function dms=rad2dms(rad)>> rad=abs(rad);>> jiaodu=rad*180.0/pi;>> % l=fix(a)>> % b=(a-l)*60.0>> % m=fix(b)>> % a=l+m/100.0+(b-m)*0.006>> % if(rad<0)>> % dms=-a;>> % else>> % dms=a;>> % end>> degree = fix(jiaodu);>> mimute = fix((jiaodu-degree)*60);>> second = ((jiaodu-degree)*60-mimute)*60;>> dms = degree+mimute/100+second/10000;>> if(rad<0)dms=-dms;elsedms=dms;endreturn>> function [x2,y2] = xy_direct(x1,y1,distance, azimuth)>>x2=x1+distance.*cos(azimuth*pi/180);>>y2=y1+distance.*sin(azimuth*pi/180);>>end2) 对文件data1.txt中数据进行坐标正算，并将已知点和计算点坐标按照格式存贮在文件data2.txt中，data1.txt格式为： x1 y1 距离方位角（dd.mmss）data2.txt格式为：x1 y1 x2 y2>> [filename,pathname]=uigetfile;>> file=[pathname,filename];>> data=importdata(file);>> %[x1,y1]=data.data(:,[1,2]);>> azimuth=dms2degree(data.data(:,4));>> distance=data.data(:,3);>> %[x2,y2]=xy_direct(x1,y1,distance,azimuth);>>[x2,y2]=xy_direct(data.data(:,1),data.data(:,2),distance,azimuth); >> [filename_out,pathname_out]=uiputfile;>> fileout=[pathname_out,filename_out];>> fid=fopen(fileout,'wt');>> fprintf(fid,'x1 y1 x2 y2\n');>> fprintf(fid,'%8.2f %8.2f %8.2f %8.2f\n',[data.data(:,1:2),x2,y2]); >> fclose('all')ans =2、编写坐标反算程序1）建立以xy_inv命名的函数，函数输入输出格式为[distance, azimuth] = xy_inv(x1,y1, x2,y2)>> function [distance, azimuth] = xy_inv(x1,y1, x2,y2)>> delt_x=x2-x1;>> delt_y=y2-y1;>> [m,x]=size(delt_x);>> azimuth=zeros(0,m);>> for i=1:mazimuth_temp=atan2(abs(delt_y(i)),abs(delt_x(i)));if delt_x(i)>0&&delt_y(i)>0azimuth(i)=azimuth_temp;elseif delt_x(i)>0&&delt_y(i)<0azimuth(i)=2*pi-azimuth_temp;elseif delt_x(i)<0&&delt_y(i)>0azimuth(i)=pi-azimuth_temp;else delt_x(i)<0&&delt_y(i)<0;azimuth(i)=pi+azimuth_temp;endend>> azimuth=rad2dms(azimuth)>> distance=sqrt((x2-x1).^2+(y2-y1).^2);>> %fprintf('两点间距离：%8.3f ；方位角为：%8.3f',distance,azimuth);2) 对文件data2.txt中数据进行坐标反算，并将计算结果按照格式存贮在文件data3.txt中，Data3.txt格式为： x1 y1 x2 y2 距离方位角（dd.mmss）>> [filename,pathname] = uigetfile;>>file = [pathname, filename];>>data=importdata(file);>> [distance, azimuth] =xy_inv(data.data(:,1),data.data(:,2),data.data(:,3),data.data(:,4)); >> [filename_out,pathname_out] = uiputfile;>>fileout = [pathname_out, filename_out];>>fid = fopen(fileout,'wt');>>fprintf(fid,' x1 y1 x2 y2 距离方位角（dd.mmss）\n'); >>fprintf(fid,'%8.2f %8.2f %8.2f %8.2f %8.2f %8.4f\n',[dat a.data(:,1:4),distance,azimuth']');>>fclose('all');3、可能用到的函数开根号，sqrt(x)sin(rad)、cos(rad)、atan2(y,x)，find。

高斯投影反算实习报告

一、实习背景高斯投影是一种广泛应用的地图投影方法，它将地球表面的经纬度坐标转换为平面直角坐标。

高斯投影在测绘、地理信息系统、地图编制等领域有着重要的应用。

为了更好地掌握高斯投影的相关知识，提高自己的实践能力，我们进行了高斯投影反算的实习。

二、实习目的1. 理解高斯投影的基本原理和方法；2. 掌握高斯投影反算的计算步骤；3. 提高自己的实践操作能力；4. 培养团队协作精神。

三、实习内容1. 高斯投影原理高斯投影是一种等角投影，其基本原理是将地球椭球面上的经纬度坐标转换为平面直角坐标。

高斯投影具有以下特点：（1）等角投影：保持地球椭球面上任意两点间的夹角不变；（2）等积投影：保持地球椭球面上任意两块区域的面积比不变；（3）高斯-克吕格投影：以中央子午线和赤道为基准线，将地球椭球面投影到平面上。

2. 高斯投影反算步骤高斯投影反算是指将平面直角坐标转换为地球椭球面上的经纬度坐标。

其计算步骤如下：（1）计算投影面大地坐标（φ，λ）：根据给定的平面直角坐标（X，Y），利用高斯投影公式计算投影面大地坐标（φ，λ）；（2）计算大地坐标（φ，λ）：根据投影面大地坐标（φ，λ）和投影带参数，计算大地坐标（φ，λ）；（3）计算经纬度坐标（B，L）：根据大地坐标（φ，λ）和椭球参数，计算经纬度坐标（B，L）。

3. 实习过程在实习过程中，我们首先学习了高斯投影的基本原理和方法，了解了高斯投影在地图编制、地理信息系统等领域的应用。

然后，我们通过查阅相关资料，掌握了高斯投影反算的计算步骤。

在实践操作环节，我们使用高斯投影软件，对给定的平面直角坐标进行反算，得到对应的经纬度坐标。

在操作过程中，我们遇到了一些问题，如坐标转换误差、投影带参数设置等。

通过查阅资料、请教老师，我们解决了这些问题，最终完成了实习任务。

四、实习总结通过本次高斯投影反算实习，我们取得了以下成果：1. 掌握了高斯投影的基本原理和方法；2. 熟悉了高斯投影反算的计算步骤；3. 提高了实践操作能力；4. 培养了团队协作精神。

高斯投影坐标正反算编程报告

高斯投影坐标正反算编程报告1.编程思想高斯投影坐标正反算的编程需要涉及大量公式。

为了使程序更清晰，各模块的数据重用性更强，本文采用结构化编程思想。

程序由四大块组成。

大地测量家庭作业。

Cpp文件用于存储main（）函数，是整个程序的入口。

尝试通过结构化编程简化main（）函数。

myfunction.h和myfunction.cpp用于存放计算过程中进行角度弧度换算时所要用到的一些自定的转换函数。

正蒜。

H和zhengsuan CPP用于存储zhengsuan类。

在正算类中，声明了用于高斯投影坐标正演计算的所有变量。

成员变量的初始化和正向计算在类的构造函数中执行。

通过get函数获得相应的阳性结果。

fansuan.h和fansuan.cpp用于存放fansuan类，类似于zhengsuan类，fansuan类中声明了高斯投影坐标反算所要用到的所有变量，在类的构造函数中进行成员变量的初始化及反算计算。

通过get函数获得相应的反算结果。

2.计算模型高斯投影正算公式十、十、nn232244？？？？sinbcosb？Lsimbcosb（5？t？9？4？）l2？？？224??? 4n5246？？sinbcosb（61？58t？t）l720？？？6y？？nn3223？？cosb？Lcosb（1？t？？）L6.3n5242225？？cosb（5？18吨？14？58吨）l120？？？五tf2mfnf5f高斯投影反算公式B男朋友？？tfy？2tf24mfn3f？5.3t2f224？？2f？9? ftfy？720mfn46y61？90t2？45泰夫？？yy32l？？1.2t2f？？f3nfcosbf6nfcosbf？？？y524222?5?28t?24t?6??8?fffftf5120nfcosbf?3.程序框图一开始输入b，l求定带号n，中央纬度l0，纬度差l正算按照实用公式计算x,y换算为国家统一坐标x，y输出x,y输入国家统一坐标x,y由y取定带号n，并换算出x，y求出中央经线l0反算按照实用公式计算b，ll=l0+l求出大地经度l输出b，l结束4.计算结果25.附录：程序代码/////主函数入口大地测量家庭作业。

高斯投影正算与反算的理论方法与实

高斯投影正算与反算的理论方法与实现代码高斯投影是正形投影的一种，同一坐标系中的高斯投影换带计算公式是根据正形投影原理推导出的两个高斯坐标系间的显函数式。

在同一大地坐标系中(例如1954北京坐标系或1980西安坐标系)，如果两个高斯坐标系只是主子午线的经度不同，那么显函数式前的系数可以根据坐标系使用的椭球元素和主子午线经度唯一确定。

但如果两个高斯坐标系除了主子午线的经度不同以外，还存在其他线性系，则将线性变换公式代入换带计算的显函数式中，仍然可以得到严密的坐标变换公式。

此时显函数式前的系数等价于使用两个坐标系主子午线的经度和线性变换参数联合求解得到的，可以唯一确定。

//6度带宽54北京坐标系//高斯投影由大地坐标(Unit:Metres)反算经纬度(Unit:DD)void GaussProjInvCal(double X, double Y, double *longitude, double *latitude){int ProjNo; int ZoneWide; ////带宽double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval;double e1,e2,f,a, ee, NN, T,C, M, D,R,u,fai, iPI;iPI = 0.0174532925199433; ////3.1415926535898/180.0;a = 6378245.0; f = 1.0/298.3; //54年北京坐标系参数////a=6378140.0; f=1/298.257; //80年西安坐标系参数ZoneWide = 6; ////6度带宽ProjNo = (int)(X/1000000L) ; //查找带号longitude0 = (ProjNo-1) * ZoneWide + ZoneWide / 2;longitude0 = longitude0 * iPI ; //中央经线X0 = ProjNo*1000000L+500000L;Y0 = 0;xval = X-X0; yval = Y-Y0; //带内大地坐标e2 = 2*f-f*f;e1 = (1.0-sqrt(1-e2))/(1.0+sqrt(1-e2));ee = e2/(1-e2);M = yval;u = M/(a*(1-e2/4-3*e2*e2/64-5*e2*e2*e2/256));fai = u+(3*e1/2-27*e1*e1*e1/32)*sin(2*u)+(21*e1*e1/16-55*e1*e1*e1*e1/32)*sin( 4*u)+(151*e1*e1*e1/96)*sin(6*u)+(1097*e1*e1*e1*e1/512)*sin(8*u);C = ee*cos(fai)*cos(fai);T = tan(fai)*tan(fai);NN = a/sqrt(1.0-e2*sin(fai)*sin(fai));R = a*(1-e2)/sqrt((1-e2*sin(fai)*sin(fai))*(1-e2*sin(fai)*sin(fai))*(1-e2*sin(fai)*sin(fai)));D = xval/NN;//计算经度(Longitude) 纬度(Latitude)longitude1 =longitude0+(D-(1+2*T+C)*D*D*D/6+(5-2*C+28*T-3*C*C+8*ee+24*T*T)*D*D*D*D*D/120)/cos(fai);latitude1 = fai -(NN*tan(fai)/R)*(D*D/2-(5+3*T+10*C-4*C*C-9*ee)*D*D*D*D/24+(61+90*T+298*C+45*T*T-256*ee-3*C*C)*D*D*D*D*D*D/720);//转换为度DD*longitude = longitude1 / iPI;*latitude = latitude1 / iPI;}//高斯投影由经纬度(Unit:DD)正算平面坐标(含带号，Unit:Metres)void GaussProjCal(double longitude, double latitude, double *X, double *Y) {int ProjNo=0; int ZoneWide; ////带宽double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval; double a,f, e2,ee, NN, T,C,A, M, iPI;iPI = 0.0174532925199433; ////3.1415926535898/180.0;ZoneWide = 6; ////6度带宽a=6378245.0; f=1.0/298.3; //54年北京坐标系参数////a=6378140.0; f=1/298.257; //80年西安坐标系参数ProjNo = (int)(longitude / ZoneWide) ;longitude0 = ProjNo * ZoneWide + ZoneWide / 2;longitude0 = longitude0 * iPI ;latitude0=0;longitude1 = longitude * iPI ; //经度转换为弧度latitude1 = latitude * iPI ; //纬度转换为弧度e2=2*f-f*f;ee=e2*(1.0-e2);NN=a/sqrt(1.0-e2*sin(latitude1)*sin(latitude1));T=tan(latitude1)*tan(latitude1);C=ee*cos(latitude1)*cos(latitude1);A=(longitude1-longitude0)*cos(latitude1);M=a*((1-e2/4-3*e2*e2/64-5*e2*e2*e2/256)*latitude1-(3*e2/8+3*e2*e2/32+45*e2 *e2*e2/1024)*sin(2*latitude1)+(15*e2*e2/256+45*e2*e2*e2/1024)*sin(4*latitude1)-(3 5*e2*e2*e2/3072)*sin(6*latitude1));xval = NN*(A+(1-T+C)*A*A*A/6+(5-18*T+T*T+72*C-58*ee)*A*A*A*A*A/120); yval = M+NN*tan(latitude1)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24+(61-58*T+T*T+600*C-330*ee)*A*A*A*A*A*A/720);X0 = 1000000L*(ProjNo+1)+500000L;Y0 = 0;xval = xval+X0; yval = yval+Y0;*X = xval;*Y = yval;}NN卯酉圈曲率半径，测量学里面用N表示M为子午线弧长，测量学里用大X表示fai为底点纬度，由子午弧长反算公式得到，测量学里用Bf表示R为底点所对的曲率半径，测量学里用Nf表示。

高斯投影反算matlab

高斯投影反算matlab高斯投影反算是地理信息系统中的重要计算方法，通过已知的高斯投影坐标，计算出对应的地理经纬度坐标。

在Matlab中，我们可以使用一些内置函数来实现高斯投影反算，方便地进行空间坐标转换和分析。

高斯投影反算的本质是将平面坐标系转换为大地坐标系，从而实现平面上的测量和分析。

在地理信息系统中，常用的高斯投影反算方法有高斯-克吕格投影反算和高斯-克吕格-赤道投影反算。

下面我们将分别介绍这两种方法在Matlab中的实现。

我们来介绍高斯-克吕格投影反算。

在Matlab中，可以使用geotiffinfo函数读取高斯投影的tiff格式地图文件，并获取地图的投影信息。

然后，使用mapshow函数显示地图，并使用ginput 函数获取鼠标点击点的像素坐标。

接下来，可以使用imtransform 函数将像素坐标转换为投影坐标。

最后，使用projinv函数将投影坐标转换为地理经纬度坐标。

通过这些步骤，我们可以实现高斯-克吕格投影反算。

接下来，我们介绍高斯-克吕格-赤道投影反算。

在Matlab中，可以使用utmgeoid函数将地理经纬度坐标转换为UTM投影坐标。

然后，使用utmzone函数获取UTM投影的区域。

接着，使用utm2ell函数将UTM投影坐标转换为椭球坐标。

最后，使用ell2xyz函数将椭球坐标转换为空间直角坐标。

通过这些步骤，我们可以实现高斯-克吕格-赤道投影反算。

在实际应用中，高斯投影反算常常用于地图测绘、导航定位和空间分析等领域。

例如，在地图测绘中，可以利用高斯投影反算将地理经纬度坐标转换为平面坐标，从而实现地图的绘制和更新。

在导航定位中，可以利用高斯投影反算将平面坐标转换为地理经纬度坐标，从而实现定位和导航功能。

在空间分析中，可以利用高斯投影反算将不同投影坐标之间进行转换，从而实现空间数据的整合和分析。

高斯投影反算是地理信息系统中常用的计算方法之一，通过将平面坐标系转换为大地坐标系，实现了地理经纬度坐标和投影坐标之间的转换。

matlab大地测量空间直角坐标转换程序设计

实验日期： 6.1 空间直角坐标转换
一、实验目的不同坐标系下空间直角坐标转换是测量工程中经常会遇到的问题，本实验通过编写空间直角坐标七参数转换的程序，并对格式化文件数据进行计算，验证程序，从而掌握空间直角转换的基本原理和方法。二、实验内容： 1、空间直角坐标转换七参数模型（布尔莎模型）布尔莎 -沃尔夫模型(在我国通常被简称为布尔莎模型 )又被称为七参数转换 (7-Parameter Transformation)或七参数赫尔墨特变换(7-Parameter Helmert Transformation) ，其数学公式为：
0 1 R1 ( X ) 0 cos X 0 sin X
cos Y 0 sin X R2 ( Y ) 0 cos X sin Y
0 sin Y 1 0 0 cos Y
cos Z R3 ( Z ) sin Z 0
七参数计算函数 TtransParam7_Comp 如下：
%由空间直角坐标求解七参数 x0,y0,z0,a1,a2,a3,a4 function Param7 = TtransParam7_Comp (x1,y1,z1,x2,y2,z2)
A=[x1,y1,z1,x2,y2,z2]; [m,n]=size(A); B=zeros(3*m,7);%创建一个 3m 行， 7 列的零矩阵 B L=zeros(3*m,1);%创建一个 3m 行， 1 列的零矩阵 L for i=1:m
-2-
B(3*i-2:3*i,:)=[1 0 0 x1(i) 0 -z1(i) y1(i); 0 1 0 y1(i) z1(i) 0 -x1(i); 0 0 1 z1(i) -y1(i) x1(i) 0] %for 循环组成 B 矩阵 L(3*i-2:3*i,:)=[x2(i),y2(i),z2(i)] %for 循环组成 L 矩阵 end gX=(inv(B'*B))*B'*L %计算误差方程的系数 %由 WGS-84 空间直角坐标和 BJ-54 空间直角坐标列方程组， %根据最小二乘法 VTPV=min 的原则，可列出法方程为求得布尔莎七参数 %根据求得的七参数通过反算，可将 WGS-84 空间直角坐标转化为 BJ-54 空间直角坐标 m=gX(4)-1; ex=(gX(5)/gX(4))*180/pi*3600; ey=(gX(6)/gX(4))*180/pi*3600; ez=(gX(7)/gX(4))*180/pi*3600; Param7=[gX(1),gX(2),gX(3),ex,ey,ez,m];%组成布尔莎七参数矩阵 End

测绘程序设计(VS2008)实验报告--高斯投影正反算

《测绘程序设计（）》上机实验报告（Visual C++.Net）班级：学号：姓名：序号：二零一一年五月实验7 常用测量程序设计1.实验目的：1.1 巩固类的创建与使用;1.2 掌握数组参数的传递；1.3 掌握常用测绘程序设计的技巧。

1.42.实验内容：编写高斯投影正反算程序。

3.设计思路：这次的实验目的是实现高斯正反算。

需要考虑投影方式即分带的方式，又要考虑椭球参数的类型，所以我添加了两个函数来完成此功能。

分别是int SetProjectType(int m)和void SetParameter(int m,double &a,double &b)。

4.界面设计：界面设计很简单，具体见运行结果。

5.主要代码：文件名：GaussProjectDlg.cpp代码：const double PI=4*atan(1.0);//获得分带方式返回中央子午线经度int CGaussProjectDlg::SetProjectType(int m){UpdateData(TRUE);int n; //记录分带带号double L; //经度L=iDegreeL+iMinL/60+dSecondL/3600;if(m==1) //6度带{n=int(L/6)+1;L0=6*n-3;}else if(m==2) //3度带{n=int((L+1.5)/3);L0=3*n;}else if(m==3) //自主分带L0=L0;return L0;}//获取椭球参数void CGaussProjectDlg::SetParameter(int m,double &a,double &b) {if(m==1) //克拉索夫斯基椭球{a=6378245.0;b=6356863.0187730473;//e=sqrt(0.006693421622966);}else if(m==2) //1975国际协议椭球{a=6378140.0;b=6356755.2881575287;//e=sqrt(0.006694384999588);}else if(m==3) //WGS-84椭球{a=6378137.0;b=6356752.3142;//e=sqrt(0.0066943799013);}}void CGaussProjectDlg::OnBnClickedButtonpositivecal(){// TODO: 在此添加控件通知处理程序代码UpdateData(TRUE);double N;double t;double Eta;double X;double A0,A2,A4,A6,A8;double RadB;double Rou;Rou=180*3600/PI;double a,b,e1,e2; //椭球参数SetParameter(iParameterType,a,b);e1=sqrt(a*a-b*b)/a;e2=sqrt(a*a-b*b)/b;double l;L0=SetProjectType(iProjectType);double L;L=iDegreeL+double(iMinL)/60+dSecondL/3600;l=(L-L0)*3600;RadB=(iDegreeB+double(iMinB)/60+dSecondB/3600)*PI/180;N=a/sqrt(1-e1*e1*sin(RadB)*sin(RadB));t=tan(RadB);Eta=e2*cos(RadB);A0=1+3.0/4*e1*e1+45.0/64*pow(e1,4)+350.0/512*pow(e1,6)+11025.0/16384*pow(e1,8);A2=-1.0/2*(3.0/4*e1*e1+60.0/64*pow(e1,4)+525.0/512*pow(e1,6)+17640.0/16384*pow(e1,8));A4=1.0/4*(15.0/64*pow(e1,4)+210.0/512*pow(e1,6)+8820.0/16384*pow(e1,8));A6=-1.0/6*(35.0/512*pow(e1,6)+2520.0/16384*pow(e1,8));A8=1.0/8*(315.0/16384*pow(e1,8));X=a*(1-e1*e1)*(A0*RadB+A2*sin(2*RadB)+A4*sin(4*RadB)+A6*sin(6*RadB)+A8*sin(8*RadB));x=X+N/(2*Rou*Rou)*sin(RadB)*cos(RadB)*l*l+N/(24*pow(Rou,4))*sin(RadB)*pow(cos(RadB),3)*(5-t*t+9*Eta*Eta+4*pow(Eta,4))*pow(l,4)+ N/(720*pow(Rou,6))*sin(RadB)*pow(cos(RadB),5)*(61-58*t*t+pow(t,4))*pow(l,6);y=N/Rou*cos(RadB)*l+N/(6*pow(Rou,3))*pow(cos(RadB),3)*(1-t*t+Eta*Eta)*pow(l,3)+N/(120*pow(Rou,5))*pow(cos(RadB),5)*(5-18*t*t+pow(t,4)+14*Eta*Eta-58*Eta*Eta*t*t)*pow( l,5);UpdateData(FALSE);}void CGaussProjectDlg::OnBnClickedButtonantical(){// TODO: 在此添加控件通知处理程序代码UpdateData(TRUE);double t_f;double Eta_f;double B_f;double N_f;double M_f;double X=x;double B0;double K0,K2,K4,K6;double a,b,e1,e2; //椭球参数SetParameter(iParameterType,a,b);e1=sqrt(a*a-b*b)/a;e2=sqrt(a*a-b*b)/b;double A0;A0=1+3.0/4*e1*e1+45.0/64*pow(e1,4)+350.0/512*pow(e1,6)+11025.0/16384*pow(e1,8);B0=X/(a*(1-e1*e1)*A0);K0=1.0/2*(3.0/4*e1*e1+45.0/64*pow(e1,4)+350.0/512*pow(e1,6)+11025.0/16384*pow(e1,8));K2=-1.0/3*(63.0/64*pow(e1,4)+1108.0/512*pow(e1,6)+58239.0/16384*pow(e1,8));K4=1.0/3*(604.0/512*pow(e1,6)+68484.0/16384*pow(e1,8));K6=-1.0/3*(26328.0/16384*pow(e1,8));B_f=B0+sin(2*B0)*(K0+sin(B0)*sin(B0)*(K4+K6*sin(B0)*sin(B0)));t_f=tan(B_f);Eta_f=e2*cos(B_f);N_f=a/sqrt(1-e1*e1*sin(B_f)*sin(B_f));M_f=N_f/(1+e2*e2*cos(B_f)*cos(B_f));double B;B=B_f-t_f/(2*M_f*N_f)*y*y+t_f/(24*M_f*pow(N_f,3))*(5+3*t_f*t_f+Eta_f*Eta_f-9*Eta_f*Eta_f*t_f*t_f)*pow(y,4)- t_f/(720*M_f*pow(N_f,5))*(61+90*t_f*t_f+45*pow(t_f,4))*pow(y,6);double l;l=1.0/(N_f*cos(B_f))*y-1.0/(6*pow(N_f,3)*cos(B_f))*(1+2*t_f*t_f+Eta_f*Eta_f)*pow(y,3)+1.0/(120*pow(N_f,5)*cos(B_f))*(5+28*t_f*t_f+24*pow(t_f,4)+6*Eta_f*Eta_f+8*Eta_f*Eta_f* t_f*t_f)*pow(y,5);//将B转化为度分秒的形式double dDegB;dDegB=B*180/PI;iDegreeB=int(dDegB);iMinB=int((dDegB-iDegreeB)*60);dSecondB=((dDegB-iDegreeB)*60-iMinB)*60;double dDegL;dDegL=l*180/PI+L0;iDegreeL=int(dDegL);iMinL=int((dDegL-iDegreeL)*60);dSecondL=((dDegL-iDegreeL)*60-iMinL)*60;UpdateData(FALSE);}6.运行结果：实验的运行结果如下图所示：正算：反算：7.实验总结这次实验是实现高斯投影的正反算，方法很多，实现并不复杂，但是计算公式复杂，变量繁多，稍有不慎，就会造成计算错误。

高斯投影正反算实习

大地测量学编程实习报告姓名：鲁尼学号：10 班级：曼联编程思想:这个投影是由德国数学家、物理学家、天文学家高斯于19 世纪20 年代拟定，后经德国大地测量学家克吕格于1912 年对投影公式加以补充，故称为高斯-克吕格投影。

即等角横切椭圆柱投影。

假想用一个圆柱横切于地球椭球体的某一经线上，这条与圆柱面相切的经线，称中央经线。

以中央经线为投影的对称轴，将东西各3°或1°30′的两条子午线所夹经差6°或3°的带状地区按数学法则、投影法则投影到圆柱面上，再展开成平面，即高斯-克吕格投影，简称高斯投影。

这个狭长的带状的经纬线网叫做高斯-克吕格投影带。

高斯投影正算公式就是由大地坐标（L，B）求解高斯平面坐标（x，y），而高斯投影反算公式则是由高斯平面坐标（x，y）求解大地坐标（L，B）。

现行的高斯投影用表都是采用克拉索夫斯基椭球参数，这次编程计算就是采用这种椭球参数，并采用实用公式按6度分带投影。

编程环境是在VC下，采用C++语言编写。

程序主要分为两部分，第一部分是高斯正反算函数，第二部分是主函数。

高斯正反算函数，参考书上175页的电算公式，正算时先将度数换算成秒，再定带号n，求中央经线l0，经度差l''，然后根据实用公式计算高斯平面坐标x，y。

最后计算国家统一坐标的x，y，再将其输出。

计算和数据模型：正算是指：由大地坐标（L,B）求得高斯平面坐标（x,y）的过程。

反算是指：由高斯平面坐标（x,y）求得大地坐标（L,B）的过程。

正算：高斯投影必须满足的三个条件：（1），中央子午线投影后为直线。

（2），中央子午线投影后长度不变。

（3），投影具有正性性质，即正性投影条件。

由第一个条件可知，中央子午线东西两侧的投影必然对称于中央子午线。

设在托球面上有P1 ，P2，且对称于中央子午线。

其大地坐标为（l,B）,(-l,B)则投影后的平面坐标一定为P1·（x,y）,P2·（x,-y）.由第二个条件可知，位于中央子午线上的点，投影后的纵坐标x应该等于投影前从赤道量至该点的子午弧长。

大地测量MATLAB实验坐标正反算程序设计

实验三坐标正反算一、实验目的编写坐标正反算程序，并对格式化文件数据进行计算，验证程序。

二、实验内容及步骤：1、编写坐标正算程序要求：1）建立以xy_direct命名的函数，函数输入输出格式为[x2,y2] = xy_direct(x1,y1,distance, azimuth)2) 对文件data1.txt中数据进行坐标正算，并将已知点和计算点坐标按照格式存贮在文件data2.txt中，data1.txt格式为：x1 y1 距离方位角（dd.mmss）data2.txt格式为：x1 y1 x2 y2编写程序共三个模块，分三个m文件组成，各m文件内容如下：1）主运行程序clear;clc;[filename, pathname]= uigetfile('*.*'); % 文件查找窗口file=fullfile(pathname, filename); % 合并路径文件名A=importdata(file);dms=A.data(:,4); %提取A.data中所有行和第4列azimuth=dms2rad(dms); %度分秒转化成弧度[x2,y2]=xy_direct(A.data(:,1),A.data(:,2),A.data(:,3),azimuth(:,1));B=[A.data(:,1),A.data(:,2),x2,y2]; %输出得到B矩阵fid=fopen('fileout','wt'); %写入文件路径fprintf(fid,' x(m) y(m) x2(m) y2(m)\n' ); fprintf(fid,'%f\t%f\t%f\t%f\n',B(:,1:4)');fclose(fid);2）xy_direct程序如下：function [x2,y2] = xy_direct(x1,y1,distance, azimuth)%求x2,y2函数x2=x1+distance.*cos(azimuth);y2=y1+distance.*sin(azimuth);end3）度分秒转弧度函数function azimuth=dms2rad(dms)%度分秒转弧度函数dms_abs=abs(dms);dgree=fix(dms_abs);min_tem=(dms_abs-dgree)*100;min=fix(min_tem);second=(min_tem-min)*100;azimuth=(dgree+min/60+second/3600)*pi/180;if dms<0azimuth=-azimuth;elseazimuth=azimuth;endreturn运行程序所得data2内容如下：2、编写坐标反算程序要求：1）建立以xy_inv命名的函数，函数输入输出格式为[distance, azimuth] = xy_inv(x1,y1, x2,y2)2) 对文件data2.txt中数据进行坐标反算，并将计算结果按照格式存贮在文件data3.txt中，Data3.txt格式为：x1 y1 x2 y2 距离方位角（dd.mmss）坐标反算程序共三个模块，分三个m文件组成，各m文件内容如下：1）主运行程序：clear;clc;[filename, pathname]= uigetfile('*.*'); % 文件查找窗口file=fullfile(pathname, filename); % 合并路径文件名A=importdata(file);for i=1:4[distance(i),azimuth(i)]=xy_inv(A.data(i,1),A.data(i,2),A.data(i,3),A.data(i,4));%提取data数据;dms(i)=rad1dms(azimuth(i)); %rad-->dmsendB=[A.data(:,1),A.data(:,2),A.data(:,3),A.data(:,4),distance',dms'];%输出得到B矩fid=fopen('data3.txt','wt'); %写入文件路径fprintf(fid,'x1 y1 x2 y2 距离方位角（dd.mmss）\n' );[m,n]=size(B); %输入矩阵m,n取行列号for i=1:mfor j=1:nif j==nfprintf(fid,'%f\n',B(i,j));elsefprintf(fid,'%f\t',B(i,j));endendendfclose(fid);2）xy_inv程序如下：function [distance, azimuth] = xy_inv(x1,y1, x2,y2)distance=sqrt((x2-x1)^2+(y2-y1)^2);if (x2-x1)~=0azimuth=atan(abs((y2-y1)/(x2-x1)));%求弧度值%判断象限if (x2-x1)>0&(y2-y1)>0azimuth=azimuth;elseif (x2-x1)>0&(y2-y1)==0azimuth=0;elseif (x2-x1)<0&(y2-y1)>0azimuth=pi-azimuth;elseif (x2-x1)<0&(y2-y1)==0azimuth=pi;elseif (x2-x1)<0&(y2-y1)<0azimuth=pi+azimuth;elseif (x2-x1)>0&(y2-y1)<0azimuth=2*pi-azimuth;endelseif(y2-y1)>0azimuth=pi/2;elseazimuth=3*pi/2;endend3）弧度转度分秒函数function dms=rad1dms(azimuth)%弧度转度分秒dgree_tem=azimuth*180/pi;dgree=fix(dgree_tem);min_tem=((dgree_tem-dgree)*60);min=fix(min_tem);second=((min_tem-min)*60);dms=dgree+min/100+second/10000;end运行程序所得data3内容如下：。

MATLAB高斯正反算

cout<<"【反算】"<<endl;
cout<<"请输入国家统一坐标 X Y。例如 3378627.1819 20243953.4517"<<endl;
cin>>myX>>myY;
FansuanmyFansuan1(myX,myY);
myFansuan1.printLocation();
printf("Radian B=%f L=%f \n",myZhengsuan1.getrB(),myZhengsuan1.getrL());
? myZhengsuan1.printLocation();?
}?
voidfansuan()
{
doublemyX,myY;
{
printf("正算得国家统一坐标为： X= %8.8f Y=%8.8f \n",X,Y);
}
doubleZhengsuan::getrB()
{
returnrB;
}
doubleZhengsuan::getrL()
{
returnrL;
}
///反算类
Fansuan.h
Fansuan.h: interface for the Fansuan class.
#if !defined(AFX_FANSUAN_H__5B5E4668_EE81_463F_8D65_FFE2AEACE745__INCLUDED_)
#define AFX_FANSUAN_H__5B5E4668_EE81_463F_8D65_FFE2AEACE745__INCLUDED_

高斯投影坐标正反算及程序设计

高斯投影坐标正反算及程序设计目录1研究背景和意义 (1)2正形投影的特性及其公式 (2)2.1正形投影特性 (2)2.2正形投影的一般条件 (2)3高斯坐标正反算公式 (6)3.1高斯投影正算 (6)3.2高斯坐标反算 (7)4使用matlab编写高斯坐标正反算函数 (10)4.1高斯坐标正算函数编写 (10)4.2高斯坐标反算函数编写 (13)5高斯坐标正反算程序的界面设计 (17)5.1单点换算模块的编写 (17)5.2批量换算模块的编写 (23)总结与讨论 (33)致谢 (34)参考文献 (34)附录A (35)1研究背景和意义高斯投影即高斯-克吕格投影这个投影是在1920年左右被拟定的，拟定者是德国人高斯，一位著名的数学家、物理学家、天文学家[1]，之后有又有一名叫克吕格的测量学家在1912 年对投影公式加以补充，所以这个投影的全称是高斯-克吕格投影，又名"等角横切椭圆柱投影”，是地球椭球面和平面间正形投影的一种。

即高斯投影是正形投影的一种，所以它具有正形投影的特点，即角度不变性、图形相似性以及在某点方向上长度比的同一性[2]。

而为了使投影区域的长度变形不致过大，采用的措施就是分带。

这样既保证了长度变形不致过大，又可以在不同的投影带里用一样的数学公式来进行各种大地问题的计算，且带与带间的互相换算也能用相同的公式和方法进行[3]。

正是因为这种特点高斯投影被广泛应用到了测绘工作中。

虽然高斯投影作用非常广泛但是因为其复杂繁琐的计算过程往往需要通过程序设计来达到简化计算，提高效率目的。

本文所使用的MATLAB 软件编程技术广泛应用于测量数据处理领域，MATLAB是matrix&laboratory两个词的组合，意为矩阵工厂（矩阵实验室）。

是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境[4]。

它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，能够开发出界面友好、使用方便的图形界面[5]，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。

高斯投影正反算实习报告

高斯投影正反算实习报告程序代码#include <stdio.h>#include <iostream.h>#include <math.h>#include<iomanip.h>class Angle{public:double degree,cent,second,Hudu,seconds;//构造函数Angle(double _degree,double _cent,double _second){degree=_degree;cent=_cent;second=_second;seconds=(degree*3600+cent*60+second);Hudu=(degree*3600+cent*60+second)*2*3.1415926535/(360*60*60);}Angle(){}SToD(){second=seconds-int(seconds/60)*60;degree=int((seconds-second)/3600);cent=int((seconds-second-degree*3600)/60);}~Angle(){}};void main(){//正算过程//定义变量Angle B(51,38,43.9023),L(126,2,13.1360),_B,_L;double L0,l,N,a0,a4,a6,a3,a5,x,y,rou;rou=(360*60*60)/(2*3.1415926535);if(int(L.degree)%6!=0){L0=6*(int(L.degree)/6+1)-3;}elseL0=6*(int(L.degree)/6)-3;l=L.Hudu-L0*3600*2*3.1415926535/(360*60*60);N=6399698.902-(21562.267-(108.973-0.612*cos(B.Hudu)*cos(B.Hudu))*cos(B.Hudu)*co s(B.Hudu))*cos(B.Hudu)*cos(B.Hudu);a0=32140.404-(135.3302-(0.7092-0.0040*cos(B.Hudu)*cos(B.Hudu))*cos(B.Hudu)*cos(B .Hudu))*cos(B.Hudu)*cos(B.Hudu);a4=(0.25+0.00252*cos(B.Hudu)*cos(B.Hudu))*cos(B.Hudu)*cos(B.Hudu)-0.04166;a6=(0.166*cos(B.Hudu)*cos(B.Hudu)-0.084)*cos(B.Hudu)*cos(B.Hudu);a3=(0.3333333+0.001123*cos(B.Hudu)*cos(B.Hudu))*cos(B.Hudu)*cos(B.Hudu)-0.1666 667;a5=0.0083-(0.1667-(0.1968+0.0040*cos(B.Hudu)*cos(B.Hudu))*cos(B.Hudu)*cos(B.Hud u))*cos(B.Hudu)*cos(B.Hudu);x=6367558.4969*B.Hudu-(a0-(0.5+(a4+a6*l*l)*l*l)*l*l*N)*sin(B.Hudu)*cos(B.Hudu);y=(1+(a3+a5*l*l)*l*l)*l*N*cos(B.Hudu);//结果输出cout<<"x="<<x<<endl;cout<<"y="<<y<<endl;// 反算过程double b,Bf,Nf,Z,b2,b3,b4,b5,l1;b=(x/6367558.4969);Bf=b+(50221746+(293622+(2350+22*cos(b)*cos(b))*cos(b)*cos(b))*cos(b)*cos(b))*0.000 0000001*sin(b)*cos(b);Nf=6399698.902-(21562.267-(108.973-0.612*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf))*cos(Bf)*c os(Bf);Z=y/(Nf*cos(Bf));b2=(0.5+0.003369*cos(Bf)*cos(Bf))*sin(Bf)*cos(Bf);b3=0.333333-(0.166667-0.001123*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b4=0.25+(0.16161+0.00562*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b5=0.2-(0.1667-0.0088*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);_B.seconds=Bf*rou-(1-(b4-0.12*Z*Z)*Z*Z)*Z*Z*b2*rou;l1=(1-(b3-b5*Z*Z)*Z*Z)*Z*rou;_L.seconds=L0*3600+l1;_B.SToD();_L.SToD();//结果输出cout<<"_B="<<_B.degree<<"du"<<_B.cent<<"fen"<<setprecision(8)<<_B.second<<" miao"<<endl;cout<<"_L="<<_L.degree<<"du"<<_L.cent<<"fen"<<setprecision(8)<<_L.second<<" miao"<<endl;计算结果x=5.72816e+006y=-205080_B=51du38fen43.902358miao_L=126du2fen13.135972miaoPress any key to continue程序设计流程图。

高斯投影坐标正反算编程报告

高斯投影坐标正反算编程报告The Standardization Office was revised on the afternoon of December 13, 2020高斯投影坐标正反算编程报告1. 编程思想进行高斯投影坐标正反算的编程需要牵涉到大量的公式，为了使程序条理更清楚，各块的数据复用性更强，这里采取了结构化的编程思想。

程序由四大块组成。

文件用于存放main()函数，是整个程序的入口。

通过结构化的编程尽力使main()函数变得简单。

和用于存放计算过程中进行角度弧度换算时所要用到的一些自定的转换函数。

和用于存放Zhengsuan 类，在Zhengsuan 类中声明了高斯投影坐标正算所要用到的所有变量，在类的构造函数中进行成员变量的初始化及正算计算。

通过get 函数获得相应的正算结果。

和用于存放Fansuan 类，类似于Zhengsuan 类，Fansuan 类中声明了高斯投影坐标反算所要用到的所有变量，在类的构造函数中进行成员变量的初始化及反算计算。

通过get 函数获得相应的反算结果。

2. 计算模型高斯投影正算公式6425644223422)5861(cos sin 720)495(cos 24cos sin 2lt t B B N lt B simB N l B B N X x ''+-''+''++-''+''⋅''+=ρηηρρ5222425532233)5814185(cos 120)1(cos 6cos lt t t B N lt B N l B N y ''-++-''+''+-''+''⋅''=ηηρηρρ高斯投影反算公式()()()()22242552233642542222328624285cos 12021cos 6cos 459061720935242f f f f f ff ff ff f f ff ff ff f f f ff f ff f f t t t B N y t B N y B N y l y t t y NM t yt tNM t y N M t B B ηηηηη+++++++-=++--+++-=3. 程序框图4. 计算结果5.附录：程序代码517"<<endl;cin>>myX>>myY;Fansuan myFansuan1(myX,myY);();}////////////////////////////////////////////////////////////////////////#include ""//////////////////////////////////////////////////////////////////////// Construction/Destruction//////////////////////////////////////////////////////////////////////Fansuan::Fansuan(){}Fansuan::Fansuan(double X,double Y){this->X=X;this->Y=Y;//初始化x,yN=(int)(Y/1000000);//取出带号L0=6*N-3;//初始化该带号的中央经线经度x=X;y=Y-1000000*N-500000;beta=x/;//初始化beta，弧度单位betaSecond=beta*rouSecond;//初始化beta，秒单位betaDegree=betaSecond/3600;//初始化beta,整度数单位Bf=beta+(+(293622+(2350+22*cos(beta)*cos(beta))*cos(beta)*cos(beta))*cos(b eta)*cos(beta))*(1e-10)*sin(beta)*cos(beta);//初始化Bf,弧度单位BfSecond=Bf*rouSecond;//初始化Bf，秒单位BfDegree=BfSecond/3600;//初始化Bf，整度数单位Nf= Z=y/(Nf*cos(Bf));b2=+*cos(Bf)*cos(Bf))*sin(Bf)*cos(Bf);b4=++*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b3= b5= Bsecond=BfSecond-(1-*Z*Z)*Z*Z)*Z*Z*b2*rouSecond;//计算大地经度B，单位为秒B=Bsecond/3600;//用整度数表示Blsecond=(1-(b3-b5*Z*Z)*Z*Z)*Z*rouSecond;//计算经度差l，单位为秒l=lsecond/3600;//用整度数表示lL=L0+l;//计算大地经度L}double Fansuan::getB(){return B;}double Fansuan::getL(){return L;}void Fansuan::printLocation(){printf("反算得大地坐标为：大地纬度B= %f°大地经度L=%f°\n",B,L);}Fansuan::~Fansuan(){}。

高斯投影正反算编程一.高斯投影正反算基本公式

高斯投影正反算编程一．高斯投影正反算基本公式（1）高斯正算基本公式（2）高斯反算基本公式以上主要通过大地测量学基础课程得到，这不进行详细的推导，只是列出基本公式指导编程的进行。

二．编程的基本方法和流程图（1）编程的基本方法高斯投影正反算基本上运用了所有的编程基本语句，本文中是利用C++语言进行基本的设计。

高斯正算中对椭球参数和带宽的选择主要运用了选择语句。

而高斯反算中除了选择语句的应用，在利用迭代算法求底点纬度还应用了循环语句。

编程中还应特别注意相关的度分秒和弧度之间的相互转换，这是极其重要的。

（2）相关流程图1）正算2）反算三．编程的相关代码（1）正算# include "stdio.h"# include "stdlib.h"# include "math.h"# include "assert.h"#define pi (4*atan(1.0))int i;struct jin{double B;double L;double L0;};struct jin g[100];main(int argc, double *argv[]){FILE *r=fopen("a.txt","r");assert(r!=NULL);FILE *w=fopen("b.txt","w");assert(r!=NULL);int i=0;while(fscanf(r,"%lf %lf %lf",&g[i].B,&g[i].L,&g[i].L0)!=EOF){double a,b;int zuobiao;printf("\n请输入坐标系：北京54=1，西安80=2，WGS84=3：");scanf("%d",&zuobiao);getchar();if(zuobiao==1){a=6378245;b=6356863.0187730473;}if(zuobiao==2){a=6378140;b=6356755.2881575287;}if(zuobiao==3){a=6378137;b=6356752.3142;} //选择坐标系//double f=(a-b)/a;double e,e2;e=sqrt(2*f-f*f);e2=sqrt((a/b)*(a/b)-1);//求椭球的第一，第二曲率//double m0,m2,m4,m6,m8;double a0,a2,a4,a6,a8;m0=a*(1-e*e);m2=3*e*e*m0/2;m4=5*e*e*m2/4;m6=7*e*e*m4/6;m8=9*e*e*m6/8;a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;a2=m2/2+m4/2+15*m6/32+7*m8/16;a4=m4/8+3*m6/16+7*m8/32;a6=m6/32+m8/16;a8=m8/128;double Bmiao,Lmiao, L0miao;Bmiao=(int)(g[i].B)*3600.0+(int)((g[i].B-(int)(g[i].B))*100.0)*60.0+( g[i].B*100-(int)(g[i].B*100))*100.0;Lmiao=(int)(g[i].L)*3600.0+(int)((g[i].L-(int)(g[i].L))*100.0)*60.0+(g [i].L*100-(int)(g[i].L*100))*100.0;L0miao=(int)(g[i].L0)*3600.0+(int)((g[i].L0-(int)(g[i].L0))*100.0)*60. 0+(g[i].L0*100-(int)(g[i].L0*100))*100.0;double db;db=pi/180.0/3600.0;double B1,L1,l;B1=Bmiao*db;L1= Lmiao*db;l=L1-L0miao*db;//角度转化为弧度//double T=tan(B1)*tan(B1);double n=e2*e2*cos(B1)*cos(B1);double A=l*cos(B1);double X,x,y;X=a0*(B1)-a2*sin(2*B1)/2+a4*sin(4*B1)/4-a6*sin(6*B1)/6+a8*sin(8*B1)/8;//求弧长//double N=a/sqrt(1-e*e*sin(B1)*sin(B1));int Zonewide;int Zonenumber;printf("\n请输入带宽：3度带或6度带Zonewide=");scanf("%d",&Zonewide);getchar();if(Zonewide==3){Zonenumber=(int)((g[i].L-Zonewide/2)/Zonewide+1);}else if(Zonewide==6){Zonenumber=(int)g[i].L/Zonewide+1;}else{printf("错误");exit(0);}//选择带宽//doubleFE=Zonenumber*1000000+500000;//改写为国家通用坐标//y=FE+N*A+A*A*A*N*(1-T*T+n*n)/6+A*A*A*A*A*N*(5-18*T*T+T *T*T*T+14*n*n-58*n*n*T*T)/120;x=X+tan(B1)*N*A*A/2+tan(B1)*N*A*A*A*A*(5-T*T+9*n*n+4*n*n *n*n)/24+tan(B1)*N*A*A*A*A*A*A*(61-58*T*T+T*T*T*T)/720;printf("\n所选坐标系的转换结果：x=%lf y=%lf\n",x,y);fprintf(w,"%lf %lf\n",x,y);//输出结果到文本文件//}fclose(r);fclose(w);system("pause");return 0;}（2）反算# include "stdio.h"# include "stdlib.h"# include "math.h"# include "assert.h"#define pi (4*atan(1.0))double X,Y,B1,B2,B3,F,t;double m0,m2,m4,m6,m8;double a0,a2,a4,a6,a8,a1,b1;double BB,LL,Bf;double e,e1;int d,m,s,i,zuobiao;double sort(double,double);struct jin{double x;double y;double L0;};struct jin g[100];//x,y,L0为输入量：x,y坐标和中央子午线经度// main(int argc, double *argv[]){FILE *r=fopen("c.txt","r");assert(r!=NULL);FILE *w=fopen("d.txt","w");assert(r!=NULL);int i=0;while(fscanf(r,"%lf %lf %lf",&g[i].x,&g[i].y,&g[i].L0)!=EOF)//文件为空，无法打开//{double a1=6378245.0000000000;//克拉索夫斯基椭球参数//double b1=6356863.0187730473;double a75=6378140.0000000000;//1975国际椭球参数//double b75=6356755.2881575287;double a84=6378137.0000000000;//WGS-84系椭球参数//double b84=6356752.3142000000;double M,N;//mouyou圈曲率半径，子午圈曲率半径//double t,n;double A,B,C;double BB,LL,Bf,LL0,BB0;double a,b;printf("\n选择参考椭球：1=克拉索夫斯基椭球，2=1975国际椭球，3=WGS-84系椭球:");scanf("%d",&zuobiao);getchar();if(zuobiao==1){a=a1;b=b1;}if(zuobiao==2){a=a75;b=b75;}if(zuobiao==3){a=a84;b=b84;}//选择参考椭球，求解第一偏心率e,第二偏心率e1// Bf=sort(a,b);//调用求解底点纬度的函数//double q=sqrt(1-e*e*sin(Bf)*sin(Bf));double G=cos(Bf);M=a*(1-e*e)/(q*q*q);N=a/q;double H,I;A=g[i].y/N;H=A*A*A;I=A*A*A*A*A;t=tan(Bf);n=e1*cos(Bf);B=t*t;C=n*n;BB0=Bf-g[i].y*t*A/(2*M)+g[i].y*t*H/(24*M)*(5+3*B+C-9*B*C)-g[i] .y*t*I/(720*M)*(61+90*B+45*B*B);LL0=g[i].L0*pi/180.0+A/G-H/(6*G)*(1.0+2*B+C)+I/(120*G)*(5.0+28 *B+24*B*B+6*C+8*B*C);//利用公式求解经纬度//int Bdu,Bfen,Ldu,Lfen;double Bmiao,Lmiao;Ldu=int(LL0/pi*180);Lfen=int((LL0/pi*180)*60-Ldu*60);Lmiao=LL0/pi*180*3600-Ldu*3600-Lfen*60;Bdu=int(BB0/pi*180);Bfen=int((BB0/pi*180)*60-Bdu*60);Bmiao=BB0/pi*180*3600-Bdu*3600-Bfen*60;//将弧度转化为角度//printf("\n所选坐标系的转换结果：%d度%d分%lf秒%d 度%d分%lf秒\n",Bdu,Bfen,Bmiao,Ldu,Lfen,Lmiao);fprintf(w,"%d°%d’%lf”%d°%d’%lf”\n",Bdu,Bfen,Bmiao,Ldu,Lfen,Lmiao);//将结果输出到文本文件//}fclose(r);fclose(w);system("pause");return 0;}double sort(double a,double b){double e,e1;e=sqrt(1-(b/a)*(b/a));e1=sqrt((a/b)*(a/b)-1);double m0,m2,m4,m6,m8;double a0,a2,a4,a6,a8;m0=a*(1-e*e);m2=3*e*e*m0/2;m4=5*e*e*m2/4;m6=7*e*e*m4/6;m8=9*e*e*m6/8;a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;a2=m2/2+m4/2+15*m6/32+7*m8/16;a4=m4/8+3*m6/16+7*m8/32;a6=m6/32+m8/16;a8=m8/128;B1=g[i].x/a0;do{F=-a2*sin(2*B1)/2+a4*sin(4*B1)/4-a6*sin(6*B1)/6+a8*sin(8*B1 )/8;B2=(g[i].x-F)/a0;B3=B1;B1=B2;} while(fabs(B3-B2)>10e-10);//利用迭代算法求解底点纬度//return B2; }。

matlab大地测量高斯投影正反算程序设计实验

F=-(a2.*sin(2.*B0))./2+(a4.*sin(4.*B0))./4-(a6.*sin(6.*B0))./6+(a8.*sin(8.*B0))./8; B=(x-F)./a0; if abs(B0-B)<10.^-6 break; end abs(B0-B) B0=B; end end 弧度转度分秒 rad2dms 函数如下： function dms=rad2dms(azimuth)%弧度转度分秒
取整数， L 为当地经度； L0=6.*M-3 ; else %3 度带： M=round(L./3); 为当地经度； L0=3.*M; end l=(L-L0).*pi/180; N=a./(sqrt(1-(e^2).*(sin(B).^2))); t=tan(B); p=e_.*cos(B); %p 表示 yita %则中央子午线经度 L0=3 × M %带号 M=round(L/3)，即对(L/3)的值四舍五入取整数，L %则中央子午线经度 L0=6 × M-3
高斯投影正反算
一、实验目的编写高斯投影正反算的程序，并对格式化文件数据进行计算，验证程序。二、实验内容： 1、高斯投影正算公式已知大地坐标（B,L）及中央子午线经度 L0，计算高斯平面坐标（ x，y），公式如下：
xX
N N sin( B) cos( B)l 2 sin( B) cos 3 ( B)(5 t 2 9 2 4 4 )l 4 2 24
720M f N
4 6 (61 90t 2 f 45t f ) y
L L0
1 1 2 3 y (1 2t 2 f f ) y 3 N f cos B f 6 N f cos B f
1 4 2 2 2 5 (5 28t 2120 N cos B f

大地电磁测深一维正反演(附matlab代码)

大地电磁测深一维正反演摘要本文推导了大地电磁测深的理论计算表达式，并以水平层状介质为例，利用推导的正演计算式在MATLAB 软件平台上进行正演，比较了不同层介质参数的视电阻率曲线。

简要介绍了阻尼最小二乘法反演的基本原理和反演迭代步骤，并对多种层介质进行了反演。

关键词大地电磁，一维正反演，阻尼最小二乘法1 引言20世纪50年代初，苏联学者吉洪诺夫和法国学者卡尼亚的经典著作奠定了大地电磁测深法（MT ）的基础。

它是利用大地仲频率范围很宽（4410~10-Hz ）广泛分布的天然变化的电磁场，进行深部地质构造研究的一种频率域电磁测深法。

由于该法不需要人工建立场源，装备轻便、成本低，且具有比人工源频率测深法更大的勘探深度，所以除主要用于研究地壳和上地幔地质构造外，也常被用来进行油气勘查、地热勘探以及地震预报等研究工作。

几十年来，由于大地电磁测深法具有以下几个优点：不受高阻屏蔽，对低阻分辨率高；不用人工供电，勘探成本低且工作方便；勘探深度范围大。

使大地电磁法在矿产勘探及普查、地壳岩石圈电性结构研究、海洋地球物理勘探、地热勘探、能源勘探、隐伏岩溶水结构、天然地震预测等都扮演着至关重要的角色。

大地电磁也存在一些缺点，比如在实际应用的过程中整理后的数据存在分散的情况；频率范围不够宽，特别是缺少高频成分，受噪音影响大信噪比低；所需观察时间长，致使野外工作效率低。

随着基础理论、技术手段、仪器设备的不断完善和发展，进一步改进和解决这些问题，才能将大地电磁法更好的应用于生产服务当中。

2 视电阻率及水平地层大地电磁测深曲线的理论计算方法 2.1大地电磁测深理论的几点假设和论证吉洪诺夫和卡尼亚提出了假设并论证了以下几点：①将场源近似地看为平面电磁波垂直入射大地。

②引入波阻抗的概念（Z=E/H ），表征地球电性分布对大地电磁场的响应。

③利用单点大地电磁场观测研究地球电性分布是可能的。

2.2视电阻率及水平地层上的理论计算表达式视电阻率概念是从均匀介质中电阻率和波阻抗关系引申出来的。

基于matlab的高斯投影正反算与相邻带坐标换算程序设计

基于matlab的⾼斯投影正反算与相邻带坐标换算程序设计基于matlab的⾼斯投影正反算与相邻带坐标换算程序设
计
作者：徐翰;周强波
作者机构：核⼯业⼆三○研究所,湖南长沙410000;核⼯业⼆三○研究所,湖南长沙410000
来源：中国⽔运（下半⽉）
ISSN：1006-7973
年：2015
卷：015
期：002
页码：72-73
页数：2
中图分类：TP391
正⽂语种：chi
关键词：⾼斯投影;正算;反算;matlab
摘要：地图投影⽅法众多,常⽤的地图投影主要是⾼斯-克吕格投影和UTM投影,UTM投影属于⾼斯投影族,故对于坐标投影的坐标正反算的学习研究以⾼斯投影正反算为主,⽂中主要利⽤matlab的强⼤软件编程功能,对常⽤的⾼斯投影正反算与相邻带坐标换算进⾏程序化设计.。