下载文档原格式
大地主题解算深度干货超精在我们生活的世界上，大地是我们所处的基本环境。
它是被压迫、被掠夺的，也是我们努力保护和改善的对象。
深度解算是一种精确的测量和计算方法，可以帮助我们更好地了解和利用大地资源。
本文将介绍大地主题解算和其带来的深度干货。
一、什么是大地主题解算？大地主题解算是一种基于卫星观测数据和大地测量学原理的技术，用于精确计算大地物体的三维位置和形态变化。
通过对不同时刻的卫星图像进行分析和计算，可以得到地面物体的精确坐标、高程等信息。
大地主题解算的核心思想是通过计算大地物体的三维位置和形态变化，来揭示地球表面的变化过程。
这种技术可以用于监测地表的沉降、地壳的变形、建筑物的结构和变化等。
它在测绘、地质灾害监测、城市规划等领域具有广泛的应用前景。
二、大地主题解算的应用领域1. 地质灾害监测地质灾害是世界各地普遍存在的问题，如地震、山体滑坡、地裂缝等。
通过大地主题解算技术，可以实时监测和预测地质灾害的发生和演变趋势。
这对于提前采取相应的防灾措施、保护人民的生命财产安全具有重要意义。
2. 城市规划随着城市化进程的不断加快，城市规划也变得愈发重要。
大地主题解算可以为城市规划提供精确的地面信息，包括土地利用、道路交通、建筑物布局等。
这有助于提高城市规划的科学性和有效性，减少资源浪费和环境污染。
3. 地表沉降监测地表沉降是地下水开采、矿井开采等人类活动造成的一个重要问题。
通过大地主题解算技术，可以实时监测地表的沉降情况，并对地下水开采和矿井开采等活动进行合理调整和管理。
这有助于减少地表沉降对城市建设和生态环境的不利影响。
4. 环境保护大地主题解算可以为环境保护提供准确的数据支持。
例如，通过对森林面积、湿地面积等进行监测和计算，可以及时发现和预防环境变化和破坏。
这对于保护生态环境和维护生态平衡具有重要意义。
三、大地主题解算的优势和挑战大地主题解算作为一种先进的测量和计算技术，具有很多优势。
首先，它可以实时获取和分析大地物体的位置和形态变化，具有高精度和高时效性。
大地主题解算一、实验目的：1.提高运用计算机语言编程开发的能力；2.加深对大地主题解算计算公式及辅助参数的理解并掌握计算步骤；3.通过编程语言实现大地主题解算。
二、工具：Windows XP Mode 环境下的Microsoft Visual C++ 6.0三、注意事项：1.计算所需变量多，容易混淆；2.正反算函数的编写；3.函数调用；4.弧度与角度之间的转化。
四、实验要求：1.提交报告，实验总结，编写代码；2.独立编程，调试运行；3.上交成果：编写思想，编写过程，问题分析，源代码，计算结果；五、编程过程实现：1.对白塞尔法大地主题解算有一定的了解，并参考教材P148-P150；2.由于参数较多，而在C语言环境下很多符号无法定义，需要符合要求的定义符号替代书本上那些无法直接在C语言环境下定义的符号来达到实现实验的目的；3.程序中采用弧度与度分秒之间转换的函数定义与调用，减轻一定的实验麻烦；4.在C语言环境下，数学函数fabs代替abs起绝对值作用，atan代替arctan起反函数作用；5.程序中尤其注意弧度与角度之间转换，在C语言环境下电脑默认为弧度。
六、源程序代码：#include<stdio.h>#include<math.h>double hudu(double,double,double); /*度分秒转换为弧度*/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;int bx2,by2,lx2,ly2,ax2,ay2;doublee2,W1,sinu1,cosu1,sinA0,coto1,sin2o1,cos2o1,sin2o,cos2o,A,B,C,r,t,o0,o,g,si nu2,q;/*以度分秒顺序输入数据*/printf("请输入大地线起点纬度度分秒\n");scanf("%lf%lf%lf",&bx,&by,&bz);printf("请输入大地线起点经度度分秒\n");scanf("%lf%lf%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)*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-sinA0*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-sin2o1))*sinA0;/*求B2*/sinu2=sinu1*cos(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=-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 && 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%d\n%d\n%lf\n",bx2,by2,bz2);printf("大地线终点经度度分秒分别为：\n%d\n%d\n%lf\n",lx2,ly2,lz2);printf("终点大地方位角度分秒分别为：\n%d\n%d\n%lf\n",ax2,ay2,az2);}/*反算*/else if(k==2){doublebz,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,si nA0,x,t1,t2,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("请输入大地线终点经度度分秒\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=sin(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)A1=fabs(A1);else if(p>0 && 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-sinA0*sinA0))*0.0000000001;t2=(28189-94*(1-sinA0*sinA0))*0.0000000001;g2=(t1*o-t2*x*sino)*sinA0;/*检验循环次数*/printf("\ng2=%lf\ng0=%lf\n",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-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)A2=PI+fabs(A2);elseA2=2*PI-fabs(A2);/*调用函数*/ax=(int)(du(A1));ay=(int)(fen(A1));az=miao(A1);ax2=(int)(du(A2));ay2=(int)(fen(A2));az2=miao(A2);printf("起点大地方位角度分秒分别为：\n%d\n%d\n%lf\n",ax,ay,az);printf("终点大地方位角度分秒分别为：\n%d\n%d\n%lf\n",ax2,ay2,az2);printf("大地线长度为：%lf\n",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;}/*弧度转换为分*/double fen(double C0){double _y,y0;_y=(int)(C0*180/PI);y0=(fabs)((int)((C0*180/PI-_y)*60));return y0;}/*弧度转换为秒*/double miao(double D0){double _z1,_z2,z0;_z1=(int)(D0*180/PI);_z2=(int)((D0*180/PI-_z1)*60);z0=(fabs)((double)(((D0*180/PI-_z1)*60-_z2)*60));return z0;}大地主题解算正算：大地主题解算反算：以上检验数据来自书本P151，P152。
用C语言实现大地主题解算
裴连磊
【期刊名称】《价值工程》
【年(卷),期】2013(32)20
【摘要】This paper expounds the direct and inverse computation thoughts of Bessel solution of geodetic problem, and uses the C language to realize the direct and inverse computation of Bessel solution of geodetic problem.% 本文阐述了白塞尔法大地主题正算和反算思想，最后用C语言实现了白塞尔大地主题解算的正反算。
【总页数】2页(P235-235,236)
【作者】裴连磊
【作者单位】新疆地矿局测绘大队，乌鲁木齐830017
【正文语种】中文
【中图分类】TP311.1
【相关文献】
1.高斯平均引数大地主题解算程序设计 [J], 田桂娥;谢露;马广涛
2.基于大地主题解算方法的无人机偏航距修正探讨 [J], 岑铭
3.高斯大地主题反解算法在海洋平台上的应用 [J], 刘洋
4.新型大地坐标系中的大地主题解算 [J], 施一民;范业明;朱紫阳
5.大地主题常微分方程组解算的数值方法r——以MathCAD为工具 [J], 刘大海;申自立
因版权原因，仅展示原文概要，查看原文内容请购买。
大地主题的数值解法范业明1,王解先1,2,刘慧芹1(1 同济大学测量与国土信息工程系,上海 200092;2 现代工程测量国家测绘局重点实验室,上海 200092)摘要:本文基于椭球面上大地线的微分方程,将法截弧方位角与大地线方位角之间的关系作为初始条件,通过用数值方法求解大地线的微分方程,进行大地主题的正反解。
并以实际数据验证了其正确性与可行性。
本法均采用封闭公式计算,精度高,公式简单,特别适用于计算机解算。
关键词:大地线;微分方程;数值解法中图分类号:P226+1文献标识码:BAbstract :Based on differential equation of geodesic on the surface of ellipsoid,the problem of direct and inverse solution of geodetic is solved through numerical method.The relationship of normal arc and geodesic is deduced and introduced as the initial condition for the differential equation.Numerical experiments are given,and the validity and feasibility of the method proposed in this paper have been proved.Besides,all the formulas are close,simple and easy to be realized by computer.Key words :geodesic;differential equation;numerical method 收稿日期:2006 02 15;修订日期:2006 05 10作者简介:范业明(1982-),男(汉族),辽宁沈阳人,硕士研究生.0 引言一直以来由于计算工具的限制,大地主题解算一般都采用级数展开形式,如短距离的高斯平均引数法,长距离的贝塞尔-赫尔默特方法等等。
1.垂线偏差：地面一点上的重力向量g和相应椭球面上的法线向量n之间的夹角定义为该点的垂线偏差。
2.参考椭球：具有确定参数（长半径a和扁率α），经过局部定位和定向，同某一地区大地水准面最佳拟合的地球椭球，叫参考椭球。
3.大地线：椭球面上两点间的最短程曲线叫做大地线。
4.力高：水准面在纬度45度处的正常高。
5.大地主题解算：已知某些大地元素推求另一些大地元素的计算工作叫大地主题解算。
6.大地主题正算：已知P１点的大地坐标(L１，Ｂ１），Ｐ１至P２的大地线长S及其大地方位角，计算P２点的大地坐标(L２，Ｂ２）和大地线S在P２点的反方位角Ａ２１，这类问题叫做大地主题正算。
7.大地基准：是指能够最佳拟合地球形状的地球椭球的参数及椭球定位和定向8.高斯投影：横轴椭圆柱等角投影（假象有一个椭圆柱横套在地球椭球体外，并与某一条子午线相切，椭球柱的中心轴通过椭球体中心，然后用一定投影方法，将中央子午线两侧各一定范围内的地区投影到椭圆柱上，再将此柱面展开成投影面）。
9.大地测量学：是指在一定的时间与空间参考系中，测量和描绘地球形状及其重力场并监测其变化，为人类活动提供关于地球的空间信息的一门科学。
10.理论闭合差：由水准面不平行而引起的水准环线闭合差，称为理论闭合差。
11.地心坐标系：地心坐标系是在大地体内建立的O-XYZ坐标系。
原点O设在大地体的质量中心，用相互垂直的X，Y，Z三个轴来表示，X轴与首子午面与赤道面的交线重合，向东为正。
Z轴与地球旋转轴重合，向北为正。
Y 轴与XOZ平面垂直构成右手系。
12.高斯投影正、反算公式进行换带计算的步骤。
这种方法的实质是把椭球面上的大地坐标作为过度坐标。
首先把某投影带内有关点的平面坐标（x，y）1利用高斯投影反算公式换算成椭球面上的大地坐标（B，l），进而得到L=L0+l，然后再由大地坐标(B,l)，利用投影正算公式换算成相邻带的平面坐标（x，y）2在计算时，要根据第2带的中央子午线来计算经差l，亦即此时l=L-L0。
大地测量实验报告实验名称：白塞尔大地主题解算（正算和反算）实验目的：1.通过编写白塞尔大地主题电算程序进一步掌握白塞尔法解算大地主题的基本思想。
2.熟练掌握将椭球面上的大地元素按照白塞尔投影条件投影到辅助球面上，继而在球面上进行大地主题解算，最后再将球面上的计算结果换算到椭球面上的基本方法和步骤。
3.学会掌握计算机编程的基本能力。
实验环境：Microsoft Visual C++注意事项：1.在编写程序的过程当中要注意代码的前后统一和重复。
2.注意数值类型的转换和度分秒的换算。
实验步骤：正算：1.计算起点的规划纬度2.计算辅助函数值3.计算系数A,B,C及d,e.4计算球面长度5.计算经度差改正数6.计算终点大地坐标及大地方位角。
反算:1.进行计算辅助函数值2.用逐次趋近法同时计算起点大地方位角、球面长度及经差。
3.计算系数A,B,C及大地线长度S.4.计算反方位角及确定符号。
程序源代码：正算：#include<stdio.h>#include<math.h>#define ee 0.006694384999588#define I 3.141592653double F(double,double,double);void main(void){double A1,B1,L1,S,A2,B2,L2; double x1,x2,x3,y1,y2,y3,z1,z2,z3; double W1,sinu1,sinu2,cosu1,sinA0;doublecota1,cos2a1,sin2a1,cosA0A0;double A,B,C,d,e,a0,a1,m;double n,a,Q,R;printf("请输入数据B1= "); scanf("%lf %lf %lf",&x1,&x2,&x3);B1=F(x1,x2,x3);printf("请输入数据L1= "); scanf("%lf %lf %lf",&y1,&y2,&y3);L1=F(y1,y2,y3);printf("请输入A1= ");scanf("%lf %lf %lf",&z1,&z2,&z3);A1=F(z1,z2,z3);printf("请输入S= ");scanf("%lf",&S);printf("B1=%f\n",B1);printf("L1=%f\n",L1);printf("A1=%f\n",A1);printf("S=%f\n",S);/*计算起点的规划纬度*/W1=sqrt(1-ee*sin(B1)*sin(B1));sinu1=sin(B1)*sqrt(1-ee)/W1;cosu1=cos(B1)/W1;printf("W1=%f\n",W1);printf("sinu1=%f\n",sinu1);printf("cosu1=%f\n",cosu1);/*计算辅助函数值*/sinA0=cosu1*sin(A1);cota1=cosu1*cos(A1)/sinu1;sin2a1=2*cota1/(cota1*cota1+1);cos2a1=(cota1*cota1-1)/(cota1*cota1+1);printf("sinA0=%f\n",sinA0);printf("cota1=%f\n",cota1);printf("sin2a1=%f\n",sin2a1);printf("cos2a1=%f\n",cos2a1);/*计算系数ABC及de*/cosA0A0=1-sinA0*sinA0;A=6356755.288+(10710.341-(13.534*cosA0A0))*cosA0A0;B=(5355.171-9.023*cosA0A0)*cosA0A0;C=(2.256*(cosA0A0))*cosA0A0+0.006;d=691.46768-(0.58143-0.00144*cosA0A0)*cosA0A0;e=(0.2907-cosA0A0*0.0010)*cosA0A0;printf("cosA0A0=%f\n",cosA0A0);printf("A=%f\n",A);printf("B=%f\n",B);printf("C=%f\n",C);printf("d=%f\n",d);printf("e=%f\n",e);/*计算球面长度*/a0=(S-(B+C*cos2a1)*sin2a1)/A;m=sin2a1*cos(2*a0)+cos2a1*sin(2*a0);n=(cos2a1)*(cos(2*a0))-(sin2a1)*(sin(2*a0));a=a0+((B+5*C*n))*m/A;printf("a0=%f\n",a0);printf("m=%f\n",m);printf("n=%f\n",n);printf("a=%f\n",a);/*计算经度差改正数*/Q=(d*a+(e*(m-sin2a1))/3600/180*I)*sinA0;printf("Q=%f\n",Q);/*计算终点大地坐标及大地方位角*/sinu2=sinu1*cos(a)+cosu1*cos(A1)*sin(a);B2=180*atan(sinu2/((sqrt(1-ee))*(sq rt(1-sinu2*sinu2))))/I;R=180*atan(sin(A1)*sin(a)/(cosu1*co s(a)-sinu1*sin(a)*cos(A1)))/I;printf("sinu2=%f\n",sinu2);printf("B2=%f\n",B2);printf("R=%f\n",R*180/I);/*确定R的值*/if(sin(A1)>0 && tan(R)>0)R=abs(R);else if(sin(A1)>0 && tan(R)<0)R=I-abs(R);else if(sin(A1)<0 && tan(R)<0)R=-abs(R);elseR=abs(R)-I;/*确定L2A2的值*/L2=(L1*180/I+R-(Q/206265*180/I));A2=atan(cosu1*sin(A1)/(cosu1*cos(a) *cos(A1)-sinu1*sin(a)));if(sin(A1)<0&&tan(A2)>0)A2=(fabs(A2))*180/I;else if(sin(A1)<0&&tan(A2)<0)A2=(I-fabs(A2))*180/I;else if(sin(A1)>0&&tan(A2)>0)A2=(I+fabs(A2))*180/I;elseA2=(2*I-fabs(A2))*180/I;printf("A2=%3f\n B2=%3f\nL2=%3f\n",A2,B2,L2);}double F(double a2,double b2,doublec2){double d2;d2=(double)(a2+1.0*b2/60+1.0*c2/3600);d2=(d2/180)*I;return (d2);}注：A1，B1，L1，S分别为大地线起点的大地方位角、纬度、经度、大地线长;B2,L2,A2为大地线终点纬度、经度及方位角。
大地主题解算(正算)代码与白塞尔大地主题解算大地主题解算(正算)代码：根据经纬度和方向角以及距离计算另外一点坐标新建模块->拷贝下面的大地主题(正算)代码，调用方法示例：起点经度：116.235(度)终点纬度：37.435(度)方向角：50（度）长度：500（米）终点经纬度（"经度,纬度"）=Computation(37.435,116.235,50,500)Const Pi = 3.1415926535898Private a, b, c, alpha, e, e2, W, V As Double'a 长轴半径'b 短轴'c 极曲率半径'alpha 扁率'e 第一偏心率'e2 第二偏心率'W 第一基本纬度函数'V 第二基本纬度函数Private B1, L1, B2, L2 As Double'B1 点1的纬度'L1 点1的经度'B2 点1的纬度'L2 点2的经度Private S As Double '''''大地线长度Private A1, A2 As Double'A1 点1到点2的方位角'A2 点2到点1的方位角Function Computation(STARTLAT, STARTLONG, ANGLE1, DISTANCE As Double) As String B1 = STARTLATL1 = STARTLONGA1 = ANGLE1S = DISTANCEa = 6378245b = 6356752.3142c = a ^ 2 / balpha = (a - b) / ae = Sqr(a ^ 2 - b ^ 2) / ae2 = Sqr(a ^ 2 - b ^ 2) / bB1 = rad(B1)L1 = rad(L1)A1 = rad(A1)W = Sqr(1 - e ^ 2 * (Sin(B1) ^ 2))V = W * (a / b)Dim W1 As DoubleE1 = e ''''第一偏心率'// 计算起点的归化纬度W1 = W ''Sqr(1 - e1 * e1 * Sin(B1 ) * Sin(B1 )) sinu1 = Sin(B1) * Sqr(1 - E1 * E1) / W1cosu1 = Cos(B1) / W1'// 计算辅助函数值sinA0 = cosu1 * Sin(A1)cotq1 = cosu1 * Cos(A1)sin2q1 = 2 * cotq1 / (cotq1 ^ 2 + 1)cos2q1 = (cotq1 ^ 2 - 1) / (cotq1 ^ 2 + 1)'// 计算系数AA,BB,CC及AAlpha, BBeta的值。
cos2A0 = 1 - sinA0 ^ 2e2 = Sqr(a ^ 2 - b ^ 2) / bk2 = e2 * e2 * cos2A0Dim aa, BB, CC, EE22, AAlpha, BBeta As Doubleaa = b * (1 + k2 / 4 - 3 * k2 * k2 / 64 + 5 * k2 * k2 * k2 / 256)BB = b * (k2 / 8 - k2 * k2 / 32 + 15 * k2 * k2 * k2 / 1024)CC = b * (k2 * k2 / 128 - 3 * k2 * k2 * k2 / 512) e2 = E1 * E1AAlpha = (e2 / 2 + e2 * e2 / 8 + e2 * e2 * e2 / 16) - (e2 * e2 / 16 + e2 * e2 * e2 / 16) * cos2A0 + (3 * e2 * e2 * e2 / 128) * cos2A0 * cos2A0BBeta = (e2 * e2 / 32 + e2 * e2 * e2 / 32) *cos2A0 - (e2 * e2 * e2 / 64) * cos2A0 * cos2A0'// 计算球面长度q0 = (S - (BB + CC * cos2q1) * sin2q1) / aasin2q1q0 = sin2q1 * Cos(2 * q0) + cos2q1 * Sin(2 * q0)cos2q1q0 = cos2q1 * Cos(2 * q0) - sin2q1 * Sin(2 * q0)q = q0 + (BB + 5 * CC * cos2q1q0) * sin2q1q0 / aa'// 计算经度差改正数theta = (AAlpha * q + beta * (sin2q1q0 - sin2q1)) * sinA0'// 计算终点大地坐标及大地方位角sinu2 = sinu1 * Cos(q) + cosu1 * Cos(A1) * Sin(q)B2 = Atn(sinu2 / (Sqr(1 - E1 * E1) * Sqr(1 - sinu2 * sinu2))) * / Pilamuda = Atn(Sin(A1) * Sin(q) / (cosu1 * Cos(q) - sinu1 * Sin(q) * Cos(A1))) * / PiIf (Sin(A1) > 0) ThenIf (Sin(A1) * Sin(q) / (cosu1 * Cos(q) - sinu1 * Sin(q) * Cos(A1)) > 0) Then lamuda = Abs(lamuda)Elselamuda = - Abs(lamuda)End IfElseIf (Sin(A1) * Sin(q) / (cosu1 * Cos(q) - sinu1 * Sin(q) * Cos(A1)) > 0) Then lamuda = Abs(lamuda) -Elselamuda = -Abs(lamuda)End IfEnd IfL2 = L1 * / Pi + lamuda - theta * / PiA2 = Atn(cosu1 * Sin(A1) / (cosu1 * Cos(q) * Cos(A1) - sinu1 * Sin(q))) * / PiIf (Sin(A1) > 0) ThenIf (cosu1 * Sin(A1) / (cosu1 * Cos(q) * Cos(A1) - sinu1 * Sin(q)) > 0) ThenA2 = + Abs(A2)ElseA2 = 360 - Abs(A2)End IfElseIf (cosu1 * Sin(A1) / (cosu1 * Cos(q) * Cos(A1) - sinu1 * Sin(q)) > 0) ThenA2 = Abs(A2)ElseA2 = - Abs(A2)End IfEnd IfComputation = L2 & "," & B2End FunctionPrivate Function rad(ByValangle_d As Double) As Doublerad = angle_d * Pi /End Function白塞尔大地主题解算一、正算代码：#include<stdio.h>#include<math.h>#define ee 0.006693421622966#define I 3.141592653double F(double,double,double);void main(void){double A1,B1,L1,S,A2,B2,L2;double x1,x2,x3,y1,y2,y3,z1,z2,z3;double W1,sinu1,sinu2,cosu1,sinA0;double cota1,cos2a1,sin2a1,cosA0A0;double A,B,C,d,e,a0,a1,m;doublen,a,Q,R;printf("请输入数据B1= ");scanf("%lf %lf %lf",&x1,&x2,&x3);B1=F(x1,x2,x3);printf("请输入数据L1= ");scanf("%lf %lf %lf",&y1,&y2,&y3);L1=F(y1,y2,y3);printf("请输入A1= ");scanf("%lf %lf %lf",&z1,&z2,&z3);A1=F(z1,z2,z3);printf("请输入S= ");scanf("%lf",&S);printf("B1=%.9f\n",B1);printf("L1=%.9f\n",L1);printf("A1=%.9f\n",A1);printf("S=%f\n",S);/*计算起点的规划纬度*/W1=sqrt(1-ee*sin(B1)*sin(B1));sinu1=sin(B1)*sqrt(1-ee)/W1;cosu1=cos(B1)/W1;printf("W1=%.9f\n",W1);printf("sinu1=%.9f\n",sinu1);printf("cosu1=%.9f\n",cosu1);/*计算辅助函数值*/sinA0=cosu1*sin(A1);cota1=cosu1*cos(A1)/sinu1;sin2a1=2*cota1/(cota1*cota1+1);cos2a1=(cota1*cota1-1)/(cota1*cota1+1);printf("sinA0=%.9f\n",sinA0);printf("cota1=%.9f\n",cota1);printf("sin2a1=%.9f\n",sin2a1);printf("cos2a1=%.9f\n",cos2a1);/*计算系数ABC及de*/cosA0A0=1-sinA0*sinA0;A=6356755.288+(10710.341-(13.534*cosA0A0))*cosA0A0;B=(5355.171-9.*cosA0A0)*cosA0A0;C=(2.256*(cosA0A0))*cosA0A0+0.006;d=691.46768-(0.58143-0.00144*cosA0A0)*cosA0A0;e=(0.2907-cosA0A0*0.0010)*cosA0A0;printf("cosA0A0=%.9f\n",cosA0A0);printf("A=%.3f\n",A);printf("B=%.6f\n",B);printf("C=%.9f\n",C);printf("d=%.7f\n",d);printf("e=%.9f\n",e);/*计算球面长度*/a0=(S-(B+C*cos2a1)*sin2a1)/A;m=sin2a1*cos(2*a0)+cos2a1*sin(2*a0);n=(cos2a1)*(cos(2*a0))-(sin2a1)*(sin(2*a0));a=a0+((B+5*C*n))*m/A;printf("a0=%.9f\n",a0);printf("m=%.9f\n",m);printf("n=%.9f\n",n);printf("a=%.9f\n",a);/*计算经度差改正数*/Q=(d*a+(e*(m-sin2a1))/3600/*I)*sinA0;printf("Q=%.7f\n",Q);/*计算终点大地坐标及大地方位角*/sinu2=sinu1*cos(a)+cosu1*cos(A1)*sin(a);B2=*atan(sinu2/((sqrt(1-ee))*(sqrt(1-sinu2*sinu2))))/I; R=*atan(sin(A1)*sin(a)/(cosu1*cos(a)-sinu1*sin(a)*cos(A1)))/I;printf("sinu2=%.9f\n",sinu2);printf("B2=%f\n",B2);printf("R=%f\n",R*/I);/*确定R的值*/if(sin(A1)>0 && tan(R)>0)R=abs(R);else if(sin(A1)>0 && tan(R)<0)R=I-abs(R);else if(sin(A1)<0 && tan(R)<0)R=-abs(R);elseR=abs(R)-I;/*确定L2A2的值*/L2=(L1*/I+R-(Q/206265*/I));A2=atan(cosu1*sin(A1)/(cosu1*cos(a)*cos(A1)-sinu1*sin(a)));if(sin(A1)<0&&tan(A2)>0)A2=(fabs(A2))*/I;else if(sin(A1)<0&&tan(A2)<0)A2=(I-fabs(A2))*/I;else if(sin(A1)>0&&tan(A2)>0)A2=(I+fabs(A2))*/I;elseA2=(2*I-fabs(A2))*/I;printf("A2=%3f\n B2=%3f\n L2=%3f\n",A2,B2,L2);}double F(double a2,double b2,double c2){double d2;d2=(double)(a2+1.0*b2/60+1.0*c2/3600);d2=(d2/)*I;return (d2);}运行结果：二、反算代码：#include<stdio.h>#include<math.h>#define ee 0.006693421622966#define I 3.141592653double F(double,double,double);void main(void){double B1,B2,L1,L2,A1,A2,S,Y;double W1,W2,L,Q,R,A,B,C;doublex,y,z,p,q;double x1,x2,x3,y1,y2,y3,z1,z2,z3,w1,w2,w3;double a1,a2,b1,b2,m,n;double sinp,cosp,sinu1,sinu2,cosu1,cosu2,sinA0,cosA0;Q=0;q=0;printf("请输入起始数据B1= ");scanf("%lf %lf %lf",&x1,&x2,&x3);B1=F(x1,x2,x3);printf("请输入起始数据L1= ");scanf("%lf %lf %lf",&y1,&y2,&y3);L1=F(y1,y2,y3);printf("请输入起始数B2= ");scanf("%lf %lf %lf",&z1,&z2,&z3);B2=F(z1,z2,z3);printf("请输入起始数据L2= ");scanf("%lf %lf %lf",&w1,&w2,&w3);L2=F(w1,w2,w3);printf("B1=%f\n",B1);printf("L1=%.9f\n",L1);printf("B2=%.9f\n",B2);printf("L2=%.9f\n",L2);/*辅助计算*/W1=sqrt(1-ee*sin(B1)*sin(B1));W2=sqrt(1-ee*sin(B2)*sin(B2));sinu1=sin(B1)*sqrt(1-ee)/W1;sinu2=sin(B2)*(sqrt(1-ee))/W2;cosu1=cos(B1)/W1;cosu2=cos(B2)/W2;L=L2-L1;a1=sinu1*sinu2;a2=cosu1*cosu2;b1=cosu1*sinu2;b2=sinu1*cosu2;printf("W1=%.9f\n",W1);printf("W2=%.9f\n",W2);printf("sinu1=%.9f\n",sinu1);printf("sinu2=%.9f\n",sinu2);printf("cosu1=%.9f\n",cosu2);printf("L=%.9f\n",L);printf("a1=%.9f\n",a1);printf("a2=%.9f\n",a2);printf("b1=%.9f\n",b1);printf("b2=%.9f\n",b2);/*用逐次趋近法同时计算起点大地方位角、球面长度及经差R*/ do{R=L+Q;x=cosu2*sin(R);y=b1-b2*cos(R);A1=atan(x/y);if(x>0&&y>0)A1=fabs(A1);else if(x>0&&y<0)A1=I-fabs(A1);else if(x<0&&y<0)A1=I+fabs(A1);elseA1=2*I-fabs(A1);sinp=x*sin(A1)+y*cos(A1);cosp=a1+a2*cos(R);p=atan(sinp/cosp);if(cosp>0)p=fabs(p);sinA0=cosu1*sin(A1);cosA0=sqrt(1-sinA0*sinA0);p=I-fabs(p);z=2*a1-cosA0*cosA0*cos(p);m=(33523299-(28189-70*cosA0*cosA0))*(1e-10);n=(28189-94*cosA0*cosA0)*(1e-10);Q=q;q=(m*p-m*z*sin(p))*sinA0;}while(fabs(206265*(q-Q))>(1e-4));printf("Q=%f\n",Q);/*计算系数ABC及大地线长度*/A=6356863.020+(10708.949-13.474*cosA0*cosA0)*cosA0*cosA0;B=10708.938-17.956*cosA0*cosA0;C=4.487;Y=(cosA0*cosA0*cosA0*cosA0-2*z*z)*cos(p);S=A*p+(B*z+C*Y)*sinp;/*计算反方位角*/A2=atan((cosu1*sin(R))/(b1*cos(R)-b2));A2=A2*/I;A1=A1*/I;if(A1>0)A2=abs(A2);elseA2=-abs(A2);printf("A1=%3f\n A2=%3f\n S=%3f\n",A1,A2,S);}double F(double a,doubleb,double c){double d;if(a<0)d=(double)(a-1.0*b/60-1.0*c/3600);elsed=(double)(a+1.0*b/60+1.0*c/3600);d=(d/)*I;return d; } 运行结果：。
