大地主题解算程序设计报告

合集下载

大地测量学白塞尔大地主题解算实习报告(模板)1

大地测量上机实习报告
一、实习时间：第9周星期三
二、实习地点：系办三楼机房
三、实习内容：白塞尔大地主题解算
四、实习目的：（1）尝试编写白塞尔主题解算的程序
（2）通过实例检验解算的结果
五、实习步骤：
1、先了解书本上关于白塞尔大地主题解算的步骤；
2、双击打开白塞尔大地主题解算的程序，如下：
3、在“选择椭球”中选择“克氏椭球”，“选择正反算”中选择“正算”；
4、然后在输入框中输入相应的数据，如下：
5、接下来单击菜单栏中的“数据正确”，记得正算结果，如下所示：
6、重复第3步，选择反算，并输入相应数据：
7、点击“数据正确”，即得反算结果：
六、实习总结：通过白塞尔大地主题结算方法的上机操作，懂得了白塞尔法结算大地主题是将椭球面上的大地元素按照白塞尔投影条件投影到辅助球面上，继而在球面上进行大地主题解算，再将球面上的计算结果换算到椭球面上的原理，并熟悉白塞尔大地主题解算的程序。

84坐标系大地主题解算

84坐标系大地主题解算84坐标系大地主题解算是地理信息系统中常见的一种坐标系解算方法，用于将地球表面上的点的经纬度坐标转换为平面坐标，以便进行地图绘制、测量和定位等工作。

本文将介绍84坐标系大地主题解算的基本原理、步骤和应用。

首先，84坐标系是一种地理坐标系，也称为WGS84坐标系，是全球定位系统（GPS）所采用的坐标系。

它使用经度和纬度来标识地球上的点，其中经度表示点位于东西方向上的位置，纬度表示点位于南北方向上的位置。

大地主题解算是将经纬度坐标转换为平面坐标的过程，需要考虑地球的椭球形状和地球表面的曲率。

大地主题解算的基本原理是根据地球的椭球形状和曲率，利用数学公式将经纬度坐标转换为平面坐标。

这个过程涉及到椭球的参数、大地主题的计算和坐标转换公式的推导。

其中，椭球的参数包括椭球的长半轴、短半轴和扁率，大地主题的计算包括椭球面上两点之间的距离、方位角和高度差等，坐标转换公式包括正算和反算两种。

大地主题解算的步骤主要包括坐标系的选择、坐标的输入、大地主题的计算和坐标的输出。

在选择坐标系时，需要确定使用的地理坐标系和投影坐标系，以便进行坐标的转换。

在输入坐标时，需要输入待转换的经纬度坐标，计算大地主题时，需要根据输入的坐标和椭球参数计算出相应的平面坐标。

最后，在输出坐标时，需要将计算得到的平面坐标转换为地理坐标，以便进行地图绘制和定位等应用。

大地主题解算在地理信息系统中具有广泛的应用，包括地图绘制、地理定位、导航引导、地质勘探和地球物理勘测等领域。

通过大地主题解算，可以实现地球表面点的精确定位和坐标转换，为地球科学研究和地理信息服务提供了重要的技术支持。

综上所述，84坐标系大地主题解算是地理信息系统中的重要技术，它通过将地球的经纬度坐标转换为平面坐标，实现地球表面点的精确定位和坐标转换。

大地主题解算的基本原理、步骤和应用对于地理信息系统的研究和应用具有重要意义，有助于提高地图绘制和地理定位的精度和效率。

高斯大地主题解算大地测量程序

double oL=L2-L1;
double Wm=sqrt(1-(e2)*sin(Bm)*sin(Bm));
double Vm=sqrt(1-(ee2)*cos(Bm)*cos(Bm));
double Mm=a*(1-e2)/(Wm*Wm*Wm);
double Nm=a/Wm;
double V=oB*Nm/(Vm*Vm)+oB*oL*oL*Nm/(24*Vm*Vm)*cos(Bm)*cos(Bm)*(2+3*tm*tm+2*nm*nm)+oB*oB*oB*Nm/(8*Vm*Vm*Vm*Vm*Vm*Vm)*(nm*nm-tm*tm*nm*nm);
double oAm=tm*cos(Bm)*oL+oB*oB*oL*1/(24*Vm*Vm*Vm*Vm)*cos(Bm)*tm*(2+7*nm*nm+9*tm*tm*nm*nm+5*nm*nm*nm*nm)+oL*oL*oL*1/24*cos(Bm)*cos(Bm)*cos(Bm)*tm*(2+tm*tm+2*nm*nm);
{
double HUDU;
if(du>0)
{
HUDU=(miao/60.0+fen)/60.0+du;
}
else
{
HUDU=-((miao/60.0+fen)/60.0+abs(du));
}//绝对值是防止输入的是负值角度
return HUDU;
}
Bm=B1+0.5*oB;
Am=A1+0.5*oA;
i++;
}while(i<4);//迭代计算3次

