%白塞爾大地主題反算
clcclear;%清理命令窗口和工作區間
format?long?g;%規定輸出數據為長整型
a=6378245;b=6356863.0187730473;%克氏橢球
%a=6378140;b=6356755.2881575287;%1975國際橢球
c=a^2/b;e1=sqrt（a^2-b^2）/b;
raw=1;%計算相關橢球參數
B1=（47+46/60+52.6470/3600）*pi/180;
L1=（35+49/60+36.3300/3600）*pi/180;
B2=（48+04/60+09.6384/3600）*pi/180;
L2=（36+14/60+45.0505/3600）*pi/180;%輸入已知值（將角度值變成弧度值）
dL=L2-L1;
dB=B2-B1;
Bm=1/2*（B1+B2）;
Vm=sqrt（1+（e1*cos（Bm））^2）;
Nm=c/Vm;
tm=tan（Bm）;
gm=sqrt（e1^2*（cos（Bm））^2）;%計算大地線中點M處相關參數
r01=（Nm/raw）*cos（Bm）;
r21=（Nm*cos（Bm）/（24*raw^3*Vm^4））*（1+gm^2-9*（gm^2）*（tm^2）*（gm^4））;
r03=（-Nm/（24*raw^3））*（cos（Bm）^3）*tm^2;
S10=Nm/（raw*Vm^2）;
S12=（Nm/（24*（raw^3）*（Vm^2）））*（cos（Bm））^2*（+2+3*tm^2+2*gm^2）;
S30=Nm/（8*（raw^3）*（Vm^6））*（gm^2-（tm^2）*（gm^2））;%公式（4-227）
U=r01*dL+r21*dB^2*dL+r03*dL^3;
V=S10*dB+S12*dB*dL^2+S30*dB^3;%公式（4-226）
t01=tm*cos（Bm）;
t21=（1/（24*raw^2*Vm^4））*cos（Bm）*tm*（2+7*gm^2+9*tm^2*gm^2+5*gm^4）;
t03=（1/（24*raw^2））*（（cos（Bm））^3）*tm*（2+tm^2+2*gm^2）;%公式（4-229）
dA=t01*dL+t21*（dB^2）*dL+t03*dL^3;%公式（4-228）
Am=atan（U/V）;%公式（4-230）
S=U/sin（Am）;%計算大地線長度
A1=Am-1/2*dA;
A2=Am+1/2*dA;
if?（A2>pi）
????A2=A2-pi;
else
????A2=A2+pi;
end
%計算方位角
SA1A2%輸出的A1，A2為弧度值

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件???????1293??2016-11-26?15:40??Bessel?fansuan.m

?????文件???????1104??2016-11-26?15:34??Bessel?zhengsuan.m

?????文件??????????0??2018-05-03?08:44??CUST.txt

-----------?---------??----------?-----??----

?????????????????2397????????????????????3