白塞尔大地主题解算

白塞尔大地主题解算班级测绘1112学号1120115230姓名陈晓旭一、正算代码：#includestdio.h#includemath.h#defineee0.006693421622966#defineI3.141592653doubleF(double,double,double);voidmain(void){doubleA1,B1,L1,S,A2,B2,L2;doublex1,x2,x3,y1,y2,y3,z1,z2,z3;doubleW1,sinu1,sinu2,cosu1,sinA0;doublecota1,cos2a1,sin2a1,cosA0A0;doubleA,B,C,d,e,a0,a1,m;doublen,a,Q,R;printf(请输入数据B1=);scanf(%lf%lf%lf,&x1,&x2,&x3);B1=F(x1,x2,x3);printf(请输入数据L1=);scanf(%lf%lf%lf,&y1,&y2,&y3);L1=F(y1,y2,y3);printf(请输入A1=);scanf(%lf%lf%lf,&z1,&z2,&z3);A1=F(z1,z2,z3);printf(请输入S=);scanf(%lf,&S);printf(B1=%.9f\n,B1);printf(L1=%.9f\n,L1);printf(A1=%.9f\n,A1);printf(S=%f\n,S);/*计算起点的规划纬度*/W1=sqrt(1-ee*sin(B1)*sin(B1));sinu1=sin(B1)*sqrt(1-ee)/W1;cosu1=cos(B1)/W1;printf(W1=%.9f\n,W1);printf(sinu1=%.9f\n,sinu1);printf(cosu1=%.9f\n,cosu1);/*计算辅助函数值*/sinA0=cosu1*sin(A1);cota1=cosu1*cos(A1)/sinu1;sin2a1=2*cota1/(cota1*cota1+1);cos2a1=(cota1*cota1-1)/(cota1*cota1+1);printf(sinA0=%.9f\n,sinA0);printf(cota1=%.9f\n,cota1);printf(sin2a1=%.9f\n,sin2a1);printf(cos2a1=%.9f\n,cos2a1);/*计算系数ABC及de*/cosA0A0=1-sinA0*sinA0;A=6356755.288+(10710.341-(13.534*cosA0A0))*cosA0A0;B=(5355.171-9.023*cosA0A0)*cosA0A0;C=(2.256*(cosA0A0))*cosA0A0+0.006;d=691.46768-(0.58143-0.00144*cosA0A0)*cosA0A0;e=(0.2907-cosA0A0*0.0010)*cosA0A0;printf(cosA0A0=%.9f\n,cosA0A0);printf(A=%.3f\n,A);printf(B=%.6f\n,B);printf(C=%.9f\n,C);printf(d=%.7f\n,d);printf(e=%.9f\n,e);/*计算球面长度*/a0=(S-(B+C*cos2a1)*sin2a1)/A;m=sin2a1*cos(2*a0)+cos2a1*sin(2*a0);n=(cos2a1)*(cos(2*a0))-(sin2a1)*(sin(2*a0));a=a0+((B+5*C*n))*m/A;printf(a0=%.9f\n,a0);printf(m=%.9f\n,m);printf(n=%.9f\n,n);printf(a=%.9f\n,a);/*计算经度差改正数*/Q=(d*a+(e*(m-sin2a1))/3600/180*I)*sinA0;printf(Q=%.7f\n,Q);/*计算终点大地坐标及大地方位角*/sinu2=sinu1*cos(a)+cosu1*cos(A1)*sin(a);B2=180*atan(sinu2/((sqrt(1-ee))*(sqrt(1-sinu2*sinu2))))/I;R=180*atan(sin(A1)*sin(a)/(cosu1*cos(a)-sinu1*sin(a)*cos(A1)))/I;printf(sinu2=%.9f\n,sinu2);printf(B2=%f\n,B2);printf(R=%f\n,R*180/I);/*确定R的值*/if(sin(A1)0&&tan(R)0)R=abs(R);elseif(sin(A1)0&&tan(R)0)R=I-abs(R);elseif(sin(A1)0&&tan(R)0)R=-abs(R);elseR=abs(R)-I;/*确定L2A2的值*/L2=(L1*180/I+R-(Q/206265*180/I));A2=atan(cosu1*sin(A1)/(cosu1*cos(a)*cos(A1)-sinu1*sin(a)));if(sin(A1)0&&tan(A2)0)A2=(fabs(A2))*180/I;elseif(sin(A1)0&&tan(A2)0)A2=(I-fabs(A2))*180/I;elseif(sin(A1)0&&tan(A2)0)A2=(I+fabs(A2))*180/I;elseA2=(2*I-fabs(A2))*180/I;printf(A2=%3f\nB2=%3f\nL2=%3f\n,A2,B2,L2);}doubleF(doublea2,doubleb2,doublec2){doubled2;d2=(double)(a2+1.0*b2/60+1.0*c2/3600);d2=(d2/180)*I;return(d2);}运行结果：二、反算代码：#includestdio.h#includemath.h#defineee0.006693421622966#defineI3.141592653doubleF(double,double,double);voidmain(void){doubleB1,B2,L1,L2,A1,A2,S,Y;doubleW1,W2,L,Q,R,A,B,C;doublex,y,z,p,q;doublex1,x2,x3,y1,y2,y3,z1,z2,z3,w1,w2,w3;doublea1,a2,b1,b2,m,n;doublesinp,cosp,sinu1,sinu2,cosu1,cosu2,sinA0,cosA0;Q=0;q=0;printf(请输入起始数据B1=);scanf(%lf%lf%lf,&x1,&x2,&x3);B1=F(x1,x2,x3);printf(请输入起始数据L1=);scanf(%lf%lf%lf,&y1,&y2,&y3);L1=F(y1,y2,y3);printf(请输入起始数B2=);scanf(%lf%lf%lf,&z1,&z2,&z3);B2=F(z1,z2,z3);printf(请输入起始数据L2=);scanf(%lf%lf%lf,&w1,&w2,&w3);L2=F(w1,w2,w3);printf(B1=%f\n,B1);printf(L1=%.9f\n,L1);printf(B2=%.9f\n,B2);printf(L2=%.9f\n,L2);/*辅助计算*/W1=sqrt(1-ee*sin(B1)*sin(B1));W2=sqrt(1-ee*sin(B2)*sin(B2));sinu1=sin(B1)*sqrt(1-ee)/W1;sinu2=sin(B2)*(sqrt(1-ee))/W2;cosu1=cos(B1)/W1;cosu2=cos(B2)/W2;L=L2-L1;a1=sinu1*sinu2;a2=cosu1*cosu2;b1=cosu1*sinu2;b2=sinu1*cosu2;printf(W1=%.9f\n,W1);printf(W2=%.9f\n,W2);printf(sinu1=%.9f\n,sinu1);printf(sinu2=%.9f\n,sinu2);printf(cosu1=%.9f\n,cosu2);printf(L=%.9f\n,L);printf(a1=%.9f\n,a1);printf(a2=%.9f\n,a2);printf(b1=%.9f\n,b1);printf(b2=%.9f\n,b2);/*用逐次趋近法同时计算起点大地方位角、球面长度及经差R*/do{R=L+Q;x=cosu2*sin(R);y=b1-b2*cos(R);A1=atan(x/y);if(x0&&y0)A1=fabs(A1);elseif(x0&&y0)A1=I-fabs(A1);elseif(x0&&y0)A1=I+fabs(A1);elseA1=2*I-fabs(A1);sinp=x*sin(A1)+y*cos(A1);cosp=a1+a2*cos(R);p=atan(sinp/cosp);if(cosp0)p=fabs(p);sinA0=cosu1*sin(A1);cosA0=sqrt(1-sinA0*sinA0);p=I-fabs(p);z=2*a1-cosA0*cosA0*cos(p);m=(33523299-(28189-70*cosA0*cosA0))*(1e-10);n=(28189-94*cosA0*cosA0)*(1e-10);Q=q;q=(m*p-m*z*sin(p))*sinA0;}while(fabs(206265*(q-Q))(1e-4));printf(Q=%f\n,Q);/*计算系数ABC及大地线长度*/A=6356863.020+(10708.949-13.474*cosA0*cosA0)*cosA0*cosA0;B=10708.938-17.956*cosA0*cosA0;C=4.487;Y=(cosA0*cosA0*cosA0*cosA0-2*z*z)*cos(p);S=A*p+(B*z+C*Y)*sinp;/*计算反方位角*/A2=atan((cosu1*sin(R))/(b1*cos(R)-b2));A2=A2*180/I;A1=A1*180/I;if(A10)A2=abs(A2);elseA2=-abs(A2);printf(A1=%3f\nA2=%3f\nS=%3f\n,A1,A2,S);}doubleF(doublea,doubleb,doublec){doubled;if(a0)d=(double)(a-1.0*b/60-1.0*c/3600);elsed=(double)(a+1.0*b/60+1.0*c/3600);d=(d/180)*I;returnd;}运行结果：