




版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领
文档简介
1、大地主题解算一、 实验目的:1. 提高运用计算机语言编程开发的能力;2. 加深对大地主题解算计算公式及辅助参数的理解并掌握计算步骤;3. 通过编程语言实现大地主题解算。二、 工具:Windows XP Mode 环境下的Microsoft Visual C+ 6.0三、 注意事项:1. 计算所需变量多,容易混淆;2. 正反算函数的编写;3. 函数调用;4. 弧度与角度之间的转化。四、 实验要求:1. 提交报告,实验总结,编写代码;2. 独立编程,调试运行;3. 上交成果:编写思想,编写过程,问题分析,源代码,计算结果;五、 编程过程实现:1. 对白塞尔法大地主题解算有一定的了解,并参考教材P1
2、48-P150;2. 由于参数较多,而在C语言环境下很多符号无法定义,需要符合要求的定义符号替代书本上那些无法直接在C语言环境下定义的符号来达到实现实验的目的;3. 程序中采用弧度与度分秒之间转换的函数定义与调用,减轻一定的实验麻烦;4. 在C语言环境下,数学函数fabs代替abs起绝对值作用,atan代替arctan起反函数作用;5. 程序中尤其注意弧度与角度之间转换,在C语言环境下电脑默认为弧度。六、源程序代码:#include<stdio.h>#include<math.h>double hudu(double,double,double); /*度分秒转换为弧度
3、*/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) double bz,lz,az,S,bz2,lz2,az2,B1,L1,A1,B2,L2,A2,bx,by,lx,ly,ax,ay
4、;intbx2,by2,lx2,ly2,ax2,ay2;double e2,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%lf",&lx,&ly,&am
5、p;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)*sin(B1);sinu1=sin(B1)*(sqrt(1-e
6、2)/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-sinA0*sinA0)*(1-sinA0*sinA0)+0.006
7、; 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-sin2o1)*sinA0;/*求B2*/sinu2=sinu1*c
8、os(o)+cosu1*cos(A1)*sin(o);B2=atan(sinu2/(sqrt(1-e2)*sqrt(1-sinu2*sinu2); /*求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 if(sin(A1)<0 && tan(q)<0)q=-fa
9、bs(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 && tan(A2)>0)A2=PI+fabs(A2);elseA2=2*P
10、I-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%dn%lfn",lx2,ly2,lz2);printf("终点大地方
11、位角度分秒分别为:n%dn%dn%lfn",ax2,ay2,az2);/*反算*/else if(k=2) double 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;double e2,W1,W2,sinu1,sinu2,cosu1,cosu2,L,a1,a2,b1,b2,g,g2,g0,r,p,q,sino,coso,o,sinA0,x,t1,t2,A,B,C,y;/*以度分秒顺序输入数据*/printf("请输入大地线起点纬度度分秒n&quo
12、t;);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("请输入大地线终点经度度分秒n");scanf("%lf%lf%lf",&
13、amp;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=sin(B2)*sqrt(1-e2)/W2;cosu1=cos(B1)/W1;cosu2=cos(B2)/W2;L=L
14、2-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) A1=fabs(A1);else if(p>0 && q<0)A1=PI-fabs(A1);else if(p<0 && q<0)A1=
15、PI+fabs(A1);else A1=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-sinA0*sinA0)*0.0000000001;t2=(28189-94*(1-sinA0*sinA0)*0.000000000
16、1;g2=(t1*o-t2*x*sino)*sinA0;/*检验循环次数*/printf("ng2=%lfng0=%lfn",g2,g0);if(g2<=g0) break;else r=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-sinA0*sinA0)*(1-sinA0*sinA0)-2*x*x)*coso;S=A*o+(B*x+C*y)*sino;/*求解A2
17、*/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)A2=PI+fabs(A2);else A2=2*PI-fabs(A2);/*调用函数*/ ax=(int)(du(A1);ay=(int)(fen(A1);az=miao(A1);ax2=(int)(du(A2);ay2=(int)(fen(A2);az2=m
18、iao(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,double b0,double c0) double A0;A0=(a0+b0/60+c0/3600)*PI/180; return A0;/*弧度转换为度*/double du(double B0)double x0;x0=(int)(B0*180/PI);return x0;/*弧度转换为分*/doub
温馨提示
- 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
- 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
- 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
- 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
- 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
- 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
- 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。
最新文档
- 2025届北京市对外经贸大学附属中学物理高一第二学期期末调研试题含解析
- 山东省邹城第一中学2025年物理高一第二学期期末达标检测试题含解析
- 车辆租赁服务合同协议条款确认书
- 2025至2030巨细胞病毒感染行业市场深度研究与战略咨询分析报告
- 数字智慧方案5813丨智慧工厂实现自动化与效率的全新方案
- 谅解赔偿协议书范本
- 养老机构护理员培训课件
- 新媒体投资协议书范本
- 地下室置换协议书范本
- 消渴病的中医护理查房
- 《别墅设计任务书》word版
- EN 4644-001-2017(高清正版)
- BICC协议介绍
- 公铁联运物流园区及配套项目建议书写作模板
- 预应力混凝土简支T形梁桥毕业论文
- 变频器变频altivar71说明书
- 第三章_同步发电机励磁自动调节
- 农村饮用水工程监理规划
- WBS BOM操作手册
- 铁路文物保护管理暂行办法
- 博士后基金申请的经验及体会PPT课件
评论
0/150
提交评论