大地坐标与空间直角坐标的转换程序代码

#include "stdio.h"

#include "math.h"

#include "stdlib.h"

#include "iostream"

#define PI 3.1415926535897323

double a,b,c,e2,ep2;

int main()

{

int m,n,t;

double RAD(double d,double f,double m);

void RBD(double hd);

void BLH_XYZ();

void XYZ_BLH();

void B_ZS();

void B_FS();

void GUS_ZS();

void GUS_FS();

printf(" 大地测量学\n");

sp1:printf("请选择功能:\n");

printf("1.大地坐标系到大地空间直角坐标的转换\n");

printf("2.大地空间直角坐标到大地坐标系的转换\n");

printf("3.贝塞尔大地问题正算\n");

printf("4.贝塞尔大地问题反算\n");

printf("5.高斯投影正算\n");

printf("6.高斯投影反算\n");

printf("0.退出程序\n");

scanf("%d",&m);

if(m==0)exit(0);

sp2:printf("请选择椭球参数(输入椭球序号):\n");

printf("1.克拉索夫斯基椭球参数\n");

printf("2.IUGG_1975椭球参数\n");

printf("3.CGCS_2000椭球参数\n");

printf("0.其他椭球参数(自行输入)\n");

scanf("%d",&n);

switch(n)

{

case

1:a=6378245.0;b=6356863.0188;c=6399698.9018;e2=0.00669342162297;ep2=0.0067385254146 8;break;

case

2:a=6378140.0;b=6356755.2882;c=6399596.6520;e2=0.00669438499959;ep2=0.0067395018194 7;break;

case

3:a=6378137.0;b=6356752.3141;c=6399593.6259;e2=0.00669438002290;ep2=0.0067394967754 7;break;

case 0:{

printf("请输入椭球参数：\n");

printf("长半径a=");scanf("%lf",&a);

printf("短半径b=");scanf("%lf",&b);

c=a*a/b;

ep2=(a*a-b*b)/(b*b);

e2=(a*a-b*b)/(a*a);

break;

}

default:printf("\n\n输入错误!\n请重新输入!\n\n");goto sp2 ;

}

while(1)

{

switch(m)

{

case 1:BLH_XYZ();break;

case 2:XYZ_BLH();break;

case 3:B_ZS();break;

case 4:B_FS();break;

case 5:GUS_ZS();break;

case 6:GUS_FS();break;

default:printf("\n\n输入错误!\n请重新输入!\n\n");goto sp1 ;

}

printf("是否继续进行此功能计算? \n\n");

printf("（若继续进行此功能计算，则输入1；\n 若选择其他功能进行计算，则输入2；\n 若退出，则输入0. ）\n");

scanf("%d",&t);

switch(t)

{

case 1:break;

case 2:goto sp1;

case 0:exit(0);

}

}

}

double RAD(double d,double f,double m)

{

double e;

double sign=(d<0.0)?-1.0:1.0;

if(d==0)

{

sign=(f<0.0)?-1.0:1.0;

if(f==0)

{

sign=(m<0.0)?-1.0:1.0;

}

}

if(d<0)

d=d*(-1.0);

if(f<0)

f=f*(-1.0);

if(m<0)

m=m*(-1.0);

e=sign*(d*3600+f*60+m)*PI/(3600*180);

return e;

}

void RBD(double hd)

{

int t;

int d,f;

double m;

double sign=(hd<0.0)?-1.0:1.0;

if(hd<0)

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;

printf("%d'%d'%lf'\n",d,f,m);

}

void BLH_XYZ()

{

double B,L,H,N,W;

double d,f,m;

double X,Y,Z;

printf(" 请输入大地坐标（输入格式为角度(例如：30'40'50')）：\n");

printf(" 大地经度L=");

scanf("%lf'%lf'%lf'",&d,&f,&m);

L=RAD(d,f,m);

printf(" 大地纬度B=");

scanf("%lf'%lf'%lf'",&d,&f,&m);

B=RAD(d,f,m);

printf(" 大地高H=");

scanf("%lf",&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);

printf("\n\n 转换后得到大地空间直角坐标为：\n\n");

printf("X=%lf\nY=%lf\nZ=%lf\n\n",X,Y,Z);

}

void XYZ_BLH()

{

double B,L,H,N,W;

double X,Y,Z;

double tgB0,tgB1;

printf(" 请输入大地空间直角坐标：\n");

printf(" X=");

scanf("%lf",&X);

printf(" Y=");

scanf("%lf",&Y);

printf(" Z=");

scanf("%lf",&Z);

printf("\n\n 转换后得到大地坐标为：\n\n");

L=atan(Y/X);

printf(" 大地经度为：L=");

RBD(L);

printf("\n");

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);

printf(" 大地纬度为：B=");

RBD(B);

printf("\n");

W=sqrt(1-e2*sin(B)*sin(B));

N=a/W;

H=sqrt(X*X+Y*Y)/cos(B)-N;

printf(" 大地高为：H=%lf\n\n",H);

}

void B_ZS()