大地主题的数值解法

大地主题的数值解法范业明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;

大地主题解算

deltaL = s * Math.Sin(Am) * (1 + s * s * (Math.Sin(Am) * Math.Sin(Am) *
tm * tm - Math.Cos(Am) * Math.Cos(Am) * (1 + gm * gm - 9 * tm * tm * gm * gm)) / (24 * Nm * Nm))
+ 7 * gm * gm + 9 * tm * tm * gm * gm + 5 * Math.Pow(gm, 4)) + Math.Sin(Am) * Math.Sin(Am) * (2
+ tm * tm + 2 * gm * gm)) / (24 * Nm * Nm)) / Nm;
Vm = Math.Sqrt(1 + e12 * Math.Pow(Math.Cos(Bm), 2));
deltaB = Vm * Vm * s * Math.Cos(Am) * (1 + s * s * (Math.Sin(Am) * Math.Sin(Am)
deltaB = B0 = s * Math.Cos(A1) / M;
deltaL = L0 = s * Math.Sin(A1) / (Math.Cos(B1) * N);
deltaA = A0 = L0 * Math.Sin(B1);
B2 = B1 + deltaB;
L2 = L1 + deltaL;
A2 = A1 + deltaA+Math .PI ;
Console.WriteLine("计算结果：");
deltaA = A0 = L0 * Math.Sin(B1);

程序设计样例

1
sin
)
第 7 页 /共 29 页
大地测量课程设计
sin A1符号
-
sin A2符号
+
A2 =
|A2|
-
180|A2|
+ +
180 | A2 |
+ -
360 | A2|
其中，|| 、|A2|是第一象限角。若A2 0, A2 A2 360;若A2 360, A2 A2 360
5.3 程序界面
计算差值小于规定限差，停止迭代。 Step4: 计算经差改正数
L
sin( ) cos(21 ) sin(2 ) cos(41 2 )sin A0
Step5: 计算终点大地坐标及坐标方位角
sin A1符号
tan 符号
sin 2 sin 1 cos cos 1 cos A1 sin
第 5 页 /共 29 页
大地测量课程设计
最终迭代到两次的δ值之差小于给定的允许值。
其中 A,B,C,及α、β、γ的计算如下:
cos2 A0 1 s in2 A0 , k 2 e'2 cos2 A0
A
(1
k2 4
7k 4 64
15k 6 256
)
/
b
B
k2 (
4
k4 8
37k 6 512
qp
cos b1
2 b2
sin cos
A1
arctan(
p) q
P 符号 q 符号
A1
+
+
-
+
-
-
|A1|
180°-|A1|
180°+|A1|

白塞尔大地主题解算

ﻩdoubleａ1,ａ2,b1，b２,p，q,x，y;
ﻩdoubleW1，W2,sinu1，sinu２,cosｕ１,ｃｏsｕ2,sinＡ０,coｔｓｉgｍａ1,sｉn２siｇma1,cｏｓ2sｉgｍa1，ｓiｇma0,sｉn2,cos2,ｓigｍa，siｎs,cｏｓｓ,deｌta0,ｄｅｌta,lamｄa；
ﻩC=(２.２38*(1－sinA０*ｓiｎA0))*(1－sinA0*sinA0）+0．００6;
ﻩaｆa=691．4６768-（0.5８14３-0.00144＊(1-siｎＡ０*ｓｉnA0）)*(1－ｓiｎＡ0*ｓｉnＡ０)；
ﻩbｅta=(0.２９0７-１.０E－３*(１-ｓinＡ０＊ｓinA０）)*(1-sｉnA０*sｉｎA0);
ﻩB１＝(Ｂ1０＋(flｏat）B11/６0+B12／36００)*pi/180;
L1=(L10+（fｌoat）L11/60＋L12/3600)*pi/180;
ﻩﻩA1=(A１０+(floaｔ)A１1/60+A12／３６00)＊pi/180；
W1=sｑｒt(1-ｅ*e*sin(B1)*sin（B1)）;//计算起点规划纬度
ﻩpｉ=４＊aｔan(1)；
pｒiｎtf("白塞尔大地主题正算请输入1\n白塞尔大地主题反算请输入2\ｎ");
ﻩｓcａnf（"%d",&k);
ｉf(k＝=1)
{
ﻩpｒintf(＂请输入大地线起点纬度B经度L,大地方位角A及大地线长度S:＼ｎ"）；
ﻩscａnf(＂％ｄ%d%lf%d%d％lf％d%d%lｆ%lf",&Ｂ１0，＆B11,&B1２,＆L1０，&L１1,&L12,&Ａ10,&A1１,&A１2,&S）;

84坐标系大地主题解算

84坐标系大地主题解算
84坐标系是一种常用的大地坐标系，用于地理测量和定位。

大
地主题解算是指通过观测数据和数学模型，计算出地球上某一点的
大地坐标。

大地主题解算涉及到以下几个方面的内容：
1. 大地椭球模型，大地主题解算基于大地椭球模型，即将地球
近似看作一个椭球体。

常用的大地椭球模型有WGS84、CGCS2000等。

这些模型定义了地球的几何形状和参数，如椭球体的长半轴、扁率等。

2. 观测数据，大地主题解算需要借助于观测数据，包括卫星导
航系统（如GPS）的测量数据、测角仪、水准仪等测量设备的观测
结果。

这些数据包括卫星的位置、接收器的观测值等。

3. 大地测量基本理论，大地主题解算基于大地测量的基本理论，如三角测量、水准测量、重力测量等。

这些理论提供了测量数据的
处理方法和数学模型，用于计算点的坐标和高程。

4. 大地主题解算方法，大地主题解算采用不同的方法，如最小二乘法、无约束最小二乘法、最大似然估计等。

这些方法根据具体的观测数据和模型，选择合适的数学模型和计算方法，以求得最优的解。

5. 大地坐标系统，大地主题解算的结果是点的大地坐标，通常包括经度、纬度和大地高。

经度表示点在东西方向上的位置，纬度表示点在南北方向上的位置，大地高表示点相对于参考椭球体的高度。

综上所述，大地主题解算是一种通过观测数据和数学模型计算地球上某一点的大地坐标的方法。

它涉及到大地椭球模型、观测数据、大地测量基本理论、解算方法和大地坐标系统等方面的内容。

通过大地主题解算，可以实现地理测量和定位的精确计算和分析。

分段累加法大地主题解算与高斯投影C语言程序

分段累加法大地主题解算与高斯投影C语言程序大地主题解算与高斯投影正反算程序设计测绘一班 XX [1**********]96一、大地主题解算：程序源码：#include#include#define PI 3.[1**********]98int main(void){double bx,by,bz,B1,lx,ly,lz,L1,ax,ay,az,A1,S;double dB,dL,dA;double e2=0.[**************],a=6378245.0000000000;double M1,N1,W1,N2,W2;double B2,L2,A2;int Bx,By,Lx,Ly,Ax,Ay;double Bz,Lz,Az;printf("请输入大地线起点纬度\n ° ′ ″\n");scanf("%lf%lf%lf",&bx,&by,&bz);printf("请输入大地线起点经度\n ° ′ ″\n");scanf("%lf%lf%lf",&lx,&ly,&lz);printf("请输入大地方位角\n ° ′ ″\n");scanf("%lf%lf%lf",&ax,&ay,&az);printf("请输入大地线长度\n");scanf("%lf",&S);double W,M,N,C;B1 = (bx + by / 60 + bz / 3600) * PI / 180; L1 = (lx + ly / 60 + lz / 3600) * PI / 180; A1 = (ax + ay / 60 + az / 3600) * PI / 180; W = sqrt(1 - e2 * (sin(B1) * sin(B1)));M = a * (1 - e2) / (W * W * W);N = a / W;C = N * cos(B1) * sin(A1);int n ,i; double s;n = (int)(S / 4000 + 1);s = S / n;for (i = 1; i{dB = cos(A1) * s / M;dL = sin(A1) * s / N / cos(B1);B2 = B1 + dB;L2 = L1 + dL;W = sqrt(1 - e2 * (sin(B2) * sin(B2)));M = a * (1 - e2) / (W * W * W);N = a / W;A2 = asin(C / N / cos(B2));B1 = B2; L1 = L2; A1 = A2;}Bx=B2/PI*180;By=(B2/PI*180-Bx)*60;Bz=fabs(((B2/PI*180-Bx)*60-By)*60);Lx=L2/PI*180;Ly=(L2/PI*180-Lx)*60;Lz=fabs(((L2/PI*180-Lx)*60-Ly)*60);Ax=A2/PI*180;Ay=(A2/PI*180-Ax)*60;Az=fabs(((A2/PI*180-Ax)*60-Ay)*60);printf("大地线终点纬度为%d°%d′%f″\n",Bx,By,Bz);printf("大地线终点经度为%d°%d′%f″\n",Lx,Ly,Lz);printf("终点大地方位角为%d°%d′%f″\n",Ax,Ay,Az);return 0;程序运行截图：高斯正算程序源码：#include#include#define PI 3.[**************]2int main(void){printf("高斯投影算\n\n");int o;printf("请选择采用的椭球参数\n1.克拉索夫斯基椭球体 2.1975年国际椭球体3.WGS-84椭球体 4.CGCS2000\n");scanf("%d",&o);double a = 0,e12,e2;if(o==1){a=6378245.0000000000; e2=0.[**************]; e12=0.[**************]; }else if(o==2){a=6378140.0000000000; e2=0.[**************]; e12=0.[**************]; }else if(o==3){a=6378137.0000000000; e2=0.[1**********]13; e12=0.[1**********]227; }else if(o==4){a=6378137.0;e2=0.[1**********]290; e12=0.[1**********]548; }int k;printf("\n执行高斯投影算法正算\n\n");double bx,by,bz,lx,ly,lz;double m0, m2, m4, m6, m8, a0, a2, a4, a6, a8, N, B, L, X, Bb, cosB ,sinB , ρ, η2, t, l, d, x, y; printf("请输入大地纬度\n ° ′ ″\n");scanf("%lf%lf%lf",&bx,&by,&bz);printf("请输入大地经度\n ° ′ ″\n");scanf("%lf%lf%lf",&lx,&ly,&lz);B=bx+by/60+bz/3600;L=lx+ly/60+lz/3600;Bb = B * PI / 180;m0 = a * (1 - e2);m2 = 3 * e2 * m0 / 2;m4 = 5 * e2 * m2 / 4;m6 = 7 * e2 * m4 / 6;m8 = 9 * e2 * m6 / 8;a0 = m0 + (m2 / 2) + ((3 * m4) / 8) + ((5 * m6) / 16) + ((35 * m8) / 128); a2 = m2 / 2 + m4 / 2 + 15 * m6 / 32 + 7 * m8 / 16;a4 = m4 / 8 + 3 * m6 / 16 + 7 * m8 / 32;a6 = m6 / 32 + m8 / 16;a8 = m8 / 128;ρ = 180 * 3600 / PI;cosB = cos(Bb);sinB = sin(Bb);η2 = e12 * cosB * cosB;t = tan(Bb);d = (int)(L / 6) + 1;l = abs(L - (6 * d - 3)) * 3600;N = a / sqrt(1 - e2 * sinB * sinB);X = a0 * Bb - sinB * cosB * ((a2 - a4 + a6) + (2 * a4 - 16 * a6 / 3) * sinB * sinB + 16 *a6 * sinB * sinB * sinB * sinB / 3);x = X + N * sinB * cosB * l * l / (2 * ρ * ρ) + N * sinB * cosB * cosB * cosB * (5 - t * t+ 9 * η2) * l * l * l * l / (24 * ρ * ρ * ρ * ρ) + N * sinB * cosB * cosB * cosB * cosB * cosB * (61 - 58 * t * t + t * t * t * t) * l * l * l * l * l * l / (720 * ρ * ρ * ρ * ρ * ρ * ρ);y = N * cosB * l / ρ + N * cosB * cosB * cosB * (1 - t * t + η2) * l * l * l / (6 * ρ * ρ *ρ)+ N * cosB * cosB * cosB * cosB * cosB * (5 - 18 * t * t + t * t * t * t) * l * l * l * l * l/ (120 * ρ * ρ * ρ * ρ * ρ);printf("x为%.4lfm\n",x);printf("y为%.4lfm\n",y);}程序截图：高斯反算程序源码：#include#include#define PI 3.[1**********]98int main(void){printf("高斯投影反算\n\n");int o;printf("请选择采用的椭球参数\n1.克拉索夫斯基椭球体 2.1975年国际椭球体3.WGS-84椭球体 4.CGCS2000\n");scanf("%d",&o);double a,e12,e2;if(o==1){a=6378245.0000000000;e2=0.[**************];e12=0.[**************];}else if(o==2){a=6378140.0000000000;e2=0.[**************];e12=0.[**************];}else if(o==3){a=6378137.0000000000;e2=0.[1**********]13;e12=0.[1**********]227;}else if(o==4){a=6378137.0;e2=0.[1**********]290;e12=0.[1**********]548;}double x,y;doublem0,m2,m4,m6,m8,a0,a2,a4,a6,a8,Bf,Bf1,B,Mf,W,Nf,tf,ηf2,n2,n4,n6,n8,l; printf("plaese enter x,y:\n");scanf("%lf %lf",&x,&y);m0=a* (1-e2);m2=3*e2*m0/2;m4=5*e2*m2/4;m6=7*e2*m4/6;m8=9*e2*m6/8;a0=m0+(m2/2)+((3*m4)/8)+((5*m6)/16)+((35*m8)/128);a2=m2/2+m4/2+15*m6/32+7*m8/16;a4=m4/8+3*m6/16+7*m8/32;a6=m6/32+m8/16;a8=m8/128;Bf1=0;Bf=x/a0;while(abs(Bf-Bf1)>1E-10){Bf1=Bf;double sinB=sin(Bf1);double cosB=cos(Bf1);double sin2B=sinB*cosB*2;double sin4B=sin2B*(1-2*sinB*sinB)*2;double sin6B=sin2B*sqrt(1-(sin4B*sin4B)+sin4B*sqrt(1-sin2B*sin2B));Bf=(x-(-a2/2.0*sin2B+a4/4.0*sin4B-a6/6.0*sin6B))/a0;}double sinBf=sin(Bf);double cosBf=cos(Bf);ηf2=e12*cosBf*cosBf;tf=tan(Bf);W=sqrt(1-e2*sinBf*sinBf);Mf=a*(1-e2)/(W*W*W);Nf=a/W;B=Bf-tf*y*y/(2*Mf*Nf)+tf*(5+3*tf*tf+ηf2-9*ηf2*tf*tf)*y*y*y*y/(24*Mf*Nf*Nf*Nf)-tf*(61+90*tf*tf+45*tf*tf*tf*tf)*y*y*y*y*y*y/(720*Mf*Nf*Nf*Nf*Nf*Nf);l=y/(Nf*cosBf)-(1+2*tf*tf+ηf2)*y*y*y/(6*Nf*Nf*Nf*cosBf)+(5+28*tf*tf+24*tf*tf*tf*tf+6*ηf2+8*ηf2*tf*tf)*y*y*y*y*y/(120*Nf*Nf*Nf*Nf*Nf*cosBf);int Bx,By,lx,ly;double Bz,lz;Bx=(B/PI*180);By=((B/PI*180-Bx)*60);Bz=abs(((B/PI*180-Bx)*60-By)*60);By=abs(By);lx=(l/PI*180);ly=((l/PI*180-lx)*60);lz=abs(((l/PI*180-lx)*60-ly)*60);ly=abs(ly);printf("大地纬度为%d°%d′%lf″\n",Bx,By,Bz); printf("经度差为%d°%d′%lf″\n",lx,ly,lz); return 0;}程序截图：。

大地测量学编程实习报告--大地主题解算

大地测量学编程实习报告——大地主题解算学院：测绘学院专业：测绘工程班级：2014级2班学号：*************姓名：**一、实习目的1. 提高运用计算机语言编程开发的能力；2. 加深对大地主题解算计算公式及辅助参数的理解并掌握计算步骤；3. 通过编程语言实现大地主题解算。

二、注意事项1. 计算所需变量多，容易混淆；2. 正反算函数的编写；3. 函数调用；4. 弧度与角度之间的转化。

三、程序框图正算反算四、源代码调用的函数：public double Mjisuan(double B)//M计算{double M;M = a * (1 - e2) *Math.Pow((1 - e2 * Math.Sin(B) * Math.Sin(B)),-1.5);return M;}public double Njisuan(double B)//N计算{double N;N = a * Math.Pow((1 - e2 * Math.Sin(B) * Math.Sin(B)), -0.5);return N;}public double angletorad(double angle)//角度转弧度{double rad;rad = Math.PI * angle / 180;return rad;}public double Amcalculate(double b,double l,double SsinAm,double ScosAm)// Am计算{double c,T,Am;c = Math.Abs(ScosAm / SsinAm);if (Math.Abs(b) > Math.Abs(l))T = Math.Atan(Math.Abs(SsinAm / ScosAm));elseT = Math.PI * 0.25 + Math.Atan((1 - c) / (1 + c));if (b > 0 && l >= 0)Am = T;else if (b < 0 && l >= 0)Am = Math.PI - T;else if (b <= 0 && l < 0)Am = Math.PI + T;else if (b == 0 && l > 0)Am = Math.PI * 0.5;else Am = 2 * Math.PI - T;return Am;}正算private void button1_Click(object sender, EventArgs e){double B1 = Convert.ToDouble(textBox1.Text) + Convert.ToDouble(textBox8.Text) / 60 + Convert.ToDouble(textBox9.Text) / 3600;double L1 = Convert.ToDouble(textBox2.Text) + Convert.ToDouble(textBox11.Text) / 60 + Convert.ToDouble(textBox10.Text) / 3600;double A12 = Convert.ToDouble(textBox4.Text) +Convert.ToDouble(textBox12.Text) / 60 + Convert.ToDouble(textBox13.Text) / 3600;double S = Convert.ToDouble(textBox3.Text);double Am,Bm,Lm,A23,B23,L23,B01,L01,A01,M,N,B02,L02,A02,t,n2,p1,p2,p3,Vm2;B1 = angletorad(B1);L1 = angletorad(L1);A12 = angletorad(A12);M=Mjisuan(B1);N=Njisuan(B1);B01 = S * Math.Cos(A12) / M;L01 = S * Math.Sin(A12)/ (N * Math.Cos(B1));A01= S *Math.Sin(A12)*Math.Tan(B1)/N;Am = A12 + 0.5 * A01;Bm = B1 + 0.5 * B01;M = Mjisuan(Bm);N = Njisuan(Bm);t=Math.Tan(Bm);n2=Math.Cos(Bm)*Math.Cos(Bm)*e2/(1-e2);Vm2 = 1 + n2;B02 = (Vm2/ N) * S * Math.Cos(Am) * (1 + S * S / (24 * N * N) * ((Math.Sin(Am) * Math.Sin(Am) * (2 + 3 * t * t + 3 * n2*t*t) + 3 * Math.Cos(Am) * Math.Cos(Am) * n2 * (t * t - 1 - n2 - 4 * n2 * t * t))));L02=(1/(N*Math.Cos(Bm)))*S*Math.Sin(Am)*(1+S*S/(24*N*N)*(t*t*Math.Sin(Am)*Math.Sin(Am)-Math.Cos(Am)*Math.Cos(Am)*(1+n2-9*n2*t*t+n2*n2)));A02=(1/N)*S*Math.Sin(Am)*t*(1+S*S/(24*N*N)*(Math.Cos(Am)*Math.Cos(Am)*(2+7*n2+9*n2*t*t+5*n2*n2)+Math.Sin(Am)*Math.Sin(Am)*(2+t*t+2*n2)));while (Math.Abs(B02 - B01) > 0.000000001 && Math.Abs(L02 - L01) > 0.000000001 && Math.Abs(A02 - A01) > 0.000000001){A01 = A02;B01 = B02;L01 = L02;Am = A12 + 0.5 * A02;Bm = B1 + 0.5 * B02;Lm = L1 + 0.5 * L02;M = Mjisuan(Bm);N = Njisuan(Bm);t = Math.Tan(Bm);n2 = Math.Cos(Bm) * Math.Cos(Bm) * e2 / (1 - e2);B02 = (Vm2/ N) * S * Math.Cos(Am) * (1 + S * S / (24 * N * N) * ((Math.Sin(Am) * Math.Sin(Am) * (2 + 3 * t * t + 3 * n2 * t * t) + 3 * Math.Cos(Am) * Math.Cos(Am) * n2 * (t * t - 1 - n2 - 4 * n2 * t * t))));L02 = (1 / (N * Math.Cos(Bm))) * S * Math.Sin(Am) * (1 + S * S / (24 * N * N) * (t * t * Math.Sin(Am) * Math.Sin(Am) - Math.Cos(Am) * Math.Cos(Am) * (1 + n2 - 9 * n2 * t * t + n2 * n2)));A02 = (1 / N) * S * Math.Sin(Am) * t * (1 + S * S / (24 * N * N) * (Math.Cos(Am) * Math.Cos(Am) * (2 + 7 * n2 + 9 * n2 * t * t + 5 * n2 * n2) + Math.Sin(Am) * Math.Sin(Am) * (2 + t * t + 2 * n2)));}B02 = B1 + B02;L02 = L1 + L02;A02 = A12 + A02+Math.PI;p1 = 180.0 * A02 / Math.PI;int A21 = (int)p1;int A22 = (int)((p1 - A21) * 60);A23 = (p1 - A21 - A22 / 60.0) * 3600;p2 = 180.0 * B02 / Math.PI;int B21 = (int)p2;int B22 = (int)((p2 - B21) * 60);B23 = (p2 - B21 - B22 / 60.0) * 3600;p3 = 180.0 * L02 / Math.PI;int L21 = (int)p3;int L22 = (int)((p3 - L21) * 60);L23 = (p3 - L21 - L22 / 60.0) * 3600;if (A21 > 360)A21 = A21 - 180;textBox5.Text = B21.ToString();textBox21.Text = B22.ToString();textBox20.Text = B23.ToString();textBox6.Text = L21.ToString();textBox18.Text = L22.ToString();textBox16.Text = L23.ToString();textBox19.Text = A21.ToString();textBox17.Text = A22.ToString();textBox15.Text = A23.ToString();}反算private void button2_Click(object sender, EventArgs e){double B1 = Convert.ToDouble(textBox1.Text) + Convert.ToDouble(textBox8.Text) / 60 + Convert.ToDouble(textBox9.Text) / 3600;double L1 = Convert.ToDouble(textBox2.Text) + Convert.ToDouble(textBox11.Text) / 60 + Convert.ToDouble(textBox10.Text) / 3600;double B2 = Convert.ToDouble(textBox5.Text) + Convert.ToDouble(textBox21.Text) / 60 + Convert.ToDouble(textBox20.Text) / 3600;double L2 = Convert.ToDouble(textBox6.Text) + Convert.ToDouble(textBox18.Text) / 60 + Convert.ToDouble(textBox16.Text) / 3600;B1 = angletorad(B1);L1 = angletorad(L1);B2 = angletorad(B2);L2 = angletorad(L2);double SsinAm, ScosAm, r01, r21, r03, S10, S12, S30, M, N, Am, Bm, t, n2, Vm2, t01, t21, t03, dA, dL, dB, S, A12, A21;double p1, p2, A123, A213;dL = L2 - L1;dB = B2 - B1;Bm = 0.5 * (B1 + B2);M = Mjisuan(Bm);N = Njisuan(Bm);t = Math.Tan(Bm);n2 = Math.Cos(Bm) * Math.Cos(Bm) * e2 / (1 - e2);Vm2 = 1 + n2;r01 = N * Math.Cos(Bm);r21 = N * Math.Cos(Bm) / (24 * Vm2 * Vm2) * (1 + n2 - 9 * n2 * t * t + n2 * n2); r03 = -N / 24 * Math.Cos(Bm) * Math.Cos(Bm) * Math.Cos(Bm) * t;S10 = N / Vm2;S12 = -N / (24 * Vm2) * Math.Cos(Bm) * Math.Cos(Bm) * (2 + 3 * t * t + 2 * n2); S30 = N / (8 * Vm2 * Vm2 * Vm2) * (n2 - t * t * n2 + n2 * n2);t01 = t * Math.Cos(Bm);t21 = 1 / (24 * Vm2 * Vm2) * Math.Cos(Bm) * t * (2 + 7 * n2 + 9 * t * t * n2 + 5 * n2 * n2);t03 = 1 / (24) * Math.Cos(Bm) * Math.Cos(Bm) * Math.Cos(Bm) * t * (2 + t * t + 2 * n2 * t*t);SsinAm = r01 * dL + r21 * dB * dB * dL + r03 * dL * dL * dL;ScosAm = S10 * dB + S12 * dB * dL * dL + S30 * dB * dB * dB;dA = t01 * dL + t21 * dB * dB * dL + t03 * dL * dL * dL;Am = Amcalculate(B2 - B1, L2 - L1, SsinAm, ScosAm);S = SsinAm / Math.Sin(Am);A12 = Am - 0.5 * dA;A21 = Am + 0.5 * dA + Math.PI;p1 = 180.0 * A12 / Math.PI;int A121 = (int)p1;int A122 = (int)((p1 - A121) * 60);A123 = (p1 - A121 - A122 / 60.0) * 3600;p2 = 180.0 * A21 / Math.PI;int A211 = (int)p2;int A212 = (int)((p2 - A211) * 60);A213 = (p2 - A211 - A212 / 60.0) * 3600;if (A121 > 360)A121 = A121 - 180;textBox4.Text = A121.ToString();textBox12.Text = A122.ToString();textBox13.Text = A123.ToString();textBox19.Text = A211.ToString();textBox17.Text = A212.ToString();textBox15.Text = A213.ToString();textBox3.Text = S.ToString();}五、运算结果正算反算（可对比结果，几乎一样）六、实习总结此次实习我收获颇多，我不仅对大地主题解算的公式更加的了解，而且对c#语言的掌握也更加熟悉了。

大地主题解算程序设计报告

}
/*数据错误*/
else
printf("数据错误，请重新输入\n");
}
/*度分秒转换为弧度*/
double hudu(double x,double y,double z)
{
double A0;
A0=(x+y/60+z/3600)*PI/180;
return A0;
}
/*弧度转换为度*/
double du(double B0)
printf("大地线终点经度度分秒分别为：\n%d\n%d\n%lf\n",ex,ey,ez);
printf("终点大地方位角度分秒分别为：\n%d\n%d\n%lf\n",fx,fy,fz);
}
/*大地主题反算*/
else if(k==2)
{
double fx,fy,fz,gx,gy,gz,hx,hy,hz,ix,iy,iz,jz,kz,B1,B2,L1,L2,S,A1,A2;
printf("请输入大地线起点经度度分秒\n");
scanf("%lf%lf%lf",&bx,&by,&bz);
printf("请输入大地方位角度分秒\n");
scanf("%lf%lf%lf",&cx,&cy,&cz);
printf("请输入大地线长度\n");
scanf("%lf",&S);
/*调用函数*/
W2=sqrt(1-e2*sin(B2)*sin(B2));
sinu1=sin(B1)*sqrt(1-e2)/W1;

高斯平均引数大地主题解算程序设计

高斯平均引数大地主题解算程序设计田桂娥;谢露;马广涛【摘要】介绍了大地主题解算基本方法,以Visual Studio 2010作为开发平台,采用C#语言编写,设计实现了一套基于高斯平均引数的大地主题解算实用程序,指出该程序功能强大,实现了单点和批量大地主题解算,能方便的将计算结果保存在文本,且精度满足大地测量及相应工程的需求.【期刊名称】《山西建筑》【年(卷),期】2016(042)027【总页数】2页(P192-193)【关键词】大地主题解算;高斯平均引数法;Visual Studio 2010;批量解算【作者】田桂娥;谢露;马广涛【作者单位】华北理工大学,河北唐山063000;中铁十一局集团第一工程有限公司,湖北襄阳441104;河北省制图院,河北石家庄050000【正文语种】中文【中图分类】TU198在天文大地测量中，为了获得点的大地坐标，需要在椭球面上进行控制点间的坐标解算。

椭球面上两点间的大地坐标(大地经度、大地纬度)、大地线、大地方位角称为大地元素，已知一些大地元素，推求另一些大地元素，通常称为大地主题解算。

大地主题解算包含大地主题正算和大地主题反算两种，已知椭球面上一点P1的大地坐标(L1，B1)，P1到P2点的大地线长度S及其大地方位角A12，计算P2点的大地坐标(L2，B2)和大地线长度S在P2点的反方位角A21，称为大地主题正解；反之，已知P1和P2点的大地坐标(L1，B1)和(L2，B2)，计算P1至P2的大地线长度S及其正反方位角A12和A21，称为大地主题反解。

由于椭球计算的复杂性，带来大地主题解算的复杂性，有的需要进行迭代计算逐步趋近，给人工计算带来极大困难。

随着计算机技术的飞速发展，计算机在大地主题解算上的应用也得到了快速的发展，迭代计算已经不再是难题，而且，可以根据精度的需要而自行确定迭代次数，极大的提高了计算效率。

同时，随着大地主题解算在空间技术领域的广泛运用，大地主题解算已经成为一项重要的研究工作。

贝塞尔大地主题正反算及其编程

存档日期：存档编号：江苏师范大学科文学院本科生毕业设计（论文）论文题目：贝塞尔大地主题正反算及程序设计*名：**系别：环境与测绘系专业：测绘工程年级、学号：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 贝塞尔法解算大地问题的基本思想基本思想：将椭球面上的大地元素按照贝塞尔投影条件投影到辅助球面上，继而在球面上进行大地问题解算，最后再将球面上的计算结果换算到椭球面上。