大地主题解算(正算)代码

大地主题解算(正算)代码：根据经纬度和方向角以及距离计算另外一点坐标新建模块-拷贝下面的大地主题(正算)代码，调用方法示例：起点经度：116.235(度)终点纬度：37.435(度)方向角：50（度）长度：500（米）终点经纬度（经度,纬度）=Computation(37.435,116.235,50,500)ConstPi=3.1415926535898Privatea,b,c,alpha,e,e2,W,VAsDouble'a长轴半径'b短轴'c极曲率半径'alpha扁率'e第一偏心率'e2第二偏心率'W第一基本纬度函数'V第二基本纬度函数PrivateB1,L1,B2,L2AsDouble'B1点1的纬度'L1点1的经度'B2点1的纬度'L2点2的经度PrivateSAsDouble'''''大地线长度PrivateA1,A2AsDouble'A1点1到点2的方位角'A2点2到点1的方位角FunctionComputation(STARTLAT,STARTLONG,ANGLE1,DISTANCEAsDouble)AsStringB1=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)DimW1AsDoubleE1=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*cos2A0Dimaa,BB,CC,EE22,AAlpha,BBetaAsDoubleaa=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)))*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-Abs(lamuda)EndIfElseIf(Sin(A1)*Sin(q)/(cosu1*Cos(q)-sinu1*Sin(q)*Cos(A1))0)Thenlamuda=Abs(lamuda)-180Elselamuda=-Abs(lamuda)EndIfEndIfL2=L1*180/Pi+lamuda-theta*180/PiA2=Atn(cosu1*Sin(A1)/(cosu1*Cos(q)*Cos(A1)-sinu1*Sin(q)))*180/PiIf(Sin(A1)0)ThenIf(cosu1*Sin(A1)/(cosu1*Cos(q)*Cos(A1)-sinu1*Sin(q))0)ThenA2=180+Abs(A2)ElseA2=360-Abs(A2)EndIfElseIf(cosu1*Sin(A1)/(cosu1*Cos(q)*Cos(A1)-sinu1*Sin(q))0)ThenA2=Abs(A2)ElseA2=180-Abs(A2)EndIfEndIfComputation=L2&,&B2EndFunctionPrivateFunctionrad(ByValangle_dAsDouble)AsDoublerad=angle_d*Pi/180EndFunction