高斯投影坐标正反算编程报告
高斯投影正反算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)截图.。
高斯投影正反算python代码
高斯投影正反算1. 什么是高斯投影高斯投影是一种常用的地图投影方法，它将地球表面的经纬度坐标转换为平面坐标，常用于地理信息系统（GIS）和测绘工程中。
高斯投影分为正算和反算两个过程。
•正算：将经纬度坐标转换为平面坐标。
•反算：将平面坐标转换为经纬度坐标。
2. 高斯投影正算2.1 原理高斯投影正算的原理是根据椭球体上某一点处的曲率半径、子午线弧长和东西方向上的距离，计算该点在平面上的x、y坐标。
2.2 具体步骤高斯投影正算的具体步骤如下：1.根据给定的椭球体参数（长半轴a、短半轴b），计算椭球体第一偏心率e。
2.根据给定的中央子午线经度λ0，计算λ - λ0 的差值Δλ。
3.计算曲率半径N和子午线弧长A0。
4.根据给定的纬度φ和经度λ，计算Δφ和Δλ。
5.计算子午线弧长A1、A2、A3和A4。
6.计算平面坐标x和y。
2.3 Python代码实现下面是使用Python实现高斯投影正算的示例代码：import math# 输入参数a = 6378137.0 # 长半轴b = 6356752.314245 # 短半轴e = math.sqrt(1 - (b/a)**2) # 第一偏心率λ0 = math.radians(120) # 中央子午线经度，单位为弧度# 输入经纬度坐标φ = math.radians(30) # 纬度，单位为弧度λ = math.radians(121) # 经度，单位为弧度# 计算ΔλΔλ = λ - λ0# 计算曲率半径N和子午线弧长A0N = a / math.sqrt(1 - e**2 * math.sin(φ)**2)A0 = a * (1 - e**2) / (1 - e**2 * math.sin(φ)**2)**1.5# 计算Δφ和ΔλΔφ = φ - φ0# 计算子午线弧长A1、A2、A3和A4A1 = A0 + N * math.tan(φ) / 2 * Δλ**2 * math.cos(φ)A2 = A0 + N * math.tan(φ) / 24 * (5 - math.tan(φ)**2 + 9 * e2 * math.cos(φ)**2 + 4 * e2**2 * math.cos(φ)**4) * Δλ**4 * math.cos(φ)A3 = A0 + N * math.tan(φ) / 720 * (61 - 58 * math.tan(φ)**2 + math.tan(φ)** 4) * Δλ**6 * math.cos(φ)A4 = A0 + N * math.tan(φ) / 40320 * (1385 - 3111*math.tan(φ)**2 + 543*math.tan(φ)**4 - math.tan(φ)**6) \* Δλ**8 \* math.cos(phi)# 计算平面坐标x和yx = A1 + A2 + A3 + A4y = N / math.cos(phi) \(Δλ - Δλ**3/6*(1+math.tan(phi))**2/N/A0^2 \+ Δλ^5/120*(5+28*math.tan(phi)^2+24*math.tan(phi)^4)*N/A0^4/N/A0^3)# 输出结果print("平面坐标(x, y):", x, y)3. 高斯投影反算3.1 原理高斯投影反算的原理是根据平面坐标和中央子午线经度，计算对应的经纬度坐标。
(完整word版)高斯投影正反算 代码
L=l+L01; ///
反算就出
L B
l=L-L02;
B=Bf-(1-(b4-0.12*Z*Z)*Z*Z)*Z*Z*b2;
t=tan(B);
n2=e2*cos(B)*cos(B);
V=sqrt(1+n2);
cB2=pow(cos(B),2);
N=6399698.902-(21562.267-(108.973-0.612*cB2)*cB2)*cB2; //N=c/V;
l=L-111*3600/P // l=((m%6)*3600+n*60+h)/P;
t=tan(B);
n2=e2*cos(B)*cos(B);
V=sqrt(1+n2);
cB2=pow(cos(B),2);
N=6399698.902-(21562.267-(108.973-0.612*cB2)*cB2)*cB2; // N=c/V;
a0=32140.404-(135.3302-(0.7092-0.004*cB2)*cB2)*cB2;
a4=(0.25+0.00252*cB2)*cB2-0.04166;
a6=(0.166*cB2-0.084)*cB2;
a3=(0.3333333+0.001123*cB2)*cB2-0.1666667;
#include "stdafx.h"
#include "iostream.h"
#include "math.h"
#include "stdio.h"
#define P 206264.806247096355
高斯投影正反算c代码
高斯投影正反算c代码 Coca-cola standardization office【ZZ5AB-ZZSYT-ZZ2C-ZZ682T-ZZT18】高斯投影正反算程序设计一．程序设计流程本程序的设计思路如下：（1），程序采用VS08版本作为开发平台，并采用C#语言作为开发语言，设计为WindowsForm窗体程序形式。
（2），程序主要的算法来自于教材。
但是本程序为了更加实用，添加了更多的解算基准，包括：WGS-84，国际椭球1975，克氏椭球，和2000国家大地坐标系。
（3），程序为了更方便的读取数据和输出数据，故需要自己定义了固定的数据输入格式和数据输出格式或形式，请老师注意查看。
二．代码using System;using ;using ;using ;using ;using Gauss{public partial class Form1 : Form{double b = (a * a * (1 - ee * ee));double c = a * a / b;double epp = ((a * a - b * b) / b / b);CRDGEODETIC pcrdGeo;CRDCARTESIAN pcrdCar;double midlong = 0;//求X，Y和带号= ;ytext = ;string temp = (0, 2);num = (temp);ytext = (0, 2);= (ytext) - 500000;try{tt = }catch{("Choose 3/6 error!");return;}if ("3度带") == 0){midlong = num * 3 * pai / 180;}if ("6度带") == 0){midlong = (6 * num - 3) * pai / 180; }b = (a * a * (1 - ee * ee));c = a * a / b;epp = (a * a - b * b) / b;double m0, m2, m4, m6, m8;double a0, a2, a4, a6, a8;m0 = a * (1 - ee * ee);m2 = / * m0 * ee * ee;m4 = / * ee * ee * m2;m6 = / * ee * ee * m4;m8 = / * ee * ee * m6;a0 = m0 + m2 / + / * m4 + / * m6 + / * m8;a2 = m2 / 2 + m4 / 2 + / * m6 + / * m8;a4 = m4 / + / * m6 + / * m8;a6 = m6 / + m8 / ;a8 = m8 / ;double Bf, B;Bf = / a0;B = ;while (Bf - B) > 1E-10){B = Bf;double sb = (B);double cb = (B);double s2b = sb * cb * 2;double s4b = s2b * (1 - 2 * sb * sb) * 2;double s6b = s2b * (1 - s4b * s4b) + s4b * (1 - s2b * s2b); Bf = - (-a2 / * s2b + a4 / * s4b - a6 / * s6b)) / a0; }double itaf, tf, Vf, Nf;itaf = epp * (Bf);tf = (Bf);Vf = (1 + epp * epp * (Bf) * (Bf));Nf = c / Vf;double ynf = / Nf;= Bf - / * Vf * Vf * tf * (ynf * ynf - / * (ynf, 4) * (5 + 3 * tf * tf + itaf * itaf - 9 * (itaf * tf, 2)) +/ * (61 + 90 * tf * tf + 45 * (tf, 4)) * (ynf, 6));= (ynf / (Bf) - (1 + 2 * tf * tf + itaf * itaf) * (ynf, 3) / / (Bf) +(5 + 28 * tf * tf + 24 * (tf, 4) + 6 * itaf * itaf + 8 * (itaf * tf, 2)) * (ynf, 5) / / (Bf));= + midlong;// = ;= "Results:\nLatitude: " + + "\nLongtitude: " + ;}private void label13_Click(object sender, EventArgs e){}}}三．程序运行结果分析通过选取书上的具体实例进行测试，本程序的精度大体满足要求，一般正算的精度在米和米之间，反算的精度在秒左右。
高斯投影正算与反算的理论方法与实
高斯投影正算与反算的理论方法与实现代码高斯投影是正形投影的一种，同一坐标系中的高斯投影换带计算公式是根据正形投影原理推导出的两个高斯坐标系间的显函数式。
在同一大地坐标系中(例如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表示。
高斯正反算实验报告(王震阳20094176)
实验2：高斯正反算以及换带程序的实现姓名：王震阳学号：20094176班级：测绘09-1班指导老师：周志易完成时间：2011.08.11实验目的和要求1．了解高斯正反算的基本思想。
2．完成高斯正反算点算程序的实现过程。
2实验环境和工具通过Windows xp 系统运行Vc 6.0 软件。
3实验结果3.1程序代码1．正反算程序char cb[100],cl[100],cl0[100],cx[100],cy[100];char cbdu[100],cbfen[100],cbmiao[100];char cldu[100],clfen[100],clmiao[100];char cl0du[100],cl0fen[100],cl0miao[100];double db=0,dl=0,dl0=0,dx=0,dy=0;double dbdu=0,dbfen=0,dbmiao=0;double dldu=0,dlfen=0,dlmiao=0;double dl0du=0,dl0fen=0,dl0miao=0;double p2=206264.8062;char cfr[100];char cfb[100],cfl[100];double dfb=0,dfl=0;double dfbdu=0,dfbfen=0,dfbmiao=0;double dfldu=0,dflfen=0,dflmiao;char cfbdu[100],cfbfen[100],cfbmiao[100];char cfldu[100],cflfen[100],cflmiao[100];m_b.GetWindowText(cb,100);for(int i=0;cb[i]!=' ';i++){if(cb[i]!=' ')cbdu[i]=cb[i];}cbdu[i]='\0';dbdu=atof((LPCTSTR)cbdu);/////for(int j=0;cb[j+strlen(cbdu)+1]!=' ';j++){cbfen[j]=cb[j+strlen(cbdu)+1];}cbfen[j]='\0';dbfen=atof((LPCTSTR)cbfen);/////for(int k=0;cb[k+strlen(cbdu)+strlen(cbfen)+2]!=' ';k++){cbmiao[k]=cb[k+strlen(cbdu)+strlen(cbfen)+2];}cbmiao[k]='\0';dbmiao=atof((LPCTSTR)cbmiao);////db=(dbdu*3600+dbfen*60+dbmiao)/206264.8062;//////纬度的弧度制////////////////////////////////////////////////////////////////////////////////获取经度L的度分秒、、、、//// m_l.GetWindowText(cl,100);for( i=0;cl[i]!=' ';i++){if(cl[i]!=' ')cldu[i]=cl[i];}cldu[i]='\0';dldu=atof((LPCTSTR)cldu);/////for( j=0;cl[j+strlen(cldu)+1]!=' ';j++){clfen[j]=cl[j+strlen(cldu)+1];}clfen[j]='\0';dlfen=atof((LPCTSTR)clfen);/////for( k=0;cl[k+strlen(cldu)+strlen(clfen)+2]!=' ';k++){clmiao[k]=cl[k+strlen(cldu)+strlen(clfen)+2];}clmiao[k]='\0';dlmiao=atof((LPCTSTR)clmiao);/////dl=(dldu*3600+dlfen*60+dlmiao)/206264.8062;//////经度的弧度制//////////////////////////////////////////////////////////////////////////获取经度L0的度分秒、、、、/////// m_l0.GetWindowText(cl0,100);double dl1=atof((LPCTSTR)cl0);/////dl0=dl1*3600/206264.8062;//////中央子午线的弧度制///////////////////////////////////////////////////////////////////////////////////////////////char czx[100],czy[100],czr[100];double dzx=0.,dzy=0.,dzr=0.;double l=dl-dl0;doublen=6399698.902-(21562.267-(108.973-0.612*cos(db)*cos(db))*cos(db)*cos(db))*cos(db)*cos(db);a0=32140.404-(135.3302-(0.7092-0.0040*cos(db)*cos(db))*cos(db)*cos(db))*cos(db)*cos(db); double a4=(0.25+0.00252*cos(db)*cos(db))*cos(db)*cos(db)-0.04166;double a6=(0.166*cos(db)*cos(db)-0.084)*cos(db)*cos(db);double a3=(0.3333333+0.001123*cos(db)*cos(db))*cos(db)*cos(db)-0.1666667;double a5=0.0083-(0.1667-(0.1968+0.004*cos(db)*cos(db))*cos(db)*cos(db))*cos(db)*cos(db); dzx=6367558.4969*db-(a0-(0.5+(a4+a6*l*l)*l*l)*l*l*n)*sin(db)*cos(db);dzy=(1.+(a3+a5*l*l)*l*l)*l*n*cos(db);_gcvt(dzx,14,czx);m_zx=czx;_gcvt(dzy,14,czy);m_zy=czy;///////////double zr1=0.33333+0.00674*cos(db)*cos(db);double zr2=(0.2*cos(db)*cos(db)-0.0067)*l*l;double dzrdu,dzrfen,dzrmiao;char czrdu[100],czrfen[100],czrmiao[100],mczr1[100];CString mczr2,mczr3,mczr4,mczr5,mczr6,mczr7,mczr8;dzr=l*sin(db)*(1+(zr1+zr2)*l*l*cos(db)*cos(db));dzrdu=(int)(dzr*p2/3600.0);dzrfen=(int)((dzr*p2-dzrdu*3600)/60.0);dzrmiao=(dzr*p2-dzrdu*3600-dzrfen*60);_gcvt(dzrdu,5,czrdu);_gcvt(dzrfen,5,czrfen);_gcvt(dzrmiao,4,czrmiao);mczr6=' ';mczr2=czrfen;mczr3=czrdu;mczr4=czrmiao;mczr5=mczr3+mczr6+mczr2+mczr6+mczr4;m_zr=mczr5;///////////////////////////反算程序//////CString cfdu1,cffen1,cfmiao1,cff1;CString cfdu2,cffen2,cfmiao2,cff2;CString mcfr4,mcfr2,mcfr3,mcfr5,mcfr6;m_x.GetWindowText(cx,100);m_y.GetWindowText(cy,100);dx=atof((LPCTSTR)cx);dy=atof((LPCTSTR)cy);double bter=dx/6367558.4969;double bf1=(293622+(2350+22*cos(bter)*cos(bter))*cos(bter)*cos(bter))*cos(bter)*cos(bter); doublebf=bter+((50221746.0+bf1)*sin(bter)*cos(bter))/(10.0*10.0*10.0*10.0*10.0*10.0*10.0*10.0*10.0*double nf1=21562.267-(108.973-0.612*cos(bf)*cos(bf))*cos(bf)*cos(bf); double nf=6399698.902-nf1*cos(bf)*cos(bf);double z=dy/(nf*cos(bf));double b2=(0.5+0.003369*cos(bf)*cos(bf))*sin(bf)*cos(bf);double b3=0.333333-(0.166667-0.001123*cos(bf)*cos(bf))*cos(bf)*cos(bf); double b4=0.25+(0.16161+0.00562*cos(bf)*cos(bf))*cos(bf)*cos(bf); double b5=0.2-(0.1667-0.0088*cos(bf)*cos(bf))*cos(bf)*cos(bf);dfb=bf-(1-(b4-0.12*z*z)*z*z)*z*z*b2;dfl=dl0+(1-(b3-b5*z*z)*z*z)*z;dfbdu=(int)(dfb*p2/3600);dfbfen=(int)((dfb*p2-dfbdu*3600)/60);dfbmiao=dfb*p2-dfbdu*3600-dfbfen*60;_gcvt(dfbdu,15,cfbdu);cfdu1=cfbdu;_gcvt(dfbfen,15,cfbfen);cffen1=cfbfen;_gcvt(dfbmiao,15,cfbmiao);cfmiao1=cfbmiao;cff1=cfdu1+mczr6+cffen1+mczr6+cfmiao1;m_fb=cff1;dfldu=(int)(dfl*p2/3600);dflfen=(int)((dfl*p2-dfldu*3600)/60);dflmiao=dfl*p2-dfldu*3600-dflfen*60;_gcvt(dfldu,10,cfldu);cfdu2=cfldu;_gcvt(dflfen,10,cflfen);cffen2=cflfen;_gcvt(dflmiao,10,cflmiao);cfmiao2=cflmiao;cff2=cfdu2+mczr6+cffen2+mczr6+cfmiao2;m_fl=cff2;//////计算反算子午收敛角du///double dfr;double l1=dl0-dfl;double fr1=0.33333+0.00674*cos(dfb)*cos(dfb);double fr2=(0.2*cos(dfb)*cos(dfb)-0.0067)*l1*l1;double dfrdu,dfrfen,dfrmiao;char cfrdu[100],cfrfen[100],cfrmiao[100],mcfr1[100];dfr=l1*sin(dfb)*(1+(fr1+fr2)*l1*l1*cos(dfb)*cos(dfb));dfrdu=(int)(dfr*p2/3600.0);dfrfen=(int)((dfr*p2-dfrdu*3600)/60.0);dfrmiao=(dfr*p2-dfrdu*3600-dfrfen*60);_gcvt(dfrdu,5,cfrdu);_gcvt(dfrfen,5,cfrfen);_gcvt(dfrmiao,4,cfrmiao);mcfr6=' ';mcfr2=czrfen;mcfr3=czrdu;mcfr4=czrmiao;mcfr5=mcfr3+mcfr6+mcfr2+mcfr6+mcfr4;m_fr=mcfr5;UpdateData(false);}///////////////////////////////////换带程序char cx[100],cy[100];char cl0[100];char cl0n[100];double dx=0.,dy=0.;double dl0=0.,dl0n=0.;m_x1.GetWindowText(cx,100);m_y1.GetWindowText(cy,100);dx=atof((LPCTSTR)cx);dy=atof((LPCTSTR)cy);m_l0.GetWindowText(cl0,100);dl0=atof((LPCTSTR)cl0);m_l0n.GetWindowText(cl0n,100);dl0n=atof((LPCTSTR)cl0n);double bter=dx/6367558.4969;double bf1=(293622+(2350+22*cos(bter)*cos(bter))*cos(bter)*cos(bter))*cos(bter)*cos(bter);double bf=bter+((50221746.0+bf1)*sin(bter)*cos(bter))/(10.0*10.0*10.0*10.0*10.0*10.0*10.0*10.0*10.0*10.0); double nf1=21562.267-(108.973-0.612*cos(bf)*cos(bf))*cos(bf)*cos(bf);double nf=6399698.902-nf1*cos(bf)*cos(bf);double z=dy/(nf*cos(bf));double b2=(0.5+0.003369*cos(bf)*cos(bf))*sin(bf)*cos(bf);double b3=0.333333-(0.166667-0.001123*cos(bf)*cos(bf))*cos(bf)*cos(bf);double b4=0.25+(0.16161+0.00562*cos(bf)*cos(bf))*cos(bf)*cos(bf);double b5=0.2-(0.1667-0.0088*cos(bf)*cos(bf))*cos(bf)*cos(bf);double dfb=0,dfl=0.;dfb=bf-(1-(b4-0.12*z*z)*z*z)*z*z*b2;CString cfdu1,cffen1,cfmiao1,cff1;CString mczr6=' ';double dfbdu=(int)(dfb*p2/3600);double dfbfen=(int)((dfb*p2-dfbdu*3600)/60);double dfbmiao=dfb*p2-dfbdu*3600-dfbfen*60;char cfbdu[100],cfbfen[100],cfbmiao[100];_gcvt(dfbdu,15,cfbdu);cfdu1=cfbdu;_gcvt(dfbfen,15,cfbfen);cffen1=cfbfen;_gcvt(dfbmiao,15,cfbmiao);cfmiao1=cfbmiao;cff1=cfdu1+mczr6+cffen1+mczr6+cfmiao1;m_b=cff1;/////////dl0=dl0*3600/206264.8062;dfl=dl0+(1-(b3-b5*z*z)*z*z)*z;double dfldu=(int)(dfl*p2/3600);double dflfen=(int)((dfl*p2-dfldu*3600)/60);double dflmiao=dfl*p2-dfldu*3600-dflfen*60;char cfldu[100],cflfen[100],cflmiao[100];CString cfdu2,cffen2,cfmiao2,cff2;_gcvt(dfldu,100,cfldu);cfdu2=cfldu;_gcvt(dflfen,100,cflfen);cffen2=cflfen;_gcvt(dflmiao,100,cflmiao);cfmiao2=cflmiao;cff2=cfdu2+mczr6+cffen2+mczr6+cfmiao2;m_l=cff2;//////UpdateData(false);double l=dfl-dl0n*3600/p2;double n=6399698.902-(21562.267-(108.973-0.612*cos(dfb)*cos(dfb))*cos(dfb)*cos(dfb))*cos(dfb)*cos(dfb); double a0=32140.404-(135.3302-(0.7092-0.0040*cos(dfb)*cos(dfb))*cos(dfb)*cos(dfb))*cos(dfb)*cos(dfb); double a4=(0.25+0.00252*cos(dfb)*cos(dfb))*cos(dfb)*cos(dfb)-0.04166;double a6=(0.166*cos(dfb)*cos(dfb)-0.084)*cos(dfb)*cos(dfb);double a3=(0.3333333+0.001123*cos(dfb)*cos(dfb))*cos(dfb)*cos(dfb)-0.1666667;double a5=0.0083-(0.1667-(0.1968+0.004*cos(dfb)*cos(dfb))*cos(dfb)*cos(dfb))*cos(dfb)*cos(dfb)char cx2[100];char cy2[100];double dx2=6367558.4969*dfb-(a0-(0.5+(a4+a6*l*l)*l*l)*l*l*n)*sin(dfb)*cos(dfb);double dy2=(1.+(a3+a5*l*l)*l*l)*l*n*cos(dfb);_gcvt(dx2,14,cx2);m_x2=cx2;_gcvt(dy2,14,cy2);m_y2=cy2;UpdateData(false);《计算机图形学》实验报告//////double zr1=0.33333+0.00674*cos(dfb)*cos(dfb);double zr2=(0.2*cos(dfb)*cos(dfb)-0.0067)*l*l;double dzrdu,dzrfen,dzrmiao;char czrdu[100],czrfen[100],czrmiao[100],mczr1[100];CString mczr2,mczr3,mczr4,mczr5,mczr7,mczr8;double dzr=l*sin(dfb)*(1+(zr1+zr2)*l*l*cos(dfb)*cos(dfb));dzrdu=(int)(dzr*p2/3600.0);dzrfen=(int)((dzr*p2-dzrdu*3600)/60.0);dzrmiao=(dzr*p2-dzrdu*3600-dzrfen*60);_gcvt(dzrdu,5,czrdu);_gcvt(dzrfen,5,czrfen);_gcvt(dzrmiao,4,czrmiao);mczr6=' ';mczr2=czrfen;mczr3=czrdu;mczr4=czrmiao;mczr5=mczr3+mczr6+mczr2+mczr6+mczr4;m_r=mczr5;UpdateData(false);3.2运行结果分析得到预期结果，经过和答案比对完全正确，程序正确。
测绘程序设计(VS2008)实验报告--高斯投影正反算
《测绘程序设计（）》上机实验报告（Visual C++.Net）班级：学号：姓名：序号：二零一一年五月实验7 常用测量程序设计1.实验目的：1.1 巩固类的创建与使用;1.2 掌握数组参数的传递；1.3 掌握常用测绘程序设计的技巧。
1.42.实验内容：编写高斯投影正反算程序。
3.设计思路：这次的实验目的是实现高斯正反算。
需要考虑投影方式即分带的方式，又要考虑椭球参数的类型，所以我添加了两个函数来完成此功能。
分别是int SetProjectType(int m)和void SetParameter(int m,double &a,double &b)。
4.界面设计：界面设计很简单，具体见运行结果。
5.主要代码：文件名：GaussProjectDlg.cpp代码：const double PI=4*atan(1.0);//获得分带方式返回中央子午线经度int CGaussProjectDlg::SetProjectType(int m){UpdateData(TRUE);int n; //记录分带带号double L; //经度L=iDegreeL+iMinL/60+dSecondL/3600;if(m==1) //6度带{n=int(L/6)+1;L0=6*n-3;}else if(m==2) //3度带{n=int((L+1.5)/3);L0=3*n;}else if(m==3) //自主分带L0=L0;return L0;}//获取椭球参数void CGaussProjectDlg::SetParameter(int m,double &a,double &b) {if(m==1) //克拉索夫斯基椭球{a=6378245.0;b=6356863.0187730473;//e=sqrt(0.006693421622966);}else if(m==2) //1975国际协议椭球{a=6378140.0;b=6356755.2881575287;//e=sqrt(0.006694384999588);}else if(m==3) //WGS-84椭球{a=6378137.0;b=6356752.3142;//e=sqrt(0.0066943799013);}}void CGaussProjectDlg::OnBnClickedButtonpositivecal(){// TODO: 在此添加控件通知处理程序代码UpdateData(TRUE);double N;double t;double Eta;double X;double A0,A2,A4,A6,A8;double RadB;double Rou;Rou=180*3600/PI;double a,b,e1,e2; //椭球参数SetParameter(iParameterType,a,b);e1=sqrt(a*a-b*b)/a;e2=sqrt(a*a-b*b)/b;double l;L0=SetProjectType(iProjectType);double L;L=iDegreeL+double(iMinL)/60+dSecondL/3600;l=(L-L0)*3600;RadB=(iDegreeB+double(iMinB)/60+dSecondB/3600)*PI/180;N=a/sqrt(1-e1*e1*sin(RadB)*sin(RadB));t=tan(RadB);Eta=e2*cos(RadB);A0=1+3.0/4*e1*e1+45.0/64*pow(e1,4)+350.0/512*pow(e1,6)+11025.0/16384*pow(e1,8);A2=-1.0/2*(3.0/4*e1*e1+60.0/64*pow(e1,4)+525.0/512*pow(e1,6)+17640.0/16384*pow(e1,8));A4=1.0/4*(15.0/64*pow(e1,4)+210.0/512*pow(e1,6)+8820.0/16384*pow(e1,8));A6=-1.0/6*(35.0/512*pow(e1,6)+2520.0/16384*pow(e1,8));A8=1.0/8*(315.0/16384*pow(e1,8));X=a*(1-e1*e1)*(A0*RadB+A2*sin(2*RadB)+A4*sin(4*RadB)+A6*sin(6*RadB)+A8*sin(8*RadB));x=X+N/(2*Rou*Rou)*sin(RadB)*cos(RadB)*l*l+N/(24*pow(Rou,4))*sin(RadB)*pow(cos(RadB),3)*(5-t*t+9*Eta*Eta+4*pow(Eta,4))*pow(l,4)+ N/(720*pow(Rou,6))*sin(RadB)*pow(cos(RadB),5)*(61-58*t*t+pow(t,4))*pow(l,6);y=N/Rou*cos(RadB)*l+N/(6*pow(Rou,3))*pow(cos(RadB),3)*(1-t*t+Eta*Eta)*pow(l,3)+N/(120*pow(Rou,5))*pow(cos(RadB),5)*(5-18*t*t+pow(t,4)+14*Eta*Eta-58*Eta*Eta*t*t)*pow( l,5);UpdateData(FALSE);}void CGaussProjectDlg::OnBnClickedButtonantical(){// TODO: 在此添加控件通知处理程序代码UpdateData(TRUE);double t_f;double Eta_f;double B_f;double N_f;double M_f;double X=x;double B0;double K0,K2,K4,K6;double a,b,e1,e2; //椭球参数SetParameter(iParameterType,a,b);e1=sqrt(a*a-b*b)/a;e2=sqrt(a*a-b*b)/b;double A0;A0=1+3.0/4*e1*e1+45.0/64*pow(e1,4)+350.0/512*pow(e1,6)+11025.0/16384*pow(e1,8);B0=X/(a*(1-e1*e1)*A0);K0=1.0/2*(3.0/4*e1*e1+45.0/64*pow(e1,4)+350.0/512*pow(e1,6)+11025.0/16384*pow(e1,8));K2=-1.0/3*(63.0/64*pow(e1,4)+1108.0/512*pow(e1,6)+58239.0/16384*pow(e1,8));K4=1.0/3*(604.0/512*pow(e1,6)+68484.0/16384*pow(e1,8));K6=-1.0/3*(26328.0/16384*pow(e1,8));B_f=B0+sin(2*B0)*(K0+sin(B0)*sin(B0)*(K4+K6*sin(B0)*sin(B0)));t_f=tan(B_f);Eta_f=e2*cos(B_f);N_f=a/sqrt(1-e1*e1*sin(B_f)*sin(B_f));M_f=N_f/(1+e2*e2*cos(B_f)*cos(B_f));double B;B=B_f-t_f/(2*M_f*N_f)*y*y+t_f/(24*M_f*pow(N_f,3))*(5+3*t_f*t_f+Eta_f*Eta_f-9*Eta_f*Eta_f*t_f*t_f)*pow(y,4)- t_f/(720*M_f*pow(N_f,5))*(61+90*t_f*t_f+45*pow(t_f,4))*pow(y,6);double l;l=1.0/(N_f*cos(B_f))*y-1.0/(6*pow(N_f,3)*cos(B_f))*(1+2*t_f*t_f+Eta_f*Eta_f)*pow(y,3)+1.0/(120*pow(N_f,5)*cos(B_f))*(5+28*t_f*t_f+24*pow(t_f,4)+6*Eta_f*Eta_f+8*Eta_f*Eta_f* t_f*t_f)*pow(y,5);//将B转化为度分秒的形式double dDegB;dDegB=B*180/PI;iDegreeB=int(dDegB);iMinB=int((dDegB-iDegreeB)*60);dSecondB=((dDegB-iDegreeB)*60-iMinB)*60;double dDegL;dDegL=l*180/PI+L0;iDegreeL=int(dDegL);iMinL=int((dDegL-iDegreeL)*60);dSecondL=((dDegL-iDegreeL)*60-iMinL)*60;UpdateData(FALSE);}6.运行结果：实验的运行结果如下图所示：正算：反算：7.实验总结这次实验是实现高斯投影的正反算，方法很多，实现并不复杂，但是计算公式复杂，变量繁多，稍有不慎，就会造成计算错误。
高斯投影正反算实习
大地测量学编程实习报告姓名：鲁尼学号：10 班级：曼联编程思想:这个投影是由德国数学家、物理学家、天文学家高斯于19 世纪20 年代拟定，后经德国大地测量学家克吕格于1912 年对投影公式加以补充，故称为高斯-克吕格投影。
即等角横切椭圆柱投影。
假想用一个圆柱横切于地球椭球体的某一经线上，这条与圆柱面相切的经线，称中央经线。
以中央经线为投影的对称轴，将东西各3°或1°30′的两条子午线所夹经差6°或3°的带状地区按数学法则、投影法则投影到圆柱面上，再展开成平面，即高斯-克吕格投影，简称高斯投影。
这个狭长的带状的经纬线网叫做高斯-克吕格投影带。
高斯投影正算公式就是由大地坐标（L，B）求解高斯平面坐标（x，y），而高斯投影反算公式则是由高斯平面坐标（x，y）求解大地坐标（L，B）。
现行的高斯投影用表都是采用克拉索夫斯基椭球参数，这次编程计算就是采用这种椭球参数，并采用实用公式按6度分带投影。
编程环境是在VC下，采用C++语言编写。
程序主要分为两部分，第一部分是高斯正反算函数，第二部分是主函数。
高斯正反算函数，参考书上175页的电算公式，正算时先将度数换算成秒，再定带号n，求中央经线l0，经度差l''，然后根据实用公式计算高斯平面坐标x，y。
最后计算国家统一坐标的x，y，再将其输出。
计算和数据模型：正算是指：由大地坐标（L,B）求得高斯平面坐标（x,y）的过程。
反算是指：由高斯平面坐标（x,y）求得大地坐标（L,B）的过程。
正算：高斯投影必须满足的三个条件：（1），中央子午线投影后为直线。
（2），中央子午线投影后长度不变。
（3），投影具有正性性质，即正性投影条件。
由第一个条件可知，中央子午线东西两侧的投影必然对称于中央子午线。
设在托球面上有P1 ，P2，且对称于中央子午线。
其大地坐标为（l,B）,(-l,B)则投影后的平面坐标一定为P1·（x,y）,P2·（x,-y）.由第二个条件可知，位于中央子午线上的点，投影后的纵坐标x应该等于投影前从赤道量至该点的子午弧长。
matlab大地测量高斯投影正反算程序设计实验
N sin( B) cos 5 ( B)(61 58t 2 t 4 270 2 330 2 t 2 )l 6 720 N N y N cos( B)l cos 3 ( B)(1 t 2 2 )l 3 cos 5 ( B)(5 18t 2 t 4 14 2 58 2 t 2 )l 5 6 120
-2-
degree = fix(jiaodu); mimute = fix((jiaodu-degree)*100); second = (jiaodu-degree-mimute/100)*10000; degree = degree+mimute/60+second/3600; end 度转度分秒函数如下： function dms=degree3dms(du) degree=fix(du); min=fix((du-degree)*60); second=(((du-degree)*60-min)*60); dms=degree+min/100+second/10000; end 度分秒转弧度 dms2rad 函数如下： function rad=dms2rad(n) deg=fix(n); min_tem=(n-deg)*100; min=fix(min_tem); sec=(min_tem-min)*100; %度取所给数 n 的整数部分 %去掉整数部分后小数点后移两位 %分取整 %秒是小数点再向后移两位的数字 %度－>度分秒 (dd.mmss)
-1-
%计算高斯平面坐标（x，y） x=X+(N.*(sin(B)).*(cos(B)).*(l.^2))./2+(N.*(sin(B)).*((cos(B)).^3).*(5-((t).^2)+9.* (p.^2)+4.*(p.^4)).*(l.^4))./24 ... +(N.*sin(B).*(cos(B)).^5.*(61-58.*(t.^2)+(t.^4)+270.*(p.^2)-330.*(p.^2).*(t.^2)).*(l .^6))./720; y=N.*cos(B).*l+(N.*(cos(B)).^3.*(1-t.^2+p.^2).*(l.^3))./6+(N.*((cos(B)).^5).*(5-18 .*(t.^2)+ ... t.^4+14.*(p.^2)-58.*(p.^2).*(t.^2)).*(l.^5))./120; L0=degree3dms(L0); end latitude2meridian 函数如下： function x = latitude2meridian(B,a,e) %由纬度 B 求对应的子午线弧长 x，计算公式 m0=a*(1-e^2); m2=(3*(e^2)*m0)/2; m4=(5*(e^2)*m2)/4; m6=(7*(e^2)*m4)/6; m8=(9*(e^2)*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; x=a0*B-(a2*sin(2*B))/2+(a4*sin(4*B))/4-(a6*sin(6*B))/6+(a8*sin(8*B))/8; end 度分秒转度函数如下： function degree = dms2degree(jiaodu) %度分秒(dd.mmss) －>度
高斯投影正反算实习报告
高斯投影正反算实习报告程序代码#include <stdio.h>#include <iostream.h>#include <math.h>#include<iomanip.h>class Angle{public:double degree,cent,second,Hudu,seconds;//构造函数Angle(double _degree,double _cent,double _second){degree=_degree;cent=_cent;second=_second;seconds=(degree*3600+cent*60+second);Hudu=(degree*3600+cent*60+second)*2*3.1415926535/(360*60*60);}Angle(){}SToD(){second=seconds-int(seconds/60)*60;degree=int((seconds-second)/3600);cent=int((seconds-second-degree*3600)/60);}~Angle(){}};void main(){//正算过程//定义变量Angle B(51,38,43.9023),L(126,2,13.1360),_B,_L;double L0,l,N,a0,a4,a6,a3,a5,x,y,rou;rou=(360*60*60)/(2*3.1415926535);if(int(L.degree)%6!=0){L0=6*(int(L.degree)/6+1)-3;}elseL0=6*(int(L.degree)/6)-3;l=L.Hudu-L0*3600*2*3.1415926535/(360*60*60);N=6399698.902-(21562.267-(108.973-0.612*cos(B.Hudu)*cos(B.Hudu))*cos(B.Hudu)*co s(B.Hudu))*cos(B.Hudu)*cos(B.Hudu);a0=32140.404-(135.3302-(0.7092-0.0040*cos(B.Hudu)*cos(B.Hudu))*cos(B.Hudu)*cos(B .Hudu))*cos(B.Hudu)*cos(B.Hudu);a4=(0.25+0.00252*cos(B.Hudu)*cos(B.Hudu))*cos(B.Hudu)*cos(B.Hudu)-0.04166;a6=(0.166*cos(B.Hudu)*cos(B.Hudu)-0.084)*cos(B.Hudu)*cos(B.Hudu);a3=(0.3333333+0.001123*cos(B.Hudu)*cos(B.Hudu))*cos(B.Hudu)*cos(B.Hudu)-0.1666 667;a5=0.0083-(0.1667-(0.1968+0.0040*cos(B.Hudu)*cos(B.Hudu))*cos(B.Hudu)*cos(B.Hud u))*cos(B.Hudu)*cos(B.Hudu);x=6367558.4969*B.Hudu-(a0-(0.5+(a4+a6*l*l)*l*l)*l*l*N)*sin(B.Hudu)*cos(B.Hudu);y=(1+(a3+a5*l*l)*l*l)*l*N*cos(B.Hudu);//结果输出cout<<"x="<<x<<endl;cout<<"y="<<y<<endl;// 反算过程double b,Bf,Nf,Z,b2,b3,b4,b5,l1;b=(x/6367558.4969);Bf=b+(50221746+(293622+(2350+22*cos(b)*cos(b))*cos(b)*cos(b))*cos(b)*cos(b))*0.000 0000001*sin(b)*cos(b);Nf=6399698.902-(21562.267-(108.973-0.612*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf))*cos(Bf)*c os(Bf);Z=y/(Nf*cos(Bf));b2=(0.5+0.003369*cos(Bf)*cos(Bf))*sin(Bf)*cos(Bf);b3=0.333333-(0.166667-0.001123*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b4=0.25+(0.16161+0.00562*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b5=0.2-(0.1667-0.0088*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);_B.seconds=Bf*rou-(1-(b4-0.12*Z*Z)*Z*Z)*Z*Z*b2*rou;l1=(1-(b3-b5*Z*Z)*Z*Z)*Z*rou;_L.seconds=L0*3600+l1;_B.SToD();_L.SToD();//结果输出cout<<"_B="<<_B.degree<<"du"<<_B.cent<<"fen"<<setprecision(8)<<_B.second<<" miao"<<endl;cout<<"_L="<<_L.degree<<"du"<<_L.cent<<"fen"<<setprecision(8)<<_L.second<<" miao"<<endl;计算结果x=5.72816e+006y=-205080_B=51du38fen43.902358miao_L=126du2fen13.135972miaoPress any key to continue程序设计流程图。
高斯投影正算与反算的理论方法与实
高斯投影正算与反算的理论方法与实现代码高斯投影是正形投影的一种，同一坐标系中的高斯投影换带计算公式是根据正形投影原理推导出的两个高斯坐标系间的显函数式。
在同一大地坐标系中(例如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)-(35 *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表示。
高斯投影坐标正反算编程报告
高斯投影坐标正反算编程报告The Standardization Office was revised on the afternoon of December 13, 2020高斯投影坐标正反算编程报告1. 编程思想进行高斯投影坐标正反算的编程需要牵涉到大量的公式，为了使程序条理更清楚，各块的数据复用性更强，这里采取了结构化的编程思想。
程序由四大块组成。
文件用于存放main()函数，是整个程序的入口。
通过结构化的编程尽力使main()函数变得简单。
和用于存放计算过程中进行角度弧度换算时所要用到的一些自定的转换函数。
和用于存放Zhengsuan 类，在Zhengsuan 类中声明了高斯投影坐标正算所要用到的所有变量，在类的构造函数中进行成员变量的初始化及正算计算。
通过get 函数获得相应的正算结果。
和用于存放Fansuan 类，类似于Zhengsuan 类，Fansuan 类中声明了高斯投影坐标反算所要用到的所有变量，在类的构造函数中进行成员变量的初始化及反算计算。
通过get 函数获得相应的反算结果。
2. 计算模型高斯投影正算公式6425644223422)5861(cos sin 720)495(cos 24cos sin 2lt t B B N lt B simB N l B B N X x ''+-''+''++-''+''⋅''+=ρηηρρ5222425532233)5814185(cos 120)1(cos 6cos lt t t B N lt B N l B N y ''-++-''+''+-''+''⋅''=ηηρηρρ高斯投影反算公式()()()()22242552233642542222328624285cos 12021cos 6cos 459061720935242f f f f f ff ff ff f f ff ff ff f f f ff f ff f f t t t B N y t B N y B N y l y t t y NM t yt tNM t y N M t B B ηηηηη+++++++-=++--+++-=3. 程序框图4. 计算结果5.附录：程序代码517"<<endl;cin>>myX>>myY;Fansuan myFansuan1(myX,myY);();}////////////////////////////////////////////////////////////////////////#include ""//////////////////////////////////////////////////////////////////////// Construction/Destruction//////////////////////////////////////////////////////////////////////Fansuan::Fansuan(){}Fansuan::Fansuan(double X,double Y){this->X=X;this->Y=Y;//初始化x,yN=(int)(Y/1000000);//取出带号L0=6*N-3;//初始化该带号的中央经线经度x=X;y=Y-1000000*N-500000;beta=x/;//初始化beta，弧度单位betaSecond=beta*rouSecond;//初始化beta，秒单位betaDegree=betaSecond/3600;//初始化beta,整度数单位Bf=beta+(+(293622+(2350+22*cos(beta)*cos(beta))*cos(beta)*cos(beta))*cos(b eta)*cos(beta))*(1e-10)*sin(beta)*cos(beta);//初始化Bf,弧度单位BfSecond=Bf*rouSecond;//初始化Bf，秒单位BfDegree=BfSecond/3600;//初始化Bf，整度数单位Nf= Z=y/(Nf*cos(Bf));b2=+*cos(Bf)*cos(Bf))*sin(Bf)*cos(Bf);b4=++*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b3= b5= Bsecond=BfSecond-(1-*Z*Z)*Z*Z)*Z*Z*b2*rouSecond;//计算大地经度B，单位为秒B=Bsecond/3600;//用整度数表示Blsecond=(1-(b3-b5*Z*Z)*Z*Z)*Z*rouSecond;//计算经度差l，单位为秒l=lsecond/3600;//用整度数表示lL=L0+l;//计算大地经度L}double Fansuan::getB(){return B;}double Fansuan::getL(){return L;}void Fansuan::printLocation(){printf("反算得大地坐标为：大地纬度B= %f°大地经度L=%f°\n",B,L);}Fansuan::~Fansuan(){}。
高斯投影正反算C#代码
高斯投影是正形投影的一种，同一坐标系中的高斯投影换带计算公式是根据正形投影原理推导出的两个高斯坐标系间的显函数式。
在同一大地坐标系中(例如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)*si n(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*lati tude1)-(35*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表。
高斯正反算实验报告
实验报告高斯正反算*名：**学号：********班级：测绘10-1班指导老师：***目录一、实验目的-----------------2二、实验内容及步骤--------2三、程序代码-----------------4四、流程图--------------------24五、运算结果-----------------26六、实验感想-----------------29一、实验目的1、了解高斯正反算的基本思想。
2.学会编写高斯正反算程序，加深理解。
二、实验内容及步骤高斯投影正算公式是由大地坐标（L，B）求解高斯平面坐标（x，y），而高斯投影反算公式则是由高斯平面坐标（x，y）求解大地坐标（L，B）。
现行的高斯投影用表都是采用克拉索夫斯基椭球参数，这次编程计算不仅采用这种椭球参数，还可以选择IAG椭球进行计算。
编程环境是在VC下，采用C++语言编写。
程序主要分为两部分，第一部分是高斯正反算函数，第二部分是主函数。
三、程序代码1.高斯投影正算// mydlg1.cpp : implementation file#include "stdafx.h"#include "高斯正反算.h"#include "mydlg1.h"#include"math.h"#define P 206264.8062471#ifdef _DEBUG#define new DEBUG_NEW#undef THIS_FILEstatic char THIS_FILE[] = __FILE__;#endifCmydlg1::Cmydlg1(CWnd* pParent /*=NULL*/) : CDialog(Cmydlg1::IDD, pParent){//{{AFX_DATA_INIT(Cmydlg1)m_num1 = 0;m_num2 = 0;m_num3 = 0.0;m_num4 = 0;m_num5 = 0;m_num6 = 0.0;m_num7 = 0;m_num8 = 0.0;m_num9 = 0.0;//}}AFX_DATA_INIT}void Cmydlg1::DoDataExchange(CDataExchange* pDX) {CDialog::DoDataExchange(pDX);//{{AFX_DATA_MAP(Cmydlg1)DDX_Text(pDX, IDC_EDIT1, m_num1);DDX_Text(pDX, IDC_EDIT2, m_num2);DDX_Text(pDX, IDC_EDIT3, m_num3);DDX_Text(pDX, IDC_EDIT4, m_num4);DDX_Text(pDX, IDC_EDIT5, m_num5);DDX_Text(pDX, IDC_EDIT6, m_num6);DDX_Text(pDX, IDC_EDIT7, m_num7);DDX_Text(pDX, IDC_EDIT8, m_num8);DDX_Text(pDX, IDC_EDIT9, m_num9);//}}AFX_DATA_MAP}BEGIN_MESSAGE_MAP(Cmydlg1, CDialog) //{{AFX_MSG_MAP(Cmydlg1)ON_BN_CLICKED(IDC_RADIO1, OnRadio1)ON_BN_CLICKED(IDC_RADIO2, OnRadio2)ON_BN_CLICKED(IDC_BUTTON2, OnButton2)ON_BN_CLICKED(IDC_BUTTON1, OnButton1)ON_BN_CLICKED(IDC_BUTTON3, OnButton3)//}}AFX_MSG_MAPEND_MESSAGE_MAP()/////////////////////////////////////////////////////////// //////////////////// Cmydlg1 message handlersvoid Cmydlg1::OnRadio1(){// TODO: Add your control notification handler code here}void Cmydlg1::OnRadio2(){// TODO: Add your control notification handler code here}void Cmydlg1::OnButton2(){// TODO: Add your control notification handler code here CEdit* pedt1 = (CEdit*)GetDlgItem(IDC_EDIT1);pedt1->SetWindowText("");CEdit* pedt2 = (CEdit*)GetDlgItem(IDC_EDIT2);pedt2->SetWindowText("");CEdit* pedt3 = (CEdit*)GetDlgItem(IDC_EDIT3);pedt3->SetWindowText("");CEdit* pedt4 = (CEdit*)GetDlgItem(IDC_EDIT4); pedt4->SetWindowText("");CEdit* pedt5= (CEdit*)GetDlgItem(IDC_EDIT5); pedt5->SetWindowText("");CEdit* pedt6 = (CEdit*)GetDlgItem(IDC_EDIT6); pedt6->SetWindowText("");CEdit* pedt7 = (CEdit*)GetDlgItem(IDC_EDIT7); pedt7->SetWindowText("");CEdit* pedt8 = (CEdit*)GetDlgItem(IDC_EDIT8); pedt8->SetWindowText("");CEdit* pedt9 = (CEdit*)GetDlgItem(IDC_EDIT9); pedt9->SetWindowText("");}void Cmydlg1::OnCancel(){// TODO: Add extra cleanup hereCDialog::OnCancel();}void Cmydlg1::OnButton1(){// TODO: Add your control notification handler code here double N,a0,a3,a4,a5,a6,B,L,l,x,y;UpdateData();B=m_num1*3600+m_num2*60+m_num3;L=m_num4*3600+m_num5*60+m_num6;B=B/P;l=(L-m_num7*3600)/P;N=6399698.902-(21562.267-(108.973-0.612*cos(B)*cos(B))*cos( B)*cos(B))*cos(B)*cos(B);a0=32140.404-(135.3302-(0.7092-0.0040*cos(B)*cos(B))*cos(B) *cos(B))*cos(B)*cos(B);a4=(0.25+0.00252*cos(B)*cos(B))*cos(B)*cos(B)-0.04166;a6=(0.166*cos(B)*cos(B)-0.084)*cos(B)*cos(B);a3=(0.3333333+0.001123*cos(B)*cos(B))*cos(B)*cos(B)-0.16666 67;a5=0.0083-(0.1667-(0.1968+0.0040*cos(B)*cos(B))*cos(B)*cos( B))*cos(B)*cos(B);x=6367558.4969*B-(a0-(0.5+(a4+a6*l*l)*l*l)*l*l*N)*sin(B)*co s(B);y=(1+(a3+a5*l*l)*l*l)*l*N*cos(B);m_num8=x;m_num9=y;UpdateData(FALSE);}void Cmydlg1::OnButton3(){// TODO: Add your control notification handler code here double N,a0,a3,a4,a5,a6,B,L,l,x,y;UpdateData();B=m_num1*3600+m_num2*60+m_num3;L=m_num4*3600+m_num5*60+m_num6;B=B/P;l=(L-m_num7*3600)/P;N=6399596.652-(21565.045-(108.996-0.603*cos(B)*cos(B))*cos( B)*cos(B))*cos(B)*cos(B);a0=32144.5189-(135.3646-(0.7034-0.0041*cos(B)*cos(B))*cos(B )*cos(B))*cos(B)*cos(B);a4=(0.25+0.00253*cos(B)*cos(B))*cos(B)*cos(B)-0.04167;a6=(0.167*cos(B)*cos(B)-0.083)*cos(B)*cos(B);a3=(0.3333333+0.001123*cos(B)*cos(B))*cos(B)*cos(B)-0.16666 67;a5=0.00878-(0.1702-(0.20382*cos(B)*cos(B)))*cos(B)*cos(B);x=6367452.1328*B-(a0-(0.5+(a4+a6*l*l)*l*l)*l*l*N)*sin(B)*co s(B);y=(1+(a3+a5*l*l)*l*l)*l*N*cos(B);m_num8=x;m_num9=y;UpdateData(FALSE);}2.高斯投影反算// mydlg2.cpp : implementation file//#include "stdafx.h"#include "高斯正反算.h"#include "mydlg2.h"#include "math.h"#define P 206264.8062471#define PI 3.1415926535898#ifdef _DEBUG#define new DEBUG_NEW#undef THIS_FILEstatic char THIS_FILE[] = __FILE__;#endif/////////////////////////////////////////////////////////// //////////////////// Cmydlg2 dialogCmydlg2::Cmydlg2(CWnd* pParent /*=NULL*/): CDialog(Cmydlg2::IDD, pParent){//{{AFX_DATA_INIT(Cmydlg2)m_num1 = 0.0;m_num2 = 0.0;m_num4 = 0;m_num5 = 0;m_num6 = 0;m_num7 = 0.0;m_num8 = 0;m_num9 = 0;m_num10 = 0.0;//}}AFX_DATA_INIT}void Cmydlg2::DoDataExchange(CDataExchange* pDX) {CDialog::DoDataExchange(pDX);//{{AFX_DATA_MAP(Cmydlg2)DDX_Text(pDX, IDC_EDIT1, m_num1);DDX_Text(pDX, IDC_EDIT2, m_num2);DDX_Text(pDX, IDC_EDIT4, m_num4);DDX_Text(pDX, IDC_EDIT5, m_num5);DDX_Text(pDX, IDC_EDIT6, m_num6);DDX_Text(pDX, IDC_EDIT7, m_num7);DDX_Text(pDX, IDC_EDIT8, m_num8);DDX_Text(pDX, IDC_EDIT9, m_num9);DDX_Text(pDX, IDC_EDIT10, m_num10);//}}AFX_DATA_MAP}BEGIN_MESSAGE_MAP(Cmydlg2, CDialog) //{{AFX_MSG_MAP(Cmydlg2)ON_BN_CLICKED(IDC_RADIO2, OnRadio2)ON_BN_CLICKED(IDC_RADIO1, OnRadio1)ON_BN_CLICKED(IDC_BUTTON1, OnButton1)ON_BN_CLICKED(IDC_BUTTON2, OnButton2)ON_BN_CLICKED(IDC_BUTTON4, OnButton4)//}}AFX_MSG_MAPEND_MESSAGE_MAP()/////////////////////////////////////////////////////////// //////////////////// Cmydlg2 message handlersvoid Cmydlg2::OnRadio2(){// TODO: Add your control notification handler code here}void Cmydlg2::OnRadio1(){// TODO: Add your control notification handler code here }void Cmydlg2::OnButton1(){// TODO: Add your control notification handler code here CEdit* pedt1 = (CEdit*)GetDlgItem(IDC_EDIT1);pedt1->SetWindowText("");CEdit* pedt2 = (CEdit*)GetDlgItem(IDC_EDIT2);pedt2->SetWindowText("");CEdit* pedt4= (CEdit*)GetDlgItem(IDC_EDIT4);pedt4->SetWindowText("");CEdit* pedt5= (CEdit*)GetDlgItem(IDC_EDIT5);pedt5->SetWindowText("");CEdit* pedt6 = (CEdit*)GetDlgItem(IDC_EDIT6); pedt6->SetWindowText("");CEdit* pedt7 = (CEdit*)GetDlgItem(IDC_EDIT7); pedt7->SetWindowText("");CEdit* pedt8 = (CEdit*)GetDlgItem(IDC_EDIT8); pedt8->SetWindowText("");CEdit* pedt9 = (CEdit*)GetDlgItem(IDC_EDIT9); pedt9->SetWindowText("");CEdit* pedt10 = (CEdit*)GetDlgItem(IDC_EDIT10); pedt10->SetWindowText("");}void Cmydlg2::OnCancel(){// TODO: Add extra cleanup hereCDialog::OnCancel();}void Cmydlg2::OnButton2(){// TODO: Add your control notification handler code here double B,L,l,Bf,Nf,b2,b3,b4,b5,b,Z;UpdateData();b=(m_num1/6367558.4969)*P;b=(b/3600)*(PI/180);Bf=b+(50221746+(293622+(2350+22*cos(b)*cos(b))*cos(b)*cos(b ))*cos(b)*cos(b))*pow(10,-10)*sin(b)*cos(b)*P/3600*(PI/180) ;Nf=6399698.902-(21562.267-(108.973-0.612*cos(Bf)*cos(Bf))*c os(Bf)*cos(Bf))*cos(Bf)*cos(Bf);Z=m_num2/(Nf*cos(Bf));b2=(0.5+0.003369*cos(Bf)*cos(Bf))*sin(Bf)*cos(Bf);b3=0.333333-(0.1666667-0.001123*cos(Bf)*cos(Bf))*cos(Bf)*co s(Bf);b4=0.25+(0.16161+0.00562*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf); b5=0.2-(0.16667-0.0088*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);Bf=Bf*(180/PI);B=Bf-(1-(b4-0.12*Z*Z)*Z*Z)*Z*Z*b2*P/3600;l=(1-(b3-b5*Z*Z)*Z*Z)*Z*P;L=m_num4+l/3600;m_num5=int(B);m_num6=int((B-m_num5)*60);m_num7=((B-m_num5)*60-m_num6)*60;m_num8=int(L);m_num9=int((L-m_num8)*60);m_num10=((L-m_num8)*60-m_num9)*60;UpdateData(FALSE);}void Cmydlg2::OnButton4(){// TODO: Add your control notification handler code here double B,L,l,Bf,Nf,b2,b3,b4,b5,b,Z;UpdateData();b=m_num1/6367558.4969*P;b=(b/3600)*(PI/180);Bf=b+(50221746+(293622+(2350+22*cos(b)*cos(b))*cos(b)*cos(b ))*cos(b)*cos(b))*pow(10,-10)*sin(b)*cos(b)*P/3600*(PI/180) ;Nf=6399698.902-(21562.267-(108.973-0.612*cos(Bf)*cos(Bf))*c os(Bf)*cos(Bf))*cos(Bf)*cos(Bf);Z=m_num2/(Nf*cos(Bf));b2=(0.5+0.00336975*cos(Bf)*cos(Bf))*sin(Bf)*cos(Bf);b3=0.333333-(0.1666667-0.001123*cos(Bf)*cos(Bf))*cos(Bf)*co s(Bf);b4=0.25+(0.161612+0.005617*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf) ;b5=0.2-(0.16667-0.00878*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);Bf=Bf*(180/PI);////////B=Bf-(1-(b4-0.147*Z*Z)*Z*Z)*Z*Z*b2*P/3600;l=(1-(b3-b5*Z*Z)*Z*Z)*Z*P;L=m_num4+l/3600;m_num5=int(B);m_num6=int((B-m_num5)*60);m_num7=((B-m_num5)*60-m_num6)*60; m_num8=int(L);m_num9=int((L-m_num8)*60);m_num10=((L-m_num8)*60-m_num9)*60; UpdateData(FALSE);}四、流程图1、高斯投影正算投影2.高斯投影反算五、运算结果1.主界面2.高斯投影正算（克氏椭球）3.高斯投影正算（IAG椭球）4.高斯投影反算（克氏椭球）5. 高斯投影反算（IAG椭球）六、实验感想我在这次实验编程中主要完成了在两种椭球下的高斯正反算。
高斯投影正反算编程一.高斯投影正反算基本公式
高斯投影正反算编程一．高斯投影正反算基本公式（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; }。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
高斯投影坐标正反算编程报告1. 编程思想进行高斯投影坐标正反算的编程需要牵涉到大量的公式，为了使程序条理更清楚，各块的数据复用性更强，这里采取了结构化的编程思想。
程序由四大块组成。
GeodesyHomework.cpp 文件用于存放main()函数，是整个程序的入口。
通过结构化的编程尽力使main()函数变得简单。
MyFunction.h 和MyFunction.cpp 用于存放计算过程中进行角度弧度换算时所要用到的一些自定的转换函数。
Zhengsuan.h 和Zhengsuan.cpp 用于存放Zhengsuan 类，在Zhengsuan 类中声明了高斯投影坐标正算所要用到的所有变量，在类的构造函数中进行成员变量的初始化及正算计算。
通过get 函数获得相应的正算结果。
Fansuan.h 和Fansuan.cpp 用于存放Fansuan 类，类似于Zhengsuan 类，Fansuan 类中声明了高斯投影坐标反算所要用到的所有变量，在类的构造函数中进行成员变量的初始化及反算计算。
通过get 函数获得相应的反算结果。
2. 计算模型高斯投影正算公式6425644223422)5861(cos sin 720)495(cos 24cos sin 2lt t B B N lt B simB N l B B N X x ''+-''+''++-''+''⋅''+=ρηηρρ5222425532233)5814185(cos 120)1(cos 6cos lt t t B N lt B N l B N y ''-++-''+''+-''+''⋅''=ηηρηρρ高斯投影反算公式()()()()22242552233642542222328624285cos 12021cos 6cos 459061720935242f f f f f ff f f ff f f ff ff ff f f f ff f ff f f t t t B N y t B N y B N y l y t t y NM t yt tNM t y N M t B B ηηηηη+++++++-=++--+++-=3. 程序框图4. 计算结果开始输入B ，L求定带号N ，中央纬度L 0，纬度差l按照实用公式计算x,y换算为国家统一坐标X ，Y输出X,Y输入国家统一坐标X,Y由Y 取定带号N ，并换算出x ，y求出中央经线L 0按照实用公式计算B ，lL=L0+l 求出大地经度L输出B ，L结束正算反算5.附录：程序代码/////主函数入口GeodesyHomework.cpp#include "MyFunction.h"#include "Zhengsuan.h"#include "Fansuan.h"#include <iostream>using namespace std;void fansuan();void zhengsuan();void main(){zhengsuan();fansuan();printf("/n over!");}void zhengsuan(){double myB,myL;cout<<"【正算】"<<endl;cout<<"请输入大地纬度B"<<endl;myB=angleToDegree();cout<<"请输入大地经度L"<<endl;myL=angleToDegree();Zhengsuan myZhengsuan1(myB,myL);printf("Radian B=%f L=%f \n",myZhengsuan1.getrB(),myZhengsuan1.getrL());myZhengsuan1.printLocation();}void fansuan(){double myX,myY;cout<<"【反算】"<<endl;cout<<"请输入国家统一坐标X Y。
例如3378627.1819 20243953.4517"<<endl; cin>>myX>>myY;Fansuan myFansuan1(myX,myY);myFansuan1.printLocation();}///自定功能函数库MyFunction.h#define PI 3.1415926#include <iostream>using namespace std;double angleToDegree(int du,int fen,float miao);double angleToDegree();//将度分秒换算为度double degreeToRadian(double degree);double degreeToRadian();//将角度换算为弧度MyFunction.cpp#include "MyFunction.h"double angleToDegree(int du,int fen,float miao){double result=0;result=miao/3600.0+fen/60.0+du;return result;}double angleToDegree(){int du,fen;float miao;double result;cout<<"请输入度分秒。
例如：30 20 00"<<endl;cin>>du>>fen>>miao;result=angleToDegree(du,fen,miao);return result;}double degreeToRadian(double degree){double result=0;result=degree/57.295779513082321;return result;}double degreeToRadian(){double result,degree;degree=angleToDegree();result=degreeToRadian(degree);return result;}///正算类Zhengsuan.h// Zhengsuan.h: interface for the Zhengsuan class.////////////////////////////////////////////////////////////////////////#if !defined(AFX_ZHENGSUAN_H__2655EA28_E810_44A3_8F14_56421A7B4466__INCL UDED_)#defineAFX_ZHENGSUAN_H__2655EA28_E810_44A3_8F14_56421A7B4466__INCLUDED_#if _MSC_VER > 1000#pragma once#endif // _MSC_VER > 1000#define rouSecond 206264.806247096355#include "MyFunction.h"#include <iostream>#include <math.h>using namespace std;class Zhengsuan{public:Zhengsuan();Zhengsuan(double fB,double fL);double getX();double getY();double getrB();double getrL();void printLocation();virtual ~Zhengsuan();private:double x;double y;//大地坐标double X;double Y;//国家统一坐标double B;double rB;int Bsecond;double L;double rL;//输入的大地纬度B，大地经度L，rB,rL为对应弧度表示值，Bsecond为换算成秒数值int n;//带号ndouble L0;//中央经线纬度L0double LDot;//纬度差L-L0int LDotSecond;//换算成秒的纬度差double l;double N;double a0;double a3;double a4;double a5;double a6;//七个计算参数};#endif// !defined(AFX_ZHENGSUAN_H__2655EA28_E810_44A3_8F14_56421A7B4466__INCLUDED _)Zhengsuan.cpp// Zhengsuan.cpp: implementation of the Zhengsuan class.////////////////////////////////////////////////////////////////////////#include "Zhengsuan.h"//////////////////////////////////////////////////////////////////////// Construction/Destruction//////////////////////////////////////////////////////////////////////Zhengsuan::Zhengsuan(){}Zhengsuan::Zhengsuan(double fB,double fL){B=fB;rB=degreeToRadian(fB);L=fL;rL=degreeToRadian(fL);Bsecond=B*3600;//初始化大地经度L，大地纬度B，Bsecond,按弧度的大地纬度rBn=(int)(L/6+1);//初始化带号nL0=6*n-3;//中央经线经度，角度单位LDot=L-L0;//经度差LDotSecond=LDot*3600;l=(LDot)*3600/rouSecond;//计算参数lN=6399698.902-(21562.267-(108.973-0.612*cos(rB)*cos(rB))*cos(rB)*cos(rB))*cos(rB) *cos(rB);//计算参数Na0=32140.404-(135.3302-(0.7092-0.004*cos(rB)*cos(rB))*cos(rB)*cos(rB))*cos(rB)*cos (rB);//计算参数a0a4=(0.25+0.00252*cos(rB)*cos(rB))*cos(rB)*cos(rB)-0.04166;//计算参数a4a6=(0.166*cos(rB)*cos(rB)-0.084)*cos(rB)*cos(rB);//计算参数a6a3=(0.3333333+0.001123*cos(rB)*cos(rB))*cos(rB)*cos(rB)-0.1666667;//计算参数a3 a5=0.0083-(0.1667-(0.1968+0.004*cos(rB)*cos(rB))*cos(rB)*cos(rB))*cos(rB)*cos(rB);//计算参数a5x=6367558.4969*Bsecond/rouSecond-(a0-(0.5+(a4+a6*l*l)*l*l)*l*l*N)*sin(rB)*cos(rB); //正算xy=(1+(a3+a5*l*l)*l*l)*l*N*cos(rB);//正算yX=x;Y=n*1000000+y+500000;//国家统一坐标}Zhengsuan::~Zhengsuan(){}double Zhengsuan::getX(){return X;}double Zhengsuan::getY(){return Y;}void Zhengsuan::printLocation(){printf("正算得国家统一坐标为：X= %8.8f Y=%8.8f \n",X,Y);}double Zhengsuan::getrB(){return rB;}double Zhengsuan::getrL(){return rL;}///反算类Fansuan.h// Fansuan.h: interface for the Fansuan class.////////////////////////////////////////////////////////////////////////#if !defined(AFX_FANSUAN_H__5B5E4668_EE81_463F_8D65_FFE2AEACE745__INCLUD ED_)#define AFX_FANSUAN_H__5B5E4668_EE81_463F_8D65_FFE2AEACE745__INCLUDED_#if _MSC_VER > 1000#pragma once#endif // _MSC_VER > 1000#define rouSecond 206264.806247096355#include <math.h>#include "MyFunction.h"#include <iostream>using namespace std;class Fansuan{public:Fansuan();Fansuan(double X,double Y);double getB();double getL();void printLocation();virtual ~Fansuan();private:double x;double y;//高斯投影坐标double X;double Y;int N;//国家统一坐标,N为带号double B,Bsecond;double L;//最后反算得到B、Ldouble L0;//中央经线经度double l,lsecond;//L=L0+l，L0=6*N-3double Bf,BfSecond,BfDegree;double beta,betaSecond,betaDegree;double Z;double Nf;double b2;double b3;double b4;double b5;//计算的8个参数};#endif// !defined(AFX_FANSUAN_H__5B5E4668_EE81_463F_8D65_FFE2AEACE745__INCLUDED_)Fansuan.cpp// Fansuan.cpp: implementation of the Fansuan class.////////////////////////////////////////////////////////////////////////#include "Fansuan.h"//////////////////////////////////////////////////////////////////////// Construction/Destruction//////////////////////////////////////////////////////////////////////Fansuan::Fansuan(){}Fansuan::Fansuan(double X,double Y){this->X=X;this->Y=Y;//初始化x,yN=(int)(Y/1000000);//取出带号L0=6*N-3;//初始化该带号的中央经线经度x=X;y=Y-1000000*N-500000;beta=x/6367558.4969;//初始化beta，弧度单位betaSecond=beta*rouSecond;//初始化beta，秒单位betaDegree=betaSecond/3600;//初始化beta,整度数单位Bf=beta+(50221746+(293622+(2350+22*cos(beta)*cos(beta))*cos(beta)*cos(beta))*co s(beta)*cos(beta))*(1e-10)*sin(beta)*cos(beta);//初始化Bf,弧度单位BfSecond=Bf*rouSecond;//初始化Bf，秒单位BfDegree=BfSecond/3600;//初始化Bf，整度数单位Nf=6399698.902-(21562.267-(108.973-0.612*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf))*cos(Bf) *cos(Bf);Z=y/(Nf*cos(Bf));b2=(0.5+0.003369*cos(Bf)*cos(Bf))*sin(Bf)*cos(Bf);b4=0.25+(0.16161+0.00562*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b3=0.333333-(0.166667-0.001123*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b5=0.2-(0.166667-0.0088*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);Bsecond=BfSecond-(1-(b4-0.12*Z*Z)*Z*Z)*Z*Z*b2*rouSecond;//计算大地经度B，单位为秒B=Bsecond/3600;//用整度数表示Blsecond=(1-(b3-b5*Z*Z)*Z*Z)*Z*rouSecond;//计算经度差l，单位为秒l=lsecond/3600;//用整度数表示lL=L0+l;//计算大地经度L}double Fansuan::getB(){return B;}double Fansuan::getL(){return L;}void Fansuan::printLocation(){printf("反算得大地坐标为：大地纬度B= %f°大地经度L=%f°\n",B,L);}Fansuan::~Fansuan(){}11。