{

double L1,B1,A1,s,d,f,mi;

double u1,u2,m,M,k2,alfa,bt,r,kp2,alfap,btp,rp;

double sgm0,sgm1,lmd,lmd1,lmd2,A2,B2,l,L2;

printf("请输入已知点的大地坐标（输入格式为角度(例如：30'40'50')，下同）：\nL1=");

scanf("%lf'%lf'%lf'",&d,&f,&mi);

L1=RAD(d,f,mi);

printf("\nB1=");

scanf("%lf'%lf'%lf'",&d,&f,&mi);

B1=RAD(d,f,mi);

printf("请输入大地方位角：\nA1=");

scanf("%lf'%lf'%lf'",&d,&f,&mi);

A1=RAD(d,f,mi);

printf("请输入该点至另一点的大地线长：\ns=");

scanf("%lf",&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)>2.8*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:A2+PI;

u2=atan(-cos(A2)*tan(M+sgm0));

lmd1=atan(sin(u1)*tan(A1));

lmd1=(lmd1>0)?lmd1:lmd1+PI;

lmd1=(m

lmd2=atan(sin(u2)*tan(A2));

lmd2=(lmd2>0)?lmd2:lmd2+PI;

lmd2=(mPI)?lmd2: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*sg m0));

L2=L1+l;

printf("\n\n得到另一点的大地坐标和大地线在该点的大地方位角为：\n\n");

printf("L2=");

RBD(L2);printf("\n");

printf("B2=");

RBD(B2);printf("\n");

printf("A2=");

RBD(A2);

printf("\n");

}

void B_FS()

{

double L1,B1,L2,B2,s,A1,A2,du,f,mi,m0,m,M;

double l,u1,u2,alfa,bt,r,lmd0,dit_lmd,lmd,sgm,dit_sgm,sgm0,sgm1,alfap,btp,rp,k2,kp2;

printf("请输入第一个点大地坐标（输入格式为角度(例如：30'40'50')，下同）：\n大地经度L1=");

scanf("%lf'%lf'%lf'",&du,&f,&mi);

L1=RAD(du,f,mi);

printf("大地纬度B1=");

scanf("%lf'%lf'%lf'",&du,&f,&mi);

B1=RAD(du,f,mi);

printf("\n请输入第二个点大地坐标：\n大地经度：L2=");

scanf("%lf'%lf'%lf'",&du,&f,&mi);

L2=RAD(du,f,mi);

printf("大地纬度：B2=");

scanf("%lf'%lf'%lf'",&du,&f,&mi);

B2=RAD(du,f,mi);

l=L2-L1;

u1=atan(sqrt(1-e2)*tan(B1));

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=0.003351831*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>0)?A1:A1+PI;

M=atan(sin(u1)*tan(A1)/sin(m));

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;

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;

sgm0=acos(sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(l));

sgm1=acos(sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(l+sin(m)*(alfap*sgm0+btp*sin(sgm0)*cos( 2*M+sgm0))));

while(fabs(sgm0-sgm1)>1*PI/180*pow(10,-8))

{

sgm0=sgm1;

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)));

A1=(A1>0)?A1:A1+PI;

A1=(m>0)?A1:A1+PI;

A2=atan(sin(lmd)/(sin(u2)*cos(lmd)-tan(u1)*cos(u2)));

A2=(A2>0)?A2:A2+PI;

A2=(m<0)?A2:A2+PI;

printf("\n\n得到两点间大地线长S和大地正反方位角A1、A2如下：\n\n");

printf("s=%lf\n",s);

printf("A1=");

RBD(A1);printf("\n");

printf("A2=");

RBD(A2);printf("\n");

}

void GUS_ZS()

{

double B,L,x3,x6,y3,y6,Y3,Y6,du,f,mi,X,N,n,t;

double At,Bt,Ct,Dt,m3,m6,l3,l6,W,L03,L06;

int DH3,DH6;

printf("请输入大地坐标（输入格式为角度(例如：30'40'50')）：\n大地经度L=");

scanf("%lf'%lf'%lf'",&du,&f,&mi);

L=RAD(du,f,mi);

printf("\n大地纬度B=");

scanf("%lf'%lf'%lf'",&du,&f,&mi);

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);

W=sqrt(1-e2*sin(B)*sin(B));

N=a/W;

n=sqrt(ep2)*cos(B);

t=tan(B);

DH3=(L-(1.5*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*m 3*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*m 6*m6*m6/120);

Y3=DH3*1000000+500000+y3;

Y6=DH6*1000000+500000+y6;

printf("\n\n 得到的高斯平面坐标为:\n\n");

printf(" 对于3度带:\n 纵坐标x=%.3lf\n 横坐标y=%.3lf(通用坐标

Y=%.3lf)\n\n",x3,y3,Y3);

printf(" 对于6度带:\n 纵坐标x=%.3lf\n 横坐标y=%.3lf(通用坐标

Y=%.3lf)\n\n",x6,y6,Y6);

}

void GUS_FS()

{

double x,y,Y,B,B0,B1,Bf,Vf,tf,Nf,nf,L,At,Bt,Ct,Dt,L3,L6;

long DH;

printf(" 请输入高斯平面坐标：\n\n");

printf(" 纵坐标X=");

scanf("%lf",&x);printf("\n");

printf(" 自然坐标y=");

scanf("%lf",&y);printf("\n");

printf(" 通用坐标Y=");

scanf("%lf",&Y);printf("\n");

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;

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*t f*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;

printf("\n\n 得到的大地坐标为:\n\n");

printf(" 大地纬度B=");

RBD(B);printf("\n");

printf(" 若为6度带，大地经度L=");

RBD(L6);printf("\n");

printf(" 若为3度带，大地经度L=");

RBD(L3);printf("\n");

}