大地主题解算深度干货超精.docx
- 文档编号:17156085
- 上传时间:2023-07-22
- 格式:DOCX
- 页数:15
- 大小:101.25KB
大地主题解算深度干货超精.docx
《大地主题解算深度干货超精.docx》由会员分享,可在线阅读,更多相关《大地主题解算深度干货超精.docx(15页珍藏版)》请在冰点文库上搜索。
大地主题解算深度干货超精
大地主题解算(正算)代码与白塞尔大地主题解算
大地主题解算(正算)代码:
根据经纬度和方向角以及距离计算另外一点坐标
新建模块->拷贝下面的大地主题(正算)代码,调用方法示例:
起点经度:
116.235(度)
终点纬度:
37.435(度)
方向角:
50(度)
长度:
500(米)
终点经纬度("经度,纬度")=Computation(37.435,116.235,50,500)
ConstPi=3.1415926535898
Privatea,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)AsString
B1=STARTLAT
L1=STARTLONG
A1=ANGLE1
S=DISTANCE
a=6378245
b=6356752.3142
c=a^2/b
alpha=(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)
DimW1AsDouble
E1=e''''第一偏心率
'//计算起点的归化纬度
W1=W''Sqr(1-e1*e1*Sin(B1)*Sin(B1))
sinu1=Sin(B1)*Sqr(1-E1*E1)/W1
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
Dimaa,BB,CC,EE22,AAlpha,BBetaAsDouble
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/512)
e2=E1*E1
AAlpha=(e2/2+e2*e2/8+e2*e2*e2/16)-(e2*e2/16+e2*e2*e2/16)*cos2A0+(3*e2*e2*e2/128)*cos2A0*cos2A0
BBeta=(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*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)))*/Pi
lamuda=Atn(Sin(A1)*Sin(q)/(cosu1*Cos(q)-sinu1*Sin(q)*Cos(A1)))*/Pi
If(Sin(A1)>0)Then
If(Sin(A1)*Sin(q)/(cosu1*Cos(q)-sinu1*Sin(q)*Cos(A1))>0)Then
lamuda=Abs(lamuda)
Else
lamuda=-Abs(lamuda)
EndIf
Else
If(Sin(A1)*Sin(q)/(cosu1*Cos(q)-sinu1*Sin(q)*Cos(A1))>0)Then
lamuda=Abs(lamuda)-
Else
lamuda=-Abs(lamuda)
EndIf
EndIf
L2=L1*/Pi+lamuda-theta*/Pi
A2=Atn(cosu1*Sin(A1)/(cosu1*Cos(q)*Cos(A1)-sinu1*Sin(q)))*/Pi
If(Sin(A1)>0)Then
If(cosu1*Sin(A1)/(cosu1*Cos(q)*Cos(A1)-sinu1*Sin(q))>0)Then
A2=+Abs(A2)
Else
A2=360-Abs(A2)
EndIf
Else
If(cosu1*Sin(A1)/(cosu1*Cos(q)*Cos(A1)-sinu1*Sin(q))>0)Then
A2=Abs(A2)
Else
A2=-Abs(A2)
EndIf
EndIf
Computation=L2&","&B2
EndFunction
PrivateFunctionrad(ByValangle_dAsDouble)AsDouble
rad=angle_d*Pi/
EndFunction
白塞尔大地主题解算
一、正算
代码:
#include
#include
#defineee0.006693421622966
#defineI3.141592653
doubleF(double,double,double);
voidmain(void)
{
doubleA1,B1,L1,S,A2,B2,L2;
doublex1,x2,x3,y1,y2,y3,z1,z2,z3;
doubleW1,sinu1,sinu2,cosu1,sinA0;
doublecota1,cos2a1,sin2a1,cosA0A0;
doubleA,B,C,d,e,a0,a1,m;
doublen,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=%.9f\n",B1);
printf("L1=%.9f\n",L1);
printf("A1=%.9f\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=%.9f\n",W1);
printf("sinu1=%.9f\n",sinu1);
printf("cosu1=%.9f\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=%.9f\n",sinA0);
printf("cota1=%.9f\n",cota1);
printf("sin2a1=%.9f\n",sin2a1);
printf("cos2a1=%.9f\n",cos2a1);
/*计算系数ABC及de*/
cosA0A0=1-sinA0*sinA0;
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.58143-0.00144*cosA0A0)*cosA0A0;
e=(0.2907-cosA0A0*0.0010)*cosA0A0;
printf("cosA0A0=%.9f\n",cosA0A0);
printf("A=%.3f\n",A);
printf("B=%.6f\n",B);
printf("C=%.9f\n",C);
printf("d=%.7f\n",d);
printf("e=%.9f\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=%.9f\n",a0);
printf("m=%.9f\n",m);
printf("n=%.9f\n",n);
printf("a=%.9f\n",a);
/*计算经度差改正数*/
Q=(d*a+(e*(m-sin2a1))/3600/*I)*sinA0;
printf("Q=%.7f\n",Q);
/*计算终点大地坐标及大地方位角*/
sinu2=sinu1*cos(a)+cosu1*cos(A1)*sin(a);
B2=*atan(sinu2/((sqrt(1-ee))*(sqrt(1-sinu2*sinu2))))/I;
R=*atan(sin(A1)*sin(a)/(cosu1*cos(a)-sinu1*sin(a)*cos(A1)))/I;
printf("sinu2=%.9f\n",sinu2);
printf("B2=%f\n",B2);
printf("R=%f\n",R*/I);
/*确定R的值*/
if(sin(A1)>0&&tan(R)>0)
R=abs(R);
elseif(sin(A1)>0&&tan(R)<0)
R=I-abs(R);
elseif(sin(A1)<0&&tan(R)<0)
R=-abs(R);
else
R=abs(R)-I;
/*确定L2A2的值*/
L2=(L1*/I+R-(Q/206265*/I));
A2=atan(cosu1*sin(A1)/(cosu1*cos(a)*cos(A1)-sinu1*sin(a)));
if(sin(A1)<0&&tan(A2)>0)
A2=(fabs(A2))*/I;
elseif(sin(A1)<0&&tan(A2)<0)
A2=(I-fabs(A2))*/I;
elseif(sin(A1)>0&&tan(A2)>0)
A2=(I+fabs(A2))*/I;
else
A2=(2*I-fabs(A2))*/I;
printf("A2=%3f\nB2=%3f\nL2=%3f\n",A2,B2,L2);
}
doubleF(doublea2,doubleb2,doublec2)
{
doubled2;
d2=(double)(a2+1.0*b2/60+1.0*c2/3600);
d2=(d2/)*I;
return(d2);
}
运行结果:
二、反算
代码:
#include
#include
#defineee0.006693421622966
#defineI3.141592653
doubleF(double,double,double);
voidmain(void)
{
doubleB1,B2,L1,L2,A1,A2,S,Y;
doubleW1,W2,L,Q,R,A,B,C;
doublex,y,z,p,q;
doublex1,x2,x3,y1,y2,y3,z1,z2,z3,w1,w2,w3;
doublea1,a2,b1,b2,m,n;
doublesinp,cosp,sinu1,sinu2,cosu1,cosu2,sinA0,cosA0;
Q=0;q=0;
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("请输入起始数B2=");
scanf("%lf%lf%lf",&z1,&z2,&z3);
B2=F(z1,z2,z3);
printf("请输入起始数据L2=");
scanf("%lf%lf%lf",&w1,&w2,&w3);
L2=F(w1,w2,w3);
printf("B1=%f\n",B1);
printf("L1=%.9f\n",L1);
printf("B2=%.9f\n",B2);
printf("L2=%.9f\n",L2);
/*辅助计算*/
W1=sqrt(1-ee*sin(B1)*sin(B1));
W2=sqrt(1-ee*sin(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=%.9f\n",W1);
printf("W2=%.9f\n",W2);
printf("sinu1=%.9f\n",sinu1);
printf("sinu2=%.9f\n",sinu2);
printf("cosu1=%.9f\n",cosu2);
printf("L=%.9f\n",L);
printf("a1=%.9f\n",a1);
printf("a2=%.9f\n",a2);
printf("b1=%.9f\n",b1);
printf("b2=%.9f\n",b2);
/*用逐次趋近法同时计算起点大地方位角、球面长度及经差R*/
do
{
R=L+Q;
x=cosu2*sin(R);
y=b1-b2*cos(R);
A1=atan(x/y);
if(x>0&&y>0)
A1=fabs(A1);
elseif(x>0&&y<0)
A1=I-fabs(A1);
elseif(x<0&&y<0)
A1=I+fabs(A1);
else
A1=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*cosA0)*(1e-10);
Q=q;
q=(m*p-m*z*sin(p))*sinA0;
}while(fabs(206265*(q-Q))>(1e-4));
printf("Q=%f\n",Q);
/*计算系数ABC及大地线长度*/
A=6356863.020+(10708.949-13.474*cosA0*cosA0)*cosA0*cosA0;
B=10708.938-17.956*cosA0*cosA0;
C=4.487;
Y=(cosA0*cosA0*cosA0*cosA0-2*z*z)*cos(p);
S=A*p+(B*z+C*Y)*sinp;
/*计算反方位角*/
A2=atan((cosu1*sin(R))/(b1*cos(R)-b2));
A2=A2*/I;
A1=A1*/I;
if(A1>0)
A2=abs(A2);
else
A2=-abs(A2);
printf("A1=%3f\nA2=%3f\nS=%3f\n",A1,A2,S);
}
doubleF(doublea,doubleb,doublec)
{
doubled;
if(a<0)
d=(double)(a-1.0*b/60-1.0*c/3600);
else
d=(double)(a+1.0*b/60+1.0*c/3600);
d=(d/)*I;
returnd;}
运行结果:
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 地主 题解 深度 干货