大地主题解算

1、精品文档大地主题解算一、实验目的：1. 提高运用计算机语言编程开发的能力；2. 加深对大地主题解算计算公式及辅助参数的理解并掌握计算步骤；3. 通过编程语言实现大地主题解算。二、工具：Windows XP Mode 环境下的 Microsoft Visual C+ 6.0 三、 注意事项：1. 计算所需变量多，容易混淆；2. 正反算函数的编写；3. 函数调用；4. 弧度与角度之间的转化。四、 实验要求：1. 提交报告，实验总结，编写代码；2. 独立编程，调试运行；3. 上交成果：编写思想，编写过程，问题分析，源代码，计算结果；五、 编程过程实现：1.对白塞尔法大地主题解算有一定的了解，并参考教

2、材P148-P150；2. 由于参数较多，而在 C语言环境下很多符号无法定义，需要符合要求的定义符号替代书本上那些无法直接在 C语言环境下定义的符号来达到实现实验的目的；3. 程序中采用弧度与度分秒之间转换的函数定义与调用，减轻一定的实验麻烦；4. 在 C 语言环境下，数学函数 fabs 代替 abs 起绝对值作用， atan 代替 arctan 起反函数作用；5. 程序中尤其注意弧度与角度之间转换，在 C语言环境下电脑默认为弧度。六、源程序代码：#include<stdio.h>#include<math.h>double hudu(double,double,dou

3、ble); /*度分秒转换为弧度 */double du(double);/*弧度转换为度 */double fen(double);/*弧度转换为分 */double miao(double);/*弧度转换为秒 */#define PI 3.1415926void main (void)int k;printf(" 请选择执行正算或者反算，若执行正算，请输入1；若执行反算，请输入 2。 n");scanf("%d",&k);/*正算 */if(k=1)。1欢迎下载精品文档doublebz,lz,az,S,bz2,lz2,az2,B1,L1,A1

