白塞尔大地主题解算
用C语言实现大地主题解算
4
6
4
( 争 寺+
‘ ) c o 5 4 A ’
裴连磊 P E I L i a n — l e i
( 新疆 地 矿 局 测 绘 大 队 , 乌鲁木齐 8 3 0 0 1 7 ) ( X i n j i a n g G e o l o g y a n d Mi n e r a l B u r e a u , S u r v e y i n g a n d Ma p p i n g B r i g a d e , U r u mq i 8 3 0 0 1 7 , C h i n a )
1 - 5计 算经 度 差 改 正 数
本 文的程序 实现 是采用 的 白塞 尔大地 主题 解 算的方 法, 根 据 白塞 尔大地主题解 算的 方法 , 分析得 出适 合计 算 机编程 的大地主题正反 算的具体 实现步骤。 1 白塞尔法大地主题正算思想 已知量 : 大 地线 起点 的纬 度 B , 经度 L , 大地 方位 角 A 1 及大地线长度 S 。 求解量 : 大地线终点 的纬度 B : , 经度 L 2
c 番 一 番 + . .
c _ b ( 薷 4 一 6 一)
2
=
c 0 s u 2 =
L a 1 = s i n u l s i n u 2 I a 2 _ c o s u l c 0 s u 2 b=
C O S Ul s i n u 2 , b 2 = s i n ul C O S U 2
之 值 。 即下 式 :
求解 量 : 大地 线长度 S及起、 终点 处的大地 方位角 A
及 A' 。
A : b ( 1 一 鲁十 击k 6 _ -
2 . 1辅 助计算 w :
大地测量学白塞尔大地主题解算实习报告(模板)1
大地测量上机实习报告
一、实习时间:第9周星期三
二、实习地点:系办三楼机房
三、实习内容:白塞尔大地主题解算
四、实习目的:(1)尝试编写白塞尔主题解算的程序
(2)通过实例检验解算的结果
五、实习步骤:
1、先了解书本上关于白塞尔大地主题解算的步骤;
2、双击打开白塞尔大地主题解算的程序,如下:
3、在“选择椭球”中选择“克氏椭球”,“选择正反算”中选择“正算”;
4、然后在输入框中输入相应的数据,如下:
5、接下来单击菜单栏中的“数据正确”,记得正算结果,如下所示:
6、重复第3步,选择反算,并输入相应数据:
7、点击“数据正确”,即得反算结果:
六、实习总结:通过白塞尔大地主题结算方法的上机操作,懂得了白塞尔法结算大地主题是将椭球面上的大地元素按照白塞尔投影条件投影到辅助球面上,继而在球面上进行大地主题解算,再将球面上的计算结果换算到椭球面上的原理,并熟悉白塞尔大地主题解算的程序。
白塞尔大地主题算法
1.球面上大地主题正解方法
此时已知量:φ 1,α1及σ;要求量:φ 2,α2及λ. 首先按(i)式计算 ,继而用下式计算φ 2球面上坐标关系式
如图13—4所示,在椭球 面极三角形PP1P2中.用B. L.S及A分别表示大地线上某 点的大地坐标,大地线长及其 大地方位角。在球面极三角形 P'P1'P2'中,与之相应,用φ , λ,σ及α分别表示球面大圆弧 上相应点的坐标,弧长及方位 角。
M-子午圈曲率半径 N-卯酉圈曲率半径
M=N=1
为简化计算,白塞尔提出如下三个投影条件: 1 . 椭球面大地线投影到球面上为大圆弧; 2.大地线和大圆弧上相应点的方位角相等; 3.球面上任意一点的纬度等于椭球面上相应 点的归化纬度.
u-归化纬度,
r-平行圈半径,a-长半轴
(三)白塞尔微分方程的积分
(一)在球面上进行大地主题解算
如图13—3示,在球面上 有两点P1和P2,其中P1点的大 地纬度φ 1,大地经度λ1,P2点 大地纬度φ 2,大地经度λ2;P1 和P2点间的大圆弧长为σ, P1P2的方位角为α1,其反方位 角为α2,球面上大地主题正算 是已知φ 1,α1,σ,要求φ 2, α2及经差λ;反算问题是已知 φ 1,φ 2及经差λ,要求σ,α1及 α2。
(四)白塞尔法大地主题正算步骤
(五)白塞尔法大地主题反算步骤
白塞尔大地主题解算方法
白塞尔法解算大地主题的基本思想是将椭球面 上的大地元素按照白塞尔投影条件投影到辅助球 面上,继而在球面上进行大地主题解算,最后再 待球面上的计算结果换算到椭球面上。由此可见 ,这种方法的关键问题是找出椭球面上的大地元 素与球面上相应元素之间的关系式。同时也要解 决在球面上进行大地主题解算的方法。
白塞尔大地主题解算的基本思想
白塞尔大地主题解算的基本思想
首先,白塞尔大地主题解算的基本思想之一是建立地壳变形的力学模型。
地壳变形是地球表面的一项重要现象,是由于地质作用和地球内部力
学过程的结果。
地壳变形的力学模型是研究地壳变形的重要方法。
常见的
地壳变形的力学模型有弹性模型、弹塑性模型和拟弹性模型等。
这些模型
可以描述地壳的变形过程,通过对地壳的变形过程进行建模,可以更好地
理解地壳变形的机制和动力学过程。
其次,白塞尔大地主题解算的思想之一是研究现今地球表面的动力学
过程。
地球表面的动力学过程包括板块运动、地震活动、火山喷发等。
这
些过程在地球的长期演化中起着重要的作用。
通过对现今地球表面动力学
过程的研究,可以揭示地球内部的结构和动力学机制,进而更好地理解地
壳变形的成因和发展变化的规律。
第三,白塞尔大地主题解算的基本思想之一是研究地壳变形的成因。
地壳变形的成因是地壳运动的基本原因,也是地球科学研究的一个重要问题。
地壳变形的成因包括构造运动、地壳应力状态的改变、地震活动等。
通过研究地壳变形的成因,可以更好地了解地壳变形的机制和规律,进而
为地震预测和地壳运动的控制提供科学依据。
总之,白塞尔大地主题解算的基本思想是通过建立地壳变形的力学模型,研究现今地球表面的动力学过程,探究地壳变形的成因,揭示地壳变
形与地球动力学过程之间的相互关系,以推动地壳运动的理论和应用研究。
这一思想对于研究地壳运动的机制和规律,了解地壳变形的成因和动力学
过程,提高地震预测和地壳运动控制的水平具有重要意义。
大地主题解算方法综述
测绘科学 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, 或相反的运算, 需要进行迭代。同时还要预先 计算辅助
白塞尔大地主题正反算
地球椭球 CGCS2000 6378137 6356752.314 6399593.626 1/298.257222101 0.00669438 0.006739497 第一偏心率e2= 第一偏心率e‘2= sinu1= cotσ1= A0= B= β= cos2(σ1+σ0)= sinu2= 大地方位角A21= 第一偏心率e2= u1= u2= cos(2σ‘1+Δσ‘)= ω2 = ω5 = s2= s5= 大地方位角A21= u1= u2= cos(2σ‘1+Δσ‘)= ω2 = ω5 = s2= s5=
2
克拉索夫斯基椭球 BJ54 6378245 6356863.019 6399698.902 1/298.3 0.006693422 0.006738525 BJ54 35°00′00.000000″ 35°00′00.000000″ 0.820055451 1.9591E-09 0.002206928 0.241484809 0.002515577 0.002514192 34°59′59.544529″ BJ54 35°00′00.000000″ 35°00′00.000000″ 0.002508527 4.64788E-09 -4.15588E-09 0.003065877 9.19076E-10 16000.00241 34°59′59.544529″ 35°10′30.957713″ 0.002508527 4.64788E-09 -4.15588E-09 0.003065877 9.19076E-10
IUGG--1979椭球 WGS84 6378137 6356752.314 6399593.626 1/298.257223563 0.00669438 0.006739497 6356863.019 89°59′59.999859″ 16000 0.820055451 -1 6360368.854 0.003351407 -0.005031132 6.90795E-06 35°10′30.957713″ 6378245 34°59′59.544529″ 35°10′30.957713″ 0.820055451 0.003358976 9.24261E-07 0.999451116 -0.000548959 90°00′00.000141″ 35°00′00.000000″ 35°00′00.000000″ 0.820055451 0.003358976 9.24261E-07 0.999451116 -0.000548959
大地测量学复习资料#(精选.)
1.垂线偏差:地面一点上的重力向量g和相应椭球面上的法线向量n之间的夹角定义为该点的垂线偏差。
2.参考椭球:具有确定参数(长半径a和扁率α),经过局部定位和定向,同某一地区大地水准面最佳拟合的地球椭球,叫参考椭球。
3.大地线:椭球面上两点间的最短程曲线叫做大地线。
4.力高:水准面在纬度45度处的正常高。
5.大地主题解算:已知某些大地元素推求另一些大地元素的计算工作叫大地主题解算。
6.大地主题正算:已知P1点的大地坐标(L1,B1),P1至P2的大地线长S及其大地方位角,计算P2点的大地坐标(L2,B2)和大地线S在P2点的反方位角A21,这类问题叫做大地主题正算。
7.大地基准:是指能够最佳拟合地球形状的地球椭球的参数及椭球定位和定向8.高斯投影:横轴椭圆柱等角投影(假象有一个椭圆柱横套在地球椭球体外,并与某一条子午线相切,椭球柱的中心轴通过椭球体中心,然后用一定投影方法,将中央子午线两侧各一定范围内的地区投影到椭圆柱上,再将此柱面展开成投影面)。
9.大地测量学:是指在一定的时间与空间参考系中,测量和描绘地球形状及其重力场并监测其变化,为人类活动提供关于地球的空间信息的一门科学。
10.理论闭合差:由水准面不平行而引起的水准环线闭合差,称为理论闭合差。
11.地心坐标系:地心坐标系是在大地体内建立的O-XYZ坐标系。
原点O设在大地体的质量中心,用相互垂直的X,Y,Z三个轴来表示,X轴与首子午面与赤道面的交线重合,向东为正。
Z轴与地球旋转轴重合,向北为正。
Y 轴与XOZ平面垂直构成右手系。
12.高斯投影正、反算公式进行换带计算的步骤。
这种方法的实质是把椭球面上的大地坐标作为过度坐标。
首先把某投影带内有关点的平面坐标(x,y)1利用高斯投影反算公式换算成椭球面上的大地坐标(B,l),进而得到L=L0+l,然后再由大地坐标(B,l),利用投影正算公式换算成相邻带的平面坐标(x,y)2在计算时,要根据第2带的中央子午线来计算经差l,亦即此时l=L-L0。
白塞尔大地主题解算(正算和反算)
大地测量实验报告实验名称:白塞尔大地主题解算(正算和反算)实验目的:1.通过编写白塞尔大地主题电算程序进一步掌握白塞尔法解算大地主题的基本思想。
2.熟练掌握将椭球面上的大地元素按照白塞尔投影条件投影到辅助球面上,继而在球面上进行大地主题解算,最后再将球面上的计算结果换算到椭球面上的基本方法和步骤。
3.学会掌握计算机编程的基本能力。
实验环境:Microsoft Visual C++注意事项:1.在编写程序的过程当中要注意代码的前后统一和重复。
2.注意数值类型的转换和度分秒的换算。
实验步骤:正算:1.计算起点的规划纬度2.计算辅助函数值3.计算系数A,B,C及d,e.4计算球面长度5.计算经度差改正数6.计算终点大地坐标及大地方位角。
反算:1.进行计算辅助函数值2.用逐次趋近法同时计算起点大地方位角、球面长度及经差。
3.计算系数A,B,C及大地线长度S.4.计算反方位角及确定符号。
程序源代码:正算:#include<stdio.h>#include<math.h>#define ee 0.006694384999588#define I 3.141592653double F(double,double,double);void main(void){double A1,B1,L1,S,A2,B2,L2; double x1,x2,x3,y1,y2,y3,z1,z2,z3; double W1,sinu1,sinu2,cosu1,sinA0;doublecota1,cos2a1,sin2a1,cosA0A0;double A,B,C,d,e,a0,a1,m;double n,a,Q,R;printf("请输入数据B1= "); scanf("%lf %lf %lf",&x1,&x2,&x3);B1=F(x1,x2,x3);printf("请输入数据L1= "); scanf("%lf %lf %lf",&y1,&y2,&y3);L1=F(y1,y2,y3);printf("请输入A1= ");scanf("%lf %lf %lf",&z1,&z2,&z3);A1=F(z1,z2,z3);printf("请输入S= ");scanf("%lf",&S);printf("B1=%f\n",B1);printf("L1=%f\n",L1);printf("A1=%f\n",A1);printf("S=%f\n",S);/*计算起点的规划纬度*/W1=sqrt(1-ee*sin(B1)*sin(B1));sinu1=sin(B1)*sqrt(1-ee)/W1;cosu1=cos(B1)/W1;printf("W1=%f\n",W1);printf("sinu1=%f\n",sinu1);printf("cosu1=%f\n",cosu1);/*计算辅助函数值*/sinA0=cosu1*sin(A1);cota1=cosu1*cos(A1)/sinu1;sin2a1=2*cota1/(cota1*cota1+1);cos2a1=(cota1*cota1-1)/(cota1*cota1+1);printf("sinA0=%f\n",sinA0);printf("cota1=%f\n",cota1);printf("sin2a1=%f\n",sin2a1);printf("cos2a1=%f\n",cos2a1);/*计算系数ABC及de*/cosA0A0=1-sinA0*sinA0;A=6356755.288+(10710.341-(13.534*cosA0A0))*cosA0A0;B=(5355.171-9.023*cosA0A0)*cosA0A0;C=(2.256*(cosA0A0))*cosA0A0+0.006;d=691.46768-(0.58143-0.00144*cosA0A0)*cosA0A0;e=(0.2907-cosA0A0*0.0010)*cosA0A0;printf("cosA0A0=%f\n",cosA0A0);printf("A=%f\n",A);printf("B=%f\n",B);printf("C=%f\n",C);printf("d=%f\n",d);printf("e=%f\n",e);/*计算球面长度*/a0=(S-(B+C*cos2a1)*sin2a1)/A;m=sin2a1*cos(2*a0)+cos2a1*sin(2*a0);n=(cos2a1)*(cos(2*a0))-(sin2a1)*(sin(2*a0));a=a0+((B+5*C*n))*m/A;printf("a0=%f\n",a0);printf("m=%f\n",m);printf("n=%f\n",n);printf("a=%f\n",a);/*计算经度差改正数*/Q=(d*a+(e*(m-sin2a1))/3600/180*I)*sinA0;printf("Q=%f\n",Q);/*计算终点大地坐标及大地方位角*/sinu2=sinu1*cos(a)+cosu1*cos(A1)*sin(a);B2=180*atan(sinu2/((sqrt(1-ee))*(sq rt(1-sinu2*sinu2))))/I;R=180*atan(sin(A1)*sin(a)/(cosu1*co s(a)-sinu1*sin(a)*cos(A1)))/I;printf("sinu2=%f\n",sinu2);printf("B2=%f\n",B2);printf("R=%f\n",R*180/I);/*确定R的值*/if(sin(A1)>0 && tan(R)>0)R=abs(R);else if(sin(A1)>0 && tan(R)<0)R=I-abs(R);else if(sin(A1)<0 && tan(R)<0)R=-abs(R);elseR=abs(R)-I;/*确定L2A2的值*/L2=(L1*180/I+R-(Q/206265*180/I));A2=atan(cosu1*sin(A1)/(cosu1*cos(a) *cos(A1)-sinu1*sin(a)));if(sin(A1)<0&&tan(A2)>0)A2=(fabs(A2))*180/I;else if(sin(A1)<0&&tan(A2)<0)A2=(I-fabs(A2))*180/I;else if(sin(A1)>0&&tan(A2)>0)A2=(I+fabs(A2))*180/I;elseA2=(2*I-fabs(A2))*180/I;printf("A2=%3f\n B2=%3f\nL2=%3f\n",A2,B2,L2);}double F(double a2,double b2,doublec2){double d2;d2=(double)(a2+1.0*b2/60+1.0*c2/3600);d2=(d2/180)*I;return (d2);}注:A1,B1,L1,S分别为大地线起点的大地方位角、纬度、经度、大地线长;B2,L2,A2为大地线终点纬度、经度及方位角。
白塞尔大地主题解算
白塞尔大地主题解算方向:学号:姓名:一.基本思路:基本思想:将椭球面上的大地元素按照白塞尔投影条件投影到辅助球面上,继而在球面上进行大地主题解算,最后在将球面上的计算结果换算到椭球面上。
其关键问题是找出椭球面上的大地元素与球面上相应元素之间的关系式,同时解决在球面上进行大地主题解算的方法。
正算流程:1.计算起点的归化纬度2.计算辅助函数值,解球面三角形3.按公式计算相关系数A,B,C 以及α,β4.计算球面长度5.计算纬度差改正数6.计算终点大地坐标及大地方位角011122S B C A{sin (cos )}σσσ=-+10101022222sin ()sin sin cos cos σσσσσσ+=+10101022222cos ()cos cos sin sin σσσσσσ+=-001101522B C A[cos ()]sin ()σσσσσσ=++++010122L A sin [(sin ()sin )]λδασβσσσ-==++-2111u u u A sin sin cos cos cos sin σσ=+2222222222222222222222111111111e u B u B W W u e B u u B B arctan e u e u sin sin cos cos tan tan sin sin tan cos -cos ⎧-==⎪⎪⎨⎪=-⎪⎩⎡⎤⎢⎥==--⎢⎥-⎣⎦1111A arctan u u A sin sin []cos cos sin sin cos σλσσ=-21L L λδ=+-112111u A A arctan u A u cos sin cos cos cos sin sin σσ⎡⎤=⎢⎥-⎣⎦反算流程:1.辅助计算2.用逐次趋近法同时计算起点大地方位角、球面长度及经差,第一次趋近时,取δ=0。
计算下式,重复上述计算过程2.3.计算大地线长度S4.计算反方位角二.已知数据序号B1(DD.MMSS)L1 (DD.MMSS)A12(DD.MMSS)S12(m)1 41.01356874 130.10122676 1.4943 8000L λδ=+211212u pA u u u u qsin cos tan cos sin sin cos cos λλ==-2121p p u q b b A arctanqsin cos cos λλ==-=11p A q A sin cos tan cos σσ+=11p A q A sin sin cos σ=+12a a cos cos σλ=+arctan sin cos σσσ⎛⎫=⎪⎝⎭011A u A sin cos sin =111u A tan tan sec σ=21+σσσ=02122L A sin [(sin sin )]λδασβσσ-==+-L λδ=+11222222S A B C B C sin (cos )sin (cos )σσσσσ=++-+1212u A b b cos sin arctan cos λλ⎡⎤=⎢⎥-⎣⎦三.源代码:#include <stdio.h>#include <math.h>#define e 0.081813334016931499 //克拉索夫斯基椭球体第一偏心率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;double a1,a2,b1,b2,p,q,x,y;doubleW1,W2,sinu1,sinu2,cosu1,cosu2,sinA0,cotsigma1,sin2sigma1,cos2sigma1,sigma0,sin2,cos2,sigm a,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,&A 11,&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=6356863.020+(10708.949-13.474*(1-sinA0*sinA0))*(1-sinA0*sinA0);B=(5354.469-8.798*(1-sinA0*sinA0))*(1-sinA0*sinA0);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);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);elselamda=pi-fabs(lamda);}else{if(tan(lamda)>0)lamda=fabs(lamda)-pi;elselamda=-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);elseA2=2*pi-fabs(A2);}else{if(tan(A2)>0)A2=fabs(A2);elseA2=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,&B 20,&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);elseA1=pi-fabs(A1);}else{if(q>0)A1=2*pi-fabs(A1);elseA1=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);elsesigma=pi-fabs(sigma);sinA0=cosu1*sin(A1);x=2*a1-(1-sinA0*sinA0)*cos(sigma);afa=(33523299-(28189-70*(1-sinA0*sinA0))*(1-sinA0*sinA0))*1.0e-10; beta=(28189-94*(1-sinA0*sinA0))*1.0e-10;delta=(afa*sigma-beta*x*sin(sigma))*sinA0;lamda=L+delta;while(fabs(delta-delta0)>4.8e-10){delta0=delta;p=cosu2*sin(lamda);q=b1-b2*cos(lamda);A1=atan(p/q);if(p>0){if(q>0)A1=fabs(A1);elseA1=pi-fabs(A1);}else{if(q>0)A1=2*pi-fabs(A1);elseA1=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);elsesigma=pi-fabs(sigma);sinA0=cosu1*sin(A1);x=2*a1-(1-sinA0*sinA0)*cos(sigma);afa=(33523299-(28189-70*(1-sinA0*sinA0))*(1-sinA0*sinA0))*1.0e-10;beta=(28189-94*(1-sinA0*sinA0))*1.0e-10;delta=(afa*sigma-beta*x*sin(sigma))*sinA0;lamda=L+delta;}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)*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);elseA2=2*pi-fabs(A2);}else{if(tan(A2)>0)A2=fabs(A2);elseA2=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);}}四.程序执行结果:。
大地测量学基础习题与思考题及答案含重点及两份武大测绘试题@
《大地测量学基础》习题与思考题一 绪论1.试述您对大地测量学的理解?2.大地测量的定义、作用与基本内容是什么?3.简述大地测量学的发展概况?大地测量学各发展阶段的主要特点有哪些?4.简述全球定位系统(GPS )、激光测卫(SLR )、 甚长基线干涉测量(VIBL )、 惯性测量系统(INS )的基本概念? 二 坐标系统与时间系统1.简述是开普勒三大行星定律? 2.什么是岁差与章动?什么是极移? 3.什么是国际协议原点 CIO?4.时间的计量包含哪两大元素?作为计量时间的方法应该具备什么条件? 5.恒星时、 世界时、 历书时与协调时是如何定义的?其关系如何? 6.什么是大地测量基准?7.什么是天球?天轴、天极、天球赤道、天球赤道面与天球子午面是如何定义的 ? 8.什么是时圈 、黄道与春分点?什么是天球坐标系的基准点与基准面? 9.如何理解大地测量坐标参考框架?10.什么是椭球的定位与定向?椭球的定向一般应该满足那些条件? 11.什么是参考椭球?什么是总地球椭球?12.什么是惯性坐标系?什么协议天球坐标系 、瞬时平天球坐标系、 瞬时真天球坐标系?13.试写出协议天球坐标系与瞬时平天球坐标系之间,瞬时平天球坐标系与瞬时真天球坐标系的转换数学关系式。
14.什么是地固坐标系、地心地固坐标系与参心地固坐标系?15.什么协议地球坐标系与瞬时地球坐标系?如何表达两者之间的关系?16.如何建立协议地球坐标系与协议天球坐标系之间的转换关系,写出其详细的数学关系式。
17.简述一点定与多点定位的基本原理。
18.什么是大地原点?大地起算数据是如何描述的?19.简述1954年北京坐标系、1980年国家大地坐标系、 新北京54坐标系的特点以及它们之间存在相互关系。
20.什么是国际地球自传服务(IERS )、国际地球参考系统(ITRS) 、国际地球参考框架(ITRF)? ITRS 的建立包含了那些大地测量技术,请加以简要说明?21. 站心坐标系如何定义的?试导出站心坐标系与地心坐标系之间的关系?22.试写出不同平面直角坐标换算、不同空间直角坐标换算的关系式?试写出上述两种坐标转换的误差方程式? 23.什么是广义大地坐标微分方程(或广义椭球变换微分方程)?该式有何作用? 三 地球重力场及地球形状的基本理论1.简述地球大气中平流层、对流层与电离层的概念。
白塞尔大地主题正解
o1 o2 − o3
据 sin A1 , tan A2 符号,给 A2 取值,化为度分秒,0.0001 秒 注: e 2 = 6.693 421623 ×10−3 ρ o = 57.295 779 513 082 32 ρ ′′ = 206 264.806 247 096
4
m
δ = m ⋅ sin A0
δ δ ′′
n1 n2 n3
δ ′′ = δ ⋅ ρ ′′
n1 = sin A1 ⋅ sin σ o n2 = cos µ1 ⋅ cos σ o n3 = sin µ1 ⋅ sin σ o ⋅ cos A1
λ = arctan
n1 n2 − n3
λ λ
据 sin A1 , tan λ 符号,给 λ 取值,并化为度分秒,0.0001 秒
4. 系数 α , β 计算
2 2 −10 α = 33 523 299 − (28 189 − 70 cos A0 ) cos A0 ×10
α
β = (14 094.3 − 46.8cos 2 A0 ) cos2 A0 × 10−10
5. 球面长度计算 g1 = C ⋅ cos 2σ 1 g 2 = ( B + g1 ) sin 2σ 1 g = S − g2
σ
σo
取小数点后 9 位
k
sin µ2 B2 B2 取有效数字 10 位 m1 m2
B2 = arctan
sin µ2 z ⋅ 1 − sin 2 µ2
将 B2 化为度分秒形式,取位至 0.0001 秒 7. 终点大地经度计算 m1 = β (h − sin 2σ 1 ) m2 = α ⋅ σ m = m1 + m2
44 797.282 6 m
4 地面观测值归算至椭球面
5.依据大地线外的其他线为基础。 连接椭球面两点的媒介除大地线之外,当然
还有其他一些有意义的线,比如弦线、法截线 等。利用弦线解决大地主题实质是三绝大地切 量问题,由电磁波测距得到法截线弧长。所以 对三边测量的大地主题而言,运用法截弧进行 解法有其优点。当然,这些解算结果还应加上 归化至大地线的改正。
39
二、电磁波测距的归算
将上式按反正弦函数展开级数,舍去五次项,得
40
4.3 大地测量主题解算概述
4.3.1 大地主题解算的一般说明
大地元素:大地经度L、大地纬度B、两点间的大地 线长度S及其正反大地方位角A12、A21。
大地主题解算:如果知道某些大地元素推求另一些 大地元素,这样的计算问题就叫大地主题解算,大地主 题解算有正解和反解。
法截面和法截线:
过椭球面上任意一点可 作一条垂直于椭球面的法线, 包含这条法线的平面叫作法 截面;法截面与椭球面的交 线叫法截线。
3
4.1.1 椭球面上的几种法截线的曲率半径
一、子午圈曲率半径
子午椭圆的一部分上取一微 分弧长DK=dS,相应地有坐标 增量dx,点n是微分弧dS的曲率 中心,于是线段Dn及Kn便是子 午圈曲率半径M。
不在同一子午圈或同一
平行圈上的两点的正反法
裁线是不重合的,它们之
间的夹角△;大地线是两
点间惟一最短线,而且位
于相对法截线之间,并靠
近正法截线。
22
二、大地线的定义和性质
在一等三角测量中,
数值可达干分之一二秒,可 见在一等或相当于一等三角 测量精度的工程三角测量中 是不容忽略的。
大地线与法截线长度之差 只有百万分之一毫米,所以 在实际计算中,这种长度差 异总是可忽略不计的。
大地测量学第四章-5大地测量主题解算
2R
2
ab sin C
A A0 , 3
B B0 , 3
bc sin A ac sin B ab sin C 2 2 2
C C0 3
化算平面角需要球面角超,而球面角超的计算又需要平面 角,因此直接用球面角计算球面角超就带有误差。 当边长不大于90km时,这种误差小于0.0005″,故可直接 用球面角代替平面角计算球面角超ε
dB d 3 B B2 B1 B ( ) M S ( 3 ) M S 3 5次项 dS 24 dS
三、高斯平均引数正算公式
(1)建立级数展开式: 同理可得:
dL d 3 L L2 L1 L ( ) M S ( 3 ) M S 3 5次项 dS 24 dS
Bm BM ,
Am AM
三、高斯平均引数正算公式
(3)求以Bm、Am为依据的导数: 经整理得:
2 Vm S2 2 2 B S cos Am {1 [sin 2 Am ( 2 3t m 2 m ) 2 Nm 24 N m 2 2 2 2 3 m cos 2 Am ( 1 m 9 m t m )]} 5次
第四章 Ⅴ大地测量主题解算
——大地主题解算思路 ——勒让德级数式 ——高斯平均引数正算公式 ——高斯平均引数反算公式
上一讲应掌握的内容 1、垂线偏差改正 垂线偏差对水平方向的影响 "u ( "sin A1 "cos A1) tan 1 2、标高差改正 e2 由照准点高度而引起的改正 "h H 2 cos2 B2 sin 2 A1
以大地线在大地坐标系中的微分方程为基础 主要特点:解算精度与距离有关,距离越长, 收敛越慢,因此只适用于较短的距离。 典型解法:高斯平均引数法
贝塞尔大地主题正反算及其编程
存档日期:存档编号:江苏师范大学科文学院本科生毕业设计(论文)论文题目:贝塞尔大地主题正反算及程序设计*名:**系别:环境与测绘系专业:测绘工程年级、学号: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 贝塞尔法解算大地问题的基本思想基本思想:将椭球面上的大地元素按照贝塞尔投影条件投影到辅助球面上,继而在球面上进行大地问题解算,最后再将球面上的计算结果换算到椭球面上。
白赛尔大地主题解算
11
Fundation of Geodesy
(1)建立级数展开式:
S
S
MP2 2 , MP1 2
B2
BM
(
dB dS
)M
S 2
1 d2B 2 ( dS 2 )
S2 4
1 d3B 6 ( dS 3 )
S3 8
B1 BM
(
dB dS
)M
S 2
1 d2B 2 ( dS 2 )M
S2 4
1 6
(
B
(
dA dS
)
dB dS
A
(
dA dS
)
dA dS
V2 c2
sin
Acos
A(1
2t
2
2
)
(4 194)
▪ 三阶导数
d3B dS 3
V5 c3
cos
A[sin2
A(1
3t 2
2
9 2t 2 )
3 2
cos2
A(1
t2
2
5η2t 2 )]
d 3L 2V 2 dS3 c2 t sec B sin Acos A
计算范例
23
Fundation of Geodesy
24
Fundation of Geodesy
25
Fundation of Geodesy
四、白塞尔大地问题解算
德国天文学家、数学家Bessel (1784~1846)
26
Fundation of Geodesy
上节知识点回顾
2 4
1 3
2 4
A" t01 L" t21B"2 L" t03 L"3
大地测量学知识点
大地测量学知识点第一篇:大地测量学知识点1.大地坐标系:地面点在参考椭圆的位置用大地经度和纬度表示,若地面的点不在椭球面上,它沿法线到椭球面的距离称为大地高2.空间大地直角坐标系:是大地坐标系相应的三维大地直角坐标系3.地心坐标系:定义大地坐标系时,如果选择的旋转椭球为总地球椭球,椭球中心就是地质中心,再定义坐标轴的指向,此时建立的大地坐标系叫做地心坐标系大地方位角:p点的子午面与过p点法线及Q点的平面所成的角度正高系统:地面上一点沿铅垂线到大地水准面的距离正常高系统:一点沿铅垂线到似水准面的距离国家水准网布设的原则:从高级到低级,从整体到局部,分为四个等级布设,逐级控制,逐级加密4.理论闭合差:在闭合的环形水准路线中,由于水准面不平行所产生的闭合差5.大地高系统:地面一点沿法线到椭球面的距离6.平面控制网的测量方法三角测量:在地面上按一定的要求选定一系列的点,他们与周围的邻近点通视,并构成相互联接的三角网状图形,称为三角网,网中各点称为三角点,在各点上可以进行水平角测量,精确观测各三角内角,另外至少精确测量一条三角形边长度D和方位角,作为网的起始边长和起始方位角,推算边长,方位角进而推算各点坐标三边测量:根据三角形的余弦公式,便可求出三角形内角,进而推算出各边的方位角和各点坐标7.国家高程基准的参考面有平均海水面,大地水准面,似大地水准面,参考椭球面1956年黄海高程系统1985年国家高程基准8.角度观测误差分析视准轴误差:视准轴不垂直于水平轴产生水平轴误差:水平轴不垂直于垂直轴产生这2个的消除误差方法为取盘左盘右读数取平均值垂直轴倾斜误差:垂直轴本身偏离铅垂线的位置,即不竖直解决的方法:观测时,气泡不得偏离一格,测回之间重新整理仪器,观测目标的垂直角大于3度,按气泡偏离的格数计算垂直轴倾斜改正9.方向观测法是在一测回内将测站上所有要观测的方向先置盘左位置,逐一照准进行观测,再盘右的位置依次观测,取盘左盘右的平均值作为各方向的观测值。
大地主题解算(深度干货-超精)
大地主题解算(正算)代码与白塞尔大地主题解算大地主题解算(正算)代码:根据经纬度和方向角以及距离计算另外一点坐标新建模块->拷贝下面的大地主题(正算)代码,调用方法示例:起点经度: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)。
白塞尔大地主题解算。
方向:学号:姓名:\一.基本思路:;基本思想:将椭球面上的大地元素按照白塞尔投影条件投影到辅助球面上,继而在球面上进行大地主题解算,最后在将球面上的计算结果换算到椭球面上。
其关键问题是找出椭球面上的大地元素与球面上相应元素之间的关系式,同时解决在球面上进行大地主题解算的方法。
正算流程:1.计算起点的归化纬度2.计算辅助函数值,解球面三角形3.按公式计算相关系数A,B,C 以及α,β4.计算球面长度#5.计算纬度差改正数6.计算终点大地坐标及大地方位角、011122S B C A{sin (cos )}σσσ=-+10101022222sin ()sin sin cos cos σσσσσσ+=+10101022222cos ()cos cos sin sin σσσσσσ+=-001101522B C A[cos ()]sin ()σσσσσσ=++++010122L A sin [(sin ()sin )]λδασβσσσ-==++-2111u u u A sin sin cos cos cos sin σσ=+22222222222222221111u B u B W u B u B B arctan e u sin cos cos tan sin tan cos ⎧==⎪⎪⎨⎪=⎪⎩⎡⎤==--1111A arctan u u A sin sin []cos cos sin sin cos σλσσ=-21L L λδ=+-反算流程:1.辅助计算2.用逐次趋近法同时计算起点大地方位角、球面长度及经差,第一次趋近时,取δ=0。
!计算下式,重复上述计算过程2.3.计算大地线长度S;4.计算反方位角二.已知数据L λδ=+211212u pA u u u u qsin cos tan cos sin sin cos cos λλ==-2121p p u q b b A arctanqsin cos cos λλ==-=11p A q A sin cos tan cos σσ+=11p A q A sin sin cos σ=+12a a cos cos σλ=+arctan sin cos σσσ⎛⎫=⎪⎝⎭011A u A sin cos sin =111u A tan tan sec σ=21+σσσ=02122L A sin [(sin sin )]λδασβσσ-==+-L λδ=+11222222S A B C B C sin (cos )sin (cos )σσσσσ=++-+…三.源代码:#include <>#include <>#define e //克拉索夫斯基椭球体第一偏心率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;double a1,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,&A1 2,&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);elselamda=pi-fabs(lamda);}else{if(tan(lamda)>0)~lamda=fabs(lamda)-pi;elselamda=-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);【elseA2=2*pi-fabs(A2);}else{if(tan(A2)>0)A2=fabs(A2);elseA2=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,&B2 1,&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);elseA1=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);elsesigma=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);elseA1=pi-fabs(A1);}else{》if(q>0)A1=2*pi-fabs(A1);elseA1=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);%elsesigma=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);elseA2=2*pi-fabs(A2);}~else{if(tan(A2)>0)A2=fabs(A2);elseA2=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);}}四.程序执行结果:&。