基于matlab的大地坐标与直角坐标间的转换精编版

matlab直角坐标和经纬度的换算摘要：I.引言- 介绍MATLAB 软件- 说明直角坐标和经纬度换算的重要性II.MATLAB 中直角坐标和经纬度的转换方法- 利用MATLAB 内置函数进行转换- 利用MATLAB 进行手动计算转换III.实际应用案例- 使用MATLAB 进行经纬度换算的实际案例- MATLAB 在地理信息系统中的应用IV.结论- 总结MATLAB 在直角坐标和经纬度换算中的作用- 展望MATLAB 在地理科学领域的未来发展正文：MATLAB 是一种功能强大的数学软件，被广泛应用于科学计算、数据分析、图像处理等领域。

在地理科学领域，MATLAB 也具有重要的应用价值，可以用于进行直角坐标和经纬度的换算。

在MATLAB 中，经纬度和直角坐标之间的转换可以通过使用MATLAB 内置的函数来实现。

具体来说，可以使用MATLAB 中的`geocode`函数将经纬度转换为直角坐标，使用`reverse`函数将直角坐标转换为经纬度。

这些函数的使用方法简单，易于操作，可以帮助用户快速完成坐标转换。

除了使用MATLAB 内置函数进行转换外，用户也可以手动计算经纬度和直角坐标之间的转换。

这种方法需要用户掌握一定的数学知识，例如地理坐标系和直角坐标系之间的转换公式。

通过手动计算，用户可以更深刻地理解坐标转换的原理，更好地掌握MATLAB 的使用方法。

在实际应用中，MATLAB 的经纬度和直角坐标换算功能被广泛应用于地理信息系统(GIS) 中。

GIS 是一种以采集、存储、管理、分析和应用地理信息为主要任务的技术系统，可以用于制作地图、分析地理数据、规划城市等方面。

MATLAB 可以与GIS 软件相结合，帮助用户更好地完成地理数据的处理和分析。

综上所述，MATLAB 在直角坐标和经纬度的换算中发挥着重要作用。

通过使用MATLAB，用户可以方便地进行地理数据的处理和分析，更好地理解和应用地理科学知识。

