大地主题解算程序设计报告

大地主题结算实验报告姓名：高成学号：20104214测绘2班一、实验目的：1.提高运用计算机语言编程开发的能力；2.加深对大地主题解算计算公式及辅助参数的理解并掌握计算步骤；3.通过编程语言实现大地主题解算。二、工具：WindowsXPMode环境下的MicrosoftVisualC++6.0三、注意事项：1.计算所需变量多，容易混淆；2.正反算函数的编写；3.函数调用；4.弧度与角度之间的转化。四、实验要求：1.提交报告，实验总结，编写代码；2.独立编程，调试运行；3.上交成果：编写思想，编写过程，问题分析，源代码，计算结果；实验源代码：#includestdio.h#includemath.hdoublehudu(double,double,double);/*度分秒转换为弧度*/doubledu(double);/*弧度转换为度*/doublefen(double);/*弧度转换为分*/doublemiao(double);/*弧度转换为秒*/#definePI3.1415926voidmain(void){intk;printf(请选择大地主题算法，若执行正算，请输入1；若执行反算，请输入2。\n);scanf(%d,&k);/*大地主题正算*/if(k==1){doubleax,ay,az,bx,by,bz,cx,cy,cz,S,dz,ez,fz,B1,B2,L1,L2,A1,A2;intdx,dy,ex,ey,fx,fy;doublee2,W1,sinu1,cosu1,sinA0,coto1,sin2o1,cos2o1,sin2o,cos2o,A,B,C,r,t,o0,o,g,sinu2,q;/*输入度分秒数据*/printf(请输入大地线起点纬度度分秒\n);scanf(%lf%lf%lf,&ax,&ay,&az);printf(请输入大地线起点经度度分秒\n);scanf(%lf%lf%lf,&bx,&by,&bz);printf(请输入大地方位角度分秒\n);scanf(%lf%lf%lf,&cx,&cy,&cz);printf(请输入大地线长度\n);scanf(%lf,&S);/*调用函数*/B1=hudu(ax,ay,az);L1=hudu(bx,by,bz);A1=hudu(cx,cy,cz);/*白塞尔大地主题解算*/e2=0.006693421622966;W1=sqrt(1-e2*sin(B1)*sin(B1));sinu1=sin(B1)*(sqrt(1-e2))/W1;cosu1=cos(B1)/W1;sinA0=cosu1*sin(A1);coto1=cosu1*cos(A1)/sinu1;sin2o1=2*coto1/(coto1*coto1+1);cos2o1=(coto1*coto1-1)/(coto1*coto1+1);A=6356863.020+(10718.949-13.474*(1-sinA0*sinA0))*(1-sinA0*sinA0);B=(5354.469-8.978*(1-sinA0*sinA0))*(1-sinA0*sinA0);C=(2.238*(1-sinA0*sinA0))*(1-sinA0*sinA0)+0.006;r=691.46768-(0.58143-0.00144*(1-sinA0*sinA0))*(1-sinA0*sinA0);t=(0.2907-0.0010*(1-sinA0*sinA0))*(1-sinA0*sinA0);o0=(S-(B+C*cos2o1)*sin2o1)/A;sin2o=sin2o1*cos(2*o0)+cos2o1*sin(2*o0);cos2o=cos2o1*cos(2*o0)-sin2o1*sin(2*o0);o=o0+(B+5*C*cos2o)*sin2o/A;g=(r*o+t*(sin2o-sin2o1))*sinA0;/*求B2*/sinu2=sinu1*cos(o)+cosu1*cos(A1)*sin(o);B2=atan(sinu2/(sqrt(1-e2)*sqrt(1-sinu2*sinu2)));/*求L2*/q=atan(sin(A1)*sin(o)/(cosu1*cos(o)-sinu1*sin(o)*cos(A1)));/*判断q*/if(sin(A1)0&&tan(q)0)q=fabs(q);elseif(sin(A1)0&&tan(q)0)q=PI-fabs(q);elseif(sin(A1)0&&tan(q)0)q=-fabs(q);elseq=fabs(q)-PI;L2=L1+q-g/3600/180*PI;/*求A2*/A2=atan(cosu1*sin(A1)/(cosu1*cos(o)*cos(A1)-sinu1*sin(o)));/*判断A2*/if(sin(A1)0&&tan(A2)0)A2=fabs(A2);elseif(sin(A1)0&&tan(A2)0)A2=PI-fabs(A2);elseif(sin(A1)0&&tan(A2)0)A2=PI+fabs(A2);elseA2=2*PI-fabs(A2);/*调用函数*/dx=(int)(du(B2));dy=(int)(fen(B2));dz=miao(B2);ex=(int)(du(L2));ey=(int)(fen(L2));ez=miao(L2);fx=(int)(du(A2));fy=(int)(fen(A2));fz=miao(A2);printf(大地线终点纬度度分秒分别为：\n%d\n%d\n%lf\n,dx,dy,dz);printf(大地线终点经度度分秒分别为：\n%d\n%d\n%lf\n,ex,ey,ez);printf(终点大地方位角度分秒分别为：\n%d\n%d\n%lf\n,fx,fy,fz);}/*大地主题反算*/elseif(k==2){doublefx,fy,fz,gx,gy,gz,hx,hy,hz,ix,iy,iz,jz,kz,B1,B2,L1,L2,S,A1,A2;intjx,jy,kx,ky;doublee2,W1,W2,sinu1,sinu2,cosu1,cosu2,L,a1,a2,b1,b2,g,g2,g0=0.0,r,p,q,sino,coso,o,sinA0,x,t1,t2,A,B,C,y;/*输入度分秒数据*/printf(请输入大地线起点纬度度分秒\n);scanf(%lf%lf%lf,&fx,&fy,&fz);printf(请输入大地线起点经度度分秒\n);scanf(%lf%lf%lf,&gx,&gy,&gz);printf(请输入大地线终点纬度度分秒\n);scanf(%lf%lf%lf,&hx,&hy,&hz);printf(请输入大地线终点经度度分秒\n);scanf(%lf%lf%lf,&ix,&iy,&iz);/*调用函数*/B1=hudu(fx,fy,fz);L1=hudu(gx,gy,gz);B2=hudu(hx,hy,hz);L2=hudu(ix,iy,iz);/*白塞尔大地主题解算*/e2=0.006693421622966;W1=sqrt(1-e2*sin(B1)*sin(B1));W2=sqrt(1-e2*sin(B2)*sin(B2));sinu1=sin(B1)*sqrt(1-e2)/W1;sinu2=sin(B2)*sqrt(1-e2)/W2;cosu1=cos(B1)/W1;cosu2=cos(B2)/W2;L=L2-L1;a1=sinu1*sinu2;a2=cosu1*cosu2;b1=cosu1*sinu2;b2=sinu1*cosu2;/*逐次趋近法求解A1*/g=0;r=L;while(1){p=cosu2*sin(r);q=b1-b2*cos(r);A1=atan(p/q);/*判断A1*/if(p0&&q0)A1=fabs(A1);elseif(p0&&q0)A1=PI-fabs(A1);elseif(p0&&q0)A1=PI+fabs(A1);elseA1=2*PI-fabs(A1);sino=p*sin(A1)+q*cos(A1);coso=a1+a2*cos(r);o=atan(sino/coso);/*判断o*/if(coso0)o=fabs(o);elseo=PI-fabs(o);sinA0=cosu1*sin(A1);x=2*a1-(1-sinA0*sinA0)*coso;t1=(33523299-(28189-70*(1-sinA0*sinA0))*(1-sinA0*sinA0))*1e-10;t2=(28189-94*(1-sinA0*sinA0))*1e-10;g2=(t1*o-t2*x*sino)*sinA0;if(fabs(g2-g0)=1e-100)break;else{r=L+g2;g0=g2;}}/*求解S*/A=6356863.020+(10708.949-13.474*(1-sinA0*sinA0))*(1-sinA0*sinA0);B=10708.938-17.956*(1-sinA0*sinA0);C=4.487;y=((1-sinA0*sinA0)*(1-sinA0*sinA0)-2*x*x)*coso;S=A*o+(B*x+C*y)*sino;/*求解A2*/A2=atan(cosu1*sin(r)/(b1*cos(r)-b2));/*判断A2*/if(p0&&q0)A2=fabs(A2);elseif(p0&&q0)A2=PI-fabs(A2);elseif(p0&&q0)A2=PI+fabs(A2);elseA2=2*PI-fabs(A2);/*调用函数*/jx=(int)(du(A1));jy=(int)(fen(A1));jz=miao(A1);kx=(int)(du(A2));ky=(int)(fen(A2));kz=miao(A2);printf(起点大地方位角度分秒分别为：\n%d\n%d\n%lf\n,jx,jy,jz);printf(终点大地方位角度分秒分别为：\n%d\n%d\n%lf\n,kx,ky,kz);printf(大地线长度为：%lf\n,S);}/*数据错误*/elseprintf(数据错误，请重新输入\n);}/*度分秒转换为弧度*/doublehudu(doublex,doubley,doublez){doubleA0;A0=(x+y/60+z/3600)*PI/180;returnA0;}/*弧度转换为度*/doubledu(doubleB0){doublex0;x0=(int)(B0*180/PI);returnx0;}/*弧度转换为分*/doublefen(doubleC0){double_y,y0;_y=(int)(C0*180/PI);y0=(fabs)((int)((C0*180/PI-_y)*60));returny0;}/*弧度转换为秒*/doublemiao(doubleD0){double_z1,_z2,z0;_z1=(int)(D0*180/PI);_z2=(int)((D0*180/PI-_z1)*60);z0=(fabs)((double)(((D0*180/PI-_z1)*60-_z2)*60));returnz0;}大地主题正算结果大地主题反算结果