白塞尔大地主题解算

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)
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); //计算辅助函数值
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);
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;
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)); //计算起点规划纬度
1.计算起点的归化纬度
2.计算辅助函数值，解球面三角形
3.按公式计算相关系数A,B,C以及α,β
4.计算球面长度
5.计算纬度差改正数
6.计算终点大地坐标及大地方位角
反算流程：
1.辅助计算
2.用逐次趋近法同时计算起点大地方位角、球面长度及经差，第一次趋近时，取δ＝０。
计算下式,重复上述计算过程2.
3.计算大地线长度S
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("白塞尔大地主题正算请输入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);
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);
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);
lamda=fabs(lamda);
else
lamda=pi-fabs(lamda);
}
else
{
if(tan(lamda)>0)
lamda=fabs(lamda)-pi;
else
lamda=-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);
else
A2=2*pi-fabs(A2);
}
else
{
if(tan(A2)>0)
A2=fabs(A2);
else
A2=pi-fabs(A2);
}
B2=B2*180*3600/pi;
L2=L2*180*3600/pi;
4.计算反方位角
二．已知数据
序号
B1（DD.MMSS)
L1 (DD.MMSS)
A12(DD.MMSS)
S12(m)
1
41.01356874
130.10122676
1.4943
8000
三．源代码：
#include <stdio.h>
#include <math.h>
#define e 0.0818499 //克拉索夫斯基椭球体第一偏心率
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);