大地坐标与空间直角坐标的转换程序代码Word格式文档下载.docx
《大地坐标与空间直角坐标的转换程序代码Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《大地坐标与空间直角坐标的转换程序代码Word格式文档下载.docx（14页珍藏版）》请在冰豆网上搜索。
,&
m）;
if（m==0）exit（0）;
sp2:
请选择椭球参数（输入椭球序号）:
1.克拉索夫斯基椭球参数\n"
椭球参数\n"
0.其他椭球参数（自行输入）\n"
n）;
switch（n）
{
case1:
a=;
b=;
c=;
e2=;
ep2=;
break;
case2:
case3:
case0:
printf（"
请输入椭球参数：
printf（"
长半径a="
scanf（"
%lf"
a）;
短半径b="
b）;
c=a*a/b;
ep2=（a*a-b*b）/（b*b）;
e2=（a*a-b*b）/（a*a）;
break;
}
default:
\n\n输入错误!
\n请重新输入!
\n\n"
gotosp2;
}
while
（1）
{
switch（m）
BLH_XYZ（）;
XYZ_BLH（）;
B_ZS（）;
case4:
B_FS（）;
case5:
GUS_ZS（）;
case6:
GUS_FS（）;
gotosp1;
是否继续进行此功能计算\n\n"
（若继续进行此功能计算，则输入1；
\n若选择其他功能进行计算，则输入2；
\n若退出，则输入0.）\n"
t）;
switch（t）
gotosp1;
exit（0）;
}
}
doubleRAD（doubled,doublef,doublem）
doublee;
doublesign=（d<
:
;
if（d==0）
sign=（f<
if（f==0）
sign=（m<
if（d<
0）
d=d*;
if（f<
f=f*;
if（m<
m=m*;
e=sign*（d*3600+f*60+m）*PI/（3600*180）;
returne;
voidRBD（doublehd）
intt;
intd,f;
doublem;
doublesign=（hd<
if（hd<
hd=fabs（hd）;
hd=hd*3600*180/PI;
t=int（hd/3600）;
d=sign*t;
hd=hd-t*3600;
f=int（hd/60）;
m=hd-f*60;
%d'
%lf'
d,f,m）;
voidBLH_XYZ（）
doubleB,L,H,N,W;
doubled,f,m;
doubleX,Y,Z;
请输入大地坐标（输入格式为角度（例如：
30'
40'
50'
））：
大地经度L="
scanf（"
d,&
f,&
L=RAD（d,f,m）;
大地纬度B="
B=RAD（d,f,m）;
大地高H="
H）;
W=sqrt（1-e2*sin（B）*sin（B））;
N=a/W;
X=（N+H）*cos（B）*cos（L）;
Y=（N+H）*cos（B）*sin（L）;
Z=（N*（1-e2）+H）*sin（B）;
\n\n转换后得到大地空间直角坐标为：
X=%lf\nY=%lf\nZ=%lf\n\n"
X,Y,Z）;
voidXYZ_BLH（）
doubletgB0,tgB1;
请输入大地空间直角坐标：
X="
X）;
Y="
Y）;
Z="
Z）;
\n\n转换后得到大地坐标为：
L=atan（Y/X）;
大地经度为：
L="
RBD（L）;
tgB0=Z/sqrt（X*X+Y*Y）;
tgB1=（1/sqrt（X*X+Y*Y））*（Z+a*e2*tgB0/sqrt（1+tgB0*tgB0-e2*tgB0*tgB0））;
while（fabs（tgB0-tgB1）>
5*pow（10,-10））
tgB0=tgB1;
tgB1=（1/sqrt（X*X+Y*Y））*（Z+a*e2*tgB0/sqrt（1+tgB0*tgB0-e2*tgB0*tgB0））;
B=atan（tgB1）;
大地纬度为：
B="
RBD（B）;
H=sqrt（X*X+Y*Y）/cos（B）-N;
大地高为：
H=%lf\n\n"
H）;
voidB_ZS（）
doubleL1,B1,A1,s,d,f,mi;
doubleu1,u2,m,M,k2,alfa,bt,r,kp2,alfap,btp,rp;
doublesgm0,sgm1,lmd,lmd1,lmd2,A2,B2,l,L2;
请输入已知点的大地坐标（输入格式为角度（例如：
），下同）：
\nL1="
mi）;
L1=RAD（d,f,mi）;
\nB1="
B1=RAD（d,f,mi）;
请输入大地方位角：
\nA1="
A1=RAD（d,f,mi）;
请输入该点至另一点的大地线长：
\ns="
s）;
u1=atan（sqrt（1-e2）*tan（B1））;
m=asin（cos（u1）*sin（A1））;
M=atan（tan（u1）/cos（A1））;
m=（m>
0）m:
m+2*PI;
M=（M>
0）M:
M+PI;
k2=ep2*cos（m）*cos（m）;
alfa=（1-k2/4+7*k2*k2/64-15*k2*k2*k2/256）/b;
bt=k2/4-k2*k2/8+37*k2*k2*k2/512;
r=k2*k2/128-k2*k2*k2/128;
sgm0=alfa*s;
sgm1=alfa*s+bt*sin（sgm0）*cos（2*M+sgm0）+r*sin（2*sgm0）*cos（4*M+2*sgm0）;
while（fabs（sgm0-sgm1）>
*PI/180*pow（10,-7））
sgm0=sgm1;
sgm1=alfa*s+bt*sin（sgm0）*cos（2*M+sgm0）+r*sin（2*sgm0）*cos（4*M+2*sgm0）;
sgm0=sgm1;
A2=atan（tan（m）/cos（M+sgm0））;
A2=（A2>
0）A2:
A2+PI;
A2=（A1>
PI）A2:
u2=atan（-cos（A2）*tan（M+sgm0））;
lmd1=atan（sin（u1）*tan（A1））;
lmd1=（lmd1>
0）lmd1:
lmd1+PI;
lmd1=（m<
PI）lmd1:
lmd2=atan（sin（u2）*tan（A2））;
lmd2=（lmd2>
0）lmd2:
lmd2+PI;
lmd2=（m<
PI）（（（M+sgm0）<
PI）lmd2:
lmd2+PI）:
（（（M+sgm0）>
lmd2+PI）;
lmd=lmd2-lmd1;
B2=atan（sqrt（1+ep2）*tan（u2））;
kp2=e2*cos（m）*cos（m）;
alfap=（e2/2+e2*e2/8+e2*e2*e2/16）-e2/16*（1+e2）*kp2+3*e2*kp2*kp2/128;
btp=e2*（1+e2）*kp2/16-e2*kp2*kp2/32;
rp=e2*kp2*kp2/256;
l=lmd-sin（m）*（alfap*sgm0+btp*sin（sgm0）*cos（2*M+sgm0）+rp*sin（2*sgm0）*cos（4*M+2*sgm0））;
L2=L1+l;
\n\n得到另一点的大地坐标和大地线在该点的大地方位角为：
L2="
RBD（L2）;
B2="
RBD（B2）;
A2="
RBD（A2）;
voidB_FS（）
doubleL1,B1,L2,B2,s,A1,A2,du,f,mi,m0,m,M;
doublel,u1,u2,alfa,bt,r,lmd0,dit_lmd,lmd,sgm,dit_sgm,sgm0,sgm1,alfap,btp,rp,k2,kp2;
请输入第一个点大地坐标（输入格式为角度（例如：
\n大地经度L1="
du,&
L1=RAD（du,f,mi）;
大地纬度B1="
B1=RAD（du,f,mi）;
\n请输入第二个点大地坐标：
\n大地经度：
L2=RAD（du,f,mi）;
大地纬度：
B2=RAD（du,f,mi）;
l=L2-L1;
u2=atan（sqrt（1-e2）*tan（B2））;
sgm0=acos（sin（u1）*sin（u2）+cos（u1）*cos（u2）*cos（l））;
m0=asin（cos（u1）*cos（u2）*sin（l）/sin（sgm0））;
dit_lmd=*sgm0*sin（m0）;
lmd0=l+dit_lmd;
dit_sgm=sin（m0）*dit_lmd;
sgm1=sgm0+dit_sgm;
m=asin（cos（u1）*cos（u2）*sin（lmd0）/sin（sgm1））;
A1=atan（sin（lmd0）/（cos（u1）*tan（u2）-sin（u1）*cos（lmd0）））;
A1=（A1>
0）A1:
A1+PI;
A1=（m>
M=atan（sin（u1）*tan（A1）/sin（m））;
sgm1=acos（sin（u1）*sin（u2）+cos（u1）*cos（u2）*cos（l+sin（m）*（alfap*sgm0+btp*sin（sgm0）*cos（2*M+sgm0））））;
1*PI/180*pow（10,-8））
sgm1=acos（sin（u1）*sin（u2）+cos（u1）*cos（u2）*cos（l+sin（m）*（alfap*sgm0+btp*sin（sgm0）*cos（2*M+sgm0））））;
sgm=sgm1;
lmd=l+sin（m）*（alfap*sgm+btp*sin（sgm）*cos（2*M+sgm））;
s=（sgm-bt*sin（sgm）*cos（2*M+sgm）-r*sin（2*sgm）*cos（4*M+2*sgm））/alfa;
A1=atan（sin（lmd）/（cos（u1）*tan（u2）-sin（u1）*cos（lmd）））;
A2=atan（sin（lmd）/（sin（u2）*cos（lmd）-tan（u1）*cos（u2）））;
A2=（m<
\n\n得到两点间大地线长S和大地正反方位角A1、A2如下：
s=%lf\n"
s）;
A1="
RBD（A1）;
voidGUS_ZS（）
doubleB,L,x3,x6,y3,y6,Y3,Y6,du,f,mi,X,N,n,t;
doubleAt,Bt,Ct,Dt,m3,m6,l3,l6,W,L03,L06;
intDH3,DH6;
请输入大地坐标（输入格式为角度（例如：
\n大地经度L="
L=RAD（du,f,mi）;
\n大地纬度B="
B=RAD（du,f,mi）;
At=1+3*e2/4+45*e2*e2/64+175*e2*e2*e2/256;
Bt=3*e2/4+15*e2*e2/16+525*e2*e2*e2/512;
Ct=15*e2*e2/64+105*e2*e2*e2/256;
Dt=35*e2*e2*e2/512;
X=a*（1-e2）*（At*B-Bt*sin（2*B）/2+Ct*sin（4*B）/4-Dt*sin（6*B）/6）;
n=sqrt（ep2）*cos（B）;
t=tan（B）;
DH3=（L-*PI/180））/（3*PI/180）+1;
DH6=L/（6*PI/180）+1;
L03=DH3*（3*PI/180）;
L06=DH6*（6*PI/180）-（3*PI/180）;
l3=L-L03;
l6=L-L06;
m3=cos（B）*l3;
m6=cos（B）*l6;
x3=X+N*t*（m3*m3/2+（5-t*t+9*n*n+4*n*n*n*n）*m3*m3*m3*m3/24+（61-58*t*t+t*t*t*t）*m3*m3*m3*m3*m3*m3/720）;
x6=X+N*t*（m6*m6/2+（5-t*t+9*n*n+4*n*n*n*n）*m6*m6*m6*m6/24+（61-58*t*t+t*t*t*t）*m6*m6*m6*m6*m6*m6/720）;
y3=N*（m3+（1-t*t+n*n）*m3*m3*m3/6+（5-18*t*t+t*t*t*t+14*n*n-58*n*n*t*t）*m3*m3*m3*m3*m3/120）;
y6=N*（m6+（1-t*t+n*n）*m6*m6*m6/6+（5-18*t*t+t*t*t*t+14*n*n-58*n*n*t*t）*m6*m6*m6*m6*m6/120）;
Y3=DH3*1000000+500000+y3;
Y6=DH6*1000000+500000+y6;
\n\n得到的高斯平面坐标为:
对于3度带:
\n纵坐标x=%.3lf\n横坐标y=%.3lf（通用坐标Y=%.3lf）\n\n"
x3,y3,Y3）;
对于6度带:
x6,y6,Y6）;
voidGUS_FS（）
doublex,y,Y,B,B0,B1,Bf,Vf,tf,Nf,nf,L,At,Bt,Ct,Dt,L3,L6;
longDH;
请输入高斯平面坐标：
纵坐标X="
x）;
自然坐标y="
y）;
通用坐标Y="
B0=x/（a*（1-e2）*At）;
B1=（x-a*（1-e2）*（-Bt*sin（2*B0）/2+Ct*sin（4*B0）/4-Dt*sin（6*B0）/6））/（a*（1-e2）*At）;
while（fabs（B1-B0）>
1*pow（10,-8））
B0=B1;
B1=（x-a*（1-e2）*（-Bt*sin（2*B0）/2+Ct*sin（4*B0）/4-Dt*sin（6*B0）/6））/（a*（1-e2）*At）;
Bf=B1;
nf=sqrt（ep2）*cos（Bf）;
tf=tan（Bf）;
Vf=sqrt（1+ep2*cos（Bf）*cos（Bf））;
Nf=c/Vf;
B=Bf-Vf*Vf*tf/2*（（y/Nf）*（y/Nf）-（5+3*tf*tf+nf*nf-9*nf*nf*tf*tf）*pow（（y/Nf）,4）/12+（61+90*tf*tf+45*tf*tf）*pow（（y/Nf）,6）/360）;
L=（（y/Nf）-（1+2*tf*tf+nf*nf）*（y/Nf）*（y/Nf）*（y/Nf）/6+（5+28*tf*tf+24*pow（tf,4）+6*nf*nf+8*nf*nf*tf*tf）*pow（（y/Nf）,5）/120）/cos（Bf）;
DH=Y/1000000;
L3=3*PI/180*double（DH）+L;
L6=6*PI/180*double（DH）-3*PI/180+L;
\n\n得到的大地坐标为:
若为6度带，大地经度L="
RBD（L6）;
若为3度带，大地经度L="
RBD（L3）;
升级会员 