白塞尔大地主题解算

白塞尔大地主题解算方向：学号：姓名：一．基本思路：基本思想：将椭球面上的大地元素按照白塞尔投影条件投影到辅助球面上，继而在球面上进行大地主题解算，最后在将球面上的计算结果换算到椭球面上。其关键问题是找出椭球面上的大地元素与球面上相应元素之间的关系式，同时解决在球面上进行大地主题解算的方法。正算流程：1.计算起点的归化纬度2.计算辅助函数值，解球面三角形3.按公式计算相关系数A,B,C以及α,β4.计算球面长度5.计算纬度差改正数6.计算终点大地坐标及大地方位角011122SBCA{sin(cos)}10101022222sin()sinsincoscos10101022222cos()coscossinsin001101522BCA[cos()]sin()010122LAsin[(sin()sin)]2111uuuAsinsincoscoscossin2222222222222222222222111111111euBuBWWueBuuBBarctaneueusinsincoscostantansinsintancos-cos1111AarctanuuAsinsin[]coscossinsincos21LL112111uAAarctanuAucossincoscoscossinsin反算流程：1.辅助计算2.用逐次趋近法同时计算起点大地方位角、球面长度及经差，第一次趋近时，取δ＝０。计算下式,重复上述计算过程2.3.计算大地线长度S4.计算反方位角二．已知数据序号B1（DD.MMSS)L1(DD.MMSS)A12(DD.MMSS)S12(m)141.01356874130.101226761.49438000L211212upAuuuuqsincostancossinsincoscos2121ppuqbbAarctanqsincoscos11pAqAsincostancos11pAqAsinsincos12aacoscosarctansincos011AuAsincossin111uAtantansec21+02122LAsin[(sinsin)]L11222222SABCBCsin(cos)sin(cos)1212uAbbcossinarctancos三．源代码：#includestdio.h#includemath.h#definee0.081813334016931499//克拉索夫斯基椭球体第一偏心率voidmain(){intk,B10,B11,L10,L11,A10,A11,B20,B21,L20,L21,A20,A21;doubleB12,L12,A12,B22,L22,A22;doubleB1,L1,A1,S,B2,L2,A2,L,pi;doubleA,B,C,afa,beta;doublea1,a2,b1,b2,p,q,x,y;doubleW1,W2,sinu1,sinu2,cosu1,cosu2,sinA0,cotsigma1,sin2sigma1,cos2sigma1,sigma0,sin2,cos2,sigma,sins,coss,delta0,delta,lamda;pi=4*atan(1);printf(白塞尔大地主题正算请输入1\n白塞尔大地主题反算请输入2\n);scanf(%d,&k);if(k==1){printf(请输入大地线起点纬度B经度L，大地方位角A及大地线长度S:\n);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=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*sinA0))*(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)+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(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);}else{if(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(A1)0){if(tan(A2)0)A2=pi+fabs(A2);elseA2=2*pi-fabs(A2);}else{if(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)A2/3600;A21=(int)A2/60-A20*60;A22=A2-A20*3600-A21*60;printf(正算得到的终点大地经度和大地纬度及A2:\n%d%d%lf\n%d%d%lf\n%d%d%lf\n,B20,B21,B22,L20,L21,L22,A20,A21,A22);}else{printf(请输入大地线起点和终点的坐标BL\n);scanf(%d%d%lf%d%d%lf%d%d%lf%d%d%lf,&B10,&B11,&B12,&L10,&L11,&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)*sqrt(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(p0){if(q0)A1=fabs(A1);elseA1=pi-fabs(A1);}else{if(q0)A1=2*pi-fabs(A1);elseA1=pi+fabs(A1);}sins=p*sin(A1)+q*cos(A1);//计算sigma的正弦值coss=a1+a2*cos(lamda);//计算sigma的余弦值sigma=atan(sins/coss);if(coss0)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=(28189-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(p0){if(q0)A1=fabs(A1);elseA1=pi-fabs(A1);}else{if(q0)A1=2*pi-fabs(A1);elseA1=pi+fabs(A1);}sins=p*sin(A1)+q*cos(A1);//计算sigma的正弦值coss=a1+a2*cos(lamda);//计算sigma的余弦值sigma=atan(sins/coss);if(coss0)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*sin