4、,B2,L2,A2,bx,by,lx,ly,ax,ay;intbx2,by2,lx2,ly2,ax2,ay2;doublee2,W1,sinu1,cosu1,sinA0,coto1,sin2o1,cos2o1,sin2o,cos2o,A,B,C,r,t,o0,o,g,sinu2,q;/* 以度分秒顺序输入数据 */printf("请输入大地线起点纬度度分秒n");scanf("%lf%lf%lf",&bx,&by,&bz);printf("请输入大地线起点经度度分秒n");scanf("%lf%lf%

5、lf",&lx,&ly,&lz);printf("请输入大地方位角度分秒n");scanf("%lf%lf%lf",&ax,&ay,&az);printf("请输入大地线长度 n");scanf("%lf",&S);/* 调用函数 */B1=hudu(bx,by,bz);L1=hudu(lx,ly,lz);A1=hudu(ax,ay,az);/* 白塞尔大地主题解算 */e2=0.006693421622966;W1=sqrt(1-e2*sin(B1

6、)*sin(B1);sinu1=sin(B1)*(sqrt(1-e2)/W1;cosu1=cos(B1)/W1;sinA0=cosu1*sin(A1);coto1=cosu1*cos(A1)/sinu1;sin2o1=2*coto1/(coto1*coto1+1);cos2o1=(coto1*coto1-1)/(coto1*coto1+1);A=6356863.020+(10718.949-13.474*(1-sinA0*sinA0)*(1-sinA0*sinA0);B=(5354.469-8.978*(1-sinA0*sinA0)*(1-sinA0*sinA0);C=(2.238*(1-si

7、nA0*sinA0)*(1-sinA0*sinA0)+0.006;r=691.46768-(0.58143-0.00144*(1-sinA0*sinA0)*(1-sinA0*sinA0); t=(0.2907-0.0010*(1-sinA0*sinA0)*(1-sinA0*sinA0);o0=(S-(B+C*cos2o1)*sin2o1)/A;sin2o=sin2o1*cos(2*o0)+cos2o1*sin(2*o0);cos2o=cos2o1*cos(2*o0)-sin2o1*sin(2*o0);o=o0+(B+5*C*cos2o)*sin2o/A;g=(r*o+t*(sin2o-sin2

8、o1)*sinA0;/* 求 B2*/sinu2=sinu1*cos(o)+cosu1*cos(A1)*sin(o);B2=atan(sinu2/(sqrt(1-e2)*sqrt(1-sinu2*sinu2);。2欢迎下载精品文档/*求 L2*/q=atan(sin(A1)*sin(o)/(cosu1*cos(o)-sinu1*sin(o)*cos(A1);/* 判断 q*/if(sin(A1)>0 && tan(q)>0)q=fabs(q);else if(sin(A1)>0 && tan(q)<0)q=PI-fabs(q);else

9、if(sin(A1)<0 && tan(q)<0)q=-fabs(q);elseq=fabs(q)-PI;L2=L1+q-g/3600/180*PI;/*求 A2*/A2=atan(cosu1*sin(A1)/(cosu1*cos(o)*cos(A1)-sinu1*sin(o);/* 判断 A2*/if(sin(A1)<0 && tan(A2)>0)A2=fabs(A2);else if(sin(A1)<0 && tan(A2)<0)A2=PI-fabs(A2);else if(sin(A1)>0 &am

10、p;& tan(A2)>0)A2=PI+fabs(A2);elseA2=2*PI-fabs(A2);/* 调用函数 */bx2=(int)(du(B2);by2=(int)(fen(B2);bz2=miao(B2);lx2=(int)(du(L2);ly2=(int)(fen(L2);lz2=miao(L2);ax2=(int)(du(A2);ay2=(int)(fen(A2);az2=miao(A2);printf("大地线终点纬度度分秒分别为：n%dn%dn%lfn",bx2,by2,bz2);printf("大地线终点经度度分秒分别为：n%dn

11、%dn%lfn",lx2,ly2,lz2);printf("终点大地方位角度分秒分别为：n%dn%dn%lfn",ax2,ay2,az2);/* 反算 */else if(k=2)double。3欢迎下载精品文档bz,lz,bz2,lz2,az,az2,B1,L1,B2,L2,S,A1,A2,bx,by,lx,ly,bx2,by2,lx2,ly2 ;int ax,ay,ax2,ay2;doublee2,W1,W2,sinu1,sinu2,cosu1,cosu2,L,a1,a2,b1,b2,g,g2,g0,r,p,q,sino,coso,o,sinA0,x,t1,t

12、2,A,B,C,y;/* 以度分秒顺序输入数据 */printf("请输入大地线起点纬度度分秒n");scanf("%lf%lf%lf",&bx,&by,&bz);printf("请输入大地线起点经度度分秒n");scanf("%lf%lf%lf",&lx,&ly,&lz);printf("请输入大地线终点纬度度分秒n");scanf("%lf%lf%lf",&bx2,&by2,&bz2);printf(

13、"请输入大地线终点经度度分秒n");scanf("%lf%lf%lf",&lx2,&ly2,&lz2);/*调用函数 */B1=hudu(bx,by,bz);L1=hudu(lx,ly,lz);B2=hudu(bx2,by2,bz2);L2=hudu(lx2,ly2,lz2);/* 白塞尔大地主题解算 */e2=0.006693421622966;W1=sqrt(1-e2*sin(B1)*sin(B1);W2=sqrt(1-e2*sin(B2)*sin(B2);sinu1=sin(B1)*sqrt(1-e2)/W1;sinu2=s

14、in(B2)*sqrt(1-e2)/W2;cosu1=cos(B1)/W1;cosu2=cos(B2)/W2;L=L2-L1;a1=sinu1*sinu2;a2=cosu1*cosu2;b1=cosu1*sinu2;b2=sinu1*cosu2;/* 逐次趋近法求解A1*/g0=-662.904266/3600*PI/180;g=0;r=L;while(1)p=cosu2*sin(r);q=b1-b2*cos(r);A1=atan(p/q);/*判断 A1*/if(p>0 && q>0)。4欢迎下载精品文档A1=fabs(A1);else if(p>0 &am

15、p;& q<0)A1=PI-fabs(A1);else if(p<0 && q<0)A1=PI+fabs(A1);elseA1=2*PI-fabs(A1);sino=p*sin(A1)+q*cos(A1);coso=a1+a2*cos(r);o=atan(sino/coso);/* 判断 o*/if(coso>0)o=fabs(o);elseo=PI-fabs(o);sinA0=cosu1*sin(A1);x=2*a1-(1-sinA0*sinA0)*coso;t1=(33523299-(28189-70*(1-sinA0*sinA0)*(1-s

16、inA0*sinA0)*0.000000000 1;t2=(28189-94*(1-sinA0*sinA0)*0.0000000001;g2=(t1*o-t2*x*sino)*sinA0;/* 检验循环次数 */printf("ng2=%lfng0=%lfn",g2,g0);if(g2<=g0)break;elser=L+g2;/* 求解 S*/A=6356863.020+(10708.949-13.474*(1-sinA0*sinA0)*(1-sinA0*sinA0); B=10708.938-17.956*(1-sinA0*sinA0);C=4.487;y=(1-

17、sinA0*sinA0)*(1-sinA0*sinA0)-2*x*x)*coso;S=A*o+(B*x+C*y)*sino;/* 求解 A2*/A2=atan(cosu1*sin(r)/(b1*cos(r)-b2);/* 判断 A2*/if(p<0 && q<0)A2=fabs(A2);else if(p<0 && q>0)A2=PI-fabs(A2);else if(p>0 && q>0)。5欢迎下载精品文档A2=PI+fabs(A2);elseA2=2*PI-fabs(A2);/* 调用函数 */ax=(i

18、nt)(du(A1);ay=(int)(fen(A1);az=miao(A1);ax2=(int)(du(A2);ay2=(int)(fen(A2);az2=miao(A2);printf("起点大地方位角度分秒分别为：n%dn%dn%lfn",ax,ay,az);printf("终点大地方位角度分秒分别为：n%dn%dn%lfn",ax2,ay2,az2);printf("大地线长度为： %lfn",S);/* 数据错误 */elseprintf("数据错误，请重新输入n");/* 度分秒转换为弧度 */double hudu(double a0,doubl