基于 MATLAB 的七参数法坐标系统转换问题分析 1张鲜妮 21, ,王磊 21,1、中国矿业大学环境与测绘学院,江苏徐州 (2210082、江苏省资源环境信息工程重点实验室,江苏徐州 (221008E-mail:摘要:GPS 测量的坐标是基于 WGS-84坐标系下的,而我国实用的测量成果大多都是基于北京 54坐标系下的。

随着 GPS 测量技术的广泛使用,由 WGS-84坐标向北京 54坐标系下坐标的转换问题一直是一个可探讨的问题, 坐标系统转换的现有模型很多, 但常用的还是经典的七参数转换模型。

随着不断的实践研究, 发现七参数在进行坐标系统转换时有一定的局限性。

本文采用 MATLAB 语言编写了七参数法坐标系统转换程序,并对七参数坐标系统转换的若干问题进行了分析讨论。

分析结果表明, 小区域范围内用正常高代替大地高对坐标转换精度影响很小; 公共点分布情况对坐标转换精度影响显著; 合适的公共点密度有利于提高坐标转换精度。

关键词:七参数法;坐标系统; MATLAB ;转换问题1. 引言随着 GPS 空间定位技术的发展, GPS 技术以其快速、精确、全天候在测量中的应用变的越来越广泛, GPS 成为建立基础控制网的首选手段 ]1[,由于 GPS 系统采用的是 WGS-84坐标系, 是一种地心坐标系, 而我国目前常用的两个坐标系 1954年北京坐标系 (以下称 BJ54 和 1980年国家大地坐标系,是一种参心坐标系,采用克拉所夫斯基椭球为参考椭球,并采用高斯克吕格投影方式进行投影, 我国的国土测量成果和在进行工程施工时大都是基于这两个坐标系下的。

所以在利用 GPS 技术进行测量过程中必然存在由 WGS-84坐标向北京 54坐标系下的转换问题。

现有的转换模型已经成熟,归纳起来主要有布尔莎 -沃尔夫模型(七参数法、莫洛登斯基 -巴代卡斯模型和范士模型 ]2[。

本文主要分析讨论是基于七参数转换模型, 分析工具是MATLAB 软件。

文档从互联网中收集，已重新修正排版，word格式支持编辑，如有帮助欢迎下载支持。

2008届毕业生毕业论文（设计）题目：基于MATLAB的坐标转换程序设计摘要本文主要阐述了基于MATLAB的坐标转换程序设计与实现的问题。

论述以MATLAB为开发平台和编程语言，设计出解决工程测量中常见的坐标转换问题的程序。

坐标转换一直是专业性强且不易解决的问题,针对目前坐标转换软件功能单一、操作不方便等问题，采用窗口、菜单、控件的操作方式，实现了所见即所得的人性化界面设计。

程序的设计主要从两个方面进行，其一保证程序有较高的转换精度，其二是友好的界面设计。

程序的运行能满足工程测量人员对坐标转换运算和坐标数据分析的需要，程序实现了不同参考椭球情况下七参数和四参数的计算过程、不同坐标系统的坐标转换和换带计算程序化。

论文还诠释了测量坐标转换的含义和内容，针对坐标转换基本模型的选用、转换参数的解算、转换计算的方法、转换计算中值得注意的问题加以研究和探讨，以便实现在测量实践和理论中各类不同坐标之间的转换计算。

关键词:坐标转换,换带,参考椭球,MATLAB,图形用户界面AbstractThis article expatiates the design and implementation of a computing program for coordinate conversion, operation of MATLAB. With programmed language, the article designs the program of solving the common coordinate conversion problems in the engineering survey, which regards MATLAB as an exploitation basis. coordinate conversion is a professional problem which is difficult to solve, to solve the existed problems , the operating modes of windows, menus and widgets are adopted. Moreover, the WYSWYG humanized program designs are realized. The program designs come from two aspects. Firstly, the powerful operation function of the program is guaranteed. Secondly, the visualization is designed. The program operation meets the needs which engineering survey personnel need to have the coordinate conversion operation and data analysis. Meanwhile, the program designs the coordinate conversion function, including coordinate conversion among different coordinate systems and between two projection zones, realizing the computation of 4 parameters as well as 7 parameters under the coordinates among different coordinate systems. Above all, the article includes the meaning and content of transformation, basic model selection of coordinates transformation, calcu1ation of transformation parameters, calculation method of transformation and problems existing in transformation. Calculation are researched and discussed in this paper in order to measure transformation calculation of different coordinate in practice and theory．Key words:Coordinate conversion, Stripe exchange, Reference ellipsoid, MATLAB,GUI目录前言..................................................... 错误!未定义书签。

基于matlab的大地坐标与直角坐标间的转换LT编写大地坐标与空间直角坐标相互转换的程序，并对格式化文件数据进行计算，验证程序。

二、实验内容：1、大地坐标向空间直角坐标换算转换公式：B h e N z LB h N y LB h N x sin ])1([sin cos )(cos cos )(2+-=+=+= （1）其中：L 为经度，B 为纬度，h 为大地高，B e aN 22sin 1-=为卯酉圈曲率半径，a b a e 22-=为第一偏心率，a 为旋转椭球长半轴，b 为短半轴。

WGS84椭球参数：长半轴 a = 6378137 扁率 f =1/298.257223563根据上式创建以geo 2xyz 命名的函数，函数输入输出格式为[x, y, z] = geo 2xyz (L, B, h)2、空间直角坐标向大地坐标换算根据式（1）推导大地坐标向空间直角坐标转换公式：N By x h yx B Ne z B x y L -+=++==cos )sin arctan()/arctan(22222 注意计算纬度时需要用到迭代，可用)arctan(22yx b az B +=作为初始值。

创建以xyz2geo 命名的函数，函数输入输出格式为[L, B, h] = xyz 2geo (x, y, z)三、实验步骤1、大地坐标向空间直角坐标换算主程序：%%大地坐标向空间直角坐标换算%函数的输入输出格式为[x,y,z]=geo2xyz(L,B,h)[filename,pathname] = uigetfile('*.txt','请选择打开的数据文件'); file = [pathname, filename];data = importdata(file);L=data.data(:,1);B=data.data(:,2);h=data.data(:,3);[x,y,z]=geo2xyz(L,B,h);A=[x,y,z];A=A';[filename_out,pathname_out] = uiputfile('*.txt','请选择要输出数据文件');fileout = [pathname_out, filename_out];fid = fopen(fileout,'wt');fprintf(fid,' x y z\n'); fprintf(fid,'%15.7f %15.7f %15.7f\n',A);close('all');函数：function [x,y,z]=geo2xyz(L,B,h)%大地坐标经纬度转换成空间直角坐标B=dms2rad(B);L=dms2rad(L);a=6378137;%a是长半轴f=1/298.257223563;%f是扁率b=a-a*f;e=sqrt(a^2-b^2)/a;N=a./(sqrt(1-e^2.*(sin(B)).^2));%N为卯酉圈半径率，e为第一偏心率x=(N+h).*cos(B).*cos(L);y=(N+h).*cos(B).*sin(L);z=(N*(1-e^2)+h).*sin(B);endfunction rad=dms2rad(jiaodu)%度分秒－>弧度（rad）degree = fix(jiaodu);mimute = fix((jiaodu-degree)*100);second = (jiaodu-degree-mimute/100)*10000;degree = degree+mimute/60+second/3600;rad=degree/180*pi;end2、空间直角坐标向大地坐标换算主程序：%% 将文件data.2.txt中的空间直角坐标系转换为大地坐标，并将计算结果按照格式存储在文件data3.txt中%data3.txt格式为：经度（ddmmss）纬度（ddmmss）大地高[filename,pathname]=uigetfile('*.txt','请选择打开的数据文件');file=[pathname,filename];data=importdata(file);x=data.data(:,1);y=data.data(:,2);z=data.data(:,3);[L,B,H]=xyz2geo(x,y,z);[filename_out,pathname_out] = uiputfile('*.txt','请选择要输出数据文件');fileout = [pathname_out, filename_out];fid = fopen(fileout,'wt');fprintf(fid,' 经度(ddmmss）纬度(ddmmss）大地高\n');fprintf(fid,'%10.7f %10.7f %7.3f\n',[L,B,H]'); fclose('all');函数：function [L,B,H]=xyz2geo(x,y,z)%将直角坐标转换为大地坐标%% 已知:WGS-84椭球参数f=1/298.257223563;%扁率f=(a-b)/aa=6378137; %长半轴b=a*(1-f); %短半轴e=(sqrt(a^2-b^2))/a; %第一偏心率%% 经度L的计算L=atan2(y,x);L=rad2dms(L);%% 纬度B的计算B0=atan2((a*z),(b.*sqrt(x.^2+y.^2))); % B初始值while 1N=a./(sqrt(1-(e^2).*(sin(B0).^2))); %卯酉圈曲率半径NB=atan2(z+N.*(e^2).*sin(B0),(sqrt(x.^2+y.^2)));if abs(B0-B)<10^-6break;endabs(B0-B)B0=B;end%% 大地高H的计算H=(sqrt(x.^2+y.^2))./cos(B)-N;B=degree2dms(B.*180/pi);%纬度BEndfunction dms= degree2dms(jiaodu)%度－>度分秒（dd.mmss）degree = fix(jiaodu);mimute = fix((jiaodu-degree)*60); second = ((jiaodu-degree)*60-mimute)*60; dms = degree+mimute/100+second/10000;。

地球坐标转换与matlab应用地球坐标转换与Matlab应用1. 引言地球坐标转换是地理信息系统中重要的一部分，它涉及到将地球上的点的位置从一种坐标系统转换成另一种坐标系统。

在时间和空间领域中，地球坐标转换具有广泛的应用，包括导航系统、地图制图、地图投影等。

而Matlab作为一款功能强大的数学软件，可以用于处理和分析地理空间数据，被广泛应用于地球坐标转换中。

本文将以地球坐标转换与Matlab应用为主题，介绍地球坐标转换的基本概念、常用的转换方法和Matlab在地球坐标转换中的应用案例。

2. 地球坐标系统简介地球坐标系统是用于描述地球上点的位置的一种坐标系统。

常用的地球坐标系统主要包括经纬度坐标系统和投影坐标系统。

经纬度坐标系统使用经度和纬度来确定地球上某个点的位置，经度表示东西方向的位置，纬度表示南北方向的位置。

投影坐标系统是经纬度坐标系统的扁平化表示，它使用投影方式将地球上的曲面投影到平面上，以方便地图制图等应用。

3. 地球坐标转换方法地球坐标转换可以通过多种方法进行，根据不同的需求和应用场景选择合适的方法进行转换。

常用的地球坐标转换方法包括经纬度与投影坐标的相互转换、不同投影坐标系之间的转换等。

3.1 经纬度与投影坐标的转换经纬度与投影坐标之间的转换是地球坐标转换中常见的任务。

其中，经纬度转投影坐标可以利用地图投影算法实现，常见的投影算法包括墨卡托投影、UTM投影等。

而投影坐标转经纬度则需要进行反投影计算，将投影坐标转换回经纬度。

3.2 不同投影坐标系的转换在地理信息系统中，常常需要将数据由一个投影坐标系转换到另一个投影坐标系，以适应不同的应用需求。

这种转换可以在Matlab中使用相关的函数进行处理，例如Matlab中的`projfwd`和`projinv`函数可以实现不同投影坐标系之间的相互转换。

4. Matlab中的地球坐标转换应用Matlab是一款功能强大的数学软件，它提供了丰富的工具箱和函数，可以用于处理地理空间数据以及进行地球坐标转换。

#include "stdio、h"#include "math、h"#include "stdlib、h"#include "iostream"#define PI 3、14323double a,b,c,e2,ep2;int main(){int m,n,t;double RAD(double d,double f,double m);void RBD(double hd);void BLH_XYZ();void XYZ_BLH();void B_ZS();void B_FS();void GUS_ZS();void GUS_FS();printf(" 大地测量学\n");sp1:printf("请选择功能:\n");printf("1、大地坐标系到大地空间直角坐标的转换\n");printf("2、大地空间直角坐标到大地坐标系的转换\n");printf("3、贝塞尔大地问题正算\n");printf("4、贝塞尔大地问题反算\n");printf("5、高斯投影正算\n");printf("6、高斯投影反算\n");printf("0、退出程序\n");scanf("%d",&m);if(m==0)exit(0);sp2:printf("请选择椭球参数(输入椭球序号):\n");printf("1、克拉索夫斯基椭球参数\n");printf("2、IUGG_1975椭球参数\n");printf("3、CGCS_2000椭球参数\n");printf("0、其她椭球参数(自行输入)\n");scanf("%d",&n);switch(n){case 1:a=6378245、0;b=6356863、0188;c=6399698、9018;e2=0、297;ep2=0、468;break;case 2:a=6378140、0;b=6356755、2882;c=6399596、6520;e2=0、959;ep2=0、947;break;case 3:a=6378137、0;b=6356752、3141;c=6399593、6259;e2=0、290;ep2=0、547;break;case 0:{printf("请输入椭球参数:\n");printf("长半径a=");scanf("%lf",&a);printf("短半径b=");scanf("%lf",&b);c=a*a/b;ep2=(a*a-b*b)/(b*b);e2=(a*a-b*b)/(a*a);break;}default:printf("\n\n输入错误!\n请重新输入!\n\n");goto sp2 ;}while(1){switch(m){case 1:BLH_XYZ();break;case 2:XYZ_BLH();break;case 3:B_ZS();break;case 4:B_FS();break;case 5:GUS_ZS();break;case 6:GUS_FS();break;default:printf("\n\n输入错误!\n请重新输入!\n\n");goto sp1 ;}printf("就是否继续进行此功能计算? \n\n");printf("( 若继续进行此功能计算,则输入1;\n 若选择其她功能进行计算,则输入2;\n 若退出, 则输入0、)\n");scanf("%d",&t);switch(t){case 1:break;case 2:goto sp1;case 0:exit(0);}}}double RAD(double d,double f,double m) {double e;double sign=(d<0、0)?-1、0:1、0;if(d==0){sign=(f<0、0)?-1、0:1、0;if(f==0){sign=(m<0、0)?-1、0:1、0;}}if(d<0)d=d*(-1、0);if(f<0)f=f*(-1、0);if(m<0)m=m*(-1、0);e=sign*(d*3600+f*60+m)*PI/(3600*180);return e;}void RBD(double hd){int t;int d,f;double m;double sign=(hd<0、0)?-1、0:1、0;if(hd<0)hd=fabs(hd);hd=hd*3600*180/PI;t=int(hd/3600);d=sign*t;hd=hd-t*3600;f=int(hd/60);m=hd-f*60;printf("%d'%d'%lf'\n",d,f,m);}void BLH_XYZ(){double B,L,H,N,W;double d,f,m;double X,Y,Z;printf(" 请输入大地坐标(输入格式为角度(例如:30'40'50')):\n");printf(" 大地经度L=");scanf("%lf'%lf'%lf'",&d,&f,&m);L=RAD(d,f,m);printf(" 大地纬度B=");scanf("%lf'%lf'%lf'",&d,&f,&m);B=RAD(d,f,m);printf(" 大地高H=");scanf("%lf",&H);W=sqrt(1-e2*sin(B)*sin(B));N=a/W;X=(N+H)*cos(B)*cos(L);Y=(N+H)*cos(B)*sin(L);Z=(N*(1-e2)+H)*sin(B);printf("\n\n 转换后得到大地空间直角坐标为:\n\n");printf("X=%lf\nY=%lf\nZ=%lf\n\n",X,Y,Z);}void XYZ_BLH(){double B,L,H,N,W;double X,Y,Z;double tgB0,tgB1;printf(" 请输入大地空间直角坐标:\n");printf(" X=");scanf("%lf",&X);printf(" Y=");scanf("%lf",&Y);printf(" Z=");scanf("%lf",&Z);printf("\n\n 转换后得到大地坐标为:\n\n");L=atan(Y/X);printf(" 大地经度为: L=");RBD(L);printf("\n");tgB0=Z/sqrt(X*X+Y*Y);tgB1=(1/sqrt(X*X+Y*Y))*(Z+a*e2*tgB0/sqrt(1+tgB0*tgB0-e2*tgB0*tgB0));while(fabs(tgB0-tgB1)>5*pow(10,-10)){tgB0=tgB1;tgB1=(1/sqrt(X*X+Y*Y))*(Z+a*e2*tgB0/sqrt(1+tgB0*tgB0-e2*tgB0*tgB0));}B=atan(tgB1);printf(" 大地纬度为:B=");RBD(B);printf("\n");W=sqrt(1-e2*sin(B)*sin(B));N=a/W;H=sqrt(X*X+Y*Y)/cos(B)-N;printf(" 大地高为:H=%lf\n\n",H);}void B_ZS(){double L1,B1,A1,s,d,f,mi;double u1,u2,m,M,k2,alfa,bt,r,kp2,alfap,btp,rp;double sgm0,sgm1,lmd,lmd1,lmd2,A2,B2,l,L2;printf("请输入已知点的大地坐标(输入格式为角度(例如:30'40'50'),下同):\nL1=");scanf("%lf'%lf'%lf'",&d,&f,&mi);L1=RAD(d,f,mi);printf("\nB1=");scanf("%lf'%lf'%lf'",&d,&f,&mi);B1=RAD(d,f,mi);printf("请输入大地方位角:\nA1=");scanf("%lf'%lf'%lf'",&d,&f,&mi);A1=RAD(d,f,mi);printf("请输入该点至另一点的大地线长:\ns=");scanf("%lf",&s);u1=atan(sqrt(1-e2)*tan(B1));m=asin(cos(u1)*sin(A1));M=atan(tan(u1)/cos(A1));m=(m>0)?m:m+2*PI;M=(M>0)?M:M+PI;k2=ep2*cos(m)*cos(m);alfa=(1-k2/4+7*k2*k2/64-15*k2*k2*k2/256)/b;bt=k2/4-k2*k2/8+37*k2*k2*k2/512;r=k2*k2/128-k2*k2*k2/128;sgm0=alfa*s;sgm1=alfa*s+bt*sin(sgm0)*cos(2*M+sgm0)+r*sin(2*sgm0)*cos(4*M+2*sgm0); while(fabs(sgm0-sgm1)>2、8*PI/180*pow(10,-7)){sgm0=sgm1;sgm1=alfa*s+bt*sin(sgm0)*cos(2*M+sgm0)+r*sin(2*sgm0)*cos(4*M+2*sgm0); }sgm0=sgm1;A2=atan(tan(m)/cos(M+sgm0));A2=(A2>0)?A2:A2+PI;A2=(A1>PI)?A2:A2+PI;u2=atan(-cos(A2)*tan(M+sgm0));lmd1=atan(sin(u1)*tan(A1));lmd1=(lmd1>0)?lmd1:lmd1+PI;lmd1=(m<PI)?lmd1:lmd1+PI;lmd2=atan(sin(u2)*tan(A2));lmd2=(lmd2>0)?lmd2:lmd2+PI;lmd2=(m<PI)?(((M+sgm0)<PI)?lmd2:lmd2+PI):(((M+sgm0)>PI)?lmd2:lmd2+PI);lmd=lmd2-lmd1;B2=atan(sqrt(1+ep2)*tan(u2));kp2=e2*cos(m)*cos(m);alfap=(e2/2+e2*e2/8+e2*e2*e2/16)-e2/16*(1+e2)*kp2+3*e2*kp2*kp2/128;btp=e2*(1+e2)*kp2/16-e2*kp2*kp2/32;rp=e2*kp2*kp2/256;l=lmd-sin(m)*(alfap*sgm0+btp*sin(sgm0)*cos(2*M+sgm0)+rp*sin(2*sgm0)*cos(4*M+2*sg m0));L2=L1+l;printf("\n\n得到另一点的大地坐标与大地线在该点的大地方位角为:\n\n");printf("L2=");RBD(L2);printf("\n");printf("B2=");RBD(B2);printf("\n");printf("A2=");RBD(A2);printf("\n");}void B_FS(){double L1,B1,L2,B2,s,A1,A2,du,f,mi,m0,m,M;double l,u1,u2,alfa,bt,r,lmd0,dit_lmd,lmd,sgm,dit_sgm,sgm0,sgm1,alfap,btp,rp,k2,kp2;printf("请输入第一个点大地坐标(输入格式为角度(例如:30'40'50'),下同):\n大地经度L1=");scanf("%lf'%lf'%lf'",&du,&f,&mi);L1=RAD(du,f,mi);printf("大地纬度B1=");scanf("%lf'%lf'%lf'",&du,&f,&mi);B1=RAD(du,f,mi);printf("\n请输入第二个点大地坐标:\n大地经度:L2=");scanf("%lf'%lf'%lf'",&du,&f,&mi);L2=RAD(du,f,mi);printf("大地纬度:B2=");scanf("%lf'%lf'%lf'",&du,&f,&mi);B2=RAD(du,f,mi);l=L2-L1;u1=atan(sqrt(1-e2)*tan(B1));u2=atan(sqrt(1-e2)*tan(B2));sgm0=acos(sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(l));m0=asin(cos(u1)*cos(u2)*sin(l)/sin(sgm0));dit_lmd=0、003351831*sgm0*sin(m0);lmd0=l+dit_lmd;dit_sgm=sin(m0)*dit_lmd;sgm1=sgm0+dit_sgm;m=asin(cos(u1)*cos(u2)*sin(lmd0)/sin(sgm1));A1=atan(sin(lmd0)/(cos(u1)*tan(u2)-sin(u1)*cos(lmd0)));A1=(A1>0)?A1:A1+PI;A1=(m>0)?A1:A1+PI;M=atan(sin(u1)*tan(A1)/sin(m));M=(M>0)?M:M+PI;k2=ep2*cos(m)*cos(m);alfa=(1-k2/4+7*k2*k2/64-15*k2*k2*k2/256)/b;bt=k2/4-k2*k2/8+37*k2*k2*k2/512;r=k2*k2/128-k2*k2*k2/128;kp2=e2*cos(m)*cos(m);alfap=(e2/2+e2*e2/8+e2*e2*e2/16)-e2/16*(1+e2)*kp2+3*e2*kp2*kp2/128;btp=e2*(1+e2)*kp2/16-e2*kp2*kp2/32;rp=e2*kp2*kp2/256;sgm0=acos(sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(l));sgm1=acos(sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(l+sin(m)*(alfap*sgm0+btp*sin(sgm0)*cos( 2*M+sgm0))));while(fabs(sgm0-sgm1)>1*PI/180*pow(10,-8)){sgm0=sgm1;sgm1=acos(sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(l+sin(m)*(alfap*sgm0+btp*sin(sgm0)*cos( 2*M+sgm0))));}sgm=sgm1;lmd=l+sin(m)*(alfap*sgm+btp*sin(sgm)*cos(2*M+sgm));s=(sgm-bt*sin(sgm)*cos(2*M+sgm)-r*sin(2*sgm)*cos(4*M+2*sgm))/alfa;A1=atan(sin(lmd)/(cos(u1)*tan(u2)-sin(u1)*cos(lmd)));A1=(A1>0)?A1:A1+PI;A1=(m>0)?A1:A1+PI;A2=atan(sin(lmd)/(sin(u2)*cos(lmd)-tan(u1)*cos(u2)));A2=(A2>0)?A2:A2+PI;A2=(m<0)?A2:A2+PI;printf("\n\n得到两点间大地线长S与大地正反方位角A1、A2如下:\n\n");printf("s=%lf\n",s);printf("A1=");RBD(A1);printf("\n");printf("A2=");RBD(A2);printf("\n");}void GUS_ZS(){double B,L,x3,x6,y3,y6,Y3,Y6,du,f,mi,X,N,n,t;double At,Bt,Ct,Dt,m3,m6,l3,l6,W,L03,L06;int DH3,DH6;printf("请输入大地坐标(输入格式为角度(例如:30'40'50')):\n大地经度L=");scanf("%lf'%lf'%lf'",&du,&f,&mi);L=RAD(du,f,mi);printf("\n大地纬度B=");scanf("%lf'%lf'%lf'",&du,&f,&mi);B=RAD(du,f,mi);At=1+3*e2/4+45*e2*e2/64+175*e2*e2*e2/256;Bt=3*e2/4+15*e2*e2/16+525*e2*e2*e2/512;Ct=15*e2*e2/64+105*e2*e2*e2/256;Dt=35*e2*e2*e2/512;X=a*(1-e2)*(At*B-Bt*sin(2*B)/2+Ct*sin(4*B)/4-Dt*sin(6*B)/6);W=sqrt(1-e2*sin(B)*sin(B));N=a/W;n=sqrt(ep2)*cos(B);t=tan(B);DH3=(L-(1、5*PI/180))/(3*PI/180)+1;DH6=L/(6*PI/180)+1;L03=DH3*(3*PI/180);L06=DH6*(6*PI/180)-(3*PI/180);l3=L-L03;l6=L-L06;m3=cos(B)*l3;m6=cos(B)*l6;x3=X+N*t*(m3*m3/2+(5-t*t+9*n*n+4*n*n*n*n)*m3*m3*m3*m3/24+(61-58*t*t+t*t*t*t)* m3*m3*m3*m3*m3*m3/720);x6=X+N*t*(m6*m6/2+(5-t*t+9*n*n+4*n*n*n*n)*m6*m6*m6*m6/24+(61-58*t*t+t*t*t*t)* m6*m6*m6*m6*m6*m6/720);y3=N*(m3+(1-t*t+n*n)*m3*m3*m3/6+(5-18*t*t+t*t*t*t+14*n*n-58*n*n*t*t)*m3*m3*m3 *m3*m3/120);y6=N*(m6+(1-t*t+n*n)*m6*m6*m6/6+(5-18*t*t+t*t*t*t+14*n*n-58*n*n*t*t)*m6*m6*m6 *m6*m6/120);Y3=DH3*1000000+500000+y3;Y6=DH6*1000000+500000+y6;printf("\n\n 得到的高斯平面坐标为:\n\n");printf(" 对于3度带:\n 纵坐标x=%、3lf\n 横坐标y=%、3lf(通用坐标Y=%、3lf)\n\n",x3,y3,Y3);printf(" 对于6度带:\n 纵坐标x=%、3lf\n 横坐标y=%、3lf(通用坐标Y=%、3lf)\n\n",x6,y6,Y6);}void GUS_FS(){double x,y,Y,B,B0,B1,Bf,Vf,tf,Nf,nf,L,At,Bt,Ct,Dt,L3,L6;long DH;printf(" 请输入高斯平面坐标:\n\n");printf(" 纵坐标X=");scanf("%lf",&x);printf("\n");printf(" 自然坐标y=");scanf("%lf",&y);printf("\n");printf(" 通用坐标Y=");scanf("%lf",&Y);printf("\n");At=1+3*e2/4+45*e2*e2/64+175*e2*e2*e2/256;Bt=3*e2/4+15*e2*e2/16+525*e2*e2*e2/512;Ct=15*e2*e2/64+105*e2*e2*e2/256;Dt=35*e2*e2*e2/512;B0=x/(a*(1-e2)*At);B1=(x-a*(1-e2)*(-Bt*sin(2*B0)/2+Ct*sin(4*B0)/4-Dt*sin(6*B0)/6))/(a*(1-e2)*At);while(fabs(B1-B0)>1*pow(10,-8)){B0=B1;B1=(x-a*(1-e2)*(-Bt*sin(2*B0)/2+Ct*sin(4*B0)/4-Dt*sin(6*B0)/6))/(a*(1-e2)*At);}Bf=B1;nf=sqrt(ep2)*cos(Bf);tf=tan(Bf);Vf=sqrt(1+ep2*cos(Bf)*cos(Bf));Nf=c/Vf;B=Bf-Vf*Vf*tf/2*((y/Nf)*(y/Nf)-(5+3*tf*tf+nf*nf-9*nf*nf*tf*tf)*pow((y/Nf),4)/12+(61+90*tf *tf+45*tf*tf)*pow((y/Nf),6)/360);L=((y/Nf)-(1+2*tf*tf+nf*nf)*(y/Nf)*(y/Nf)*(y/Nf)/6+(5+28*tf*tf+24*pow(tf,4)+6*nf*nf+8*nf *nf*tf*tf)*pow((y/Nf),5)/120)/cos(Bf);DH=Y/1000000;L3=3*PI/180*double(DH)+L;L6=6*PI/180*double(DH)-3*PI/180+L;printf("\n\n 得到的大地坐标为:\n\n");printf(" 大地纬度B=");RBD(B);printf("\n");printf(" 若为6度带,大地经度L=");RBD(L6);printf("\n");printf(" 若为3度带,大地经度L=");RBD(L3);printf("\n"); }。