贝塞尔大地主题解算分析 PPT
贝塞尔函数PPT演示课件
1
r 2 sin 2
2u
2
k 2u
0
设u(r, ,) R(r)( )()，代入原方程
''() m2() 0
1
s in
d
d
s in
d ( 2 d
m2
sin 2 ) 0
d r 2 dR (k 2r 2 2 )R 0
要使等式两边成立，则x各次幂的系数为零
(1) (c2 v2 ) C0 0 (k 0)
(c2 v2 ) 0
c v
(2) [(c 1)2 v2 ]C1 0 (k 1)
(3) [(c k)2 v2 ]Ck Ck2 0 (k 2)
将c=v代入(2)，得C1=0
k 2u
0
u(,, z) R()()Z(z)
''() m2() 0
Z''(z) 2Z(z) 0
2
d 2R
d 2
dR
d
(k 2
2 ) 2
m2
R
0
x (k 2 2) y(x) R()
贝塞尔方程
x2
0
0
0
0
(1) etdt et 1 0 0
(2) 1 (1) 1
(3) 2 (2) 2!
(4) 3(3) 3! (n 1) n!
求证： 1 2
(x) ett x1dt
令t=u2
(1)m
2(2mv) m ! (m 1 v)
白赛尔大地主题解算
dS
M
cos A c
dL
dS
sin A N cos B
V c
sec B sin B
dA
dS
tan B sin N
AV c
tan B sin
A
二阶导数：
d2B dS 2
B
( dB dS
)
dB dS
A
( dB dS
)
dA dS
V4 c2
t (3 2
10
Fundation of Geodesy
4.7.3 高斯平均引数正算公式
高斯平均引数正算公式推导 的基本思想：
首先把勒让德级数在 P１点展 开改在大地线长度中点M展开，以 使级数公式项数减少，收敛快， 精度高；其次，考虑到求定中点 M 的复杂性，将 M 点用大地线两 端点平均纬度及平均方位角相对 应的 m 点来代替，并借助迭代计 算便可顺利地实现大地主题正解。
f (BM , AM )
F (Bm
BM
Bm
,
Am
AM
Am )
+
dB ( dS )M
f f (Bm , Am ) ( Bm )(BM
f Bm ) ( Am )( AM
Am )
+
dB
dB
dB ( dS )M
( )
f (Bm , Am ) (
dS Bm
S sin Am 24 N m2
[S 2tm2
sin 2
Am
S 2 cos2 Am (1m2 9m2tm2 )]
大地主题解算
大地主题解算一、实验目的：1.提高运用计算机语言编程开发的能力；2.加深对大地主题解算计算公式及辅助参数的理解并掌握计算步骤；3.通过编程语言实现大地主题解算。
二、工具：Windows XP Mode 环境下的Microsoft Visual C++ 6.0三、注意事项：1.计算所需变量多，容易混淆；2.正反算函数的编写；3.函数调用；4.弧度与角度之间的转化。
四、实验要求：1.提交报告，实验总结，编写代码；2.独立编程，调试运行；3.上交成果：编写思想，编写过程，问题分析，源代码，计算结果；五、编程过程实现：1.对白塞尔法大地主题解算有一定的了解，并参考教材P148-P150；2.由于参数较多，而在C语言环境下很多符号无法定义，需要符合要求的定义符号替代书本上那些无法直接在C语言环境下定义的符号来达到实现实验的目的；3.程序中采用弧度与度分秒之间转换的函数定义与调用，减轻一定的实验麻烦；4.在C语言环境下，数学函数fabs代替abs起绝对值作用，atan代替arctan起反函数作用；5.程序中尤其注意弧度与角度之间转换，在C语言环境下电脑默认为弧度。
六、源程序代码：#include<stdio.h>#include<math.h>double hudu(double,double,double); /*度分秒转换为弧度*/double du(double); /*弧度转换为度*/double fen(double); /*弧度转换为分*/double miao(double); /*弧度转换为秒*/#define PI 3.1415926void main (void){int k;printf("请选择执行正算或者反算，若执行正算，请输入1；若执行反算，请输入2。
\n");scanf("%d",&k);/*正算*/if(k==1){double bz,lz,az,S,bz2,lz2,az2,B1,L1,A1,B2,L2,A2,bx,by,lx,ly,ax,ay;int bx2,by2,lx2,ly2,ax2,ay2;doublee2,W1,sinu1,cosu1,sinA0,coto1,sin2o1,cos2o1,sin2o,cos2o,A,B,C,r,t,o0,o,g,si nu2,q;/*以度分秒顺序输入数据*/printf("请输入大地线起点纬度度分秒\n");scanf("%lf%lf%lf",&bx,&by,&bz);printf("请输入大地线起点经度度分秒\n");scanf("%lf%lf%lf",&lx,&ly,&lz);printf("请输入大地方位角度分秒\n");scanf("%lf%lf%lf",&ax,&ay,&az);printf("请输入大地线长度\n");scanf("%lf",&S);/*调用函数*/B1=hudu(bx,by,bz);L1=hudu(lx,ly,lz);A1=hudu(ax,ay,az);/*白塞尔大地主题解算*/e2=0.006693421622966;W1=sqrt(1-e2*sin(B1)*sin(B1));sinu1=sin(B1)*(sqrt(1-e2))/W1;cosu1=cos(B1)/W1;sinA0=cosu1*sin(A1);coto1=cosu1*cos(A1)/sinu1;sin2o1=2*coto1/(coto1*coto1+1);cos2o1=(coto1*coto1-1)/(coto1*coto1+1);A=6356863.020+(10718.949-13.474*(1-sinA0*sinA0))*(1-sinA0*sinA0);B=(5354.469-8.978*(1-sinA0*sinA0))*(1-sinA0*sinA0);C=(2.238*(1-sinA0*sinA0))*(1-sinA0*sinA0)+0.006;r=691.46768-(0.58143-0.00144*(1-sinA0*sinA0))*(1-sinA0*sinA0);t=(0.2907-0.0010*(1-sinA0*sinA0))*(1-sinA0*sinA0);o0=(S-(B+C*cos2o1)*sin2o1)/A;sin2o=sin2o1*cos(2*o0)+cos2o1*sin(2*o0);cos2o=cos2o1*cos(2*o0)-sin2o1*sin(2*o0);o=o0+(B+5*C*cos2o)*sin2o/A;g=(r*o+t*(sin2o-sin2o1))*sinA0;/*求B2*/sinu2=sinu1*cos(o)+cosu1*cos(A1)*sin(o);B2=atan(sinu2/(sqrt(1-e2)*sqrt(1-sinu2*sinu2)));/*求L2*/q=atan(sin(A1)*sin(o)/(cosu1*cos(o)-sinu1*sin(o)*cos(A1)));/*判断q*/if(sin(A1)>0 && tan(q)>0)q=fabs(q);else if(sin(A1)>0 && tan(q)<0)q=PI-fabs(q);else if(sin(A1)<0 && tan(q)<0)q=-fabs(q);elseq=fabs(q)-PI;L2=L1+q-g/3600/180*PI;/*求A2*/A2=atan(cosu1*sin(A1)/(cosu1*cos(o)*cos(A1)-sinu1*sin(o))); /*判断A2*/if(sin(A1)<0 && tan(A2)>0)A2=fabs(A2);else if(sin(A1)<0 && tan(A2)<0)A2=PI-fabs(A2);else if(sin(A1)>0 && tan(A2)>0)A2=PI+fabs(A2);elseA2=2*PI-fabs(A2);/*调用函数*/bx2=(int)(du(B2));by2=(int)(fen(B2));bz2=miao(B2);lx2=(int)(du(L2));ly2=(int)(fen(L2));lz2=miao(L2);ax2=(int)(du(A2));ay2=(int)(fen(A2));az2=miao(A2);printf("大地线终点纬度度分秒分别为：\n%d\n%d\n%lf\n",bx2,by2,bz2);printf("大地线终点经度度分秒分别为：\n%d\n%d\n%lf\n",lx2,ly2,lz2);printf("终点大地方位角度分秒分别为：\n%d\n%d\n%lf\n",ax2,ay2,az2);}/*反算*/else if(k==2){doublebz,lz,bz2,lz2,az,az2,B1,L1,B2,L2,S,A1,A2,bx,by,lx,ly,bx2,by2,lx2,ly2;int ax,ay,ax2,ay2;doublee2,W1,W2,sinu1,sinu2,cosu1,cosu2,L,a1,a2,b1,b2,g,g2,g0,r,p,q,sino,coso,o,si nA0,x,t1,t2,A,B,C,y;/*以度分秒顺序输入数据*/printf("请输入大地线起点纬度度分秒\n");scanf("%lf%lf%lf",&bx,&by,&bz);printf("请输入大地线起点经度度分秒\n"); scanf("%lf%lf%lf",&lx,&ly,&lz);printf("请输入大地线终点纬度度分秒\n"); scanf("%lf%lf%lf",&bx2,&by2,&bz2); printf("请输入大地线终点经度度分秒\n"); scanf("%lf%lf%lf",&lx2,&ly2,&lz2);/*调用函数*/B1=hudu(bx,by,bz);L1=hudu(lx,ly,lz);B2=hudu(bx2,by2,bz2);L2=hudu(lx2,ly2,lz2);/*白塞尔大地主题解算*/e2=0.006693421622966;W1=sqrt(1-e2*sin(B1)*sin(B1));W2=sqrt(1-e2*sin(B2)*sin(B2));sinu1=sin(B1)*sqrt(1-e2)/W1;sinu2=sin(B2)*sqrt(1-e2)/W2;cosu1=cos(B1)/W1;cosu2=cos(B2)/W2;L=L2-L1;a1=sinu1*sinu2;a2=cosu1*cosu2;b1=cosu1*sinu2;b2=sinu1*cosu2;/*逐次趋近法求解A1*/g0=-662.904266/3600*PI/180; g=0;r=L;while(1){p=cosu2*sin(r);q=b1-b2*cos(r);A1=atan(p/q);/*判断A1*/if(p>0 && q>0)A1=fabs(A1);else if(p>0 && q<0)A1=PI-fabs(A1);else if(p<0 && q<0)A1=PI+fabs(A1);elseA1=2*PI-fabs(A1);sino=p*sin(A1)+q*cos(A1);coso=a1+a2*cos(r);o=atan(sino/coso);/*判断o*/if(coso>0)o=fabs(o);elseo=PI-fabs(o);sinA0=cosu1*sin(A1);x=2*a1-(1-sinA0*sinA0)*coso;t1=(33523299-(28189-70*(1-sinA0*sinA0))*(1-sinA0*sinA0))*0.0000000001;t2=(28189-94*(1-sinA0*sinA0))*0.0000000001;g2=(t1*o-t2*x*sino)*sinA0;/*检验循环次数*/printf("\ng2=%lf\ng0=%lf\n",g2,g0);if(g2<=g0)break;elser=L+g2;}/*求解S*/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)*coso; S=A*o+(B*x+C*y)*sino;/*求解A2*/A2=atan(cosu1*sin(r)/(b1*cos(r)-b2));/*判断A2*/if(p<0 && q<0)A2=fabs(A2);else if(p<0 && q>0)A2=PI-fabs(A2);else if(p>0 && q>0)A2=PI+fabs(A2);elseA2=2*PI-fabs(A2);/*调用函数*/ax=(int)(du(A1));ay=(int)(fen(A1));az=miao(A1);ax2=(int)(du(A2));ay2=(int)(fen(A2));az2=miao(A2);printf("起点大地方位角度分秒分别为：\n%d\n%d\n%lf\n",ax,ay,az);printf("终点大地方位角度分秒分别为：\n%d\n%d\n%lf\n",ax2,ay2,az2);printf("大地线长度为：%lf\n",S);}/*数据错误*/elseprintf("数据错误，请重新输入\n");}/*度分秒转换为弧度*/double hudu(double a0,double b0,double c0){double A0;A0=(a0+b0/60+c0/3600)*PI/180;return A0;}/*弧度转换为度*/double du(double B0){double x0;x0=(int)(B0*180/PI);return x0;}/*弧度转换为分*/double fen(double C0){double _y,y0;_y=(int)(C0*180/PI);y0=(fabs)((int)((C0*180/PI-_y)*60));return y0;}/*弧度转换为秒*/double miao(double D0){double _z1,_z2,z0;_z1=(int)(D0*180/PI);_z2=(int)((D0*180/PI-_z1)*60);z0=(fabs)((double)(((D0*180/PI-_z1)*60-_z2)*60));return z0;}大地主题解算正算：大地主题解算反算：以上检验数据来自书本P151，P152。
大地主题解算深度干货超精
大地主题解算(正算)代码与白塞尔大地主题解算大地主题解算(正算)代码：根据经纬度和方向角以及距离计算另外一点坐标新建模块->拷贝下面的大地主题(正算)代码，调用方法示例：起点经度：116.235(度)终点纬度：37.435(度)方向角：50（度）长度：500（米）终点经纬度（"经度,纬度"）=Computation(37.435,116.235,50,500)Const Pi = 3.14Private a, b, c, alpha, e, e2, W, V As Double'a 长轴半径'b 短轴'c 极曲率半径'alpha 扁率'e 第一偏心率'e2 第二偏心率'W 第一基本纬度函数'V 第二基本纬度函数Private B1, L1, B2, L2 As Double'B1 点1的纬度'L1 点1的经度'B2 点1的纬度'L2 点2的经度Private S As Double '''''大地线长度Private A1, A2 As Double'A1 点1到点2的方位角'A2 点2到点1的方位角Function Computation(STARTLAT, STARTLONG, ANGLE1, DISTANCE As Double) As StringB1 = STARTLATL1 = STARTLONGA1 = ANGLE1S = DISTANCEa = 6378245b = 6356752.3142c = a ^ 2 / balpha = (a - b) / ae = Sqr(a ^ 2 - b ^ 2) / ae2 = Sqr(a ^ 2 - b ^ 2) / bB1 = rad(B1)L1 = rad(L1)A1 = rad(A1)W = Sqr(1 - e ^ 2 * (Sin(B1) ^ 2))V = W * (a / b)Dim W1 As DoubleE1 = e ''''第一偏心率'// 计算起点的归化纬度W1 = W ''Sqr(1 - e1 * e1 * Sin(B1 ) * Sin(B1 )) sinu1 = Sin(B1) * Sqr(1 - E1 * E1) / W1cosu1 = Cos(B1) / W1'// 计算辅助函数值sinA0 = cosu1 * Sin(A1)cotq1 = cosu1 * Cos(A1)sin2q1 = 2 * cotq1 / (cotq1 ^ 2 + 1)cos2q1 = (cotq1 ^ 2 - 1) / (cotq1 ^ 2 + 1)'// 计算系数AA,BB,CC及AAlpha, BBeta的值。
大地主题解算方法综述
测绘科学 Sc ience o f Survey ing and M app ing
V o l132 N o14 Ju l1
大地主题解算方法综述
周振宇, 郭广礼, 贾新果
(中国矿业大学 环测学院, 江苏徐 州 221008)
【摘 要】大地主题解算是大地测量中的 重要问 题, 由于 椭球的 复杂性, 随 之产生 的解算 方法也 是多种 多样。本
和方位角差。
像这类方法 还有 很多, 就 不再一 一赘 述, 其 基本 原理
大致相同。
21 2 以白塞尔 大地投影为基础
可以近似认 为地球 椭球 的形状 与圆 球区 别不大。 在椭
球面上解算大地主 题问题 时可 借助于 球面 三角 学公式 简单
而严密的进 行。因此, 如 将椭 球面上 的大 地线长 度投 影到
在于: 解算精 度与距 离有 关, 距离 越长, 收敛 越慢, 因此
只适用于较短的距离。当公式 为四次项时 , 在中纬度 地区,
公式可用于 200km 以下的 大地主 题解算, 经纬 度计算 精度
可达 01 0001", 方位 角计 算精度 达 01001", 计算 边长 可精
确至厘 米。是 短 距离 解 算 的理 想 公式。 当距 离 小于 70km
椭球的过渡。
白塞尔 ( Besse l) 首先提出并解决了投影条件, 使这一解
法得以实现。
这类公式的特点是计算 公式展 开 e2 或 的幂级 数, 解算
精度与距离 长短无 关。因此 它既 适用于 短距 离解 算, 也适 用于长距离 解算 [ 5] 。其主 要缺 点在 于: 由 S 求 R、由 L 求 K, 或相反的运算, 需要进行迭代。同时还要预先 计算辅助
大地主题的数值解法
大地主题的数值解法范业明1,王解先1,2,刘慧芹1(1 同济大学测量与国土信息工程系,上海 200092;2 现代工程测量国家测绘局重点实验室,上海 200092)摘要:本文基于椭球面上大地线的微分方程,将法截弧方位角与大地线方位角之间的关系作为初始条件,通过用数值方法求解大地线的微分方程,进行大地主题的正反解。
并以实际数据验证了其正确性与可行性。
本法均采用封闭公式计算,精度高,公式简单,特别适用于计算机解算。
关键词:大地线;微分方程;数值解法中图分类号:P226+1文献标识码:BAbstract :Based on differential equation of geodesic on the surface of ellipsoid,the problem of direct and inverse solution of geodetic is solved through numerical method.The relationship of normal arc and geodesic is deduced and introduced as the initial condition for the differential equation.Numerical experiments are given,and the validity and feasibility of the method proposed in this paper have been proved.Besides,all the formulas are close,simple and easy to be realized by computer.Key words :geodesic;differential equation;numerical method 收稿日期:2006 02 15;修订日期:2006 05 10作者简介:范业明(1982-),男(汉族),辽宁沈阳人,硕士研究生.0 引言一直以来由于计算工具的限制,大地主题解算一般都采用级数展开形式,如短距离的高斯平均引数法,长距离的贝塞尔-赫尔默特方法等等。
Bessel大地主题解算程序
// 计 算 终 点 大 地 坐 标 及 方 位 角
B2,L2,A2
sinu2=sinu1*cos(sigma)+cosu1*cos(A1)*sin(sigma);
B2=atan(sinu2/sqrt((1-e2)*(1-sinu2*sinu2)));
//计算 B2
lambda=atan(sin(A1)*sin(sigma)/(cosu1*cos(sigma)-sinu1*sin(sigma)*cos(A1))); //求 λ
大地主题解算的意义bessel函数bessel卫星坐标解算程序基线解算besseljmatlabbesselgps基线解算原理gps基线解算南方gps解算软件
#include<stdio.h> #include<math.h> #include<stdlib.h>
double trans1() { double B1,B11,B12,B13,B111;
//计算归化纬度
double sinA0,cotsigma1,sin2sigma1,cos2sigma1; //计算辅助三角函数值 sinA0=cosu1*sin(A1); cotsigma1=cosu1*cos(A1)/sinu1; sin2sigma1=2*cotsigma1/(1+cotsigma1*cotsigma1); cos2sigma1=(cotsigma1*cotsigma1-1)/(cotsigma1*cotsigma1+1);
//计算 A,B,C 以及 α，β 的值 //P144(4-265)
//P146(4-284)
（好像有问
double sigma0,sin2sigma1sigma0,cos2sigma1sigma0,sigma; //计算球面长度 σ sigma0=(S-(B+C*cos2sigma1)*sin2sigma1)/A; sin2sigma1sigma0=sin2sigma1*cos(2*sigma0)+cos2sigma1*sin(2*sigma0); cos2sigma1sigma0=cos2sigma1*cos(2*sigma0)-sin2sigma1*sin(2*sigma0); sigma=sigma0+(B+5*C*cos2sigma1sigma0)*sin2sigma1sigma0/A;
贝塞尔大地坐标解算
贝塞尔大地坐标解算（正算)#include<stdio.h>#include<math.h>#define PI 3.1415926535897323int main(){int q,p;double h,i,j,L,B,A1,S;double a,g,c,d,e,f;void n(double b,double l,double a1,double s,double a,double d,double e,double f);double z(double d,double m,double s);void r(double dd);printf("请选择：\n1.我国1954北京坐标系\n2.我国1980西安坐标系\n3.我国2000国家大地坐标系\n");scanf("%d",&p);if(p==1){a=6378245.0;g=6356863.0188;c=6399698.9018;d=0.00335232986926;e=0.00669342162297;f=0.00673852541468;}else if(p==2){a=6378140.0;g=6356755.2882;c=6399596.6520;d=0.0033528131779;e=0.00669438499959;f=0.00673950181947;}else if(p==3){a=6378137.0;g=6356752.3141;c=6399593.6259;d=1/298.257222101;e=0.0066943800229;f=0.00673949677547;}printf("请输入注：输入格式为角度(例如：30'40'50')\nL1=");scanf("%lf'%lf'%lf'",&h,&i,&j);L=z(h,i,j);printf("请输入注：输入格式为角度(例如：30'40'50')\nB1=");scanf("%lf'%lf'%lf'",&h,&i,&j);B=z(h,i,j);printf("请输入注：输入格式为角度(例如：30'40'50')\nA1=");scanf("%lf'%lf'%lf'",&h,&i,&j);A1=z(h,i,j);printf("请输入\nS=");scanf("%lf",&S);n(B,L,A1,S,a,d,e,f);scanf("%d",&q);}double z(double d,double m,double s){double e;double sign=(d<0.0)?-1.0:1.0;if(d==0){sign=(m<0.0)?-1.0:1.0;if(m==0){sign=(s<0.0)?-1.0:1.0;}}if(d<0)d=d*(-1.0);if(m<0)m=m*(-1.0);if(s<0)s=s*(-1.0);e=sign*(d*3600+m*60+s)*PI/(3600*180);return e;}void r(double dd){int de;int d,m;double s;double sign=(dd<0.0)?-1.0:1.0;if(dd<0)dd=fabs(dd);dd=dd*3600*180/PI;de=int(dd/3600);d=sign*de;dd=dd-de*3600;m=int(dd/60);s=dd-m*60;printf("%d'%d'%lf'\n",d,m,s);}void n(double b,double l,double a1,double s,double a,double d,double e,double f) {double u1,m,M,A2,u2,ss,ss1,k,aa,bb,cc,rr,rr1,rr2,B2,L2;double aa2,bb2,cc2,k2,ll;u1=atan(sqrt(1-e)*tan(b));m=asin(cos(u1)*sin(a1));if(m<0.0)m=m+2*PI;M=atan(tan(u1)/cos(a1));if(M<0.0)M=M+PI;k=f*cos(m)*cos(m);aa=sqrt(1+f)*(1-k/4+7*k*k/64-15*k*k*k/256)/a;bb=k/4-k*k/8+37*k*k*k/512;cc=k*k/128-k*k*k/128;ss1=aa*s;ss=aa*s+bb*sin(ss1)*cos(2*M+ss1)+cc*sin(2*ss1)*cos(4*M+2*ss1);for(;fabs(ss-ss1)>(2.8*PI/180*pow(10.0,-7.0));){ss1=ss;ss=aa*s+bb*sin(ss1)*cos(2*M+ss1)+cc*sin(2*ss1)*cos(4*M+2*ss1);}A2=atan(tan(m)/cos(M+ss));if(A2<0.0)A2=A2+PI;if(a1<=PI)A2=A2+PI;u2=atan(-cos(A2)*tan(M+ss));rr1=atan(sin(u1)*tan(a1));if(rr1<0)rr1=rr1+PI;if(m>=PI)rr1=rr1+PI;rr2=atan(sin(u2)*tan(A2));if(rr1<0)rr2=rr2+PI;if(m<PI){if((M+ss)>=PI)rr2=rr2+PI;}elseif((M+ss)<=PI)rr2=rr2+PI;rr=rr2-rr1;B2=atan(sqrt(1+f)*tan(u2));k2=e*cos(m)*cos(m);aa2=(e/2+e*e/8+e*e*e/16)-e*(1+e)*k2/16+3*e*k2*k2/128;bb2=e*(1+e)*k2/16-e*k2*k2/32;cc2=e*k2*k2/256;ll=rr-sin(m)*(aa2*ss+bb2*sin(ss)*cos(2*M+ss)+cc2*sin(2*ss)*cos(4*M+2*ss));L2=l+ll;printf("L2=");r(L2);printf("\nB2=");r(B2);printf("\nA2=");r(A2);}反算#include<stdio.h>#include<math.h>#define PI 3.1415926535897323int main(){int q,p;double h,i,j,L,B,L2,B2;double a,g,c,d,e,f;void n(double b,double l,double l2,double b2,double a,double d,double e,double f);double z(double d,double m,double s);void r(double dd);printf("请选择：\n1.我国1954北京坐标系\n2.我国1980西安坐标系\n3.我国2000国家大地坐标系\n");scanf("%d",&p);if(p==1){a=6378245.0;g=6356863.0188;c=6399698.9018;d=0.00335232986926;e=0.00669342162297;f=0.00673852541468;}else if(p==2){a=6378140.0;g=6356755.2882;c=6399596.6520;d=0.0033528131779;e=0.00669438499959;f=0.00673950181947;}else if(p==3){a=6378137.0;g=6356752.3141;c=6399593.6259;d=1/298.257222101;e=0.0066943800229;f=0.00673949677547;}printf("请输入注：输入格式为角度(例如：30'40'50')\nL1=");scanf("%lf'%lf'%lf'",&h,&i,&j);L=z(h,i,j);printf("请输入注：输入格式为角度(例如：30'40'50')\nB1=");scanf("%lf'%lf'%lf'",&h,&i,&j);B=z(h,i,j);printf("请输入注：输入格式为角度(例如：30'40'50')\nL2=");scanf("%lf'%lf'%lf'",&h,&i,&j);L2=z(h,i,j);printf("请输入注：输入格式为角度(例如：30'40'50')\nB2=");scanf("%lf'%lf'%lf'",&h,&i,&j);B2=z(h,i,j);n(B,L,L2,B2,a,d,e,f);scanf("%s",&q);}double z(double d,double m,double s){double e;double sign=(d<0.0)?-1.0:1.0;if(d==0){sign=(m<0.0)?-1.0:1.0;if(m==0){sign=(s<0.0)?-1.0:1.0;}}if(d<0)d=d*(-1.0);if(m<0)m=m*(-1.0);if(s<0)s=s*(-1.0);e=sign*(d*3600+m*60+s)*PI/(3600*180);return e;}void r(double dd){int de;int d,m;double s;double sign=(dd<0.0)?-1.0:1.0;if(dd<0)dd=fabs(dd);dd=dd*3600*180/PI;de=int(dd/3600);d=sign*de;dd=dd-de*3600;m=int(dd/60);s=dd-m*60;printf("%d'%d'%lf'\n",d,m,s);}void n(double b,double l,double l2,double b2,double a,double d,double e,double f) {double u1,m,M,A2,u2,ss,ss1,k,aa,bb,cc,rr,rr0,B2,L2,m0,A10,A1,S;double aa2,bb2,cc2,k2,ll;u1=atan(sqrt(1-e)*tan(b));u2=atan(sqrt(1-e)*tan(b2));ss=acos(sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(l2-l));m0=asin(cos(u1)*cos(u2)*sin(l2-l)/sin(ss));rr0=l2-l+0.003351*ss*sin(m0);k2=e*cos(m0)*cos(m0);aa2=(e/2+e*e/8+e*e*e/16)-e*(1+e)*k2/16+3*e*k2*k2/128;ss1=ss+sin(m0)*aa2*ss*sin(m0);m=asin(cos(u1)*cos(u2)*sin(rr0)/sin(ss1));A10=atan(sin(rr0)/(cos(u1)*tan(u2)-sin(u1)*cos(rr0)));if(A10<=0.0)A10=A10+PI;if(m<=0.0)A10=A10+PI;M=atan(sin(u1)*tan(A10)/sin(m));if(M<=0.0)M=M+PI;k2=e*cos(m)*cos(m);aa2=(e/2+e*e/8+e*e*e/16)-e*(1+e)*k2/16+3*e*k2*k2/128;bb2=e*(1+e)*k2/16-e*k2*k2/32;rr=l2-l+sin(m)*(aa2*ss1+bb2*sin(ss1)*cos(2*M+ss1));ss=acos(sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(rr));A1=atan(sin(rr)/(cos(u1)*tan(u2)-sin(u1)*cos(rr)));if(A1<=0.0)A1=A1+PI;if(m<=0.0)A1=A1+PI;A2=atan(sin(rr)/(-cos(u2)*tan(u1)+sin(u2)*cos(rr)));if(A2<=0.0)A2=A2+PI;if(m>=0.0)A2=A2+PI;k=f*cos(m)*cos(m);aa=sqrt(1+f)*(1-k/4+7*k*k/64-15*k*k*k/256)/a;bb=k/4-k*k/8+37*k*k*k/512;cc=k*k/128-k*k*k/128;S=(ss-bb*sin(ss)*cos(2*M+ss)-cc*sin(2*ss)*cos(4*M+2*ss))/aa;printf("A1=");r(A1);printf("A2=");r(A2);printf("S=%lf",S);}。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
3.计算分析
1.B=0，大地线长误差
大家有疑问的，可以询问和交流
可以互相讨论下，但要小声点
3.计算分析
1.B=30，大地线长误差
3.计算分析
1.B=70，大地线长误差
3.计算分析
1.B=90，大地线长误差
3.计算分析
200km时大地方位角误差
3.计算分析
800km时大地方位角误差
3.计算分析
也适用于长距离解算。可适应10000km或更长的距离，这对 于国际联测，精密导航，远程导弹发射等都具有重要意义。
2.问题分析
1.球面三角函数计算
球面大地主题解算不同球面角距计算误差 检验方法：先反算再正算
2.问题分析
1.球面三角函数计算
球面大地主题解算不同方位角计算误差
2.问题分析
q接近0赋值1Biblioteka -20贝塞尔大地主题解算分析
目录
• 1.基本原理 • 2.问题分析 • 3.数值计算 • 4.结论及问题
1.基本原理
大地主题解算
• 主题解算分为: • 短距离(<400km) • 中距离(<1000km) • 长距离(1000km以上)
1.基本原理
1.白塞尔大地主题解算
•以白塞尔大地投影为基础： •1)按椭球面上的已知值计算球面相应值； •2)在球面上解算大地主题问题； •3)按球面上得到的数值计算椭球面上的相应数值。 • 特点：解算精度与距离长短无关，它既适用于短距离解算，