大地测量主题解算白赛尔大地主题解算

大地测量上机实习报告
一、实习时间：第9周星期三
二、实习地点：系办三楼机房
三、实习内容：白塞尔大地主题解算
四、实习目的：（1）尝试编写白塞尔主题解算的程序
（2）通过实例检验解算的结果
五、实习步骤：
1、先了解书本上关于白塞尔大地主题解算的步骤；
2、双击打开白塞尔大地主题解算的程序，如下：
3、在“选择椭球”中选择“克氏椭球”，“选择正反算”中选择“正算”；
4、然后在输入框中输入相应的数据，如下：
5、接下来单击菜单栏中的“数据正确”，记得正算结果，如下所示：
6、重复第3步，选择反算，并输入相应数据：
7、点击“数据正确”，即得反算结果：
六、实习总结：通过白塞尔大地主题结算方法的上机操作，懂得了白塞尔法结算大地主题是将椭球面上的大地元素按照白塞尔投影条件投影到辅助球面上，继而在球面上进行大地主题解算，再将球面上的计算结果换算到椭球面上的原理，并熟悉白塞尔大地主题解算的程序。

白塞尔大地主题解算的基本思想
首先，白塞尔大地主题解算的基本思想之一是建立地壳变形的力学模型。

地壳变形是地球表面的一项重要现象，是由于地质作用和地球内部力
学过程的结果。

地壳变形的力学模型是研究地壳变形的重要方法。

常见的
地壳变形的力学模型有弹性模型、弹塑性模型和拟弹性模型等。

这些模型
可以描述地壳的变形过程，通过对地壳的变形过程进行建模，可以更好地
理解地壳变形的机制和动力学过程。

其次，白塞尔大地主题解算的思想之一是研究现今地球表面的动力学
过程。

地球表面的动力学过程包括板块运动、地震活动、火山喷发等。

这
些过程在地球的长期演化中起着重要的作用。

通过对现今地球表面动力学
过程的研究，可以揭示地球内部的结构和动力学机制，进而更好地理解地
壳变形的成因和发展变化的规律。

第三，白塞尔大地主题解算的基本思想之一是研究地壳变形的成因。

地壳变形的成因是地壳运动的基本原因，也是地球科学研究的一个重要问题。

地壳变形的成因包括构造运动、地壳应力状态的改变、地震活动等。

通过研究地壳变形的成因，可以更好地了解地壳变形的机制和规律，进而
为地震预测和地壳运动的控制提供科学依据。

总之，白塞尔大地主题解算的基本思想是通过建立地壳变形的力学模型，研究现今地球表面的动力学过程，探究地壳变形的成因，揭示地壳变
形与地球动力学过程之间的相互关系，以推动地壳运动的理论和应用研究。

这一思想对于研究地壳运动的机制和规律，了解地壳变形的成因和动力学
过程，提高地震预测和地壳运动控制的水平具有重要意义。

白塞尔大地主题解算课程设计大地测量学课程设计设计题目:白塞尔大地主题解算学院: 矿业学院专业: 测绘工程班级: 测绘122学号: 1208010255学生姓名: 郑英智指导教师: 张俊2014年12月26日目录1. 白塞尔法大地主题正算步骤 12. 白塞尔法大地主题反算步骤 33. 同一平行圈弧长、子午线弧长与大地线比较大小 44. 程序代码 55. 演算示例 96. 参考文献 127. 心得体会 13教师评语 14白塞尔大地主题解算摘要:白塞尔法解算大地主题的基本思想是将椭球面上的大地元素按照白塞尔投影条件(,椭球面大地线投影到球面上为大圆弧;,大地线和大圆弧上相应点的方位角相等;,球面上任意一点的纬度等于椭球面上相应点的归化纬度。

)投影到辅助球面上，继而在球面上进行大地主题解算，最后再将球面上的计算结果换算到椭球面上。

1.白塞尔法大地主题正算步骤2已知、、、、()、S，计算、、()。

