分段累加法大地主题解算与高斯投影C语言程序

大地主题解算与高斯投影正反算程序设计测绘一班XX2014301610096一、大地主题解算：程序源码：#includestdio.h#includemath.h#definePI3.1415926535898intmain(void){doublebx,by,bz,B1,lx,ly,lz,L1,ax,ay,az,A1,S;doubledB,dL,dA;doublee2=0.006693421622966,a=6378245.0000000000;doubleM1,N1,W1,N2,W2;doubleB2,L2,A2;intBx,By,Lx,Ly,Ax,Ay;doubleBz,Lz,Az;printf(请输入大地线起点纬度\n°′″\n);scanf(%lf%lf%lf,&bx,&by,&bz);printf(请输入大地线起点经度\n°′″\n);scanf(%lf%lf%lf,&lx,&ly,&lz);printf(请输入大地方位角\n°′″\n);scanf(%lf%lf%lf,&ax,&ay,&az);printf(请输入大地线长度\n);scanf(%lf,&S);doubleW,M,N,C;B1=(bx+by/60+bz/3600)*PI/180;L1=(lx+ly/60+lz/3600)*PI/180;A1=(ax+ay/60+az/3600)*PI/180;W=sqrt(1-e2*(sin(B1)*sin(B1)));M=a*(1-e2)/(W*W*W);N=a/W;C=N*cos(B1)*sin(A1);intn,i;doubles;n=(int)(S/4000+1);s=S/n;for(i=1;i=n;i++){dB=cos(A1)*s/M;dL=sin(A1)*s/N/cos(B1);B2=B1+dB;L2=L1+dL;W=sqrt(1-e2*(sin(B2)*sin(B2)));M=a*(1-e2)/(W*W*W);N=a/W;A2=asin(C/N/cos(B2));B1=B2;L1=L2;A1=A2;}Bx=B2/PI*180;By=(B2/PI*180-Bx)*60;Bz=fabs(((B2/PI*180-Bx)*60-By)*60);Lx=L2/PI*180;Ly=(L2/PI*180-Lx)*60;Lz=fabs(((L2/PI*180-Lx)*60-Ly)*60);Ax=A2/PI*180;Ay=(A2/PI*180-Ax)*60;Az=fabs(((A2/PI*180-Ax)*60-Ay)*60);printf(大地线终点纬度为%d°%d′%f″\n,Bx,By,Bz);printf(大地线终点经度为%d°%d′%f″\n,Lx,Ly,Lz);printf(终点大地方位角为%d°%d′%f″\n,Ax,Ay,Az);return0;程序运行截图：高斯正算程序源码：#includestdio.h#includemath.h#definePI3.1415926535897932intmain(void){printf(高斯投影算\n\n);into;printf(请选择采用的椭球参数\n1.克拉索夫斯基椭球体2.1975年国际椭球体3.WGS-84椭球体4.CGCS2000\n);scanf(%d,&o);doublea=0,e12,e2;if(o==1){a=6378245.0000000000;e2=0.006693421622966;e12=0.006738525414683;}elseif(o==2){a=6378140.0000000000;e2=0.006694384999588;e12=0.006739501819473;}elseif(o==3){a=6378137.0000000000;e2=0.0066943799013;e12=0.00673949674227;}elseif(o==4){a=6378137.0;e2=0.00669438002290;e12=0.00673949677548;}intk;printf(\n执行高斯投影算法正算\n\n);doublebx,by,bz,lx,ly,lz;doublem0,m2,m4,m6,m8,a0,a2,a4,a6,a8,N,B,L,X,Bb,cosB,sinB,ρ,η2,t,l,d,x,y;printf(请输入大地纬度\n°′″\n);scanf(%lf%lf%lf,&bx,&by,&bz);printf(请输入大地经度\n°′″\n);scanf(%lf%lf%lf,&lx,&ly,&lz);B=bx+by/60+bz/3600;L=lx+ly/60+lz/3600;Bb=B*PI/180;m0=a*(1-e2);m2=3*e2*m0/2;m4=5*e2*m2/4;m6=7*e2*m4/6;m8=9*e2*m6/8;a0=m0+(m2/2)+((3*m4)/8)+((5*m6)/16)+((35*m8)/128);a2=m2/2+m4/2+15*m6/32+7*m8/16;a4=m4/8+3*m6/16+7*m8/32;a6=m6/32+m8/16;a8=m8/128;ρ=180*3600/PI;cosB=cos(Bb);sinB=sin(Bb);η2=e12*cosB*cosB;t=tan(Bb);d=(int)(L/6)+1;l=abs(L-(6*d-3))*3600;N=a/sqrt(1-e2*sinB*sinB);X=a0*Bb-sinB*cosB*((a2-a4+a6)+(2*a4-16*a6/3)*sinB*sinB+16*a6*sinB*sinB*sinB*sinB/3);x=X+N*sinB*cosB*l*l/(2*ρ*ρ)+N*sinB*cosB*cosB*cosB*(5-t*t+9*η2)*l*l*l*l/(24*ρ*ρ*ρ*ρ)+N*sinB*cosB*cosB*cosB*cosB*cosB*(61-58*t*t+t*t*t*t)*l*l*l*l*l*l/(720*ρ*ρ*ρ*ρ*ρ*ρ);y=N*cosB*l/ρ+N*cosB*cosB*cosB*(1-t*t+η2)*l*l*l/(6*ρ*ρ*ρ)+N*cosB*cosB*cosB*cosB*cosB*(5-18*t*t+t*t*t*t)*l*l*l*l*l/(120*ρ*ρ*ρ*ρ*ρ);printf(x为%.4lfm\n,x);printf(y为%.4lfm\n,y);}程序截图：高斯反算程序源码：#includestdio.h#includemath.h#definePI3.1415926535898intmain(void){printf(高斯投影反算\n\n);into;printf(请选择采用的椭球参数\n1.克拉索夫斯基椭球体2.1975年国际椭球体3.WGS-84椭球体4.CGCS2000\n);scanf(%d,&o);doublea,e12,e2;if(o==1){a=6378245.0000000000;e2=0.006693421622966;e12=0.006738525414683;}elseif(o==2){a=6378140.0000000000;e2=0.006694384999588;e12=0.006739501819473;}elseif(o==3){a=6378137.0000000000;e2=0.0066943799013;e12=0.00673949674227;}elseif(o==4){a=6378137.0;e2=0.00669438002290;e12=0.00673949677548;}doublex,y;doublem0,m2,m4,m6,m8,a0,a2,a4,a6,a8,Bf,Bf1,B,Mf,W,Nf,tf,ηf2,n2,n4,n6,n8,l;printf(plaeseenterx,y:\n);scanf(%lf%lf,&x,&y);m0=a*(1-e2);m2=3*e2*m0/2;m4=5*e2*m2/4;m6=7*e2*m4/6;m8=9*e2*m6/8;a0=m0+(m2/2)+((3*m4)/8)+((5*m6)/16)+((35*m8)/128);a2=m2/2+m4/2+15*m6/32+7*m8/16;a4=m4/8+3*m6/16+7*m8/32;a6=m6/32+m8/16;a8=m8/128;Bf1=0;Bf=x/a0;while(abs(Bf-Bf1)1E-10){Bf1=Bf;doublesinB=sin(Bf1);doublecosB=cos(Bf1);doublesin2B=sinB*cosB*2;doublesin4B=sin2B*(1-2*sinB*sinB)*2;doublesin6B=sin2B*sqrt(1-(sin4B*sin4B)+sin4B*sqrt(1-sin2B*sin2B));Bf=(x-(-a2/2.0*sin2B+a4/4.0*sin4B-a6/6.0*sin6B))/a0;}doublesinBf=sin(Bf);doublecosBf=cos(Bf);ηf2=e12*cosBf*cosBf;tf=tan(Bf);W=sqrt(1-e2*sinBf*sinBf);Mf=a*(1-e2)/(W*W*W);Nf=a/W;B=Bf-tf*y*y/(2*Mf*Nf)+tf*(5+3*tf*tf+ηf2-9*ηf2*tf*tf)*y*y*y*y/(24*Mf*Nf*Nf*Nf)-tf*(61+90*tf*tf+45*tf*tf*tf*tf)*y*y*y*y*y*y/(720*Mf*Nf*Nf*Nf*Nf*Nf);l=y/(Nf*cosBf)-(1+2*tf*tf+ηf2)*y*y*y/(6*Nf*Nf*Nf*cosBf)+(5+28*tf*tf+24*tf*tf*tf*tf+6*ηf2+8*ηf2*tf*tf)*y*y*y*y*y/(120*Nf*Nf*Nf*Nf*Nf*cosBf);intBx,By,lx,ly;doubleBz,lz;Bx=(B/PI*180);By=((B/PI*180-Bx)*60);Bz=abs(((B/PI*180-Bx)*60-By)*60);By=abs(By);lx=(l/PI*180);ly=((l/PI*180-lx)*60);lz=abs(((l/PI*180-lx)*60-ly)*60);ly=abs(ly);printf(大地纬度为%d°%d′%lf″\n,Bx,By,Bz);printf(经度差为%d°%d′%lf″\n,lx,ly,lz);return0;}程序截图：