贝塞尔大地主题解算分析 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);}。
贝塞尔大地主题正反算及其编程
存档日期：存档编号：江苏师范大学科文学院本科生毕业设计（论文）论文题目：贝塞尔大地主题正反算及程序设计*名：**系别：环境与测绘系专业：测绘工程年级、学号：08测绘、************师：***江苏师范大学科文学院教务部印制摘要在大地测量计算过程中，大地主题解算计算繁琐复杂，手工计算易于出错，而且费时费力。
随着计算机技术的高速发展，计算机计算的速度快、准确度高、计算机语言的丰富、编程可视化等优点为我们将复杂烦琐的计算过程简单、简洁、高效化带来了契机。
为了便于工程计算，本课题着眼于研究借助计算机及其编程语言MATLAB来实现大地主题解算问题。
大地主题解算方法，主要有高斯平均引数法、勒让德级数法、贝塞尔法。
前两种方法受到大地线长度的制约，随着大地线两端点的距离加大，其解算精度明显降低。
而贝塞尔法具有不受大地线长度制约的优点，解算精度最大不超过5毫米，是大地主题解算方法中解算精度最高的一种。
因此，本文就以贝塞尔法为研究对象，开发贝塞尔大地主题解算小程序。
关键词：贝塞尔大地主题正反算，程序设计A b s t r a c tIn Geodetic computation process, the solution of geodetic problem computational complexity of manual calculation, error prone, and took the time and trouble. With the rapid development of computer technology, computational speed, high accuracy, computer language, the advantages of rich programming visualization for we will complex complicated calculating process is simple, concise, efficient change brings opportunity. For the convenience of engineering calculation, this paper focus on the research of have the aid of computer and programming language MATLAB to realize the geodetic problem solving.Solution of geodetic problem method, mainly Gauss average argument method, Legendre series expansion method, Bessel method. The former two methods by geodesic length constraints, along with the line ends point distance increase, the calculation precision significantly reduced. Bessel law is not affected by the advantages of geodesic length restriction, calculation accuracy of less than 5mm, is the theme of the earth solution method of calculating precision is highest kind. Therefore, this article on Bessel law as the object of study, the development of Bessel solution of geodetic problem of small procedures.Key words：Direct and inverse solution of geodetic problem，The designing of program目录摘要 .................................................... I Abstract (II)1. 椭球面和球面上对应元素间的关系 (1)1.1 贝塞尔法解算大地问题的基本思想 (1)1.2 对应元素关系式 (1)2. 在球面上进行大地主题解算 (5)2.1 球面上大地主题正解方法 (6)2.2 球面上大地主题反解方法 (6)3. 贝塞尔微分方程的积分 (8)3.1 用于大地主题反算时的大地线长度公式 (8)3.2 用于大地主题正算时的大地线长度公式 (10)3.3 椭球面大地线端点经差与球面经差的关系式 (11)3.4 反解时，大地线长度和球面长度关系式的简化 (13)4. 贝塞尔大地主题正解算步骤 (15)4.1 计算起点的归化纬度 (15)4.2 计算辅助函数值 (15)4.3 计算系数 (15)4.4 计算球面长度 (15)4.5 计算经差改正数 (16)4.6 计算终点大地坐标及大地方位角 (16)5. 贝塞尔大地主题正解算MATLAB 程序设计 (17)5.1 正算流程 (17)5.2 界面设计及功能模块编写 (18)6. 贝赛尔大地主题反解算步骤 (23)6.1 辅助计算 (23)6.2 用逐次趋近法同时计算起点大地方位角、球面长度及经差δλ+=l ：236.3 计算系数及大地线长度S (24)6.4 计算反方位角 (24)7. 贝塞尔大地主题反解算MATLAB 程序设计 (25)7.1 反算流程 (25)7.2 界面设计及功能模块编写 (26)8 总结 (29)参考文献 (31)致谢 (29)1. 椭球面和球面上对应元素间的关系1.1 贝塞尔法解算大地问题的基本思想基本思想：将椭球面上的大地元素按照贝塞尔投影条件投影到辅助球面上，继而在球面上进行大地问题解算，最后再将球面上的计算结果换算到椭球面上。
贝塞尔函数3幻灯片
o
246
-0.5
8 10 12
9
以
(n) m
(m1,2,L)表示
J
n
(
x
)
的非负零点,
则
lim
m
(n)
(n)
m1 m
.
1.0 J 0 ( x )
0.5
J1( x )
函数以为周期振荡
o
2 4 6 8 10 12
-0.5
10
方程 Jn R 0 的解为：
R m n , m 1 ,2 ,L
由 条 件 （ 8 ） 知 D 0.
29
二、求本征值、本征函数
再 由 条 件 （ 9） 得 ,
R(b)CJ0( b)0
即 ， J 0 ( b ) 0， 由 此 可 知b是 J 0 (x )的 零 点 。
以 (0 ) m
表 示 J 0 (x )的 正 零 点 ， 有
J0(m(0)) 0
从 而 ， 得 到 方 程 （ 7 ） 在 条 件 （ 8 ） 、 （ 9 ） 下 的
由 条 件 (4） ， 得
z 0
z h
u(,0)
m 1
(C mD m)J0(b m (0))0
于是得
C m D m 0 ( m 1 , 2 ,L )( 1 1 )
再 由 条 件 （ 5） 得
u b 0 (5)
(0)
mh
(0)
mh
(0)
u(,h) (C me m 1
D me )J0(b m
的 通 解 为 P ( r ) A J n (
r ) B Y n (
r ) 26
一、建立方程 方 程 （ 7 ） 为 零 阶 贝 塞 尔 方 程 ， 其 通 解 为
剖析贝塞尔曲线.PPT
■ 显示器分辨率（Screen Resolution）
■ 概念：1、显卡在显示器上显示的像素数量。 2、显示器中显示的分辨率，
■ 格式：1、宽度×高度（800 ×600） 2、此值固定为72dpi（PC）96dpi（Mac）
7
相关概念
打印机分辨率（Print Resolution）
■ 适用于字体、商标、漫画的设计
■ 矢量图像主要适用于以曲线为主的艺术字，商标等制 作因为字体本身就是一种特殊矢量对象
■ 文件占用空间小
■ 矢量图像文件只需要记录曲线数学函数关系，数据量 比逐点记录色彩的位图图像小很多
15
常用矢量图绘制软件
■ Adobe Illustrator ■ Corel Draw ■ Macromedia Freehand ■ CreatureHouse Expression 对于矢量图像的创作我们称为Drawing，
■ 概念：打印设备产生物理点的频率，又称为输 出分辨率（Export Resolution）。用来设置 Halftone Screen（网线板，网屏）的稠密程度。
■ 单位：dpi（dot per inch 点每英寸）300、600、 1200
■ 颜色深度（Color Depth）
■ 概念：图像中每一个像素所能显示的颜色数目。 也叫位分辨率（Bit Resolution)
它包含了用颜料给图形上色的意思。
11
由贝塞尔曲线构成的矢量图像
■ Illustrator是利用贝塞尔曲线的图形软件， 只有熟练掌握贝塞尔曲线才能轻松对 Illustrator经行操作。
■ 贝塞尔曲线的特征是，用四个点控制一条 曲线的位置和曲率
贝塞尔大地主题解算分析
3.计算分析
1.B=0，大地线长误差
3.计算分析
1.B=30，大地线长误差
3.计算分析
1.B=70，大地线长误差
3.计算分析
1.B=90，大地线长误差
3.计算分析
200km时大地方位角误差
3.计算分析
800km时大地方位角误差
3.计算分析
10000km时大地方位角误差
结论及问题
• ቤተ መጻሕፍቲ ባይዱ离精度优于2cm • 角度精度优于2ms • 距离精度有提高空间
也适用于长距离解算。可适应10000km或更长的距离，这对 于国际联测，精密导航，远程导弹发射等都具有重要意义。
2.问题分析
1.球面三角函数计算
球面大地主题解算不同球面角距计算误差 检验方法：先反算再正算
2.问题分析
1.球面三角函数计算
球面大地主题解算不同方位角计算误差
2.问题分析
q接近0赋值1d-20
目录
• 1.基本原理 • 2.问题分析 • 3.数值计算 • 4.结论及问题
1.基本原理
大地主题解算
• 主题解算分为: • 短距离(<400km) • 中距离(<1000km) • 长距离(1000km以上)
1.基本原理
1.白塞尔大地主题解算
•以白塞尔大地投影为基础： •1)按椭球面上的已知值计算球面相应值； •2)在球面上解算大地主题问题； •3)按球面上得到的数值计算椭球面上的相应数值。 • 特点：解算精度与距离长短无关，它既适用于短距离解算，
大地测量主题解算
一般情况下主项趋近3次，改正项趋近2次 就可满足要求。
四、高斯平均引数反算公式
高斯平均引数反算公式可以依正算公式导出:
S
sin
Am
L" "
Nm
cos Bm
S sin Am
24
N
2 m
[S 2tm2
sin 2
Am
S 2 cos2 Am (1 m2 9m2tm2 )]
t
2 m
2m2 )]}
5次
B2 B1 B, L2 L1 L, A21 A12 A 180
Am
1 2
( A12
A21 )
A12
1 2
A
Bm
1 2
(B2
B1)
B1
1 2
B
注意：
从公式可知，欲求ΔＬ，ΔＢ及ΔＡ，必先有Bm及 Am。但由于Ｂ2和Ａ21未知，故精确值尚不知，为
2f
2
上一讲应掌握的内容
四、整周数N值解算的一般原理
有可变频率法和固定频率法两种
五、全站仪中测距新技术
• 使用高频测距技术
· 温控与动态频率校正技术
• 无棱镜测距技术
· 目标自动识别技术
六、测距的误差分析和精度表达式
D c0 C 4fn
mD2
( c0 )2 m2
4 fn 2
Nm Vm2
高斯平均引数反算公式（续）
S sin Am r01L" r21B"2L" r03L"3 S cos Am s10B" s12B"L"2 s30B"3
大地主题解算(深度干货-超精)
大地主题解算(正算)代码与白塞尔大地主题解算大地主题解算(正算)代码：根据经纬度和方向角以及距离计算另外一点坐标新建模块->拷贝下面的大地主题(正算)代码，调用方法示例：起点经度：116.235(度)终点纬度：37.435(度)方向角：50（度）长度：500（米）终点经纬度（"经度,纬度"）=Computation(37.435,116.235,50,500)Const Pi = 3.1415926535898Private 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的值。
- 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)按球面上得到的数值计算椭球面上的相应数值。 • 特点：解算精度与距离长短无关，它既适用于短距离解算，