BLABLAaAAe1112221221(1)将椭球面元素投影到球面上2utanu,1,etanB求:,由B1111A,计算辅助量和 ,01sinA,cosusinAtan,,tanusecA , 011111, ,计算球面长度，将S化为,,,S，,sin,cos(2,，,)，,sin2,cos(4,，2,) 11式中系数分别为:''''246k7k15k,,,,(1,，,，？) ,bAb464256''246,Bkk37k'',,,,(,，，？) A48512''46,Ckk'',,,,(,，？) A12812822'2k,ecosA 024635kkk ,1，,,，？A464256246kk15kB,,，，？41651246kk3C,,？128512,,,,S上式右端含有，因此需要迭代计算。

第一次迭代取近似值，第二次计算取 0,,,S，,sin,cos(2,，,)，,sin2,cos(4,，2,) 01010- 1 -,以后计算用代换,代入上式迭代计算，直到所要求的精度为止。

大地主题解算(正算)代码与白塞尔大地主题解算大地主题解算(正算)代码：根据经纬度和方向角以及距离计算另外一点坐标新建模块->拷贝下面的大地主题(正算)代码，调用方法示例：起点经度：116.235(度)终点纬度：37.435(度)方向角：50（度）长度：500（米）终点经纬度（"经度,纬度"）=Computation(37.435,116.235,50,500)Const Pi = 3.14Private 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 StringB1 = 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的值。

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为大地线终点纬度、经度及方位角。

大地测量试题集(4.2)一、判断正误(正确的在括号内打√，错误的在括号内打×)● 高斯平均引数正算公式必须采用迭代计算。

（）● 高斯平均引数反算公式必须采用迭代计算。

（）● 白塞尔法大地主题解算是在球面上进行大地主题解算。

（）● 白塞尔法大地主题解算公式的精度与距离长短无关。

（）● 高斯平均引数大地主题解算公式的精度与距离长短无关。

（）● 地图投影的实质是求平面上的点位与椭球面上的点位之间的关系。

（） ● 最大及最小长度比的方向称为主方向，它们在投影平面上是正交的。

（） ● 地图投影长度比近似等于1。

（）● 地图投影长度变形可能为正、为负或等于0。

（）● 等角投影各方向的长度比相同。

（）● 等距离投影各方向的长度比相同。

（）● 高斯投影长度比与方向无关。

（）● 某点的高斯平面坐标x ，y 分别等于在椭球面上该点至赤道和中央子午线的距离。

（）● 正反大地方位角互差180°。

（）● 子午线收敛角当经差l 不变时，随纬度增加而增大。

（）● 高斯投影的距离改化是大地线描写形与其弦线长的差别。

（）二、填空题● 大地主题解算思路之一是以大地线的 为基础进行积分运算，并用级数展开积分式。

● 高斯投影的长度比m 与 有关。

当高斯平面坐标y =0时，长度比m = 。

这是因为 。

● 在长度比的通用公式 中，满足正形条件是： =0， = 。

● 柯西-黎曼微分条件（方程）： 。

● 最大角度变形是最大方向变形的 倍。

● 设a 、b 为变形椭圆的长、短半轴，等角投影的条件是 ，等积投影的条件22222cos sin cos sin E A F A A G A m r ++=是，等距离投影的条件是。

●高斯投影中央子午线投影后为直线，其它经线投影后为向中央子午线的曲线，且长度改变；赤道投影后为一直线，其它纬线投影后为向赤道的曲线；投影后的经线与纬度。

●等量纬度与大地纬度的关系式是d q= dB。

大地主题解算(正算)代码与白塞尔大地主题解算大地主题解算(正算)代码：根据经纬度和方向角以及距离计算另外一点坐标新建模块->拷贝下面的大地主题(正算)代码，调用方法示例：起点经度：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 StringB1 = 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的值。

白塞尔大地主题解算方向：学号：姓名：一．基本思路：基本思想：将椭球面上的大地元素按照白塞尔投影条件投影到辅助球面上，继而在球面上进行大地主题解算，最后在将球面上的计算结果换算到椭球面上。

其关键问题是找出椭球面上的大地元素与球面上相应元素之间的关系式，同时解决在球面上进行大地主题解算的方法。

正算流程：1.计算起点的归化纬度2.计算辅助函数值，解球面三角形3.按公式计算相关系数A,B,C 以及α,β4.计算球面长度5.计算纬度差改正数6.计算终点大地坐标及大地方位角011122S B C A{sin (cos )}σσσ=-+10101022222sin ()sin sin cos cos σσσσσσ+=+10101022222cos ()cos cos sin sin σσσσσσ+=-001101522B C A[cos ()]sin ()σσσσσσ=++++010122L A sin [(sin ()sin )]λδασβσσσ-==++-2111u u u A sin sin cos cos cos sin σσ=+2222222222222222222222111111111e u B u B W W u e B u u B B arctan e u e u sin sin cos cos tan tan sin sin tan cos -cos ⎧-==⎪⎪⎨⎪=-⎪⎩⎡⎤⎢⎥==--⎢⎥-⎣⎦1111A arctan u u A sin sin []cos cos sin sin cos σλσσ=-21L L λδ=+-112111u A A arctan u A u cos sin cos cos cos sin sin σσ⎡⎤=⎢⎥-⎣⎦反算流程：1.辅助计算2.用逐次趋近法同时计算起点大地方位角、球面长度及经差，第一次趋近时，取δ＝０。

计算下式,重复上述计算过程2.3.计算大地线长度S4.计算反方位角二．已知数据序号B1（DD.MMSS)L1 (DD.MMSS)A12(DD.MMSS)S12(m)1 41.01356874 130.10122676 1.4943 8000L λδ=+211212u pA u u u u qsin cos tan cos sin sin cos cos λλ==-2121p p u q b b A arctanqsin cos cos λλ==-=11p A q A sin cos tan cos σσ+=11p A q A sin sin cos σ=+12a a cos cos σλ=+arctan sin cos σσσ⎛⎫=⎪⎝⎭011A u A sin cos sin =111u A tan tan sec σ=21+σσσ=02122L A sin [(sin sin )]λδασβσσ-==+-L λδ=+11222222S A B C B C sin (cos )sin (cos )σσσσσ=++-+1212u A b b cos sin arctan cos λλ⎡⎤=⎢⎥-⎣⎦三．源代码：#include <stdio.h>#include <math.h>#define e 0.081813334016931499 //克拉索夫斯基椭球体第一偏心率void main(){int k,B10,B11,L10,L11,A10,A11,B20,B21,L20,L21,A20,A21;double B12,L12,A12,B22,L22,A22;double B1,L1,A1,S,B2,L2,A2,L,pi;double A,B,C,afa,beta;double a1,a2,b1,b2,p,q,x,y;doubleW1,W2,sinu1,sinu2,cosu1,cosu2,sinA0,cotsigma1,sin2sigma1,cos2sigma1,sigma0,sin2,cos2,sigm a,sins,coss,delta0,delta,lamda;pi=4*atan(1);printf("白塞尔大地主题正算请输入1\n白塞尔大地主题反算请输入2\n");scanf("%d",&k);if(k==1){printf("请输入大地线起点纬度B经度L，大地方位角A及大地线长度S:\n");scanf("%d%d%lf%d%d%lf%d%d%lf%lf",&B10,&B11,&B12,&L10,&L11,&L12,&A10,&A 11,&A12,&S);B1=(B10+(float)B11/60+B12/3600)*pi/180;L1=(L10+(float)L11/60+L12/3600)*pi/180;A1=(A10+(float)A11/60+A12/3600)*pi/180;W1=sqrt(1-e*e*sin(B1)*sin(B1)); //计算起点规划纬度sinu1=sin(B1)*sqrt(1-e*e)/W1; //计算起点规划纬度cosu1=cos(B1)/W1; //计算起点规划纬度sinA0=cosu1*sin(A1); //计算辅助函数值cotsigma1=cosu1*cos(A1)/sinu1; //计算辅助函数值sin2sigma1=2*cotsigma1/(cotsigma1*cotsigma1+1); //计算辅助函数值cos2sigma1=(cotsigma1*cotsigma1-1)/(cotsigma1*cotsigma1+1); //计算辅助函数值A=6356863.020+(10708.949-13.474*(1-sinA0*sinA0))*(1-sinA0*sinA0);B=(5354.469-8.798*(1-sinA0*sinA0))*(1-sinA0*sinA0);C=(2.238*(1-sinA0*sinA0))*(1-sinA0*sinA0)+0.006;afa=691.46768-(0.58143-0.00144*(1-sinA0*sinA0))*(1-sinA0*sinA0);beta=(0.2907-1.0E-3*(1-sinA0*sinA0))*(1-sinA0*sinA0);sigma0=(S-(B+C*cos2sigma1)*sin2sigma1)/A;sin2=sin2sigma1*cos(2*sigma0)+cos2sigma1*sin(2*sigma0);cos2=cos2sigma1*cos(2*sigma0)-sin2sigma1*sin(2*sigma0);sigma=sigma0+(B+5*C*cos2)*sin2/A;delta=(afa*sigma+beta*(sin2-sin2sigma1))*sinA0; //计算经度差改正数delta=delta/3600*pi/180;sinu2=sinu1*cos(sigma)+cosu1*cos(A1)*sin(sigma);B2=atan(sinu2/(sqrt(1-e*e)*sqrt(1-sinu2*sinu2)));lamda=atan(sin(A1)*sin(sigma)/(cosu1*cos(sigma)-sinu1*sin(sigma)*cos(A1))); if(sin(A1)>0){if(tan(lamda)>0)lamda=fabs(lamda);elselamda=pi-fabs(lamda);}else{if(tan(lamda)>0)lamda=fabs(lamda)-pi;elselamda=-1*fabs(lamda);}L2=L1+lamda-delta;A2=atan(cosu1*sin(A1)/(cosu1*cos(sigma)*cos(A1)-sinu1*sin(sigma)));if(sin(A1)>0){if(tan(A2)>0)A2=pi+fabs(A2);elseA2=2*pi-fabs(A2);}else{if(tan(A2)>0)A2=fabs(A2);elseA2=pi-fabs(A2);}B2=B2*180*3600/pi;L2=L2*180*3600/pi;A2=A2*180*3600/pi;B20=(int)B2/3600;B21=(int)B2/60-B20*60;B22=B2-B20*3600-B21*60;L20=(int)L2/3600;L21=(int)L2/60-L20*60;L22=L2-L20*3600-L21*60;A20=(int)A2/3600;A21=(int)A2/60-A20*60;A22=A2-A20*3600-A21*60;printf("正算得到的终点大地经度和大地纬度及A2:\n%d %d %lf\n%d %d %lf\n%d %d %lf\n",B20,B21,B22,L20,L21,L22,A20,A21,A22);}else{printf("请输入大地线起点和终点的坐标BL\n");scanf("%d%d%lf%d%d%lf%d%d%lf%d%d%lf",&B10,&B11,&B12,&L10,&L11,&L12,&B 20,&B21,&B22,&L20,&L21,&L22);B1=(B10+(double)B11/60+B12/3600)*pi/180;L1=(L10+(double)L11/60+L12/3600)*pi/180;B2=(B20+(double)B21/60+B22/3600)*pi/180;L2=(L20+(double)L21/60+L22/3600)*pi/180;W1=sqrt(1-e*e*sin(B1)*sin(B1));W2=sqrt(1-e*e*sin(B2)*sin(B2));sinu1=sin(B1)*sqrt(1-e*e)/W1;sinu2=sin(B2)*sqrt(1-e*e)/W2;cosu1=cos(B1)/W1;cosu2=cos(B2)/W2;L=L2-L1;a1=sinu1*sinu2;a2=cosu1*cosu2;b1=cosu1*sinu2;b2=sinu1*cosu2;delta0=0;lamda=L+delta0;p=cosu2*sin(lamda);q=b1-b2*cos(lamda);A1=atan(p/q);if(p>0){if(q>0)A1=fabs(A1);elseA1=pi-fabs(A1);}else{if(q>0)A1=2*pi-fabs(A1);elseA1=pi+fabs(A1);}sins=p*sin(A1)+q*cos(A1); //计算sigma的正弦值coss=a1+a2*cos(lamda); //计算sigma的余弦值sigma=atan(sins/coss);if(coss>0)sigma=fabs(sigma);elsesigma=pi-fabs(sigma);sinA0=cosu1*sin(A1);x=2*a1-(1-sinA0*sinA0)*cos(sigma);afa=(33523299-(28189-70*(1-sinA0*sinA0))*(1-sinA0*sinA0))*1.0e-10; beta=(28189-94*(1-sinA0*sinA0))*1.0e-10;delta=(afa*sigma-beta*x*sin(sigma))*sinA0;lamda=L+delta;while(fabs(delta-delta0)>4.8e-10){delta0=delta;p=cosu2*sin(lamda);q=b1-b2*cos(lamda);A1=atan(p/q);if(p>0){if(q>0)A1=fabs(A1);elseA1=pi-fabs(A1);}else{if(q>0)A1=2*pi-fabs(A1);elseA1=pi+fabs(A1);}sins=p*sin(A1)+q*cos(A1); //计算sigma的正弦值coss=a1+a2*cos(lamda); //计算sigma的余弦值sigma=atan(sins/coss);if(coss>0)sigma=fabs(sigma);elsesigma=pi-fabs(sigma);sinA0=cosu1*sin(A1);x=2*a1-(1-sinA0*sinA0)*cos(sigma);afa=(33523299-(28189-70*(1-sinA0*sinA0))*(1-sinA0*sinA0))*1.0e-10;beta=(28189-94*(1-sinA0*sinA0))*1.0e-10;delta=(afa*sigma-beta*x*sin(sigma))*sinA0;lamda=L+delta;}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)*cos(sigma);S=A*sigma+(B*x+C*y)*sin(sigma);A2=atan((cosu1*sin(lamda))/(b1*cos(lamda)-b2));if(sin(A1)>0){if(tan(A2)>0)A2=pi+fabs(A2);elseA2=2*pi-fabs(A2);}else{if(tan(A2)>0)A2=fabs(A2);elseA2=pi-fabs(A2);}A1=A1*3600*180/pi;A2=A2*3600*180/pi;A10=(int)A1/3600;A11=(int)A1/60-A10*60;A12=A1-A10*3600-A11*60;A20=(int)A2/3600;A21=(int)A2/60-A20*60;A22=A2-A20*3600-A21*60;printf("反算得到的方位角A1A2及大地线长S:\n%d %d %lf\n%d %d %lf\n%lf\n",A10,A11,A12,A20,A21,A22,S);}}四．程序执行结果：。