白塞爾大地主題解算

1、 白塞爾大地主題解算 方向： 學(xué)號： 姓名：一基本思路：基本思想：將橢球面上的大地元素按照白塞爾投影條件投影到輔助球面上，繼而在球面上進行大地主題解算，最后在將球面上的計算結(jié)果換算到橢球面上。其關(guān)鍵問題是找出橢球面上的大地元素與球面上相應(yīng)元素之間的關(guān)系式，同時解決在球面上進行大地主題解算的方法。正算流程：1.計算起點的歸化緯度 2.計算輔助函數(shù)值，解球面三角形3.按公式計算相關(guān)系數(shù)A,B,C以及,4.計算球面長度5.計算緯度差改正數(shù)6.計算終點大地坐標(biāo)及大地方位角 反算流程：1.輔助計算2.用逐次趨近法同時計算起點大地方位角、球面長度及經(jīng)差，第一次趨近時，取。計算下式,重復(fù)上述計算過程2.3.

2、計算大地線長度S4.計算反方位角二已知數(shù)據(jù)序號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經(jīng)度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); /計算起點規(guī)劃緯度sinu1=

5、sin(B1)*sqrt(1-e*e)/W1; /計算起點規(guī)劃緯度cosu1=cos(B1)/W1; /計算起點規(guī)劃緯度sinA0=cosu1*sin(A1); /計算輔助函數(shù)值cotsigma1=cosu1*cos(A1)/sinu1; /計算輔助函數(shù)值 sin2sigma1=2*cotsigma1/(cotsigma1*cotsigma1+1); /計算輔助函數(shù)值 cos2sigma1=(cotsigma1*cotsigma1-1)/(cotsigma1*cotsigma1+1); /計算輔助函數(shù)值 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; /計算經(jīng)度差改正數(shù)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("正算得到的終點大地經(jīng)度和大地緯度及A2:n%d %d %lfn%d %d %lfn%d %d %lfn",B20,B21,B22,L20,L21,L22,A20,A21,A22);else printf("請輸入大地線起點和終點的坐標(biāo)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);