高斯投影正反算编程
高斯投影正反算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
高斯投影坐标正反算编程报告
高斯投影坐标正反算编程报告1.编程思想高斯投影坐标正反算的编程需要涉及大量公式。
为了使程序更清晰，各模块的数据重用性更强，本文采用结构化编程思想。
程序由四大块组成。
大地测量家庭作业。
Cpp文件用于存储main（）函数，是整个程序的入口。
尝试通过结构化编程简化main（）函数。
myfunction.h和myfunction.cpp用于存放计算过程中进行角度弧度换算时所要用到的一些自定的转换函数。
正蒜。
H和zhengsuan CPP用于存储zhengsuan类。
在正算类中，声明了用于高斯投影坐标正演计算的所有变量。
成员变量的初始化和正向计算在类的构造函数中执行。
通过get函数获得相应的阳性结果。
fansuan.h和fansuan.cpp用于存放fansuan类，类似于zhengsuan类，fansuan类中声明了高斯投影坐标反算所要用到的所有变量，在类的构造函数中进行成员变量的初始化及反算计算。
通过get函数获得相应的反算结果。
2.计算模型高斯投影正算公式十、十、nn232244？？？？sinbcosb？Lsimbcosb（5？t？9？4？）l2？？？224??? 4n5246？？sinbcosb（61？58t？t）l720？？？6y？？nn3223？？cosb？Lcosb（1？t？？）L6.3n5242225？？cosb（5？18吨？14？58吨）l120？？？五tf2mfnf5f高斯投影反算公式B男朋友？？tfy？2tf24mfn3f？5.3t2f224？？2f？9? ftfy？720mfn46y61？90t2？45泰夫？？yy32l？？1.2t2f？？f3nfcosbf6nfcosbf？？？y524222?5?28t?24t?6??8?fffftf5120nfcosbf?3.程序框图一开始输入b，l求定带号n，中央纬度l0，纬度差l正算按照实用公式计算x,y换算为国家统一坐标x，y输出x,y输入国家统一坐标x,y由y取定带号n，并换算出x，y求出中央经线l0反算按照实用公式计算b，ll=l0+l求出大地经度l输出b，l结束4.计算结果25.附录：程序代码/////主函数入口大地测量家庭作业。
高斯投影正算
高斯投影正、反算代码//高斯投影正、反算//////6度带宽 54年北京坐标系//高斯投影由经纬度(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/30 72)*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;}//高斯投影由大地坐标(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;}。
高斯投影正算与反算的理论方法与实
高斯投影正算与反算的理论方法与实现代码高斯投影是正形投影的一种，同一坐标系中的高斯投影换带计算公式是根据正形投影原理推导出的两个高斯坐标系间的显函数式。
在同一大地坐标系中(例如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表示。
高斯投影正反算c代码
高斯投影正反算程序设计一．程序设计流程本程序的设计思路如下：（1），程序采用VS08版本作为开发平台，并采用C#语言作为开发语言，设计为WindowsForm 窗体程序形式。
（2），程序主要的算法来自于教材。
但是本程序为了更加实用，添加了更多的解算基准，包括：WGS-84，国际椭球1975，克氏椭球，和2000国家大地坐标系。
（3），程序为了更方便的读取数据和输出数据，故需要自己定义了固定的数据输入格式和数据输出格式或形式，请老师注意查看。
二．代码using System;using System.Collections.Generic;using ponentModel;using System.Data;using System.Drawing;using System.Text;using System.Windows.Forms;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{tt = boBox1.Items[boBox1.SelectedIndex].ToString(); }catch{MessageBox.Show("Gauss Inverse: Choose datum error!");return;}if (pareTo("克氏椭球")==0){a = 6378245.00;ee = Math.Sqrt(0.006693421622);}if (pareTo("WGS-84") == 0){a = 6378137.00;ee = Math.Sqrt(0.00669437999013);}if (pareTo("1975国际椭球") == 0){a = 6378140.00;ee = Math.Sqrt(0.006694384999588);}if (pareTo("2000国家大地坐标系") == 0){a = 6378137.0;ee =Math.Sqrt(0.0066943802290);}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;temp = textBox1.Text.Split(' ');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;//求经度temp = textBox2.Text.Split(' ');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.dLongitude * 180 / pai);//求中央经度int num; //带号 midlong = 0; //默认值，需要制定分带try{tt = boBox3.Items[boBox3.SelectedIndex].ToString();}catch{MessageBox.Show("Choose 3/6 error!");return;}if (pareTo("3度带") == 0){num = Convert.ToInt32(deglon / 6 + 1);midlong = (6 * num - 3) / 180.0 * pai;}if (pareTo("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 * epp * Math.Cos(pcrdGeo.dLatitude) * Math.Cos(pcrdGeo.dLatitude));double M = c / Math.Pow(1 + epp * epp * Math.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) + Ncosb * Math.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{tt = boBox2.Items[boBox2.SelectedIndex].ToString(); }catch{MessageBox.Show("Gauss Inverse: Choose datum error!");return;}if (pareTo("克氏椭球") == 0){a = 6378245.00;ee = Math.Sqrt(0.006693421622);}if (pareTo("WGS-84") == 0){a = 6378137.00;ee = Math.Sqrt(0.00669437999013);}if (pareTo("1975国际椭球") == 0){a = 6378140.00;ee = Math.Sqrt(0.006694384999588);}if (pareTo("2000国家大地坐标系") == 0){a = 6378137.0;ee =Math.Sqrt(0.0066943802290);}const double pai = 3.1415926535898;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{tt = boBox4.Items[boBox4.SelectedIndex].ToString(); }catch{MessageBox.Show("Choose 3/6 error!");return;}if (pareTo("3度带") == 0){midlong = num * 3 * pai / 180;}if (pareTo("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){}}}三．程序运行结果分析通过选取书上的具体实例进行测试，本程序的精度大体满足要求，一般正算的精度在0.01米和0.001米之间，反算的精度在0.0001秒左右。
测绘程序设计(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.实验总结这次实验是实现高斯投影的正反算，方法很多，实现并不复杂，但是计算公式复杂，变量繁多，稍有不慎，就会造成计算错误。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
高斯投影正反算编程一．高斯投影正反算基本公式（1）高斯正算基本公式（2）高斯反算基本公式以上主要通过大地测量学基础课程得到，这不进行详细的推导，只是列出基本公式指导编程的进行。
二．编程的基本方法和流程图（1）编程的基本方法高斯投影正反算基本上运用了所有的编程基本语句，本文中是利用C++语言进行基本的设计。
高斯正算中对椭球参数和带宽的选择主要运用了选择语句。
而高斯反算中除了选择语句的应用，在利用迭代算法求底点纬度还应用了循环语句。
编程中还应特别注意相关的度分秒和弧度之间的相互转换，这是极其重要的。
（2）相关流程图1）正算输入大地坐标B，L 和经差L0选择带宽3/6度带计算带号计算弧长计算平面坐标x,y打印x,y计算带号计算弧长计算平面坐标x,y打印x,y开始选择椭球参数3度带6度带2）反算开始输入自然值坐标x,y和经差L0选择椭球参数利用迭代算法求解底点纬度利用公式计算B和L打印B和L三．编程的相关代码（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))*1 00.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))*10 0.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-1 8*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;}。