大地坐标与空间直角坐标转换_C程序

#include
#include

double HD(double a,double b,double c)
{b=b+c/60;
a=a+b/60;
a=a/180*3.1415926535897932384626433832795;
return a;}

void JD(double d,int *d1,int *d2,double *d3)
{d=d/3.1415926535897932384626433832795*180;
*d1=(int)d;
*d2=(int)((d-*d1)*60);
*d3=((d-*d1)*60-*d2)*60;}

void lbhxyz(double a1,double e)
{double l,b,h,x,y,z,a,b1,c,n,w;
printf("请输入L:\n");
scanf("%lf%lf%lf",&a,&b1,&c);
l=HD(a,b1,c);
printf("请输入B:\n");
scanf("%lf%lf%lf",&a,&b1,&c);
b=HD(a,b1,c);
printf("请输入H:\n");
scanf("%lf",&h);
w=sqrt((1-e*e*sin(b)*sin(b)));
n=a1/w;
x=(n+h)*cos(b)*cos(l);
y=(n+h)*cos(b)*sin(l);
z=(n*(1-e*e)+h)*sin(b);
printf("X=%f\nY=%f\nZ=%f\n",x,y,z);}

void xyzlbh(double a1,double e)
{double l,b,h,x,y,z,n,w,m,g,m1,g1,l3,b3;int l1,l2,b1,b2;
printf("请输入X:\n");
scanf("%lf",&x);
printf("请输入Y:\n");
scanf("%lf",&y);
printf("请输入Z:\n");
scanf("%lf",&z);
l=atan2(y,x);
do{m=(z+a1*e*e*g/sqrt(1+g*g-e*e*g*g))/sqrt(x*x+y*y);
m1=atan(m);
g1=atan(g);
m1=m1/3.1415926535897932384626433832795*180*3600;
g1=g1/3.1415926535897932384626433832795*180*3600;
g=m;}
while(sqrt((m1-g1)*(m1-g1))>0.0001);
b=atan(m);
w=sqrt(1-e*e*sin(b)*sin(b));
n=a1/w;
h=sqrt(x*x+y*y)/cos(b)-n;
JD(l,&l1,&l2,&l3);
JD(b,&b1,&b2,&b3);
printf("L=%d°%d'%2.6f''\nB=%d°%d'%2.6f''\nH=%5.6fm\n",l1,l2,l3,b1,b2,b3,h);}

void DY(int a)
{double a1,e,e2;int b,c;
printf("请选择：L,B,H→X,Y,Z(请输入1) 或 X,Y,Z→L,B,H(请输入2)\n");
scanf("%d",&b);
printf("请选择参考椭球面：\n克拉索夫斯基椭球(请输入1) \n或 IUGG 1975椭球(请输入2) \n或 CGCS 2000椭球(请输入3)\n");
scanf("%d",&c);
if (c==1){a1=6378245.0000,e2=0.00669342162297;}
else if (c==2){a1=6378140.0000,e2=0.00669438499959;}
else if (c==3){a1=6378137.0000,e2=0.00669438002290;};
e=sqrt(e2);
if (b==1)lbhxyz(a1,e);
else if(b==2)xyzlbh(a1,e);}

main()
{int a=1;
do{DY(a);
printf("继续使用请按1，结束使用请按0\n");
scanf("%d",&a);}
while(a);
}