大地主题解算深度干货-超精

1、主题解算（正算）代码与白塞尔主题解算主题解算（正算）代侶：根据经纬度和方向角以及距离 计算另外_点坐标新建模挟-拷贝下面的主题（正算）代码，调用方法例起点经度:116.235（度）终点纬度:37.435（度）方向角：50度长度：500米终点经纬度"经度，纬度"=putation（37.435,116.235,50,500）Const Pi = 3.98Private a, b, c, alpha, e, e2, W, V As Double'a长轴半径'b短轴'c极曲率半径'alpha扁率'e第一偏心率'e2第二偏心率W 第一

2、根本纬度函数'V 第二根本纬度函数Private B1, L1, B2, L2 As Double'B1点1的纬度11点1的经度'B2点1的纬度12点2的经度Private S As Double ”小线长度Private A1, A2 As Double'A1点1到点2的方位角'A2点2到点1的方位角Function putation(STARTLAT, STARTLONG, ANGLE1,DISTANCE As Double) As StringB1 =STARTLATL1 = STARTLONGA1 =ANGL E1S = DISTANCEa=63

3、78245b = 6356752.3142c = a2/b-可修编.alpha = (a - b) / ae = Sqr(a" 2 - b 八 2) / ae2 = Sqr(a2-b2)/bB1 = rad(B1)L1 =rad(L1)A1 = rad(A1)W = Sqr(1 - e2*(Sin(B1)2)V = W* (a/b)Dim W1 As Double臼二y第一偏心率7/计算起目的I月化纬度W1 =W HSqr(1 -e1 *e1 * Sin(B1 ) * Sin(B1 )sinul =Sin(B1)*Sqr(1 - E1 * E1 )/W1cosul = Cos(B1)

4、 / W17/ it算稱助函数值sinAO = cosul * Sin(A1)cotql = cosul * Cos(A1)sin2q1 = 2 * cotql / (cotql "2 + 1)cos2q1 = (cotql " 2 - 1) / (cotql "2 + 1) 7/ 计算系数 AA,BB,CC H AAlpha, BBeta 的值 ocos2A0 = 1 - sinAO 八 2 e2 = Sqr(a2-b2)/bk2 = e2 * e2 * cos2A0Dim aa, BB, CC, EE22, AAlpha, BBeta As Double aa

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 * e2 / 128) * co

6、s2A0 * cos2A0BBeta = (e2 * e2 / 32 + e2 * e2 * e2 / 32) * cos2A0 - (e2 * e2 * e2 / 64) * cos2A0 * cos2A07/计算球面长度qO = (S - (BB + CC * cos2q1) * sin2q1)/aasin2q1q0 = sin2q1 * Cos(2 * qO) + cos2q1 * Sin(2 *qO)cos2q1q0 = cos2q1 * Cos(2 * qO) - sin2q1 * Sin(2 * qO)q = qO + (BB + 5 * CC * cos2q1 qO) * sin

7、2q1 qO / aa7/计算经度差改正数theta = (AAlpha * q + beta * (sin2q1q0 - sin2q1) * sinAO7/廿算终点坐标及方位角sinu2 = sinul * Cos(q) + cosul * Cos(A1) * Sin(q)B2 = Atn(sinu2 / (Sqr(1 - E1 *E1)* Sqr(1 - sinu2 * sinu2) * / Pilamuda = Atn(Sin(A1) * Sin(q) / (cosul * Cos(q) - sinul * Sin(q) * Cos(A1) * / PiIf (Sin(A1) >

8、0) ThenIf (Sin(A1) * Sin(q)/(cosul *Cos(q) - sinul * Sin(q) * Cos(A1) > 0) Then lamuda = Abs(lamuda)Elselamuda = - Abs(lamuda)End IfElseIf (Sin(A1) * Sin(q)/(cosu1 * Cos(q) - sinul * Sin(q) * Cos(A1) > 0) Then lamuda = Abs(lamuda)-Elselamuda = -Abs(lamuda)End IfEnd IfL2 = L1 * / Pi + lamuda -

9、theta *PiA2 = Atn(cosu1 * Sin(A1)/(cosu1 *Cos(q) * Cos(A1) - sinul * Sin(q) * / PiIf (Sin(A1) > 0) ThenIf (cosul * Sin(A1) / (cosul *Cos(q) * Cos(A1) - sinul * Sin(q) > 0) ThenA2= + Abs(A2)ElseA2 = 360 - Abs(A2)End IfElseIf (cosul * Sin(A1) / (cosul *Cos(q) * Cos(A1)-sinul * Sin(q) > 0) The

10、nA2 = Abs(A2)ElseA2= - Abs(A2)End IfEnd Ifputation = L2 & T & B2End FunctionPrivate Function rad(ByValangle_d As Double) As Double rad = angle_d * Pi /End Function白塞尔主题解算正算代碍:#include<stdio.h>#include<math.h>#define ee 0.2966#define I 3.141592653double F(double,double,double);voi

11、d main(void)double A1,B1,L1,S,A2,B2,L2;double x1,x2,x3,y1 ,y2,y3,z1 ,z2,z3;double W1 tsinu1 ,sinu2,cosu1 tsinA0;double cotai tcos2a1 ,sin2a1,cosA0A0;double ABC,d,e,a0,a1,m;doublenaQR；printfC请输人数据B1=H);scanfC%lf %lf %lf",&x1,&x2,&x3);B1=F(x1,x2,x3);printf(请输入数据L1=H);scanfC%lf %lf %lf

12、",& y1,&y2,&y3);L1=F(y1,y2,y3)；printfCiBi A A1=H);scanff%lf %lf %lfw,&z1,&z2,&z3);A1=F(z1,z2,z3);printf(请 i A S=H);scanf(lf",&S);printf(HB1=%.9fn-,B1);printfCL1=%.9fn,L1);printfCA1=%.9fn,A1); printf('S=%fn',S);厂计算起点的W1=sqrt(1-ee*sin(B1)*sin(B1); sinul =s

13、in(B1)*sqrt(1-ee)/W1; cosu1=cos(B1)/W1;printfCW1=%.9fn,W1);printf('sinu1 =%.9fn",sinul);printf('cosu1 =%.9fn",cosu1);厂计算辅BJ函数值”/sinA0=cosu1*sin(A1);cotai =cosu1 *cos(A1 )/sinu1;Sin2a1=2*cota1/(cota1*cota1+1);cos2a1 =(cota1 *cota1-1 )/(cota1 *cota1 +1);printf("sinA0=%.9fn',

14、sinA0);printf("cota1 =%.9fn",cota1);printf('sin2a1 =%.9fn",sin2a1);printf(Bcos2a1 =%.9fn",cos2a1);广计算系» ABC K de*/cosA0A0=1 -sinAO*sinAO;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.5

15、8143-0.00144*cosA0A0)*cosA0A0;e=(0.2907-cosA0A0*0.0010)*cosA0A0;printfCcosA0A0=%.9fn",cosA0A0);printfCA=%.3fn",A);printf("B=%.6fn-,B);pnntfCC=%.9fn",C);printfCd=%.7fn",d);printf("e=%.9fn",e);广计算球面长rar/aO=(S-(B+C*cos2a1 )*sin2a1 )/A;m=sin2a1 *cos(2*a0)+cos2a1 *sin(2

16、*a0);n=(cos2a1 )*(cos(2*a0)-(sin2a1)*(sin(2*a0);a=aO+(B+5*C*n)*m/A;printf("a0=%.9fn",a0);printf(Bm=%.9fn",m);printfCn=%.9fn",n);printf('a=%.9fn",a);厂计算经度差改正数*/Q=(d*a+(e*(m-sin2a1)/3600/*l)*sinA0;printfCQ=%.7fn,Q);广计算终自坐标15方位肝/sinu2=sinu1 *cos(a)+cosu1 *cos(A1)*sin(a);B2=

17、*atan (sinu2/(sqrt(1-ee)*(sqrt(1-sinu2*sinu2)/l;R=*atan (sin(A1 )*sin(a)/(cosul *cos(a)-sinu1 *sin(a)*cos(A1)/!; printf("sinu2=%.9fn',sinu2);printf(,B2=%fn,JB2);printfCR=%fn-,R7l);厂确定R的值*/if(sin(A1)> 0&& tan(R)>0)R=abs(R);else if(sin(A1)>0 && tan(R)<0)R=l-abs(R);

18、else if(sin(A1)<0 && tan(R)<0)R=-abs(R);elseR=abs(R)-l;/*确定L2A2的值*/L2=(L1 */l+R-(Q/206265*/l);A2=atan(cosu1 *sin(A1 )/(cosu1 *cos(a)*cos(A1 )-sinu1 *sin(a); if(sin(A1)<0&&tan(A2)>0)A2=(fabs(A2)*/l;else if(sin(A1 )<0&&tan(A2)<0)A2=(l-fabs(A2)*/l;else if(sin(A

19、1 )>0&&tan(A2)>0)A2=(l+fabs(A2)*/l;elseA2=(2*l-fabs(A2)*/l;printf(BA2=%3fn B2=%3fn L2=%3fn',A2,B2丄2);double F(double a2,double b2,double c2)double d2;d2=(double)(a2+1.052/60+1.0*c2/3600);d2=(d2/)*l;return (d2);运仃结果:E: 白 Kexle>:l.lDebugmyfl.exe-S3Al=3.929899698W1=0.999137531 sinu

20、l=0.506273572 cosul-0.0G237293B sinA0=-06丄1560970 cotal=-l.200958132 sin2al=-0.983464700 coe2al=0.181099929 uusA0A0-0G252夕3180A=6363454.5858=3348.764701 C=R.89R052992 d-G91.1012731 «-0.181584350 a0=l.S71991H34 n=0.983029170 n=-0.182449314 a-1.572508226 Q=-664.6243766 sinu2=-0.6W8879161 B2=-3?.6

21、01471 R=-3631.26062702-50.438693B2=-37.601471L2=51.684618Press any key to continue反算代码：#include<stdio.h>#include<math.h>#define ee 0.2966#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,q;double x1 ,x2,x3

22、,y1 ,y2y3,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; printfCffi#人起始数据B1 = '); scanf(H%lf %lf %lf,&x1，&x2,&x3);B1=F(x1fx2,x3);printfCig输人起始数据L1=*); scanf(n%lf %lf %lf&y1，& y2,&y3);L1=F(y1,y2,y3);printfCig# 人起始釵

23、 B2="); scanf(,%lf %lf %ir,&zh&z2,&z3);B2=F(z1 ,z2,z3); printfC请输人起始数据L2J);scanf(H%lf %lf %lf：&w1,&w2,&w3); L2=F(w1 ,w2tw3);printf(® 二 fn；B1);printf(,L1=%.9fM,L1);printf(,B2=%.9fnB2);printf(,L2=%.9f,L2);/*»» itw*/W1=sqrt(1-ee*sin(B1)*sin(B1);W2=sqrt(1-ee*s

24、in(B2)*sin(B2); sinu1=sin(B1)*sqrt(1-ee)/W1;sinu2=sin(B2)*(sqrt(1 -ee)/W2; cosu1=cos(B1)AA/1;cosu2=cos(B2)/W2;L=L2-L1;a1=sinursinu2; a2=cosu1*cosu2; b1=cosursinu2;b2=sinu1 *cosu2;pnntf(,W1=%.9,W1);printf(,W2=%.9nn,lW2); printf(,sinu1=%.W,sinu1);printf("sinu2=%.9fn,1sinu2);printf(,cosu1 =%.9fnco

25、su2);printfCL=%.9fnM1L);printfCa1=%.9fnM,a1);printf(Ma2=%.9fnM,a2);printtrb1=%9fn；b1);printf(,b2=%.9nn,b2);厂用逐次tlifi法同时廿算起点方位角、球面长度及经差R"/ doR=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=kfabs(A1);else if(x<0&&y<0)A1=l+fabs(A