大地主題解算(深度干貨-超精)

1、大地主題解算（正算）代碼與白塞爾大地主題解算大地主題解算（正算）代碼：根據(jù)經(jīng)緯度和方向角以及距離計算另外一點坐標新建模塊-拷貝下面的大地主題（正算）代碼，調用方法示例：起點經(jīng)度：116.235（度 ）終點緯度：37.435（度 ）方向角：50（度）長度：500（米）終點經(jīng)緯度（"經(jīng)度,緯度 "）=Computation（37.435,116.235,50,500）Const Pi = 3.1415926535898Private a, b, c, alpha, e, e2, W, V As Double'a長軸半徑'b短軸'c極曲率半徑'al

2、pha 扁率'e第一偏心率'e2第二偏心率'W 第一基本緯度函數(shù)'V 第二基本緯度函數(shù)Private B1, L1, B2, L2 As Double'B1點1的緯度'L1點1的經(jīng)度'B2點1的緯度'L2點2的經(jīng)度Private S As Double '''''大地線長度Private A1, A2 As Double'A1點 1 到點2 的方位角'A2點 2 到點1 的方位角Function Computation(STARTLAT, STARTLONG, ANGLE1,

3、 DISTANCE As Double) As StringB1 = STARTLATL1 = STARTLONGA1 = ANGLE1S = DISTANCEa = 6378245b = 6356752.3142alpha = (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 '''

4、'第一偏心率'/ 計算起點的歸化緯度W1 = W ''Sqr(1 - e1 * e1 * Sin(B1 ) * Sin(B1 ) sinu1 = Sin(B1) * Sqr(1 - E1 * E1) / W1 cosu1 = Cos(B1) / W1 '/ 計算輔助函數(shù)值sinA0 = cosu1 * Sin(A1)cotq1 = cosu1 * Cos(A1)sin2q1 = 2 * cotql / (cotql 八 2 + 1)cos2q1 = (cotql 八 2 - 1) / (cotql 八 2 + 1)'/ 計算系數(shù) AA,BB,CC

5、M AAlpha, BBeta的值。cos2A0 = 1 - sinA0 八 2e2 = Sqr(a 八 2 - b 八 2) / bk2 = e2 * e2 * cos2A0Dim aa, BB, CC, EE22, AAlpha, BBeta As Double aa = 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

6、 / 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) * sin2

7、q1) / 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'/ 計算經(jīng)度差改正數(shù)theta = (AAlpha * q + beta * (sin2q1q0 - sin2q1) * sinA0'/ 計算終點大地坐標及大地方位角sinu2 = sinu1 * Cos(q) + cosu1 * Cos(A1) * S

