白塞尔大地主题解算

1、 白塞尔大地主题解算 方向： 学号： 姓名：一基本思路：基本思想：将椭球面上的大地元素按照白塞尔投影条件投影到辅助球面上，继而在球面上进行大地主题解算，最后在将球面上的计算结果换算到椭球面上。其关键问题是找出椭球面上的大地元素与球面上相应元素之间的关系式，同时解决在球面上进行大地主题解算的方法。正算流程：1.计算起点的归化纬度 2.计算辅助函数值，解球面三角形3.按公式计算相关系数A,B,C以及,4.计算球面长度5.计算纬度差改正数6.计算终点大地坐标及大地方位角 反算流程：1.辅助计算2.用逐次趋近法同时计算起点大地方位角、球面长度及经差，第一次趋近时，取。计算下式,重复上述计算过程2.3.

2、计算大地线长度S4.计算反方位角二已知数据序号B1（DD.MMSS)L1 (DD.MMSS)A12(DD.MMSS)S12(m)141.01356874130.101226761.49438000三源代码：#include <stdio.h>#include <math.h>void main()int k,B10,B11,L10,L11,A10,A11,B20,B21,L20,L21,A20,A21;double B12,L12,A12,B22,L22,A22;double B1,L1,A1,S,B2,L2,A2,L,pi;double A,B,C,afa,beta;

