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

1、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

2、 gus_zs();void gus_fs();printfc人地测量学 nh);sphprintf("请选择功能:n");printf("l.大地坐标系到大地空间直角坐标的转换n"> printf("2.大地空间直角坐标到大地坐标系的转换n"> printf("3.w塞尔大地问题正算n");printf("4.贝塞尔大地问题反算nn); printf("5.高斯投影正算川; printf("6.高斯投影反算|1" printf("o.退出程序 n&qu

3、ot;scanf("%d",&m);if(m=0)exit(0);sp2:printf("请选择椭球参数(输入椭球序号) printf("l.克拉索夫斯基椭球参数n"); printf("2.iugg_1975 椭球参数门")； printf("3.cgcs_2000 椭球参数|1"); printf("o.其他椭球参数自行输入)n"); scanf(,，o/od",&n);switch(n)casel:a=6378245.0;b=6356863.0188;c

4、=6399698.9018;e2=0.00669342162297;ep2=0.00673852541468;break;case2:a=6378140.0;b=6356755.2882;c=6399596.6520;e2=0.00669438499959;ep2=0.00673950181947;break;case3:a=6378137.0;b=6356752.3141;c=6399593.6259;e2=0.00669438002290;ep2=0.00673949677547;break;case 0:printf("请输入椭球参数:n"); printff 长半径

5、 a="scanf("%lf'&a); printff 短半径 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("nn 输入错误!n 请重新输入!nn"goto sp2 ;while switch(m)case l:blh_xyz();break; case 2:xyz_blh();break; case 3:b_zs();break; case 4:b_fs();break; case 5

6、:gus_zs();break; case 6:gus_fs();break;default:printf("nn 输入错误!n 请重新输入!nn;goto spl;printf("足否继续进行此功能计算？ nn-);printfc1 (若继续进行此功能计算，则输入1; n若选择其他功能进行计算,则输入2; n 若退出，则输入0.) n")；scanf(,，o/od",&t);switch(t)case l:break; case 2:goto spl; case 0:exit(0);double rad(double d,double f,do

7、uble m)double e;double sign=(d<0.0)?-1.0:1.0; if(d=o)sign=(f<o.o)?-l.o:l.o;if(f=o)sign=(m<0.0)?-l.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)?-l.0:1.0;

8、 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'%cf%lf|n",d,f,m>void blh_xyz()double b,l,h,n，w; double d,f,m; double x,y,z;printf("请输入大地坐标(输入格式为角度(例如：soo'so')： nu)；printf("大地经度 l=");scanf(m%lf%lf

9、，o/olf,/&d/&f/&m);l=rad(d,f,m);printf("人地纬度 b=");b=rad(d,f,m);printf("大地高h=m);scanf(,0/0|f'&h);w=sqrt(l-e2*sin(b)*sin(b);n=a/w;x=(n+h)*cos(b)*cos(l);y=(n+h)*cos(b)*sin(l);z=(n*(l-e2)+h)*sin(b);printf("nn转换后得到大地空间直角坐标为：w);printf("x=%lfny=%lfnz=%lfnn",

10、x,y,z);void xyz_blh()double b,l,h,n，w; double x,y,z; double tgbo,tgbl;phntfc请输入人地空叫直角坐标：n");printf("x=");scanfcw'&x);printf("y=h);scanf("%lf'&y>printff"z=");scanf(n%ir,&z);printf("nn转换后得到大地坐标为:nn");l=atan(y/x);printf("大地经度为：l=&

11、quot;);rbd(l);printf("nh);tgbo=z/sqrt(x*x+y*y);tgbl=(l/sqrt(x*x+y*y)*(z+a*e2*tgb0/sqrt(l+tgb0*tgb0-e2*tgb0*tgb0);while(fabs(tgb0-tgbl)>5*pow(10,-10)tgbo=tgbl;tgbl=(l/sqrt(x*x+y*y)*(z+a*e2*tgb0/sqrt(l+tgb0*tgb0-e2*tgb0*tgb0);b=atan(tgbl);printf("人地纬度为：b=h)；rbd(b);printf("n"w=sq

12、rt(l-e2*sin(b)*sin(b);n=a/w;h=sqrt(x*x+y*y)/cos(b)-n;printf("大地高为:h=%lfnn"/h);void b_zs()double ll,bl,al,s,d,f,mi;double ul,u2/m,mzk2/alfa/btzr/kp2/alfap,btp/rp;double sgmo,sgml,lmd,lmdl,lmd2,a2,b2,l,l2;printf("请输入己知点的大地坐标(输入格式为角度(例如：3ql40ls0l)f下同):nll="); scanf("%lf%lf%lf&q

13、uot;,&d,&f,&mi);ll=rad(d,f,mi);pnntf(,nbl=");bl=rad(d,f,mi);printf("请输入大地方位角：nal=">al=rad(d,fzmi);printf("请输入该点至另一点的大地线长:ns=h); scanf(,，o/olfh,&s);ul=atan(sqrt(l-e2)*tan(bl);m=asin(cos(ul)*sin(al);m=atan(tan(ul)/cos(al);m=(m>0)?m:m+2*pi;m=(m0)?m:m+pi;k2=ep2*

14、cos(m)*cos(m);alfa=(l-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;sgmo=alfa*s;sgml=alfa*s+bt*sin(sgmo)*cos(2*m+sgmo)+r*sin(2*sgmo)*cos(4*m+2*sgmo);while(fabs(sgm0-sgml)>2.8*pi/180*pow(10,-7)sgml=alfa*s+bt*sin(sgmo)*cos(2*m+sgmo)+r*sin(2*sgmo)*cos(4*

15、m+2*sgmo);sgmo=sgml;a2=atan(tan(m)/cos(m+sgmo);a2=(a2>0)?a2:a2+pi;a2=(a1>pi)?a2:a2+pi;u2=atan(-cos(a2)*tan(m+sgmo);lmdl=atan(sin(ul)*tan(al);lmdl=(lmdl>o)?lmdl:lmdl+pi；lmdl=(m<pi)?lmdl:lmdl+pi;lmd2=atan(sin(u2)*tan(a2);lmd2=(lmd2>0)?lmd2:lmd2+pi；lmd2=(m<pi)?(m+sgmo)<pi)?lmd2:lmd

16、2+pi):(m+sgmo)>pi)?lmd2:lmd2+pi);imd=lmd2-lmdl;b2=atan(sqrt(l+ep2)*tan(u2);kp2=e2*cos(m)*cos(m);alfap=(e2/2+e2*e2/8+e2*e2*e2/16)-e2/16*(l+e2)*kp2+3*e2*kp2*kp2/128;btp=e2*(l+e2)*kp2/16-e2*kp2*kp2/32;rp=e2*kp2*kp2/256;l=lmd-sin(m)*(alfap*sgmo+btp*sin(sgmo)*cos(2*m+sgmo)+rp*sin(2*sgmo)*cos(4*l/l+2*s

17、gmo);l2=l1+i;printf("nn得到另一点的大地呼标和大地线在该点的大地方位角为：nn"); printf("l2=");rbd(l2);printf(,n")；printf("b2="rbd(b2);printf("n"pnntf(,a2=");rbd(a2);printf("n");void b_fs()double ll,bl,l2,b2,s,al,a2,du,f,rrd,mo,m,m;double l,ul/u2/alfa/bt,rjmdo,dit_lm

18、djmd/sgm/dit_sgm,sgmo,sgml,alfap,btp,rp,k2,kp2;printf("请输入第一个点人地坐标(输入格式为角度(例如：30'40'50'),下同)：n人地经 度 ll=n);ll=rad(du,f,mi);printf("大地纬度 bl=");scanf(n%lf%lf%lf"',&du,&f,&mi>bl=rad(du,f,mi);printf("n请输入第二个点人地坐标：n人地经度：l2="); scanf(n%lf%lf%lf&q

19、uot;,&du,&f,&mi);l2=rad(du,f,mi);printf("大地纬度：b2=")；b2=rad(duj,mi);i=l2-l1;ul=atan(sqrt(l-e2)*tan(bl);u2=atan(sqrt(l-e2)*tan(b2);sgm0=acos(sin(ul)*sin(u2)+cos(ul)*cos(u2)*cos(l);mo=asin(cos(ul)*cos(u2)*sin(l)/sin(sgmo);dit_lmd=0.003351831*sgm0*sin(m0);lmdo=l+dit_lmd;dit_sgm=sin

20、(mo)*dit_lmd;sgml=sgmo+dit_sgm;m=asin(cos(ul)*cos(u2)*sin(lmdo)/sin(sgml);al=atan(sin(lmdo)/(cos(ul)*tan(u2)-sin(ul)*cos(lmdo);a1=(a1>o)?a1:a1+pi;al=(m>o)?al:al+pi;m=atan(sin(ul)*tan(al)/sin(m);m=(m>o)?m:m+pi;k2=ep2*cos(m)*cos(m);alfa=(l-k2/4+7*k2*k2/64-15*k2*k2*k2/256)/b;bt=k2/4-k2*k2/8+37

21、*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*(l+e2)*kp2+3*e2*kp2*kp2/128;btp=e2*(l+e2)*kp2/16-e2*kp2*kp2/32;rp=e2*kp2*kp2/256;sgm0=acos(sin(ul)*sin(u2)+cos(ul)*cos(u2)*cos(l);sgml=acos(sin(ul)*sin(u2)+cos(ul)*cos(u2)*cos(l+sin(m)*(alfap*sgmo+btp*

22、sin(sgmo)*cos(2*m+sgmo);while(fabs(sgm0-sgml)>l*pi/180*pow(10,-8)sgml=acos(sin(ul)*sin(u2)+cos(ul)*cos(u2)*cos(l+sin(m)*(alfap*sgmo+btp*sin(sgmo)*cos(2*m+sgmo);sgm=sgml;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;al=atan(sin(l

23、md)/(cos(ul)*tan(u2)-sin(ul)*cos(lmd);a1=(a1>o)?a1:a1+pi;al=(m>o)?al:al+pi;a2=atan(sin(lmd)/(sin(u2)*cos(lmd)-tan(ul)*cos(u2);a2=(a2>0)?a2:a2+pi;a2=(m<0)?a2:a2+pi;printf("nn得到两点fuj人地线长s和人地正反方位角al、a2如k: nn");pnntf("s=%lfn",s);printf("al=");rbd(al);printf(&quo

24、t;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*501): n大地经度！_=");l=rad(du,tmi);printf("n 大地纬度 b=");scanf(,%lf%lfo/olf,"

25、/&duz&f/&mi);b=rad(du,f,mi);at=l+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*(l-e2)*(at*b-bt*sin(2*b)/2+ct*sin(4*b)/4-dt*sin(6*b)/6);w=sqrt(l-e2*sin(b)*sin(b);n=a/w;n=sqrt(ep2)*cos(b);t=tan(b);dh3=(l-(1

26、.5*pi/180)/(3*pi/180)+l;dh6=l/(6*pi/180)+l;l03=dh3*(3*pi/180);l06=dh6*(6*pi/180)-(3*pi/180);i3=l-lo3;i6=l-lo6;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*

27、t*t+t*t*t*t)*m6*m6*m6*m6*m6*m6/720);y3=n*(m3+(l-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+(l-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;printf("nn得到的高斯平面坐标为:nn"pr

28、intf("对于3度带:n纵少标x=%.3lfn横.來标y=%.3lf(通用少标y=%.3lf)nn"/x3/y3,y3);printf("对于6度带:n纵坐称x=%.3lfn横坐标y=%.3lf(通用坐称y=o/o.3lf)nn"zx6,y6,y6);void gus_fs()double x,y,y,b,b0,bl,bf,vf,tf,nf,nfxat,bt,ct,dt,l3,l6; long dh;printf("请输入高斯平面坐标:w);printf(" 纵坐标 x=n);scanf("%lf",&x);printf("n");printf("自然坐标 y=h)；scanf"%lf.,&y;printfn");printff"通用坐标 y=");scanf(n%lfn,&y);printf