《大地主题解算(深度干货-超精)》由会员分享，可在线阅读，更多相关《大地主题解算(深度干货-超精)（20页珍藏版）》请在金锄头文库上搜索。
1、大地主题解算(正算)代码与白塞尔大地主题解算大地主题解算(正算)代码：根据经纬度和方向角以及距离计算另外一点坐标新建模块-拷贝下面的大地主题(正算)代码，调用方法示例：起点经度：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 Doublea 长轴半径b 短轴c 极曲率半径alpha 扁率e 第一偏心率e2 第二偏心率W 第一基本纬度函数
2、V 第二基本纬度函数Private B1, L1, B2, L2 As DoubleB1 点1的纬度L1 点1的经度B2 点1的纬度L2 点2的经度Private S As Double 大地线长度Private A1, A2 As DoubleA1 点1到点2的方位角A2 点2到点1的方位角Function Computation(STARTLAT, STARTLONG, ANGLE1, DISTANCE As Double) As String B1 = STARTLAT L1 = STARTLONG A1 = ANGLE1 S = DISTANCE a = 6378245 b = 635
3、6752.3142 c = a 2 / balpha = (a - b) / a e = Sqr(a 2 - b 2) / a e2 = Sqr(a 2 - b 2) / b B1 = rad(B1) L1 = rad(L1) A1 = rad(A1) W = Sqr(1 - e 2 * (Sin(B1) 2) V = W * (a / b) Dim W1 As Double E1 = e 第一偏心率 / 计算起点的归化纬度 W1 = W Sqr(1 - e1 * e1 * Sin(B1 ) * Sin(B1 ) sinu1 = Sin(B1) * Sqr(1 - E1 * E1) / W1
4、cosu1 = 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 2 e2 = Sqr(a 2 - b 2) / b k2 = e2 * e2 * cos2A0 Dim aa, BB, CC, EE22, AAlpha, BBeta As Doubleaa
5、= 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 *
6、e2 / 128) * cos2A0 * cos2A0BBeta = (e2 * e2 / 32 + e2 * e2 * e2 / 32) * cos2A0 - (e2 * e2 * e2 / 64) * cos2A0 * cos2A0 / 计算球面长度 q0 = (S - (BB + CC * cos2q1) * sin2q1) / aa sin2q1q0 = sin2q1 * Cos(2 * q0) + cos2q1 * Sin(2 * q0) cos2q1q0 = cos2q1 * Cos(2 * q0) - sin2q1 * Sin(2 * q0) q = q0 + (BB + 5 *
7、 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) * 180 / Pilamuda = Atn(Sin(A1) * Sin(q) / (cosu1 * Cos(q) - sinu1 * Sin(q)
8、 * Cos(A1) * 180 / Pi If (Sin(A1) 0) Then If (Sin(A1) * Sin(q) / (cosu1 * Cos(q) - sinu1 * Sin(q) * Cos(A1) 0) Thenlamuda = Abs(lamuda) Elselamuda = 180 - Abs(lamuda) End If Else If (Sin(A1) * Sin(q) / (cosu1 * Cos(q) - sinu1 * Sin(q) * Cos(A1) 0) Thenlamuda = Abs(lamuda) - 180 Elselamuda = -Abs(lam
9、uda) End If End If L2 = L1 * 180 / Pi + lamuda - theta * 180 / Pi A2 = Atn(cosu1 * Sin(A1) / (cosu1 * Cos(q) * Cos(A1) - sinu1 * Sin(q) * 180 / Pi If (Sin(A1) 0) Then If (cosu1 * Sin(A1) / (cosu1 * Cos(q) * Cos(A1) - sinu1 * Sin(q) 0) Then A2 = 180 + Abs(A2) Else A2 = 360 - Abs(A2) End If Else If (c
10、osu1 * Sin(A1) / (cosu1 * Cos(q) * Cos(A1) - sinu1 * Sin(q) 0) Then A2 = Abs(A2) Else A2 = 180 - Abs(A2) End If End If Computation = L2 & , & B2End FunctionPrivate Function rad(ByValangle_d As Double) As Doublerad = angle_d * Pi / 180End Function白塞尔大地主题解算一、 正算代码：#include#include#define ee 0.00669342
11、1622966#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,
12、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=%.9fn,B1);printf(L1=%.9fn,L1);printf(A1=%.9fn,A1);printf(S=%fn,S);/*计算起点的规划纬度*/W1=sqrt(1-ee*sin(B1)*sin(B1); sinu1=sin(B1)*sqrt(1-ee)/W1;cosu1=cos(B1)/W1;printf(W1=%.9fn,W1);printf(sinu1=%.9fn,sinu1);printf(cosu1=%.9fn,cosu1);/*计算辅助函数值*/sinA0=cosu1*sin(A1); cota1=cosu1*cos(A1)/sinu1;sin2a1=2*cota1/(cota1*cota1+1)