3、double a1,a2,b1,b2,p,q,x,y;double W1,W2,sinu1,sinu2,cosu1,cosu2,sinA0,cotsigma1,sin2sigma1,cos2sigma1,sigma0,sin2,cos2,sigma,sins,coss,delta0,delta,lamda;pi=4*atan(1);printf("白塞尔大地主题正算请输入1n白塞尔大地主题反算请输入2n");scanf("%d",&k);if(k=1) printf("请输入大地线起点纬度B经度L，大地方位角A及大地线长度S:n&quo

4、t;);scanf("%d%d%lf%d%d%lf%d%d%lf%lf",&B10,&B11,&B12,&L10,&L11,&L12,&A10,&A11,&A12,&S);B1=(B10+(float)B11/60+B12/3600)*pi/180; L1=(L10+(float)L11/60+L12/3600)*pi/180;A1=(A10+(float)A11/60+A12/3600)*pi/180;W1=sqrt(1-e*e*sin(B1)*sin(B1); /计算起点规划纬度sinu1=

5、sin(B1)*sqrt(1-e*e)/W1; /计算起点规划纬度cosu1=cos(B1)/W1; /计算起点规划纬度sinA0=cosu1*sin(A1); /计算辅助函数值cotsigma1=cosu1*cos(A1)/sinu1; /计算辅助函数值 sin2sigma1=2*cotsigma1/(cotsigma1*cotsigma1+1); /计算辅助函数值 cos2sigma1=(cotsigma1*cotsigma1-1)/(cotsigma1*cotsigma1+1); /计算辅助函数值 A=6356863.020+(10708.949-13.474*(1-sinA0*sinA

6、0)*(1-sinA0*sinA0);B=(5354.469-8.798*(1-sinA0*sinA0)*(1-sinA0*sinA0);C=(2.238*(1-sinA0*sinA0)*(1-sinA0*sinA0)+0.006;afa=691.46768-(0.58143-0.00144*(1-sinA0*sinA0)*(1-sinA0*sinA0);beta=(0.2907-1.0E-3*(1-sinA0*sinA0)*(1-sinA0*sinA0);sigma0=(S-(B+C*cos2sigma1)*sin2sigma1)/A;sin2=sin2sigma1*cos(2*sigma0

7、)+cos2sigma1*sin(2*sigma0);cos2=cos2sigma1*cos(2*sigma0)-sin2sigma1*sin(2*sigma0);sigma=sigma0+(B+5*C*cos2)*sin2/A;delta=(afa*sigma+beta*(sin2-sin2sigma1)*sinA0; /计算经度差改正数delta=delta/3600*pi/180;sinu2=sinu1*cos(sigma)+cosu1*cos(A1)*sin(sigma);B2=atan(sinu2/(sqrt(1-e*e)*sqrt(1-sinu2*sinu2);lamda=atan

8、(sin(A1)*sin(sigma)/(cosu1*cos(sigma)-sinu1*sin(sigma)*cos(A1);if(sin(A1)>0)if(tan(lamda)>0)lamda=fabs(lamda);elselamda=pi-fabs(lamda);elseif(tan(lamda)>0)lamda=fabs(lamda)-pi;elselamda=-1*fabs(lamda);L2=L1+lamda-delta;A2=atan(cosu1*sin(A1)/(cosu1*cos(sigma)*cos(A1)-sinu1*sin(sigma); if(sin

9、(A1)>0)if(tan(A2)>0)A2=pi+fabs(A2);elseA2=2*pi-fabs(A2);elseif(tan(A2)>0)A2=fabs(A2);elseA2=pi-fabs(A2);B2=B2*180*3600/pi;L2=L2*180*3600/pi;A2=A2*180*3600/pi;B20=(int)B2/3600;B21=(int)B2/60-B20*60;B22=B2-B20*3600-B21*60;L20=(int)L2/3600;L21=(int)L2/60-L20*60;L22=L2-L20*3600-L21*60;A20=(int)

10、A2/3600;A21=(int)A2/60-A20*60;A22=A2-A20*3600-A21*60;printf("正算得到的终点大地经度和大地纬度及A2:n%d %d %lfn%d %d %lfn%d %d %lfn",B20,B21,B22,L20,L21,L22,A20,A21,A22);else printf("请输入大地线起点和终点的坐标BLn");scanf("%d%d%lf%d%d%lf%d%d%lf%d%d%lf",&B10,&B11,&B12,&L10,&L11,&

11、;L12,&B20,&B21,&B22,&L20,&L21,&L22);B1=(B10+(double)B11/60+B12/3600)*pi/180; L1=(L10+(double)L11/60+L12/3600)*pi/180;B2=(B20+(double)B21/60+B22/3600)*pi/180; L2=(L20+(double)L21/60+L22/3600)*pi/180; W1=sqrt(1-e*e*sin(B1)*sin(B1);W2=sqrt(1-e*e*sin(B2)*sin(B2);sinu1=sin(B1)*sqr

12、t(1-e*e)/W1;sinu2=sin(B2)*sqrt(1-e*e)/W2;cosu1=cos(B1)/W1;cosu2=cos(B2)/W2;L=L2-L1; a1=sinu1*sinu2;a2=cosu1*cosu2;b1=cosu1*sinu2;b2=sinu1*cosu2;delta0=0;lamda=L+delta0;p=cosu2*sin(lamda);q=b1-b2*cos(lamda);A1=atan(p/q);if(p>0)if(q>0)A1=fabs(A1);elseA1=pi-fabs(A1);elseif(q>0)A1=2*pi-fabs(A1)

13、;elseA1=pi+fabs(A1);sins=p*sin(A1)+q*cos(A1); /计算sigma的正弦值coss=a1+a2*cos(lamda); /计算sigma的余弦值sigma=atan(sins/coss); if(coss>0)sigma=fabs(sigma);elsesigma=pi-fabs(sigma);sinA0=cosu1*sin(A1); x=2*a1-(1-sinA0*sinA0)*cos(sigma);afa=(33523299-(28189-70*(1-sinA0*sinA0)*(1-sinA0*sinA0)*1.0e-10;beta=(281

14、89-94*(1-sinA0*sinA0)*1.0e-10;delta=(afa*sigma-beta*x*sin(sigma)*sinA0;lamda=L+delta;while(fabs(delta-delta0)>4.8e-10)delta0=delta; p=cosu2*sin(lamda); q=b1-b2*cos(lamda); A1=atan(p/q); if(p>0)if(q>0)A1=fabs(A1); else A1=pi-fabs(A1); elseif(q>0)A1=2*pi-fabs(A1); else A1=pi+fabs(A1); sins

15、=p*sin(A1)+q*cos(A1); /计算sigma的正弦值 coss=a1+a2*cos(lamda); /计算sigma的余弦值 sigma=atan(sins/coss); if(coss>0)sigma=fabs(sigma); else sigma=pi-fabs(sigma); sinA0=cosu1*sin(A1); x=2*a1-(1-sinA0*sinA0)*cos(sigma); afa=(33523299-(28189-70*(1-sinA0*sinA0)*(1-sinA0*sinA0)*1.0e-10; beta=(28189-94*(1-sinA0*si

16、nA0)*1.0e-10; delta=(afa*sigma-beta*x*sin(sigma)*sinA0; lamda=L+delta;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)*cos(sigma);S=A*sigma+(B*x+C*y)*sin(sigma);A2=atan(cosu1*sin(lamda)/(b1*cos(lamda)-b2);if(sin(A1)>0)if(tan(A2)>0)A2=pi+fabs(A2);elseA2=2*pi-fabs(A2);elseif(tan(A2)>0)A2=fabs(A2);elseA2=pi-fabs(A2);