8、in(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) * Cos(A1) * 180 / PiIf (Sin(A1) > 0) ThenIf (Sin(A1) * Sin(q) / (cosu1 *Cos(q) - sinu1 * Sin(q) * Cos(A1) > 0) Thenlamuda = Abs(lamuda)Elselamuda = 180 -

9、Abs(lamuda)End IfElseIf (Sin(A1) * Sin(q) / (cosu1 *Cos(q) - sinu1 * Sin(q) * Cos(A1) > 0) Thenlamuda = Abs(lamuda) - 180Elselamuda = -Abs(lamuda)End IfEnd IfL2 = L1 * 180 / Pi + lamuda - theta180 / PiA2 = Atn(cosu1 * Sin(A1) / (cosu1 *Cos(q) * Cos(A1) - sinu1 * Sin(q) * 180 / PiIf (Sin(A1) >

10、0) ThenIf (cosu1 * Sin(A1) / (cosu1 *Cos(q) * Cos(A1) - sinu1 * Sin(q) > 0) ThenA2 = 180 + Abs(A2)ElseA2 = 360 - Abs(A2)End IfElseIf (cosu1 * Sin(A1) / (cosu1 *Cos(q) * Cos(A1) - sinu1 * Sin(q) > 0) ThenA2 = Abs(A2)ElseA2 = 180 - Abs(A2)End IfEnd IfComputation = L2 & "," & B2

11、End FunctionPrivate Function rad(ByValangle_d As Double) As Doublerad = angle_d * Pi / 180End 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,y

12、1,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(" 請輸入數(shù)據(jù)B1= ");scanf("%lf %lf %lf",&x1,&x2,&x3);B1=F(x1,x2,x3);printf(" 請輸入數(shù)據(jù)L1= ");scanf("%lf %lf %lf",&y1,&am

13、p;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=

14、%fn",S);/* 計算起點的規(guī)劃緯度*/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);/* 計算輔助函數(shù)值*/sinA0=cosu1*sin(A1); cota1=cosu1*cos(A1)/sinu1; sin2a1=2*cota1/(cota1*cota1

15、+1);cos2a1=(cota1*cota1-1)/(cota1*cota1+1);printf("sinA0=%.9fn",sinA0);printf("cota1=%.9fn",cota1);printf("sin2a1=%.9fn",sin2a1);printf("cos2a1=%.9fn",cos2a1);/*計算系數(shù)ABC及de*/cosA0A0=1-sinA0*sinA0;A=6356755.288+(10710.341-(13.534*cosA0A0)*cosA0A0;B=(5355.171-9.0

16、23*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=%.9fn",cosA0A0);printf("A=%.3fn",A);printf("B=%.6fn",B);printf("C=%.9fn",C);printf("d=%.7fn",d);print

17、f("e=%.9fn",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=%.9fn",a0);printf("m=%.9fn",m);printf("n=%.9fn",n);printf("a=%.9fn",a);/* 計算經(jīng)度差改正數(shù)*/Q=(

18、d*a+(e*(m-sin2a1)/3600/180*I)*sinA0;printf("Q=%.7fn",Q);/* 計算終點大地坐標及大地方位角*/sinu2=sinu1*cos(a)+cosu1*cos(A1)*sin(a);B2=180*atan(sinu2/(sqrt(1-ee)*(sqrt(1-sinu2*sinu2)/I;R=180*atan(sin(A1)*sin(a)/(cosu1*cos(a)-sinu1*sin(a)*cos(A1)/I; printf("sinu2=%.9fn",sinu2);printf("B2=%fn&

19、quot;,B2);printf("R=%fn",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*co

20、s(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=%3fn B2=%3fn L2=%3fn",A2,B2,L2);double F(

21、double a2,double b2,double c2)double d2;d2=(double)(a2+1.0*b2/60+1.0*c2/3600);d2=(d2/180)*I;return (d2);運行結果：'Eg 垂爾d'uxl.nXDBbu 動 rryfLc .|ftl=3.929B999fl 生二i時蒯爆狗一 M口削I削tUl0.9991375310 1m*1=0.5062735 72 cosul-0.32 37S930 3inA0- e.tllC097B cotl- 1.20S9SQ132 iri2a±- 0,9fi34647BB cusadi-B.

22、iSlft5MZ? posA0A0=0,62S3£8B H=6363454_5fi5B=J 348.76-17 Hl E-twHusayy；! d=691-1042731 6=0-18158-1350 a0=l.571991G34 n=0-983029170 n=-0.183449314 a=l-572508226Q=-664.62437Sfi=-37.6(11471R=TRmi .2(>flfi27 ft2=50.43S693B2-37-GG1471L2i .GSGISPrcas 6rl y kccont: in lac反算代碼：#include<stdio.h>

23、#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,

24、cosA0;Q=0;q=0;printf("請輸入起始數(shù)據(jù) B1=");scanf("%lf %lf %lf",&x1,&x2,&x3);B1=F(x1,x2,x3);printf("請輸入起始數(shù)據(jù) L1=");scanf("%lf %lf %lf",&y1,&y2,&y3);L1=F(y1,y2,y3);printf("請輸入起始數(shù) B2=");scanf("%lf %lf %lf",&z1,&z2,&

25、z3);B2=F(z1,z2,z3);printf(" 請輸入起始數(shù)據(jù)L2= ");scanf("%lf %lf %lf",&w1,&w2,&w3);L2=F(w1,w2,w3);printf("B1=%fn",B1);printf("L1=%.9fn",L1);printf("B2=%.9fn",B2);printf("L2=%.9fn",L2);/* 輔助計算*/W1=sqrt(1-ee*sin(B1)*sin(B1);W2=sqrt(1-ee*s

26、in(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=%.9fn",W1);printf("W2=%.9fn",W2);printf("sinu1=%.9fn",sinu1);printf("sinu2=%.9fn

27、",sinu2);printf("cosu1=%.9fn",cosu2);printf("L=%.9fn",L);printf("a1=%.9fn",a1);printf("a2=%.9fn",a2);printf("b1=%.9fn",b1);printf("b2=%.9fn",b2);/* 用逐次趨近法同時計算起點大地方位角、球面長度及經(jīng)差R*/doR=L+Q;x=cosu2*sin(R);y=b1-b2*cos(R);A1=atan(x/y);if(x>

28、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*c