白塞尔大地主题解算.docx
《白塞尔大地主题解算.docx》由会员分享，可在线阅读，更多相关《白塞尔大地主题解算.docx（13页珍藏版）》请在冰豆网上搜索。
白塞尔大地主题解算
白塞尔大地主题解算
。
方向：
学号：
姓名：
\
一．基本思路：
；
基本思想：
将椭球面上的大地元素按照白塞尔投影条件投影到辅助球面上，继而在球面上进行大地主题解算，最后在将球面上的计算结果换算到椭球面上。
其关键问题是找出椭球面上的大地元素与球面上相应元素之间的关系式，同时解决在球面上进行大地主题解算的方法。
正算流程：
1.计算起点的归化纬度
2.计算辅助函数值，解球面三角形
3.按公式计算相关系数A,B,C以及α,β
4.计算球面长度
#
5.计算纬度差改正数
6.计算终点大地坐标及大地方位角
、
反算流程：
1.辅助计算
2.用逐次趋近法同时计算起点大地方位角、球面长度及经差，第一次趋近时，取δ＝０。
!
计算下式,重复上述计算过程2.
3.计算大地线长度S
;
4.计算反方位角
二．已知数据
序号
B1（
L1
·
A12
S12（m）
1
130.
8000
…
三．源代码：
#include<>
#include<>
#definee//克拉索夫斯基椭球体第一偏心率
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=+B=C=*（1-sinA0*sinA0））*（1-sinA0*sinA0）+;
afa=beta=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）;
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;
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（p>0）
{
"
if（q>0）
A1=fabs（A1）;
else
A1=pi-fabs（A1）;
}
else
{
if（q>0）
A1=2*pi-fabs（A1）;
else
）
A1=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）;
else
sigma=pi-fabs（sigma）;
sinA0=cosu1*sin（A1）;
>
x=2*a1-（1-sinA0*sinA0）*cos（sigma）;
afa=（-（28189-70*（1-sinA0*sinA0））*（1-sinA0*sinA0））*;
beta=（28189-94*（1-sinA0*sinA0））*;
delta=（afa*sigma-beta*x*sin（sigma））*sinA0;
lamda=L+delta;
while（fabs（delta-delta0）>
{
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）;
}
else
{
》
if（q>0）
A1=2*pi-fabs（A1）;
else
A1=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）;
%
else
sigma=pi-fabs（sigma）;
sinA0=cosu1*sin（A1）;
x=2*a1-（1-sinA0*sinA0）*cos（sigma）;
afa=（-（28189-70*（1-sinA0*sinA0））*（1-sinA0*sinA0））*;
beta=（28189-94*（1-sinA0*sinA0））*;
delta=（afa*sigma-beta*x*sin（sigma））*sinA0;
lamda=L+delta;
}
A=+B=C=;
~
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）;
else
A2=2*pi-fabs（A2）;
}
~
else
{
if（tan（A2）>0）
A2=fabs（A2）;
else
A2=pi-fabs（A2）;
}
A1=A1*3600*180/pi;
A2=A2*3600*180/pi;
A10=（int）A1/3600;
A11=（int）A1/60-A10*60;
A12=A1-A10*3600-A11*60;
A20=（int）A2/3600;
A21=（int）A2/60-A20*60;
A22=A2-A20*3600-A21*60;
printf（"反算得到的方位角A1A2及大地线长S:
\n%d%d%lf\n%d%d%lf\n%lf\n",A10,A11,A12,A20,A21,A22,S）;
}
}
四．程序执行结果：
&
升级会员 