VC++ MFC高斯平均引数大地主题正反算
高斯投影正反算c代码
高斯投影正反算程序设计一.程序设计流程本程序(de)设计思路如下:(1),程序采用VS08版本作为开发平台,并采用C语言作为开发语言,设计为WindowsForm窗体程序形式.(2),程序主要(de)算法来自于教材.但是本程序为了更加实用,添加了更多(de)解算基准,包括:WGS-84,国际椭球1975,克氏椭球,和2000国家大地坐标系.(3),程序为了更方便(de)读取数据和输出数据,故需要自己定义了固定(de)数据输入格式和数据输出格式或形式,请老师注意查看.二.代码using System;using Systusing SystemponentModel;using System.Data;using System.Drawing;using System.Text;namespace Gauss{public partial class Form1 : Form{//大地坐标//Geodetic Coordinatepublic struct CRDGEODETIC{public double dLongitude;public double dLatitude;public double dHeight;}//笛卡尔坐标//Cartesian Coordinatepublic struct CRDCARTESIAN{public double x;public double y;public double z;}public Form1(){InitializeComponent();}private void button1_Click(object sender, EventArgs e) {double ee = 0;double a = 0;string tt;try{}catch{MessageBox.Show("Gauss Inverse: Choose datum error");return;}if (ttpareTo("克氏椭球")==0){a = 6378245.00;}if (ttpareTo("WGS-84") == 0){a = 6378.00;}if (ttpareTo("1975国际椭球") == 0){a = 6378140.00;ee =}if (ttpareTo("2000国家大地坐标系") == 0){a = 6378.0;}const double pai = 3.1415926;double b = Math.Sqrt(a a (1 - ee ee));double c = a a / b;double epp = Math.Sqrt((a a - b b) / b / b);CRDGEODETIC pcrdGeo;CRDCARTESIAN pcrdCar;double midlong;//求纬度string[] temp;double[] tempradius = new double[3];for (int i = 0; i < 3; i++){tempradius[i] = Convert.ToDouble(temp[i]);}pcrdGeo.dLatitude = tempradius[0] / 180.0 pai + tempradius[1] / 180.0 / 60.0 pai + tempradius[2] / 180 / 60.0 / 60 pai;//求经度for (int i = 0; i < 3; i++){tempradius[i] = Convert.ToDouble(temp[i]);}pcrdGeo.dLongitude = tempradius[0] / 180.0 pai + tempradius[1] / 180.0 / 60.0 pai + tempradius[2] / 180 / 60.0 / 60 pai;int deglon = Convert.ToInt32(pcrdGeo.dLongitude180 / pai);//求中央经度int num; //带号midlong = 0; //默认值,需要制定分带try{}catch{MessageBox.Show("Choose 3/6 error");return;}if (ttpareTo("3度带") == 0){num = Convert.ToInt32(deglon / 6 + 1);midlong = (6 num - 3) / 180.0 pai;}if (ttpareTo("6度带") == 0){num = Convert.ToInt32((deglon + 1.5) / 3);midlong = num 3 pai / 180;}double lp=pcrdGeo.dLongitude - midlong;double N = c / Math.Sqrt(1 + epp eppMath.Cos(pcrdGeo.dLatitude) Math.Cos(pcrdGeo.dLatitude));double M = c / Math.Pow(1 + epp eppMath.Cos(pcrdGeo.dLatitude) Math.Cos(pcrdGeo.dLatitude), 1.5);double ita = epp Math.Cos(pcrdGeo.dLatitude);double t = Math.Tan(pcrdGeo.dLatitude);double Nscnb = N Math.Sin(pcrdGeo.dLatitude) Math.Cos(pcrdGeo.dLatitude);double Ncosb = N Math.Cos(pcrdGeo.dLatitude);double cosb = Math.Cos(pcrdGeo.dLatitude);double X;double m0, m2, m4, m6, m8;double a0, a2, a4, a6, a8;m0 = a (1 - ee ee);m2 = 3.0 / 2.0 m0 ee ee;m4 = 5.0 / 4.0 ee ee m2;m6 = 7.0 / 6.0 ee ee m4;m8 = 9.0 / 8.0 ee ee m6;a0 = m0 + m2 / 2.0 + 3.0 / 8.0 m4 + 5.0 / 16.0 m6 + 35.0 / 128.0 m8;a2 = m2 / 2 + m4 / 2 + 15.0 / 32.0 m6 + 7.0 / 16.0 m8;a4 = m4 / 8.0 + 3.0 / 16.0 m6 + 7.0 / 32.0 m8;a6 = m6 / 32.0 + m8 / 16.0;a8 = m8 / 128.0;double B = pcrdGeo.dLatitude;double sb = Math.Sin(B);double cb = Math.Cos(B);double s2b = sb cb 2;double s4b = s2b (1 - 2 sb sb) 2;double s6b = s2b Math.Sqrt(1 - s4b s4b) + s4b Math.Sqrt(1 - s2b s2b);X = a0 B - a2 / 2.0 s2b + a4 s4b / 4.0 - a6 / 6.0 s6b;pcrdCar.x = Nscnb lp lp / 2.0 + Nscnb cosb cosb Math.Pow(lp, 4) (5 - t t + 9 ita ita + 4 Math.Pow(ita, 4)) / 24.0 + Nscnb Math.Pow(cosb, 4) Math.Pow(lp, 6) (61 - 58 t t + Math.Pow(t, 4)) / 720.0 + X;pcrdCar.y = Ncosb Math.Pow(lp, 1) + Ncosb cosb cosb (1 - t t + ita ita) / 6.0 Math.Pow(lp, 3) + NcosbMath.Pow(lp, 5) Math.Pow(cosb, 4) (5 - 18 t t+ Math.Pow(t, 4) + 14 ita ita - 58 ita ita t t) / 120.0 ;if (pcrdCar.y < 0)pcrdCar.y += 500000;richTextBox1.Text = "Results:\nX:\t" +Convert.ToString(pcrdCar.x) +"\nY:\t"+ Convert.ToString(pcrdCar.y);}private void button2_Click(object sender, EventArgs e) {double ee = 0;double a = 0;string tt;int num; //带号string ytext; //利用y值求带号和中央经线try{}catch{MessageBox.Show("Gauss Inverse: Choose datumerror");return;}if (ttpareTo("克氏椭球") == 0){a = 6378245.00;}if (ttpareTo("WGS-84") == 0){a = 6378.00;}if (ttpareTo("1975国际椭球") == 0){a = 6378140.00;}if (ttpareTo("2000国家大地坐标系") == 0){a = 6378.0;}double b = Math.Sqrt(a a (1 - ee ee));double c = a a / b;double epp = Math.Sqrt((a a - b b) / b / b); CRDGEODETIC pcrdGeo;CRDCARTESIAN pcrdCar;double midlong = 0;//求X,Y和带号pcrdCar.x = Convert.ToDouble(textBox4.Text); ytext = textBox5.Text;string temp = ytext.Substring(0, 2);num = Convert.ToInt32(temp);ytext = ytext.Remove(0, 2);pcrdCar.y = Convert.ToDouble(ytext) - 500000; try{}catch{MessageBox.Show("Choose 3/6 error");return;}if (ttpareTo("3度带") == 0){midlong = num 3 pai / 180;}if (ttpareTo("6度带") == 0){midlong = (6 num - 3) pai / 180;}b = Math.Sqrt(a a (1 - ee ee));c = a a / b;epp = Math.Sqrt(a a - b b) / b;double m0, m2, m4, m6, m8;double a0, a2, a4, a6, a8;m0 = a (1 - ee ee);m2 = 3.0 / 2.0 m0 ee ee;m4 = 5.0 / 4.0 ee ee m2;m6 = 7.0 / 6.0 ee ee m4;m8 = 9.0 / 8.0 ee ee m6;a0 = m0 + m2 / 2.0 + 3.0 / 8.0 m4 + 5.0 / 16.0 m6 + 35.0 / 128.0 m8;a2 = m2 / 2 + m4 / 2 + 15.0 / 32.0 m6 + 7.0 / 16.0 m8;a4 = m4 / 8.0 + 3.0 / 16.0 m6 + 7.0 / 32.0 m8;a6 = m6 / 32.0 + m8 / 16.0;a8 = m8 / 128.0;double Bf, B;Bf = pcrdCar.x / a0;B = 0.0;while (Math.Abs(Bf - B) > 1E-10){B = Bf;double sb = Math.Sin(B);double cb = Math.Cos(B);double s2b = sb cb 2;double s4b = s2b (1 - 2 sb sb) 2;double s6b = s2b Math.Sqrt(1 - s4b s4b) + s4b Math.Sqrt(1 - s2b s2b);Bf = (pcrdCar.x - (-a2 / 2.0 s2b + a4 / 4.0 s4b - a6 / 6.0 s6b)) / a0;}double itaf, tf, Vf, Nf;itaf = epp Math.Cos(Bf);tf = Math.Tan(Bf);Vf = Math.Sqrt(1 + epp epp Math.Cos(Bf)Math.Cos(Bf));Nf = c / Vf;double ynf = pcrdCar.y / Nf;pcrdGeo.dLatitude = Bf - 1.0 / 2.0 Vf Vf tf (ynf ynf - 1.0 / 12.0 Math.Pow(ynf, 4) (5 + 3 tf tf + itaf itaf - 9 Math.Pow(itaf tf, 2)) +1.0 / 360.0 (61 + 90 tf tf + 45 Math.Pow(tf, 4)) Math.Pow(ynf, 6));pcrdGeo.dLongitude = (ynf / Math.Cos(Bf) - (1 + 2 tf tf + itaf itaf) Math.Pow(ynf, 3) / 6.0 / Math.Cos(Bf) +(5 + 28 tf tf + 24 Math.Pow(tf, 4) + 6 itaf itaf + 8 Math.Pow(itaf tf, 2)) Math.Pow(ynf, 5) / 120.0 /Math.Cos(Bf));pcrdGeo.dLongitude = pcrdGeo.dLongitude + midlong;//pcrdGeo.dLatitude = pcrdGeo.dLatitude;richTextBox2.Text = "Results:\nLatitude: " + Convert.ToString(pcrdGeo.dLatitude) + "\nLongtitude: " +Convert.ToString(pcrdGeo.dLongitude);}private void label13_Click(object sender, EventArgs e) {}}}三.程序运行结果分析通过选取书上(de)具体实例进行测试,本程序(de)精度大体满足要求,一般正算(de)精度在0.01米和0.001米之间,反算(de)精度在0.0001秒左右.以下是程序运行(de)截图.。
高斯投影正算与反算的理论方法与实
高斯投影正算与反算的理论方法与实现代码高斯投影是正形投影的一种,同一坐标系中的高斯投影换带计算公式是根据正形投影原理推导出的两个高斯坐标系间的显函数式。
在同一大地坐标系中(例如1954北京坐标系或1980西安坐标系),如果两个高斯坐标系只是主子午线的经度不同,那么显函数式前的系数可以根据坐标系使用的椭球元素和主子午线经度唯一确定。
但如果两个高斯坐标系除了主子午线的经度不同以外,还存在其他线性系,则将线性变换公式代入换带计算的显函数式中,仍然可以得到严密的坐标变换公式。
此时显函数式前的系数等价于使用两个坐标系主子午线的经度和线性变换参数联合求解得到的,可以唯一确定。
//6度带宽54北京坐标系//高斯投影由大地坐标(Unit:Metres)反算经纬度(Unit:DD)void GaussProjInvCal(double X, double Y, double *longitude, double *latitude){int ProjNo; int ZoneWide; ////带宽double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval;double e1,e2,f,a, ee, NN, T,C, M, D,R,u,fai, iPI;iPI = 0.0174532925199433; ////3.1415926535898/180.0;a = 6378245.0; f = 1.0/298.3; //54年北京坐标系参数////a=6378140.0; f=1/298.257; //80年西安坐标系参数ZoneWide = 6; ////6度带宽ProjNo = (int)(X/1000000L) ; //查找带号longitude0 = (ProjNo-1) * ZoneWide + ZoneWide / 2;longitude0 = longitude0 * iPI ; //中央经线X0 = ProjNo*1000000L+500000L;Y0 = 0;xval = X-X0; yval = Y-Y0; //带内大地坐标e2 = 2*f-f*f;e1 = (1.0-sqrt(1-e2))/(1.0+sqrt(1-e2));ee = e2/(1-e2);M = yval;u = M/(a*(1-e2/4-3*e2*e2/64-5*e2*e2*e2/256));fai = u+(3*e1/2-27*e1*e1*e1/32)*sin(2*u)+(21*e1*e1/16-55*e1*e1*e1*e1/32)*sin( 4*u)+(151*e1*e1*e1/96)*sin(6*u)+(1097*e1*e1*e1*e1/512)*sin(8*u);C = ee*cos(fai)*cos(fai);T = tan(fai)*tan(fai);NN = a/sqrt(1.0-e2*sin(fai)*sin(fai));R = a*(1-e2)/sqrt((1-e2*sin(fai)*sin(fai))*(1-e2*sin(fai)*sin(fai))*(1-e2*sin(fai)*sin(fai)));D = xval/NN;//计算经度(Longitude) 纬度(Latitude)longitude1 =longitude0+(D-(1+2*T+C)*D*D*D/6+(5-2*C+28*T-3*C*C+8*ee+24*T*T)*D*D*D*D*D/120)/cos(fai);latitude1 = fai -(NN*tan(fai)/R)*(D*D/2-(5+3*T+10*C-4*C*C-9*ee)*D*D*D*D/24+(61+90*T+298*C+45*T*T-256*ee-3*C*C)*D*D*D*D*D*D/720);//转换为度DD*longitude = longitude1 / iPI;*latitude = latitude1 / iPI;}//高斯投影由经纬度(Unit:DD)正算平面坐标(含带号,Unit:Metres)void GaussProjCal(double longitude, double latitude, double *X, double *Y) {int ProjNo=0; int ZoneWide; ////带宽double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval; double a,f, e2,ee, NN, T,C,A, M, iPI;iPI = 0.0174532925199433; ////3.1415926535898/180.0;ZoneWide = 6; ////6度带宽a=6378245.0; f=1.0/298.3; //54年北京坐标系参数////a=6378140.0; f=1/298.257; //80年西安坐标系参数ProjNo = (int)(longitude / ZoneWide) ;longitude0 = ProjNo * ZoneWide + ZoneWide / 2;longitude0 = longitude0 * iPI ;latitude0=0;longitude1 = longitude * iPI ; //经度转换为弧度latitude1 = latitude * iPI ; //纬度转换为弧度e2=2*f-f*f;ee=e2*(1.0-e2);NN=a/sqrt(1.0-e2*sin(latitude1)*sin(latitude1));T=tan(latitude1)*tan(latitude1);C=ee*cos(latitude1)*cos(latitude1);A=(longitude1-longitude0)*cos(latitude1);M=a*((1-e2/4-3*e2*e2/64-5*e2*e2*e2/256)*latitude1-(3*e2/8+3*e2*e2/32+45*e2 *e2*e2/1024)*sin(2*latitude1)+(15*e2*e2/256+45*e2*e2*e2/1024)*sin(4*latitude1)-(3 5*e2*e2*e2/3072)*sin(6*latitude1));xval = NN*(A+(1-T+C)*A*A*A/6+(5-18*T+T*T+72*C-58*ee)*A*A*A*A*A/120); yval = M+NN*tan(latitude1)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24+(61-58*T+T*T+600*C-330*ee)*A*A*A*A*A*A/720);X0 = 1000000L*(ProjNo+1)+500000L;Y0 = 0;xval = xval+X0; yval = yval+Y0;*X = xval;*Y = yval;}NN卯酉圈曲率半径,测量学里面用N表示M为子午线弧长,测量学里用大X表示fai为底点纬度,由子午弧长反算公式得到,测量学里用Bf表示R为底点所对的曲率半径,测量学里用Nf表示。
高斯投影正反算原理
高斯投影正反算原理高斯投影是一种常用于地图制图的投影方式,也被广泛应用于其他领域的空间数据处理。
高斯投影正反算是对于已知的地球坐标系上的位置(经纬度),通过计算得到该点的平面坐标(东、北坐标),或者对于已知的平面坐标(东、北坐标),通过计算得到该点的地球坐标系上的位置(经纬度)的过程。
本文将详细介绍高斯投影正反算的原理。
一、高斯投影简介高斯投影是一种圆锥投影,其投影面在地球表面的某个经线上,也就是说,投影面是以该经线为轴的圆锥面。
经过对圆锥体的调整后,使其切于地球椭球面,在该经线上进行投影,同时保持沿该经线方向的比例尺一致,从而达到地图上各点在包括该经线的垂直面上映射的目的。
这种投影方式在某一特定区域内得到高精度的结果,因此广泛应用于地图制图。
二、高斯投影数学模型对于高斯投影正反算,需要先建立高斯投影坐标系与地球坐标系的转换模型。
1.高斯投影坐标系的建立高斯投影坐标系的建立需要确定圆锥面的基本参数,首先需要确定其所处的中央子午线,再确定该子午线上的经度为零点,并利用该经线上某一点的经度和该点的高度来确定该点所在的圆锥体。
圆锥体的底面包括所有与地球椭球面相切的圆面,通过对这些圆面进行调整,使得圆锥体转动后能够在中央子午线上进行投影。
在此基础上,可建立高斯投影坐标系,其中投影面为圆锥面,且中央子午线与投影面的交点称为该投影坐标系的中心,投影面的上端点和下端点分别对应正北方向和正南方向。
2.地球坐标系的建立地球坐标系是以地球椭球体为基础建立的,其坐标系原点确定为地球椭球体上的一个特定点。
在已知该点经纬度和高度的前提下,可确定以该点为中心的地球椭球体,并可根据它与地球坐标系之间的转换关系得到平面坐标系。
3.高斯投影坐标系与地球坐标系之间的转换关系由于高斯投影坐标系与地球坐标系存在不同的坐标体系和基准面,因此需要通过数学关系式来建立它们之间的转换关系。
(1)高斯投影坐标系转地球坐标系:已知高斯投影坐标系中任意一点的东北坐标(N,E),以及所属的中央子午线经度λ0、椭球参数a和e,则可通过以下公式求出该点的地球坐标系经纬度(φ,λ)和高度H:A0为以地球椭球体中心为原点,高斯投影坐标系中心投影坐标为(0,0)的点到椭球面的距离。
高斯投影正反算公式_新
高斯投影坐大地坐标系由大地基准面和地图投影确定,由地图投影到特定椭圆柱面后在南北两极剪开展开而成,是对地球表面的逼近,各国或地区有各自的大地基准面,我国目前主要采用的基准面为:1.WGS84基准面,为GPS基准面,17届国际大地测量协会上推荐,椭圆柱长半轴a=6378137m,短半轴b=6356752.3142451m;2.西安80坐标系,1975年国际大地测量协会上推荐,椭圆柱长半轴a=6378140m,短半轴b=6356755.2881575m;3.北京54坐标系,参照前苏联克拉索夫斯基椭球体建立,椭圆柱长半轴6度0-6度,角。
值为正,Y增加500公里;反算则是由高斯平面坐标(X,Y)求解大地坐标(L,B)。
二、计算模/***************************************本文直接依据空间立体三角函数关系得出结果。
*****/(一)正算由图表1,由方程式(1),令,可得在图表2中,,则由椭圆方程,令可知:(三、程序代doubleL=(m_L-6.0*L0//换算成弧度doublexita=atan(b*b*tan(B)/a/a/cos(L));doubledxita=0.000001;doublexi=dxita;x=0.0;doublec=a*a/b/b;while(xi<xita){x+=dxita/sqrt(c*sin(xi)*sin(xi)+cos(xi)*cos(xi));xi+=dxita;坐标*&B,do ubledoubledxi=0.000001;doublexi=dxi;doubleX=0.0;doublec=a*a/b/b;while(X<x/a){X+=dxi/sqrt(c*sin(xi)*sin(xi)+cos(xi)*cos(xi));xi+=dxi;}doubler=a/sqrt(c*sin(xi)*sin(xi)+cos(xi)*cos(xi));。
高斯正反算及空间直角坐标与大地地理坐标转换
高斯正反算及空间直角坐标与大地地理坐标转换一、实验目的与要求1.对以上理论内容的验证与应用。
2.通过学习掌握测绘软件开发过程与方法,初步具备测绘软件开发基本技能。
3.熟练掌握Visual C++编程环境的使用,了解其特点与程序开发过程,掌软件调试、测试的技术方法。
4.分析测绘程序设计技术课程中相关软件的结构和模块功能,掌握结构化程序设计方法和技术,掌握测绘数据处理问题的基本特点。
5.开发相关程序功能模块,独立完成相关问题概念结构分析、程序结构设计、模块设计、代码编写、调试、测试等工作。
二、实验安排1.实验时数12学时。
2.每实验小组可以由3~4人组成,或独立完成。
若由几个人完成程序设计,应进行合理的分工。
三、实验步骤和要点1.熟悉程序设计任务书的基本内容,调查了解软件需求状况,进行需求分析;2.进行总体设计。
根据所调查收集的资料和任务书的要求,对系统的硬件资源进行初步设计,提出硬件配置计划;进行软件总体设计,设计出软件程序功能的模块;3.根据总体设计的结果,进行详细设计,进行数据存储格式设计、算法等,写出逻辑代码;4.编写程序代码,调试运行;5.程序试运行。
最后同学们可根据自己的选题,写出软件开发设计书一份,打印程序代码和运行结果。
四实验原理高斯正反算:高斯正反算包括两部分内容:高斯正算和高斯反算。
简单的说就是大地地理坐标系坐标(B,L)与其对应的高斯平面直角坐标系坐标(x,y)之间的转换。
若已知大地地理坐标系坐标(B,L)解求对应的高斯平面直角坐标系坐标(x,y)称为高斯正算;反之,则为高斯反算。
空间直角坐标与大地地理坐标转换:地球表面可用一个椭球面表示。
设空间直角坐标系为OXYZ,当椭球的中心与空间直角坐标系原点重合,空间坐标系Z 轴与地球旋转重合(北极方向为正),X 轴正向经度为零时,就可以确定空间直角坐标系与大地地理坐标系的数学关系。
⎪⎩⎪⎨⎧+-=+=+=B H e N Z LB H N Y L B H N X sin ])1([sin cos )(cos cos )(2 式中 N 为卯酉圈曲率半径,B e a N 22sin 1-=; e 为椭球偏心率,222a b a e -=(a ,b 为椭球长半轴和短半轴)。
高斯正反算程序
高斯正反算程序
高斯正反算程序是一种用于计算地理坐标和直角坐标之间转换的程序。
它基于高斯投影原理,将地球表面上的点投影到平面上,实现地理坐标和直角坐标之间的转换。
高斯正算程序是将地理坐标(经度和纬度)转换为直角坐标(X 和Y)的程序。
其计算步骤包括:
1. 将经度和纬度转换为弧度;
2. 计算经度余弦、纬度正弦和纬度余弦;
3. 计算X和Y坐标。
高斯反算程序是将直角坐标(X和Y)转换为地理坐标(经度和纬度)的程序。
其计算步骤包括:
1. 计算X和Y坐标的平方和;
2. 计算平方差;
3. 计算经度和纬度的弧度值;
4. 将弧度值转换为度数。
以上是高斯正反算程序的基本步骤,具体的计算公式和细节可能因不同的投影方法和地区而有所不同。
在实际应用中,需要根据具体情况选择合适的投影方法和参数进行计算。
高斯投影正反算编程
高斯投影正反算编程一.高斯投影正反算基本公式(1)高斯正算基本公式(2)高斯反算基本公式以上主要通过大地测量学基础课程得到.这不进行详细的推导.只是列出基本公式指导编程的进行。
二.编程的基本方法和流程图(1)编程的基本方法高斯投影正反算基本上运用了所有的编程基本语句.本文中是利用C++语言进行基本的设计。
高斯正算中对椭球参数和带宽的选择主要运用了选择语句。
而高斯反算中除了选择语句的应用.在利用迭代算法求底点纬度还应用了循环语句。
编程中还应特别注意相关的度分秒和弧度之间的相互转换.这是极其重要的。
(2)相关流程图1)正算2)反算三.编程的相关代码(1)正算# include "stdio.h"# include "stdlib.h"# include "math.h"# include "assert.h"#define pi (4*atan(1.0))int i;struct jin{double B;double L;double L0;};struct jin g[100];main(int argc, double *argv[]) {FILE *r=fopen("a.txt","r"); assert(r!=NULL);FILE *w=fopen("b.txt","w"); assert(r!=NULL);int i=0;while(fscanf(r,"%lf %lf %lf",&g[i].B,&g[i].L,&g[i].L 0)!=EOF){double a,b;int zuobiao;printf("\n请输入坐标系:北京54=1.西安80=2.WGS84=3:");scanf("%d",&zuobiao);getchar();if(zuobiao==1){a=6378245;b=6356863.0187730473;}if(zuobiao==2){a=6378140;b=6356755.2881575287;}if(zuobiao==3){a=6378137;b=6356752.3142;} //选择坐标系//double f=(a-b)/a;double e,e2;e=sqrt(2*f-f*f);e2=sqrt((a/b)*(a/b)-1);//求椭球的第一.第二曲率//double m0,m2,m4,m6,m8;double a0,a2,a4,a6,a8;m0=a*(1-e*e);m2=3*e*e*m0/2;m4=5*e*e*m2/4;m6=7*e*e*m4/6;m8=9*e*e*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;double Bmiao,Lmiao, L0miao;Bmiao=(int)(g[i].B)*3600.0+(int)((g[i].B-(int)(g[i].B) )*100.0)*60.0+(g[i].B*100-(int)(g[i].B*100))*100.0;Lmiao=(int)(g[i].L)*3600.0+(int)((g[i].L-(int)(g[i].L) )*100.0)*60.0+(g[i].L*100-(int)(g[i].L*100))*100.0;L0miao=(int)(g[i].L0)*3600.0+(int)((g[i].L0-(int)(g[i] .L0))*100.0)*60.0+(g[i].L0*100-(int)(g[i].L0*100))*100 .0;double db;db=pi/180.0/3600.0;double B1,L1,l;B1=Bmiao*db;L1= Lmiao*db;l=L1-L0miao*db;//角度转化为弧度//double T=tan(B1)*tan(B1);double n=e2*e2*cos(B1)*cos(B1);double A=l*cos(B1);double X,x,y;X=a0*(B1)-a2*sin(2*B1)/2+a4*sin(4*B1)/4-a6*sin(6*B1)/6 +a8*sin(8*B1)/8;//求弧长//double N=a/sqrt(1-e*e*sin(B1)*sin(B1));int Zonewide;int Zonenumber;printf("\n请输入带宽:3度带或6度带Zonewide=");scanf("%d",&Zonewide);getchar();if(Zonewide==3){Zonenumber=(int)((g[i].L-Zonewide/2)/Zonewide+1);}else if(Zonewide==6){Zonenumber=(int)g[i].L/Zonewide+1;}else{printf("错误");exit(0);}//选择带宽//doubleFE=Zonenumber*1000000+500000;//改写为国家通用坐标//y=FE+N*A+A*A*A*N*(1-T*T+n*n)/6+A*A*A*A*A*N*(5-18*T*T+T *T*T*T+14*n*n-58*n*n*T*T)/120;x=X+tan(B1)*N*A*A/2+tan(B1)*N*A*A*A*A*(5-T*T+9*n*n+4*n *n*n*n)/24+tan(B1)*N*A*A*A*A*A*A*(61-58*T*T+T*T*T*T)/7 20;printf("\n所选坐标系的转换结果:x=%lf y=%lf\n",x,y);fprintf(w,"%lf %lf\n",x,y);//输出结果到文本文件//}fclose(r);fclose(w);system("pause");return 0;}(2)反算# include "stdio.h"# include "stdlib.h"# include "math.h"# include "assert.h"#define pi (4*atan(1.0))double X,Y,B1,B2,B3,F,t;double m0,m2,m4,m6,m8;double a0,a2,a4,a6,a8,a1,b1;double BB,LL,Bf;double e,e1;int d,m,s,i,zuobiao;double sort(double,double);struct jin{double x;double y;double L0;};struct jin g[100];//x,y,L0为输入量:x,y坐标和中央子午线经度//main(int argc, double *argv[]){FILE *r=fopen("c.txt","r");assert(r!=NULL);FILE *w=fopen("d.txt","w");assert(r!=NULL);int i=0;while(fscanf(r,"%lf %lf %lf",&g[i].x,&g[i].y,&g[i].L 0)!=EOF)//文件为空.无法打开//{double a1=6378245.0000000000;//克拉索夫斯基椭球参数//double b1=6356863.0187730473;double a75=6378140.0000000000;//1975国际椭球参数//double b75=6356755.2881575287;double a84=6378137.0000000000;//WGS-84系椭球参数//double b84=6356752.3142000000;double M,N;//mouyou圈曲率半径.子午圈曲率半径// double t,n;double A,B,C;double BB,LL,Bf,LL0,BB0;double a,b;printf("\n选择参考椭球:1=克拉索夫斯基椭球.2=1975国际椭球.3=WGS-84系椭球:");scanf("%d",&zuobiao);getchar();if(zuobiao==1){a=a1;b=b1;}if(zuobiao==2){a=a75;b=b75;}if(zuobiao==3){a=a84;b=b84;}//选择参考椭球.求解第一偏心率e,第二偏心率e1// Bf=sort(a,b);//调用求解底点纬度的函数//double q=sqrt(1-e*e*sin(Bf)*sin(Bf));double G=cos(Bf);M=a*(1-e*e)/(q*q*q);N=a/q;double H,I;A=g[i].y/N;H=A*A*A;I=A*A*A*A*A;t=tan(Bf);n=e1*cos(Bf);B=t*t;C=n*n;BB0=Bf-g[i].y*t*A/(2*M)+g[i].y*t*H/(24*M)*(5+3*B+C-9*B *C)-g[i].y*t*I/(720*M)*(61+90*B+45*B*B);LL0=g[i].L0*pi/180.0+A/G-H/(6*G)*(1.0+2*B+C)+I/(120*G)*(5.0+28*B+24*B*B+6*C+8*B*C);//利用公式求解经纬度// int Bdu,Bfen,Ldu,Lfen;double Bmiao,Lmiao;Ldu=int(LL0/pi*180);Lfen=int((LL0/pi*180)*60-Ldu*60);Lmiao=LL0/pi*180*3600-Ldu*3600-Lfen*60;Bdu=int(BB0/pi*180);Bfen=int((BB0/pi*180)*60-Bdu*60);Bmiao=BB0/pi*180*3600-Bdu*3600-Bfen*60;//将弧度转化为角度//printf("\n所选坐标系的转换结果:%d度%d分%lf秒 %d 度%d分%lf秒 \n",Bdu,Bfen,Bmiao,Ldu,Lfen,Lmiao);fprintf(w,"%d°%d’%lf”%d°%d’%lf”\n",Bdu,Bfen,Bmiao,Ldu,Lfen,Lmiao);//将结果输出到文本文件//}fclose(r);fclose(w);system("pause");return 0;}double sort(double a,double b){double e,e1;e=sqrt(1-(b/a)*(b/a));e1=sqrt((a/b)*(a/b)-1);double m0,m2,m4,m6,m8;double a0,a2,a4,a6,a8;m0=a*(1-e*e);m2=3*e*e*m0/2;m4=5*e*e*m2/4;m6=7*e*e*m4/6;m8=9*e*e*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;B1=g[i].x/a0;do{F=-a2*sin(2*B1)/2+a4*sin(4*B1)/4-a6*sin(6*B1)/6+a8*s in(8*B1)/8;B2=(g[i].x-F)/a0;B3=B1;B1=B2;} while(fabs(B3-B2)>10e-10);//利用迭代算法求解底点纬度//return B2;}。
高斯投影坐标反算公式
d B 0
x dx dl l y dy dl l
x dx l tan dy y l
x N sin B cos 3 N N sin B cos Bl (5 t 2 9 2 4 4 )l 3 l 6 N sin B cos 5 B (61 58t 2 t 4 )l 5 120 y cos 2 B cos 4 B 2 2 2 N cos B{1 (1 t )l (5 18t 2 t 4 )l 4 } l 2 24
121212角子午线收敛角方大地方位角坐标方位对于克拉索夫斯基椭球1根据高斯坐标确定带号计算中央子午线经度六度带24114923计算中央子午线经度2迭代法求取大地纬度迭代开始时设以后每次迭代按下式计算为止3计算反算公式中的各符号的值2428cos120掌握子午线收敛角的定义及作用
上节回顾
高斯投影坐标正算要求: 1、掌握公式中各符号的意义,公式不要求记住; 2、每个同学必须会计算,一般计算两遍。
1 sin B cos 4 Bl 5 (2 t 2 ) 15
1 3
说明:⑴ 为 l 的奇函数,而且 l 越大, 越大; ⑵ 有正负,当描写点在中央子午线以东时, 为正, 以西时, 为负; ⑶ 当 l 不变时,则 随纬度增加而增大。 (2)由平面坐标f B Bf y2 2M f N f cos B f 1 1 2 2 l y (1 2t f n f ) y 3 3 N f cos B f 6 N f cos B f
算例:带第38带 3 x 5586514 .369 m y 4374 .724 m X 5586514 .369 m Y 38504374 .724 m
大地主题解算C语言程序
printf("终点纬度%d°%d′%lf″",BX,BY,BZ);
printf("终点经度%d°%d′%lf″",LX,LY,LZ);
printf("终点角度%d°%d′%lf″",AX,AY,AZ);
}
if(k==2)
{double c3,c4,c5,c6,c7,c8,d;
printf("输入起点纬度");scanf("%lf,%lf,%lf",&bx,&by,&bz);
printf("输入起点经度");scanf("%lf,%lf,%lf",&lx,&ly,&lz);
printf("输入起点角度");scanf("%lf,%lf,%lf",&ax,&ay,&az);
c8=A21-a8*3600-b8*60;
printf("S= %lf\n",S);
printf("A12= %d°%d′%.4lf″\n",a7,b7,c7);
printf("A21= %d°%d′%.4lf″\n",a8,b8,c8);
}
}
L=L2-L1;
B=B2-B1;
Bm=(B1+B2)/2;
Vm=sqrt(1+ee*cos(Bm/(3600*180)*PI)*cos(Bm/(3600*180)*PI));
Nm=c/Vm;
tm=tan(Bm/(3600*180)*PI);
bm=b1+db/2;
程序设计样例
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|
高斯投影坐标正反算.ppt
2
4
6
8
高斯投影坐标正算(3)
dm0 = dX
dB =M
N cos B =N cos B
,
dq dB dq
M
c
m1 = N cos B
= V
cos B
子午线曲率半径
等量纬度定义式
m2
N 2
sin B cos B
m3
N b
cos3 B(1 t 2
2)
m4
N sin B cos3 24
y 2 dn4 dx
y4
N
cos M
B
(n1
3n3
y2
5n3
y4
)
2n2 y 4n4 y3
N cos B ( dn1 y dn3 y3 dn5 y5 )
M
dx
dx
dx
由恒等式两边对应系数相等,从而得待定系数的递推公式
n1
M N cos B
• 高斯平面直角坐标系: 区分为:自然坐标;国家统一坐标。(掌握两者的换算)
§4.9.2 正形投影的一般条件
一、长度比的通用公式推导
dS 2 (MdB)2 ( N cos Bdl )2
ds2 dx2 dy2
长度比平方为:
m2
ds dS
2
dx2 dy2 (MdB)2 (N cos
dn2 dx
y2 dn4 dx
y4
N
cos M
B
(n1
(完整版)高斯投影正反算
高斯投影正反算学院:资源与环境工程工程学院专业:测绘工程学号:X51414012姓名:孙超一、高斯投影概述想象有一个椭圆柱面横套在地球椭球体外面,并与某一条子午线相切,椭圆柱的中心轴通过椭球体的中心,然后用一定投影方法,将中央子午线两侧各一定经差范围内的地区投影到椭圆柱面上,再将此柱面展开即成为投影面。
高斯投影由于是正形投影,故保证了投影的角度不变性,图形的相似性以及在某点各方向上长度比的同一性。
由于采用了同样法则的分带投影,这即限制了长度变形,又保证了在不同投影带中采用相同的简便公式和数表进行变形引起的各项改正的计算,并且带与带间的互相换算也能用相同的公式和方法进行。
高斯投影的这些优点必将使它得到广泛的推广和具有国际意义。
二、高斯投影坐标正算公式1.高斯投影必须满足以下三个条件1)中央子午线投影后为直线2)中央子午线投影后长度不变3)投影具有正形性质,即正形投影条件2.高斯正算公式推导1)由第一个条件可知,由于地球椭球体是一个旋转椭球体,所以高斯投影必然有这样一个性质,即中央子午线东西两侧的投影必然对称于中央子午线。
2)由于高斯投影是换带投影,在每带内经差l是不大的,lρ是一个微小量,所以可以将 X=X (l,q ),Y=Y (l ,q )展开为经差为l 的幂级数,它可写成如下的形式X=m 0+m 2l 2+m 4l 4+…Y=m 1l+m 3l 2+m 5l 5+…式中m 0,m1,m2,…是待定系数,他们都是纬度B 的函数。
3)由第三个条件:∂y ∂l =∂x ∂q 和∂x ∂l =-∂y ∂q ,将上式分别对l 和q 求偏导2340123423401234...........x m m l m l m l m l y n n l n l n l n l =+++++=+++++可得到下式0312123403121234111,,,, 234111,,,,234dm dm dm dm n n n n dq dq dq dq dn dn dn dn m m m m dq dq dq dq ⎧====⎪⎪⎨⎪=-=-=-=-⎪⎩L L 经过计算可以得出232244524632235242225sin cos sin cos (594)224sin cos (6158)720cos cos (1) 6cos (5181458)120N N x X B B l B B t l N B B t t l N y N B l B t l N B t t t l ηηηηη=+⋅+-+++-+=⋅+-++-++-三、高斯投影坐标反算公式推导1.思路:级数展开,应用高斯投影三个条件,待定系数法求解。
高斯投影正反算编程一.高斯投影正反算基本公式
高斯投影正反算编程一.高斯投影正反算基本公式(1)高斯正算基本公式(2)高斯反算基本公式以上主要通过大地测量学基础课程得到,这不进行详细的推导,只是列出基本公式指导编程的进行。
二.编程的基本方法和流程图(1)编程的基本方法高斯投影正反算基本上运用了所有的编程基本语句,本文中是利用C++语言进行基本的设计。
高斯正算中对椭球参数和带宽的选择主要运用了选择语句。
而高斯反算中除了选择语句的应用,在利用迭代算法求底点纬度还应用了循环语句。
编程中还应特别注意相关的度分秒和弧度之间的相互转换,这是极其重要的。
(2)相关流程图1)正算2)反算三.编程的相关代码(1)正算# include "stdio.h"# include "stdlib.h"# include "math.h"# include "assert.h"#define pi (4*atan(1.0))int i;struct jin{double B;double L;double L0;};struct jin g[100];main(int argc, double *argv[]){FILE *r=fopen("a.txt","r");assert(r!=NULL);FILE *w=fopen("b.txt","w");assert(r!=NULL);int i=0;while(fscanf(r,"%lf %lf %lf",&g[i].B,&g[i].L,&g[i].L0)!=EOF){double a,b;int zuobiao;printf("\n请输入坐标系:北京54=1,西安80=2,WGS84=3:");scanf("%d",&zuobiao);getchar();if(zuobiao==1){a=6378245;b=6356863.0187730473;}if(zuobiao==2){a=6378140;b=6356755.2881575287;}if(zuobiao==3){a=6378137;b=6356752.3142;} //选择坐标系//double f=(a-b)/a;double e,e2;e=sqrt(2*f-f*f);e2=sqrt((a/b)*(a/b)-1);//求椭球的第一,第二曲率//double m0,m2,m4,m6,m8;double a0,a2,a4,a6,a8;m0=a*(1-e*e);m2=3*e*e*m0/2;m4=5*e*e*m2/4;m6=7*e*e*m4/6;m8=9*e*e*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;double Bmiao,Lmiao, L0miao;Bmiao=(int)(g[i].B)*3600.0+(int)((g[i].B-(int)(g[i].B))*100.0)*60.0+( g[i].B*100-(int)(g[i].B*100))*100.0;Lmiao=(int)(g[i].L)*3600.0+(int)((g[i].L-(int)(g[i].L))*100.0)*60.0+(g [i].L*100-(int)(g[i].L*100))*100.0;L0miao=(int)(g[i].L0)*3600.0+(int)((g[i].L0-(int)(g[i].L0))*100.0)*60. 0+(g[i].L0*100-(int)(g[i].L0*100))*100.0;double db;db=pi/180.0/3600.0;double B1,L1,l;B1=Bmiao*db;L1= Lmiao*db;l=L1-L0miao*db;//角度转化为弧度//double T=tan(B1)*tan(B1);double n=e2*e2*cos(B1)*cos(B1);double A=l*cos(B1);double X,x,y;X=a0*(B1)-a2*sin(2*B1)/2+a4*sin(4*B1)/4-a6*sin(6*B1)/6+a8*sin(8*B1)/8;//求弧长//double N=a/sqrt(1-e*e*sin(B1)*sin(B1));int Zonewide;int Zonenumber;printf("\n请输入带宽:3度带或6度带Zonewide=");scanf("%d",&Zonewide);getchar();if(Zonewide==3){Zonenumber=(int)((g[i].L-Zonewide/2)/Zonewide+1);}else if(Zonewide==6){Zonenumber=(int)g[i].L/Zonewide+1;}else{printf("错误");exit(0);}//选择带宽//doubleFE=Zonenumber*1000000+500000;//改写为国家通用坐标//y=FE+N*A+A*A*A*N*(1-T*T+n*n)/6+A*A*A*A*A*N*(5-18*T*T+T *T*T*T+14*n*n-58*n*n*T*T)/120;x=X+tan(B1)*N*A*A/2+tan(B1)*N*A*A*A*A*(5-T*T+9*n*n+4*n*n *n*n)/24+tan(B1)*N*A*A*A*A*A*A*(61-58*T*T+T*T*T*T)/720;printf("\n所选坐标系的转换结果:x=%lf y=%lf\n",x,y);fprintf(w,"%lf %lf\n",x,y);//输出结果到文本文件//}fclose(r);fclose(w);system("pause");return 0;}(2)反算# include "stdio.h"# include "stdlib.h"# include "math.h"# include "assert.h"#define pi (4*atan(1.0))double X,Y,B1,B2,B3,F,t;double m0,m2,m4,m6,m8;double a0,a2,a4,a6,a8,a1,b1;double BB,LL,Bf;double e,e1;int d,m,s,i,zuobiao;double sort(double,double);struct jin{double x;double y;double L0;};struct jin g[100];//x,y,L0为输入量:x,y坐标和中央子午线经度// main(int argc, double *argv[]){FILE *r=fopen("c.txt","r");assert(r!=NULL);FILE *w=fopen("d.txt","w");assert(r!=NULL);int i=0;while(fscanf(r,"%lf %lf %lf",&g[i].x,&g[i].y,&g[i].L0)!=EOF)//文件为空,无法打开//{double a1=6378245.0000000000;//克拉索夫斯基椭球参数//double b1=6356863.0187730473;double a75=6378140.0000000000;//1975国际椭球参数//double b75=6356755.2881575287;double a84=6378137.0000000000;//WGS-84系椭球参数//double b84=6356752.3142000000;double M,N;//mouyou圈曲率半径,子午圈曲率半径//double t,n;double A,B,C;double BB,LL,Bf,LL0,BB0;double a,b;printf("\n选择参考椭球:1=克拉索夫斯基椭球,2=1975国际椭球,3=WGS-84系椭球:");scanf("%d",&zuobiao);getchar();if(zuobiao==1){a=a1;b=b1;}if(zuobiao==2){a=a75;b=b75;}if(zuobiao==3){a=a84;b=b84;}//选择参考椭球,求解第一偏心率e,第二偏心率e1// Bf=sort(a,b);//调用求解底点纬度的函数//double q=sqrt(1-e*e*sin(Bf)*sin(Bf));double G=cos(Bf);M=a*(1-e*e)/(q*q*q);N=a/q;double H,I;A=g[i].y/N;H=A*A*A;I=A*A*A*A*A;t=tan(Bf);n=e1*cos(Bf);B=t*t;C=n*n;BB0=Bf-g[i].y*t*A/(2*M)+g[i].y*t*H/(24*M)*(5+3*B+C-9*B*C)-g[i] .y*t*I/(720*M)*(61+90*B+45*B*B);LL0=g[i].L0*pi/180.0+A/G-H/(6*G)*(1.0+2*B+C)+I/(120*G)*(5.0+28 *B+24*B*B+6*C+8*B*C);//利用公式求解经纬度//int Bdu,Bfen,Ldu,Lfen;double Bmiao,Lmiao;Ldu=int(LL0/pi*180);Lfen=int((LL0/pi*180)*60-Ldu*60);Lmiao=LL0/pi*180*3600-Ldu*3600-Lfen*60;Bdu=int(BB0/pi*180);Bfen=int((BB0/pi*180)*60-Bdu*60);Bmiao=BB0/pi*180*3600-Bdu*3600-Bfen*60;//将弧度转化为角度//printf("\n所选坐标系的转换结果:%d度%d分%lf秒%d 度%d分%lf秒\n",Bdu,Bfen,Bmiao,Ldu,Lfen,Lmiao);fprintf(w,"%d°%d’%lf”%d°%d’%lf”\n",Bdu,Bfen,Bmiao,Ldu,Lfen,Lmiao);//将结果输出到文本文件//}fclose(r);fclose(w);system("pause");return 0;}double sort(double a,double b){double e,e1;e=sqrt(1-(b/a)*(b/a));e1=sqrt((a/b)*(a/b)-1);double m0,m2,m4,m6,m8;double a0,a2,a4,a6,a8;m0=a*(1-e*e);m2=3*e*e*m0/2;m4=5*e*e*m2/4;m6=7*e*e*m4/6;m8=9*e*e*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;B1=g[i].x/a0;do{F=-a2*sin(2*B1)/2+a4*sin(4*B1)/4-a6*sin(6*B1)/6+a8*sin(8*B1 )/8;B2=(g[i].x-F)/a0;B3=B1;B1=B2;} while(fabs(B3-B2)>10e-10);//利用迭代算法求解底点纬度//return B2; }。
高斯平均引数计算大地坐标主题反解的迭代算法
高斯平均引数计算大地坐标主题反解的迭代算法近些年来,随着测量新技术的发展,对大地坐标系统的研究已经取得了很大的进步,改进了测量技术,提高了大地坐标系统的精度,但是,由于大地坐标系统是一种复杂的数学模型,大地坐标的反解受到各种因素的影响,这在一定程度上制约了大地坐标系统的发展。
考虑到这个问题,研究者提出运用高斯平均引数计算大地坐标的方法,以解决大地坐标系统反解的问题。
高斯平均引数计算大地坐标主题反解的迭代算法(MAIT)是利用最小二乘法进行拟合,进行大地坐标主题反解的有效算法。
MAIT算法通过利用多次参数迭代,改变拟合参数的统计均值,利用高斯平均引数,使得系统的反解更加准确,从而解决测量精度的问题。
MAIT算法的基础是利用大地坐标系统拟合方程,拟合关系是彩色拟合函数,该函数是由多个参数确定的,可以用于描述大地坐标系统中两点间的关系。
MAIT算法每次迭代时,通过利用大地坐标系统拟合时计算出的拟合参数均值,改变拟合参数统计均值,并利用高斯平均引数更新拟合参数,从而达到大地坐标主题反解的效果。
同时,MAIT算法还能够实现对数据的正则化处理,即利用一致性条件来校正和修正数据以适应拟合函数,进一步提高大地坐标系统的精度。
MAIT算法的主要优势在于,算法的处理过程非常简洁,可以显著缩短大地坐标主题反解的处理时间;另外,由于拟合参数比较稳定,MAIT算法可以降低大地坐标系统中反解失真的可能性。
总之,MAIT算法是一种有效的大地坐标主题反解方法,它利用高斯平均引数来解决拟合参数不稳定的问题,能够显著提高大地坐标系统的精度和处理效率。
因此,MAIT算法有着极大的发展前景,预期它会在大地坐标系统中得到更广泛的应用。
本文以《高斯平均引数计算大地坐标主题反解的迭代算法》为标题,介绍了MAIT算法,即高斯平均引数计算大地坐标主题反解的迭代算法,该算法通过多次迭代不断改变拟合参数的统计均值,利用高斯平均引数,提高了大地坐标系统的精度,这可以有效解决大地坐标系统的反解问题,有望在大地坐标系统中得到更广泛的应用。
高斯平均引数大地主题正反算
(VC++-MFC高斯平均引数大地主题正反算SouthwestJiaotong University地球科学与环境工程学院实验报告书课程名:_______________学号:___________________姓名:___________________指导老师:________________日期:_____________________、实验内容六、计算结果实验体会、目的与要求 (1)…三、计算公式整理程序代1.5.1.6.、目的与要求参考椭球面是大地测量计算的基准面。
大地坐标是椭球面上的基本坐标系,根据大地测量的观测成果(如距离与方向),从大地原点出发,逐点计算在椭球面上的大地坐标;或根据两点的大地坐标,计算它们之间的大地线长度和大地方位角,这类计算称为大地问题解算(或称为大地主题解算)。
大地问题解算的用途是多方面的,随着现代空间技术和航空航天、航海等领域的发展,大地问题解算(尤其是大地反算)有着更为重要的作用,因此需要熟练掌握其计算。
二、实验内容在《大地测量学基础》教材中,介绍了高斯平均引数法与白塞尔方法的计算过程、步骤。
鉴于此,需要熟练掌握高斯平均引数法与白塞尔方法解大地主题问题的基本方法与原理。
采用所熟悉的计算机语言编程计算。
计算时采用克拉索夫椭球参数,至少完成其中一种方法正反算,按照数据序号选取不同的已知数据,在计算结果中注明所选取的数据序号,选取其它数据作为无效数据处理。
三、计算公式整理3.1高斯平均引数正算计算公式(S< 200 km)(1)计算辅助量公式(2)计算甘气AZ的初值t= 耳=萨cm亠方A 7, j —--------- 5 - sic T]■比厶 =P S r,卫1 ■扌]I科1鸟=門旨m E\〈3)计算兀、〔八仏比胡+ ;込1・厶1 =A + 5亠匚》A(=4+〒*再次计算占%* «2却'寸$ C0S A —百肩【曲血Q +圮十城)AI^24 V;—匚2厂匕F U —1—叽一丄汛€山$血厶彳1十-3; 4/1 +兀-9川:)卜立=半-S侖心川亠磐孑[m* 4/2十了恋十沏北十加)A旳24 V fH亠品;占0十叹+2心]}(5>重复计算(3),直到计算满足丄5 出—AS. < 8AZ M-AL;?也(i-A< < f如按弧度计算可取按角度计算可取0.0001 o(6)计算虽、5、缶的最后值£; =■+■ xXS-L: N £]+ AZ:虫:■内4土1±1&护(易>1呦’取+,岛<130^-(1) 凤二丄(场+昂) 也艮二冬一场AL - L:—知二-壽■"J氏Q十九-2心24 rarctan(4) T—亠5■匚05月芹ll-C^cos J xT丸一T — T 2^-T 托2:i-C当^>0, 当ASg 当A5<0, 当卫沁为卫=dAl>0A£>0AKOAZ<0U>03.2、高斯平均引数正算计算公式(S< 200 km)(2) U = Ssin= r cl A£十V - S cos + +A^-/fl AS + f:1AS:AZ 址山式中各系数咖=^€05^心=_善血曲24AL切-市?5月*启+7即:+ /》;) 心帀心弘U2心注:这里对教材公式中相应系数进行了修正与改进©C3) A=ar^fS^^ct4⑸「也_ U sin sin 月比对■厶-;四、程序代码4.1、角度转换类的头文件:#pragma onceconst double Pi=3.141592653589793;class AngleTrans{public :AngleTrans( void );〜AngleTrans( void );double D,F,M,DFM,Rad,Ten;double trans1( double DFM), //度分秒形式的角度转换为弧度形式trans2( double Rad), //弧度形式的角度转换为度分秒形式trans3( double D);//十进制度转化为弧度};4.2、角度转换类的源文件:#include "StdAfx.h"#include "AngleTrans.h"#include <cmath>AngleTrans::AngleTrans( void ){}AngleTrans::〜AngleTrans( void ){}//度分秒转换为弧度double AngleTrans::trans1( double DFM) { D=floor(DFM);F=floor((DFM-D)*100);M=(DFM-D-F/100)*10000; Ten=D+F/60+M/3600;Rad=Ten/180*Pi;return Rad;} //十进 doub { } 制度转化为弧度 class {publitZhengFanSuanzL1,double zA12, double zS), ZhengS uanL( d double };4.3 、正反算类的源文件://弧度转换为度分秒 double AngleTrans::trans2( double Rad) { Ten=Rad/Pi*180;D=floor(Ten);F=(T en-D)*60;M=(F-floor(F))*60;F=floor(F);DFM=D+F/100+M/10000;return DFM;e AngleTrans::trans3( double D)Rad=D/180*Pi;return Rad;once ZhengFanSuan( void );~ZhengFanSuan( void );double zB1,zL1,zA12,zS,fB1,fL1,fB2,fL2;double ZhengSuanB( double zB1, doublezB1,double zL1,double zA12, double zS),ZhengSuanA( double zB1,double zL1,double zA12, double zS); FanSuanA12( double fB1,doublefL1, double fB2,double fL2), FanSuanS( double fB1, double fL1,double fB2, double fL2),FanSuanA21( double fB1, double fL1,double fB2, double fL2);#include "StdAfx.h"#include "ZhengFanSuan.h"#include "AngleTrans.h"#include <cmath>4.3 、正反算类的头文件:ZhengFanSuan::ZhengFanSuan( void ) {}ZhengFanSuan::〜ZhengFanSuan( void ){}AngleTrans _AngleTrans;const double e仁0.0066934216622966,e2=0.006738525414683,a=6378245.0000.b= 6356863.01877,temp=pow(10.0, -10); // 精度要求doub e Calc_M( double z) // 计算Mm{ double x=a*(1-pow(e1,2))/sqrt(pow(1-pow(e1,2)*pow(sin(z),2),3));return x;}doub e Calc_N( double z) // 计算Nm{ double x=a/sqrt(1-pow(e1,2)*pow(sin(z),2));return x;}doub e Calc_t( double z)//计算tm{ double x=tan(z);return x;}doub e Calc_yita( double z) // 计算yitam{ double x=pow(e2,2)*pow(cos(z),2); returnx;}//正算纬度doub e ZhengFanSuan::ZhengSuanB(double zB1,double zL1,double zA12, double zS) { double M[10000],N[10000],t[10000],B[10000],Bm[10000],L[10000],Lm[10000],A[10000],Am[10000],yita[10000];),2)* (2+7,2))* (2+3A[i]=_AngleTrans.trans3((206265/N[i]*zS*sin(Am[i])*t[i]*(1+pow(zS,2)/(24*pow(N[i],2))*(pow(co s(Am[i]pow(yita[i],2)+9*pow(yita[i],2)*pow(t[i],2)+5*pow(yita[i],4))+pow(sin(Am[i]),2)*(2+pow(t[i],2)+2*pow(yita[i],2)))))/3600);while (B[i]-B[i-1]>=temp){B[i+1]=_AngleTrans.trans3((206265/M[i]*zS*cos(Am[i])*(1+pow(zS,2)/(24*pow(N[i],2))*(pow(sin(Am[i])pow(t[i],2)+2*pow(yita[i],2))+3*pow(cos(Am[i]),2)*pow(yita[i],2)*(pow(t[i],2)-1-pow(yita[i],2)-4*pow(yita[i],2)*pow(t[i],2))))/3600);L[i+1]=_AngleTrans.trans3((206265/(N[i]*cos(Bm[i]))*zS*sin(Am[i])*(1+pow(zS,2)/(24*pow(N[i],2))*(po w(t[i],2)* pow(sin(Am[i]),2)-pow(cos(Am[i]),2)*(1+pow(yita[i],2)-9*pow(yita[i],2)*pow(t[i],2)))))/3600);A[i+1]=_AngleTrans.trans3((206265/N[i]*zS*sin(Am[i])*t[i]*(1+pow(zS,2)/(24*pow(N[i],2))*(pow(cos(Adouble .zB1=_AngleTrans.trans1(zB1), _zL1=_AngleTrans.trans1(zL1), _zA12=_AngleTrans.trans1(zA12), final; t[1]=tan(_zB1); yita[1]=pow(e2,2)*pow(cos(_zB1),2); N[1]=a/sqrt(1-pow(e1,2)*pow(sin(_zB1),2)); B[1]=M[1]=a*(1-pow(e1,2))/sqrt(pow(1-pow(e1,2)*pow(sin(_zB1),2),3)); B[0]=_AngleTrans.trans3(206265/M[1]*zS*cos(_zA12)/3600), L[0]=_AngleTrans.trans3(206265*zS*sin(_zA12)/(N[1]*cos(_zB1))/3600), A[0]=L[0]*sin(_zB1); Bm[1]=_zB1+1/2*B[0]; Lm[1]=_zL1+1/2*L[0]; Am[1]=_zA12+1/2*A[0]; int i=1; _AngleTrans.trans3((206265/M[1]*zS*cos(Am[1])*(1+pow(zS,2)/(24*pow(N[1],2))*(pow(sin(Am[1]),2)) (2+3 pow(t[1],2)+2*pow(yita[1],2))+3*pow(cos(Am[1]),2)*pow(yita[1],2)*(pow(t[1],2)-1-pow(yita[1],2)- 4*pow(yita[1],2)*pow(t[1],2))))/3600); L[i]=_AngleTrans.trans3((206265/(N[1]*cos(Bm[1]))*zS*sin(Am[i])*(1+pow(zS,2)/(24*pow(N[i],2))*(pow pow(}//正Idoub{ double M[10000],N[10000],t[10000],B[10000],Bm[10000],L[10000],Lm[10000],A[10000],Am[10000],yita[10000];double_zB1=_AngleTrans.trans1(zB1),_zL1=_AngleTrans.trans1(zL1), _zA12=_AngleTrans.trans1(zA12),final;m 【i]),2)* (2+7*pow(yita[i],2)+9*pow(yita[i],2)*pow(t[i],2)+5*pow(yita[i],4))+pow(sin(Am[i]),2)*(2+pow(t[i],2)+2* pow(yita[i],2)))))/3600); Bm[i+1]=_zB1+1/2*B[i]; Lm[i+1]=_zL1+1/2*L[i]; Am[i+1]=_zA12+1/2*A[i]; M[i+1]=Calc_M(B[i]); N[i+1]=Calc_N(N[i]); yita[i+1]=Calc_yita(yita[i]); t[i+1]=Calc_t(t[i]); final=B[i]; i++; } double Final=final+_zB1; return _AngleTrans.trans2(Final); t[1]=tan(_zB1);yita[1]=pow(e2,2)*pow(cos(_zB1),2);N[1]=a/sqrt(1-pow(e1,2)*pow(sin(_zB1),2));M[1]=a*(1-pow(e1,2))/sqrt(pow(1-pow(e1,2)*pow(sin(_zB1),2),3));B[0]=_AngleTrans.trans3(206265/M[1]*zS*cos(_zA12)/3600),L[0]=_AngleTrans.trans3(206265*zS*sin(_zA12)/(N[1]*cos(_zB1))/3600),A[0]=L[0]*sin(_zB1);Bm[1]=_zB1+1/2*B[0];Am[1]=_zA12+1/2*A[0];int i=1;B[1]=_AngleTrans.trans3((206265/M[1]*zS*cos(Am[1])*(1+pow(zS,2)/(24*pow(N[1],2))*(pow(sin(Am[1]),2))ZhengFanSuan::ZhengSuanL( double zB1,double zL1,double zA12, double zS)(2+7* (2+3* w(t[i] 2)*(2+7* (2+3*pow(t[1],2)+2*pow(yita[1],2))+3*pow(cos(Am[1]),2)*pow(yita[1],2)*(pow(t[1],2)-1-pow(yita[1],2)-4*pow(yita[1],2)*pow(t[1],2))))/3600);L[i]=_AngleTrans.trans3((206265/(N[1]*cos(Bm[1]))*zS*sin(Am[i])*(1+pow(zS,2)/(24*pow(N[i],2))*(pow (t[i],2)* pow(sin(Am[i]),2)-pow(cos(Am[i]),2)*(1+pow(yita[i],2)-9*pow(yita[i],2)*pow(t[i],2)))))/3600);A[i]=_AngleTrans.trans3((206265/N[i]*zS*sin(Am[i])*t[i]*(1+pow(zS,2)/(24*pow(N[i],2))*(pow(cos(Am[i] ),2)*pow(yita[i],2)+9*pow(yita[i],2)*pow(t[i],2)+5*pow(yita[i],4))+pow(sin(Am[i]),2)*(2+pow(t[i],2)+2*pow(yita[i],2)))))/3600);if(L[i]-L[i-1]>=temp)while (L[i]-L[i-1]>=temp)B[i+1]=_AngleTrans.trans3((206265/M[i]*zS*cos(Am[i])*(1+pow(zS,2)/(24*pow(N[i],2))*(pow(sin(Am[i]) ,2))*pow(t[i],2)+2*pow(yita[i],2))+3*pow(cos(Am[i]),2)*pow(yita[i],2)*(pow(t[i],2)-1-pow(yita[i],2)-4*pow(yita[i],2)*pow(t[i],2))))/3600);L[i+1]=_AngleTrans.trans3((206265/(N[i]*cos(Bm[i]))*zS*sin(Am[i])*(1+pow(zS,2)/(24*pow(N[i],2))*(popow(sin(Am[i]),2)-pow(cos(Am[i]),2)*(1+pow(yita[i],2)-9*pow(yita[i],2)*pow(t[i],2)))))/3600);A[i+1]=_AngleTrans.trans3((206265/N[i]*zS*sin(Am[i])*t[i]*(1+pow(zS,2)/(24*pow(N[i],2))*(pow(cos(A m[i]),2)*pow(yita[i],2)+9*pow(yita[i],2)*pow(t[i],2)+5*pow(yita[i],4))+pow(sin(Am[i]),2)*(2+pow(t[i],2)+2*pow(yita[i],2)))))/3600);Bm[i+1]=_zB1+1/2*B[i];Lm[i+1]=_zL1+1/2*L[i];Am[i+1]=_zA12+1/2*A[i];M[i+1]=Calc_M(B[i]);N[i+1]=Calc_N(N[i]);yita[i+1]=Calc_yita(yita[i]);t[i+1]=Calc_t(t[i]);final=L[i];elsefinal=L[i];double Final=final+_zL1;B[1] = (2+3 return _AngleTrans.trans2(Final);//正算大地方位角double ZhengFanSuan::ZhengSuanA( double zB1,double zL1,double zA12, double zS) { double M[100],N[100], t[100],B[100],Bm[100],L[100],Lm[100],A[100],Am[100],yita[10000];double _zB1=_AngleTrans.trans1(zB1),_zL1=_AngleTrans.trans1(zL1), _zA12=_AngleTrans.trans1(zA12), final;t[1]=tan(_zB1);yita[1]=pow(e2,2)*pow(cos(_zB1),2);N[1]=a/sqrt(1-pow(e1,2)*pow(sin(_zB1),2));M[1]=a*(1-pow(e1,2))/sqrt(pow(1-pow(e1,2)*pow(sin(_zB1),2),3));B[0]=_AngleTrans.trans3(206265/M[1]*zS*cos(_zA12)/3600),L[0]=_AngleTrans.trans3(206265*zS*sin(_zA12)/(N[1]*cos(_zB1))/3600), A[0]=L[0]*sin(_zB1);Bm[1]=_zB1+1/2*B[0];Lm[1]=_zL1+1/2*L[0];Am[1]=_zA12+1/2*A[0];int i=1;_AngleTrans.trans3((206265/M[1]*zS*cos(Am[1])*(1+pow(zS,2)/(24*pow(N[1],2))*(pow(sin(Am[1]),2)) pow(t[1],2)+2*pow(yita[1],2))+3*pow(cos(Am[1]),2)*pow(yita[1],2)*(pow(t[1],2)-1-pow(yita[1],2)-4*pow(yita[1],2)*pow(t[1],2))))/3600);L[i]=_AngleTrans.trans3((206265/(N[1]*cos(Bm[1]))*zS*sin(Am[i])*(1+pow(zS,2)/(24*pow(N[i],2)) *(pow pow(),2)*(2+7*pow(yita[i],2)+9*pow(yita[i],2)*pow(t[i],2)+5*pow(yita[i],4))+pow(sin(Am[i]),2)*(2+pow(t[i],2)+2*pow(yita[i],2)))))/3600);while (A[i]-A[i-1]>=temp){w(t[i] 2)* pow(sin(Am[i]),2)-pow(cos(Am[i]),2)*(1+pow(yita[i],2)-9*pow(yita[i],2)*pow(t[i],2)))))/3600);A[i+1]=_AngleTrans.trans3((206265/N[i]*zS*sin(Am[i])*t[i]*(1+pow(zS,2)/(24*pow(N[i],2))*(pow(cos(A m[i]),2)*(2+7*p ow(yita[i],2)+9*pow(yita[i],2)*pow(t[i],2)+5*pow(yita[i],4))+pow(sin(Am[i]),2)*(2+pow(t[i],2)+2* pow(yita[i],2)))))/3600);Bm[i+1]=_zB1+1/2*B[i];Lm[i+1]=_zL1+1/2*L[i];Am[i+1]=_zA12+1/2*A[i];M[i+1]=Calc_M(B[i]);N[i+1]=Calc_N(N[i]);yita[i+1]=Calc_yita(yita[i]);t[i+1]=Calc_t(t[i]);final=A[i];i++;double Final;if(_zA12>Pi)Final=final+_zA12-Pi;elseFinal=final+_zA12+Pi;return _AngleTrans.trans2(Final);//反算Sdoub e ZhengFanSuan::FanSuanS( double fB1, double fL1,double fB2, double fL2)double A12,A21,S;double 」B1=_AngleTrans.trans1(fB1), 」L1=_An gleTrans.trans1」B2=_AngleTrans.trans1(fB2), _f L2=_An double Bm=(_fB1+_fB2)/2,B[i+1]=_AngleTrans.trans3((206265/M[i]*zS*cos(Am[i])*(1+pow(zS,2)/(24*pow(N[i],2))*(pow(sin(Am[i]) ,2))*(2+3*pow(t[i],2)+2*pow(yita[i],2))+3*pow(cos(Am[i]),2)*pow(yita[i],2)*(pow(t[i],2)-1-pow(yita[i],2)-4*pow(yita[i],2)*pow(t[i],2))))/3600);L[i+1]=_AngleTrans.trans3((206265/(N[i]*cos(Bm[i]))*zS*sin(Am[i])*(1+pow(zS,2)/(24*pow(N[i],2))*(po_B=_fB2-_fB1, _L=_fL2-_fL1;}//反Idoub{ double T,A12,A21,S;double _flB1=_AngleTrans.trans1(fB1),」L1=_AngleTrans.trans1(fL1),」B2=_AngleTrans.trans1(fB2),」L2=_AngleTrans.trans1(fL2); double Bm=(_fB1+_fB2)/2,_B=_fB2-_fB1,_L=_fL2-_fL1;double Nm=Calc_N(Bm),tm=Calc_t(Bm), yitam=Calc_yita(Bm);double Vm=sqrt(1+pow(e2,2)*pow(cos(Bm),2)),double Nm=Calc_N(Bm),tm=Calc_t(Bm), yitam=Calc_yita(Bm);double Vm=sqrt(1+pow(e2,2)*pow(cos(Bm),2)),r01= Nm*cos(Bm),r03=(-Nm*pow(cos(Bm),3)*pow(tm,2))/24,r21= (Nm*cos(Bm)*(1+pow(yitam,2)-9*pow(yitam,2)*pow(tm,2)*pow(yitam,4))/pow(Vm,4))/24,s10=Nm/pow(Vm,2),s12=(Nm/pow(Vm,2)*pow(cos(Bm),2)*(2+3*pow(tm,2)+2*pow(yitam,2)))/24,s30=(Nm/pow(Vm,6)*(pow(yitam,2)-pow(yitam,2)*pow(tm,2)))/8,t01=tm*cos(Bm),t21 =m)*tm/(24*pow(Vm,4))*(2+7*pow(yitam,2)+9*pow(yitam,2)*pow(tm,2)+5*pow(yitam,4)),t03=(tm*pow(cos(Bm),3)*(2+pow(tm,2)+2*pow(yitam,2)))/24;double U=r01*_L+r21*pow(_B,2)*_L+r03*pow(_L,3),V=s10*_B+s12*_B*pow(_L,2)+s30*pow(_B,3),_A=t01*_L+t21*pow(_B,2)*_L+t03*pow(_L,3);double Am=atan(U/V);S=V/cos(Am);return S; r03=(-Nm*pow(cos(Bm),3)*pow(tm,2))/24,r21= (Nm*cos(Bm)*(1+pow(yitam,2)-9*pow(yitam,2)*pow(tm,2)*pow(yitam,4))/pow(Vm,4))/24,s10=Nm/pow(Vm,2), s12=(Nm/pow(Vm,2)*pow(cos(Bm),2)*(2+3*pow(tm,2)+2*pow(yitam,2)))/24,ZhengFanSuan::FanSuanA12( double fB1, double fL1,double fB2, double fL2)s30=(Nm/pow(Vm,6)*(pow(yitam,2)-pow(yitam,2)*pow(tm,2)))/8,t01=tm*cos(Bm),t21 =cos(Bm)*tm/(24*pow(Vm,4))*(2+7*pow(yitam,2)+9*pow(yitam,2)*pow(tm,2)+5*pow(yitam,4)),t03=(tm*pow(cos(Bm),3)*(2+pow(tm,2)+2*pow(yitam,2)))/24;double U=r01*_L+r21*pow(_B,2)*_L+r03*pow(_L,3),V=s10*_B+s12*_B*pow(_L,2)+s30*pow(_B,3),A=t01*_L+t21*pow(_B,2)*_L+t03*pow(_L,3);double Am=atan(U/V);A12=Am-0.5*_A;return _AngleTrans.trans2(A12);r21= (Nm*cos(Bm)*(1+pow(yitam,2)-9*pow(yitam,2)*pow(tm,2)*pow(yitam,4))/pow(Vm,4))/24,s10=Nm/pow(Vm,2),s12=(Nm/pow(Vm,2)*pow(cos(Bm),2)*(2+3*pow(tm,2)+2*pow(yitam,2)))/24,s30=(Nm/pow(Vm,6)*(pow(yitam,2)-pow(yitam,2)*pow(tm,2)))/8, t01=tm*cos(Bm),t21 =m)*tm/(24*pow(Vm,4))*(2+7*pow(yitam,2)+9*pow(yitam,2)*pow(tm,2)+5*pow(yitam,4)),double U=r01*_L+r21*pow(_B,2)*_L+r03*pow(_L,3),V=s10*_B+s12*_B*pow(_L,2)+s30*pow(_B,3),_A=t01*_L+t21*pow(_B,2)*_L+t03*pow(_L,3);double Am=atan(U/V);A12=Am-0.5*_A;}//反;doub{ double T,A12,A21,S; double _flB1=_AngleTrans.trans1(fB1), 」L1=_AngleTrans.trans1(fL1), 」B2=_AngleTrans.trans1(fB2), 」L2=_AngleTrans.trans1(fL2); double Bm=(_fB1+_fB2)/2, _B=_fB2-_fB1, _L=_fL2-_fL1; double Nm=Calc_N(Bm), tm=Calc_t(Bm), yitam=Calc_yita(Bm); double Vm=sqrt(1+pow(e2,2)*pow(cos(Bm),2)), r01= Nm*cos(Bm), r03=(-Nm*pow(cos(Bm),3)*pow(tm,2))/24,ZhengFanSuan::FanSuanA21( double fB1, double fL1,double fB2, double fL2)if(A12>Pi){ A2仁Am+0.5*_A-Pi;}else{ A21=Am+0.5*_A+Pi;}return _AngleTrans.trans2(A21);}4.4 、正算的计算按钮代码:void C大地主题高斯引数正反算Dlg::OnBnClickedButton2(){UpdateData( true );ZhengFanSuan _ZhengFanSuan;double B1=_wtof(zB1),L1=_wtof(zL1),A12=_wtof(zA12),S=_wtof(zS),B2=_ZhengFanSuan.ZhengSuanB(B1,L1,A12,S),L2=_ZhengFanSuan.ZhengSuanL(B1,L1,A12,S),A21=_ZhengFanSuan.ZhengSuanA(B1,L1,A12,S); zB2.Format(_T( "%.7f"),B2);zL2.Format(_T( "%.7f"),L2),zA21.Format(_T( "%.7f"),A21);UpdateData( false);}4.5 、反算的计算按钮代码:void C大地主题高斯引数正反算Dlg::OnBnClickedButton1()UpdateData( true );{ZhengFanSuan _ZhengFanSuan;double B1=_wtof(fB1),L1=_wtof(fL1),B2=_wtof(fB2),L2=_wtof(fL2),S=_ZhengFanSuan.FanSuanS(B1,L1,B2,L2),A12=_ZhengFanSuan.FanSuanA12(B1,L1,B2,L2),A21=_ZhengFanSuan.FanSuanA21(B1,L1,B2,L2);fS.Format(_T( "%.7f"),S);fA12.Format(_T( "%.7f"),A12);fA21.Format(_T( "%.7f"),A21);UpdateData( false);}void{ 五计算结果据4.6、清零按钮代码:void C 大地主题高斯引数正反算Dlg::OnBnClickedButton3(){ UpdateData( true );fB1.Format(_T( "%.0f"),0);fL1.Format(_T( "%.0f"),0);fB2.Format(_T( "%.0f"),0);fL2.Format(_T( "%.0f"),0);fA12.Format(_T( "%.0f"),0);fA21.Format(_T( "%.0f"),0);fS.Format(_T( "%.0f"),0);UpdateData( false); C 大地主题高斯引数正反算 Dlg::OnBnClickedButton4()UpdateData( true );zB1.Format(_T( "%.0f"),0);zL1.Format(_T( %0f"),0);zB2.Format(_T( "%.0f"),0);zL2.Format(_T( "%.0f"),0);zA12.Format(_T( "%.0f"),0);zA21.Format(_T( "%.0f"),0);zS.Format(_T( "%.0f"),0);UpdateData( false);六、实习体会实验一开始本来考虑用控制台做,后来发现用控制台做出来可视化效果不好,而且程序中要用很多cin 和cout,写代码十分麻烦,所以最后选择了用MFC做。
高斯平均引数计算大地坐标主题反解的迭代算法
S
im
1) }
当两次迭代值之差小于给定值, 即
ui - ui- 1
v i - v i- 1
( 9)
停止迭代。
将 u, v 代入式( 5) 中的第 3 式, 则可算得 A , 于是可计算出
A 12=
Am-
1 2
A
A 21=
Am+
1 2
A
并由 u , v 可以反算出
( 10)
S = u2+ v 2
A m=
L 1= 35 49 36. 330 0 L 2= 36 14 45. 050 4 利用反算程序求得
S = 44 797. 282 7 m
A 12= 44 12 13. 664
A 21= 224 30 53. 550
四、结 论
高斯平均引数反算的迭代法计算, 一般只需要 迭代 2~ 3 次, 即能满足精 度要求。该 方法理论简 单, 编程容易, 便于掌握与应用。建议用迭代法取代 传统的高斯平均引数反算公式。
一、高斯平均引数的正解公式
高斯平均引数正解公式[ 1] 为
B
=
( B 2-
B 1)
=
V
2 m
Nm
S cos A m
{ 1+
2
S2 4N
2 m
[
sin2
A m ( 2+
3
t
2 m
+
2
2 m
)
+
3
2 m
cos2
A m(
-
1+
t
2 m
-
2 m
-
4
t
2 m
2 m
大地测量高斯投影正反算
高斯投影正反算姓名:王义文班级:四班学号:2009301610135程序说明:本程序基于MFC基本对话框,由于MFC程序构建框架的代码冗长,打印出来恐怕影响阅读效率,所以下面的代码是经过剔除之后的一些核心代码,希望老师谅解!界面关联变量说明:纬度值的度、分、秒分别与B_DD,B_MM,B_SS相关联,缺省项自动设定为0。
经度值的度、分、秒分别与L_DD,L_MM,L_SS相关联,同上。
L0_DD则表示为L0的值,单位:度。
x,y分别与坐标值关联,且y为加500KM以后的值。
单位:米。
m_keshi=1是克氏椭球,=0为IAG椭球,且此处不能缺省。
m_zheng=1表示正算,=0为反算,且此处不能缺省。
核心代码说明:void C高斯投影正反算4Dlg::OnBnClickedRadio3() //此处为正算的设定,如果正算,设定输入焦{ 点在B_DD。
x,y为只读项。
m_zheng=1;((CEdit*)GetDlgItem(IDC_EDIT1))->SetFocus();((CEdit*)GetDlgItem(IDC_EDIT1))->SetReadOnly(FALSE);((CEdit*)GetDlgItem(IDC_EDIT2))->SetReadOnly(FALSE);((CEdit*)GetDlgItem(IDC_EDIT3))->SetReadOnly(FALSE);((CEdit*)GetDlgItem(IDC_EDIT4))->SetReadOnly(FALSE);((CEdit*)GetDlgItem(IDC_EDIT5))->SetReadOnly(FALSE);((CEdit*)GetDlgItem(IDC_EDIT6))->SetReadOnly(FALSE);((CEdit*)GetDlgItem(IDC_EDIT7))->SetReadOnly(FALSE);((CEdit*)GetDlgItem(IDC_EDIT8))->SetReadOnly();((CEdit*)GetDlgItem(IDC_EDIT9))->SetReadOnly();}void C高斯投影正反算4Dlg::OnBnClickedRadio4() //反算的设定:设置输入焦点在L0_DD,且B,{ L的角度值为只读。
分段累加法大地主题解算与高斯投影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;}程序截图:。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
地球科学与环境工程学院实验报告书课程名:学号:姓名:指导老师:日期:目录一、目的与要求 (1)二、实验容 (1)三、计算公式整理 (1)四、程序代码 (4)五、计算结果 (15)六、实验体会 (16)一、目的与要求参考椭球面是测量计算的基准面。
坐标是椭球面上的基本坐标系,根据测量的观测成果(如距离与方向),从原点出发,逐点计算在椭球面上的坐标;或根据两点的坐标,计算它们之间的线长度和方位角,这类计算称为问题解算(或称为主题解算)。
问题解算的用途是多方面的,随着现代空间技术和航空航天、航海等领域的发展,问题解算(尤其是反算)有着更为重要的作用,因此需要熟练掌握其计算。
二、实验容在《测量学基础》教材中,介绍了高斯平均引数法与白塞尔方法的计算过程、步骤。
鉴于此,需要熟练掌握高斯平均引数法与白塞尔方法解主题问题的基本方法与原理。
采用所熟悉的计算机语言编程计算。
计算时采用克拉索夫椭球参数,至少完成其中一种方反算,按照数据序号选取不同的已知数据,在计算结果中注明所选取的数据序号,选取其它数据作为无效数据处理。
三、计算公式整理3.1、高斯平均引数正算计算公式(S< 200 km)3.2、高斯平均引数正算计算公式(S< 200 km)四、程序代码4.1、角度转换类的头文件:#pragma onceconst double Pi=3.141592653589793;class AngleTrans{public:AngleTrans(void);~AngleTrans(void);double D,F,M,DFM,Rad,Ten;double trans1(double DFM), //度分秒形式的角度转换为弧度形式trans2(double Rad), //弧度形式的角度转换为度分秒形式trans3(double D); //十进制度转化为弧度};4.2、角度转换类的源文件:#include"StdAfx.h"#include"AngleTrans.h"#include<cmath>AngleTrans::AngleTrans(void){}AngleTrans::~AngleTrans(void){}//度分秒转换为弧度double AngleTrans::trans1(double DFM){ D=floor(DFM);F=floor((DFM-D)*100);M=(DFM-D-F/100)*10000;Ten=D+F/60+M/3600;Rad=Ten/180*Pi;return Rad;}//弧度转换为度分秒double AngleTrans::trans2(double Rad){ Ten=Rad/Pi*180;D=floor(Ten);F=(Ten-D)*60;M=(F-floor(F))*60;F=floor(F);DFM=D+F/100+M/10000;return DFM;}//十进制度转化为弧度double AngleTrans::trans3(double D){ Rad=D/180*Pi;return Rad;}4.3、正反算类的头文件:#pragma onceclass ZhengFanSuan{public:ZhengFanSuan(void);~ZhengFanSuan(void);double zB1,zL1,zA12,zS,fB1,fL1,fB2,fL2;double ZhengSuanB(double zB1,double zL1,double zA12,double zS),ZhengSuanL(double zB1,double zL1,double zA12,double zS),ZhengSuanA(double zB1,double zL1,double zA12,double zS);double FanSuanA12(double fB1,double fL1,double fB2,double fL2),FanSuanS(double fB1,double fL1,double fB2,double fL2),FanSuanA21(double fB1,double fL1,double fB2,double fL2); };4.3、正反算类的源文件:#include"StdAfx.h"#include"ZhengFanSuan.h"#include"AngleTrans.h"#include<cmath>ZhengFanSuan::ZhengFanSuan(void)}ZhengFanSuan::~ZhengFanSuan(void){}AngleTrans _AngleTrans;const double e1= 0.0066934216622966,e2=0.006738525414683,a=6378245.0000,b= 6356863.01877,temp=pow(10.0, -10);//精度要求double Calc_M(double z) //计算Mm{ double x=a*(1-pow(e1,2))/sqrt(pow(1-pow(e1,2)*pow(sin(z),2),3));return x;}double Calc_N(double z) //计算Nm{ double x=a/sqrt(1-pow(e1,2)*pow(sin(z),2));return x;}double Calc_t(double z) //计算tm{ double x=tan(z);return x;}double Calc_yita(double z) //计算yitam{ double x=pow(e2,2)*pow(cos(z),2);return x;}//正算纬度double ZhengFanSuan::ZhengSuanB(double zB1,double zL1,double zA12,double zS) { double M[10000],N[10000],t[10000],B[10000],Bm[10000],L[10000],Lm[10000],A[10000],Am[10000],yita[10000];double_zB1=_AngleTrans.trans1(zB1),_zL1=_AngleTrans.trans1(zL1),_zA12=_AngleTrans.trans1(zA12),final;t[1]=tan(_zB1);yita[1]=pow(e2,2)*pow(cos(_zB1),2);N[1]=a/sqrt(1-pow(e1,2)*pow(sin(_zB1),2));M[1]=a*(1-pow(e1,2))/sqrt(pow(1-pow(e1,2)*pow(sin(_zB1),2),3));B[0]=_AngleTrans.trans3(206265/M[1]*zS*cos(_zA12)/3600),L[0]=_AngleTrans.trans3(206265*zS*sin(_zA12)/(N[1]*cos(_zB1))/3600),A[0]=L[0]*sin(_zB1);Bm[1]=_zB1+1/2*B[0];Lm[1]=_zL1+1/2*L[0];Am[1]=_zA12+1/2*A[0];int i=1;B[1]=_AngleTrans.trans3((206265/M[1]*zS*cos(Am[1])*(1+pow(zS,2)/(24*pow(N[1],2))*(pow(sin(Am[1]),2))*(2+3*pow(t[1],2)+2*pow(yita[1],2))+3*pow(cos(Am[1]),2)*pow(yita[1],2)*(pow(t[1],2)-1-pow(yita[1],2)-4*pow(yita[1],2)*pow(t[1],2))))/3600);L[i]=_AngleTrans.trans3((206265/(N[1]*cos(Bm[1]))*zS*sin(Am[i])*(1+pow(zS,2)/(24*pow(N[i],2))*(pow(t[i] ,2)*pow(sin(Am[i]),2)-pow(cos(Am[i]),2)*(1+pow(yita[i],2)-9*pow(yita[i],2)*pow(t[i],2)))))/3600);A[i]=_AngleTrans.trans3((206265/N[i]*zS*sin(Am[i])*t[i]*(1+pow(zS,2)/(24*pow(N[i],2))*(pow(cos(Am[i]),2 )*(2+7*pow(yita[i],2)+9*pow(yita[i],2)*pow(t[i],2)+5*pow(yita[i],4))+pow(sin(Am[i]),2)*(2+pow(t[i],2)+2*pow(yita[i],2)))))/3600);while(B[i]-B[i-1]>=temp){B[i+1]=_AngleTrans.trans3((206265/M[i]*zS*cos(Am[i])*(1+pow(zS,2)/(24*pow(N[i],2))*(pow(sin(Am[i]),2) )*(2+3*pow(t[i],2)+2*pow(yita[i],2))+3*pow(cos(Am[i]),2)*pow(yita[i],2)*(pow(t[i],2)-1-pow(yita[i],2)-4*pow(yita[i],2)*pow(t[i],2))))/3600);L[i+1]=_AngleTrans.trans3((206265/(N[i]*cos(Bm[i]))*zS*sin(Am[i])*(1+pow(zS,2)/(24*pow(N[i],2))*(pow(t[ i],2)*pow(sin(Am[i]),2)-pow(cos(Am[i]),2)*(1+pow(yita[i],2)-9*pow(yita[i],2)*pow(t[i],2)))))/3600);A[i+1]=_AngleTrans.trans3((206265/N[i]*zS*sin(Am[i])*t[i]*(1+pow(zS,2)/(24*pow(N[i],2))*(pow(cos(Am[i] ),2)*(2+7*pow(yita[i],2)+9*pow(yita[i],2)*pow(t[i],2)+5*pow(yita[i],4))+pow(sin(Am[i]),2)*(2+pow(t[i],2)+2*pow(yita[i],2)))))/3600);Bm[i+1]=_zB1+1/2*B[i];Lm[i+1]=_zL1+1/2*L[i];Am[i+1]=_zA12+1/2*A[i];M[i+1]=Calc_M(B[i]);N[i+1]=Calc_N(N[i]);yita[i+1]=Calc_yita(yita[i]);t[i+1]=Calc_t(t[i]);final=B[i];i++;}double Final=final+_zB1;return _AngleTrans.trans2(Final);}//正算经度double ZhengFanSuan::ZhengSuanL(double zB1,double zL1,double zA12,double zS){ double M[10000],N[10000],t[10000],B[10000],Bm[10000],L[10000],Lm[10000],A[10000],Am[10000],yita[10000];double_zB1=_AngleTrans.trans1(zB1),_zL1=_AngleTrans.trans1(zL1),_zA12=_AngleTrans.trans1(zA12),final;t[1]=tan(_zB1);yita[1]=pow(e2,2)*pow(cos(_zB1),2);N[1]=a/sqrt(1-pow(e1,2)*pow(sin(_zB1),2));M[1]=a*(1-pow(e1,2))/sqrt(pow(1-pow(e1,2)*pow(sin(_zB1),2),3));B[0]=_AngleTrans.trans3(206265/M[1]*zS*cos(_zA12)/3600),L[0]=_AngleTrans.trans3(206265*zS*sin(_zA12)/(N[1]*cos(_zB1))/3600),A[0]=L[0]*sin(_zB1);Bm[1]=_zB1+1/2*B[0];Lm[1]=_zL1+1/2*L[0];Am[1]=_zA12+1/2*A[0];int i=1;B[1]=_AngleTrans.trans3((206265/M[1]*zS*cos(Am[1])*(1+pow(zS,2)/(24*pow(N[1],2))*(pow(sin(Am[1]),2))*(2+3*pow(t[1],2)+2*pow(yita[1],2))+3*pow(cos(Am[1]),2)*pow(yita[1],2)*(pow(t[1],2)-1-pow(yita[1],2)-4*pow(yita[1],2)*pow(t[1],2))))/3600);L[i]=_AngleTrans.trans3((206265/(N[1]*cos(Bm[1]))*zS*sin(Am[i])*(1+pow(zS,2)/(24*pow(N[i],2))*(pow(t[i],2)*pow(sin(Am[i]),2)-pow(cos(Am[i]),2)*(1+pow(yita[i],2)-9*pow(yita[i],2)*pow(t[i],2)))))/3600);A[i]=_AngleTrans.trans3((206265/N[i]*zS*sin(Am[i])*t[i]*(1+pow(zS,2)/(24*pow(N[i],2))*(pow(cos(Am[i]),2 )*(2+7*pow(yita[i],2)+9*pow(yita[i],2)*pow(t[i],2)+5*pow(yita[i],4))+pow(sin(Am[i]),2)*(2+pow(t[i],2)+2*pow(yita[i],2)))))/3600);if(L[i]-L[i-1]>=temp){ while(L[i]-L[i-1]>=temp){B[i+1]=_AngleTrans.trans3((206265/M[i]*zS*cos(Am[i])*(1+pow(zS,2)/(24*pow(N[i],2))*(pow(sin(Am[i]),2) )*(2+3*pow(t[i],2)+2*pow(yita[i],2))+3*pow(cos(Am[i]),2)*pow(yita[i],2)*(pow(t[i],2)-1-pow(yita[i],2)-4*pow(yita[i],2)*pow(t[i],2))))/3600);L[i+1]=_AngleTrans.trans3((206265/(N[i]*cos(Bm[i]))*zS*sin(Am[i])*(1+pow(zS,2)/(24*pow(N[i],2))*(pow(t[ i],2)*pow(sin(Am[i]),2)-pow(cos(Am[i]),2)*(1+pow(yita[i],2)-9*pow(yita[i],2)*pow(t[i],2)))))/3600);A[i+1]=_AngleTrans.trans3((206265/N[i]*zS*sin(Am[i])*t[i]*(1+pow(zS,2)/(24*pow(N[i],2))*(pow(cos(Am[i] ),2)*(2+7*pow(yita[i],2)+9*pow(yita[i],2)*pow(t[i],2)+5*pow(yita[i],4))+pow(sin(Am[i]),2)*(2+pow(t[i],2)+2*pow(yita[i],2)))))/3600);Bm[i+1]=_zB1+1/2*B[i];Lm[i+1]=_zL1+1/2*L[i];Am[i+1]=_zA12+1/2*A[i];M[i+1]=Calc_M(B[i]);N[i+1]=Calc_N(N[i]);yita[i+1]=Calc_yita(yita[i]);t[i+1]=Calc_t(t[i]);final=L[i];i++;}}elsefinal=L[i];double Final=final+_zL1;return _AngleTrans.trans2(Final);}//正算方位角double ZhengFanSuan::ZhengSuanA(double zB1,double zL1,double zA12,double zS){ double M[100],N[100],t[100],B[100],Bm[100],L[100],Lm[100],A[100],Am[100],yita[10000];double_zB1=_AngleTrans.trans1(zB1),_zL1=_AngleTrans.trans1(zL1),_zA12=_AngleTrans.trans1(zA12),final;t[1]=tan(_zB1);yita[1]=pow(e2,2)*pow(cos(_zB1),2);N[1]=a/sqrt(1-pow(e1,2)*pow(sin(_zB1),2));M[1]=a*(1-pow(e1,2))/sqrt(pow(1-pow(e1,2)*pow(sin(_zB1),2),3));B[0]=_AngleTrans.trans3(206265/M[1]*zS*cos(_zA12)/3600),L[0]=_AngleTrans.trans3(206265*zS*sin(_zA12)/(N[1]*cos(_zB1))/3600),A[0]=L[0]*sin(_zB1);Bm[1]=_zB1+1/2*B[0];Lm[1]=_zL1+1/2*L[0];Am[1]=_zA12+1/2*A[0];int i=1;B[1]=_AngleTrans.trans3((206265/M[1]*zS*cos(Am[1])*(1+pow(zS,2)/(24*pow(N[1],2))*(pow(sin(Am[1]),2))*(2+3*pow(t[1],2)+2*pow(yita[1],2))+3*pow(cos(Am[1]),2)*pow(yita[1],2)*(pow(t[1],2)-1-pow(yita[1],2)-4*pow(yita[1],2)*pow(t[1],2))))/3600);L[i]=_AngleTrans.trans3((206265/(N[1]*cos(Bm[1]))*zS*sin(Am[i])*(1+pow(zS,2)/(24*pow(N[i],2))*(pow(t[i] ,2)*pow(sin(Am[i]),2)-pow(cos(Am[i]),2)*(1+pow(yita[i],2)-9*pow(yita[i],2)*pow(t[i],2)))))/3600);A[i]=_AngleTrans.trans3((206265/N[i]*zS*sin(Am[i])*t[i]*(1+pow(zS,2)/(24*pow(N[i],2))*(pow(cos(Am[i]),2 )*(2+7*pow(yita[i],2)+9*pow(yita[i],2)*pow(t[i],2)+5*pow(yita[i],4))+pow(sin(Am[i]),2)*(2+pow(t[i],2)+2*pow(yita[i],2)))))/3600);while(A[i]-A[i-1]>=temp){B[i+1]=_AngleTrans.trans3((206265/M[i]*zS*cos(Am[i])*(1+pow(zS,2)/(24*pow(N[i],2))*(pow(sin(Am[i]),2) )*(2+3*pow(t[i],2)+2*pow(yita[i],2))+3*pow(cos(Am[i]),2)*pow(yita[i],2)*(pow(t[i],2)-1-pow(yita[i],2)-4*pow(yita[i],2)*pow(t[i],2))))/3600);L[i+1]=_AngleTrans.trans3((206265/(N[i]*cos(Bm[i]))*zS*sin(Am[i])*(1+pow(zS,2)/(24*pow(N[i],2))*(pow(t[ i],2)*pow(sin(Am[i]),2)-pow(cos(Am[i]),2)*(1+pow(yita[i],2)-9*pow(yita[i],2)*pow(t[i],2)))))/3600);A[i+1]=_AngleTrans.trans3((206265/N[i]*zS*sin(Am[i])*t[i]*(1+pow(zS,2)/(24*pow(N[i],2))*(pow(cos(Am[i] ),2)*(2+7*pow(yita[i],2)+9*pow(yita[i],2)*pow(t[i],2)+5*pow(yita[i],4))+pow(sin(Am[i]),2)*(2+pow(t[i],2)+2*pow(yita[i],2)))))/3600);Bm[i+1]=_zB1+1/2*B[i];Lm[i+1]=_zL1+1/2*L[i];Am[i+1]=_zA12+1/2*A[i];M[i+1]=Calc_M(B[i]);N[i+1]=Calc_N(N[i]);yita[i+1]=Calc_yita(yita[i]);t[i+1]=Calc_t(t[i]);final=A[i];i++;}double Final;if(_zA12>Pi){ Final=final+_zA12-Pi;}else{ Final=final+_zA12+Pi;}return _AngleTrans.trans2(Final);}//反算Sdouble ZhengFanSuan::FanSuanS(double fB1,double fL1,double fB2,double fL2){ double A12,A21,S;double _fB1=_AngleTrans.trans1(fB1),_fL1=_AngleTrans.trans1(fL1),_fB2=_AngleTrans.trans1(fB2),_fL2=_AngleTrans.trans1(fL2);double Bm=(_fB1+_fB2)/2,_B=_fB2-_fB1,_L=_fL2-_fL1;double Nm=Calc_N(Bm),tm=Calc_t(Bm),yitam=Calc_yita(Bm);double Vm=sqrt(1+pow(e2,2)*pow(cos(Bm),2)),r01=Nm*cos(Bm),r03=(-Nm*pow(cos(Bm),3)*pow(tm,2))/24,r21=(Nm*cos(Bm)*(1+pow(yitam,2)-9*pow(yitam,2)*pow(tm,2)*pow(yitam,4))/pow(Vm,4))/24,s10=Nm/pow(Vm,2),s12=(Nm/pow(Vm,2)*pow(cos(Bm),2)*(2+3*pow(tm,2)+2*pow(yitam,2)))/24,s30=(Nm/pow(Vm,6)*(pow(yitam,2)-pow(yitam,2)*pow(tm,2)))/8,t01=tm*cos(Bm),t21=cos(Bm)*tm/(24*pow(Vm,4))*(2+7*pow(yitam,2)+9*pow(yitam,2)*pow(tm,2)+5*pow(yitam,4)),t03=(tm*pow(cos(Bm),3)*(2+pow(tm,2)+2*pow(yitam,2)))/24;double U=r01*_L+r21*pow(_B,2)*_L+r03*pow(_L,3),V=s10*_B+s12*_B*pow(_L,2)+s30*pow(_B,3),_A=t01*_L+t21*pow(_B,2)*_L+t03*pow(_L,3);double Am=atan(U/V);S=V/cos(Am);return S;}//反算A12double ZhengFanSuan::FanSuanA12(double fB1,double fL1,double fB2,double fL2){ double T,A12,A21,S;double _fB1=_AngleTrans.trans1(fB1),_fL1=_AngleTrans.trans1(fL1),_fB2=_AngleTrans.trans1(fB2),_fL2=_AngleTrans.trans1(fL2);double Bm=(_fB1+_fB2)/2,_B=_fB2-_fB1,_L=_fL2-_fL1;double Nm=Calc_N(Bm),tm=Calc_t(Bm),yitam=Calc_yita(Bm);double Vm=sqrt(1+pow(e2,2)*pow(cos(Bm),2)),r01=Nm*cos(Bm),r03=(-Nm*pow(cos(Bm),3)*pow(tm,2))/24,r21=(Nm*cos(Bm)*(1+pow(yitam,2)-9*pow(yitam,2)*pow(tm,2)*pow(yitam,4))/pow(Vm,4))/24,s10=Nm/pow(Vm,2),s12=(Nm/pow(Vm,2)*pow(cos(Bm),2)*(2+3*pow(tm,2)+2*pow(yitam,2)))/24,s30=(Nm/pow(Vm,6)*(pow(yitam,2)-pow(yitam,2)*pow(tm,2)))/8,t01=tm*cos(Bm),t21=cos(Bm)*tm/(24*pow(Vm,4))*(2+7*pow(yitam,2)+9*pow(yitam,2)*pow(tm,2)+5*pow(yitam,4)),t03=(tm*pow(cos(Bm),3)*(2+pow(tm,2)+2*pow(yitam,2)))/24;double U=r01*_L+r21*pow(_B,2)*_L+r03*pow(_L,3),V=s10*_B+s12*_B*pow(_L,2)+s30*pow(_B,3),_A=t01*_L+t21*pow(_B,2)*_L+t03*pow(_L,3);double Am=atan(U/V);A12=Am-0.5*_A;return _AngleTrans.trans2(A12);}//反算A21double ZhengFanSuan::FanSuanA21(double fB1,double fL1,double fB2,double fL2){ double T,A12,A21,S;double _fB1=_AngleTrans.trans1(fB1),_fL1=_AngleTrans.trans1(fL1),_fB2=_AngleTrans.trans1(fB2),_fL2=_AngleTrans.trans1(fL2);double Bm=(_fB1+_fB2)/2,_B=_fB2-_fB1,_L=_fL2-_fL1;double Nm=Calc_N(Bm),tm=Calc_t(Bm),yitam=Calc_yita(Bm);double Vm=sqrt(1+pow(e2,2)*pow(cos(Bm),2)),r01=Nm*cos(Bm),r03=(-Nm*pow(cos(Bm),3)*pow(tm,2))/24,r21=(Nm*cos(Bm)*(1+pow(yitam,2)-9*pow(yitam,2)*pow(tm,2)*pow(yitam,4))/pow(Vm,4))/24,s10=Nm/pow(Vm,2),s12=(Nm/pow(Vm,2)*pow(cos(Bm),2)*(2+3*pow(tm,2)+2*pow(yitam,2)))/24,s30=(Nm/pow(Vm,6)*(pow(yitam,2)-pow(yitam,2)*pow(tm,2)))/8,t01=tm*cos(Bm),t21=cos(Bm)*tm/(24*pow(Vm,4))*(2+7*pow(yitam,2)+9*pow(yitam,2)*pow(tm,2)+5*pow(yitam,4)),t03=(tm*pow(cos(Bm),3)*(2+pow(tm,2)+2*pow(yitam,2)))/24;double U=r01*_L+r21*pow(_B,2)*_L+r03*pow(_L,3),V=s10*_B+s12*_B*pow(_L,2)+s30*pow(_B,3),_A=t01*_L+t21*pow(_B,2)*_L+t03*pow(_L,3);double Am=atan(U/V);A12=Am-0.5*_A;if(A12>Pi){ A21=Am+0.5*_A-Pi;}else{ A21=Am+0.5*_A+Pi;}return _AngleTrans.trans2(A21);}4.4、正算的计算按钮代码:void C主题高斯引数正反算Dlg::OnBnClickedButton2(){ UpdateData(true);ZhengFanSuan _ZhengFanSuan;double B1=_wtof(zB1),L1=_wtof(zL1),A12=_wtof(zA12),S=_wtof(zS),B2=_ZhengFanSuan.ZhengSuanB(B1,L1,A12,S),L2=_ZhengFanSuan.ZhengSuanL(B1,L1,A12,S),A21=_ZhengFanSuan.ZhengSuanA(B1,L1,A12,S);zB2.Format(_T("%.7f"),B2);zL2.Format(_T("%.7f"),L2),zA21.Format(_T("%.7f"),A21);UpdateData(false);}4.5、反算的计算按钮代码:void C主题高斯引数正反算Dlg::OnBnClickedButton1(){ UpdateData(true);ZhengFanSuan _ZhengFanSuan;double B1=_wtof(fB1),L1=_wtof(fL1),B2=_wtof(fB2),L2=_wtof(fL2),S=_ZhengFanSuan.FanSuanS(B1,L1,B2,L2),A12=_ZhengFanSuan.FanSuanA12(B1,L1,B2,L2),A21=_ZhengFanSuan.FanSuanA21(B1,L1,B2,L2);fS.Format(_T("%.7f"),S);fA12.Format(_T("%.7f"),A12);fA21.Format(_T("%.7f"),A21);UpdateData(false);}4.6、清零按钮代码:void C主题高斯引数正反算Dlg::OnBnClickedButton3(){ UpdateData(true);fB1.Format(_T("%.0f"),0);fL1.Format(_T("%.0f"),0);fB2.Format(_T("%.0f"),0);fL2.Format(_T("%.0f"),0);fA12.Format(_T("%.0f"),0);fA21.Format(_T("%.0f"),0);fS.Format(_T("%.0f"),0);UpdateData(false);}void C主题高斯引数正反算Dlg::OnBnClickedButton4(){ UpdateData(true);zB1.Format(_T("%.0f"),0);zL1.Format(_T("%.0f"),0);zB2.Format(_T("%.0f"),0);zL2.Format(_T("%.0f"),0);zA12.Format(_T("%.0f"),0);zA21.Format(_T("%.0f"),0);zS.Format(_T("%.0f"),0);UpdateData(false);}五、计算结果数据组号:18六、实习体会实验一开始本来考虑用控制台做,后来发现用控制台做出来可视化效果不好,而且程序中要用很多cin和cout,写代码十分麻烦,所以最后选择了用MFC做。