(完整word版)高斯投影正反算 代码
高斯投影正反算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 原理高斯投影反算的原理是根据平面坐标和中央子午线经度，计算对应的经纬度坐标。
大地测量高斯投影正反算[解说]
![大地测量高斯投影正反算[解说]](https://img.taocdn.com/s3/m/f0947b707f21af45b307e87101f69e314332faba.png)
大地测量高斯投影正反算程序代码课程：大地测量学基础姓名：林江伟学号：2008301610045班级： 0804 界面如下：输入数据计算：using System;using System.Collections.Generic; using ponentModel;using System.Data;using System.Drawing;using System.Text;using System.Windows.Forms;namespace大地{public partial class Form1 : Form {double B, L;double x, y;double X, Y;int N;double L0, l1;double p=206264.80625;public Form1(){InitializeComponent();}private void Form1_Load(object sender, EventArgs e){//自动化控件显示初始值radioButton2.Checked=true;radioButton3.Checked = true;textBox1.Focus();textBox1.Text = textBox2.Text = textBox3.Text = textBox4.Text = textBox5.Text = textBox6.Text = "0";this.richTextBox2.Text = "说明：输入的坐标需为按6°带投影且采用克氏椭球参数所得的国家统一坐标";}private void button1_Click(object sender, EventArgs e){//获取输入数据double bb1=Convert.ToDouble(this.textBox1.Text);double bb2=Convert.ToDouble(this.textBox2.Text);double bb3=Convert.ToDouble(this.textBox3.Text);double ll1=Convert.ToDouble(this.textBox4.Text);double ll2=Convert.ToDouble(this.textBox5.Text);double ll3=Convert.ToDouble(this.textBox6.Text);//检查输入格式的正确性if (bb1 >= 0 && bb1 <90 && bb2 >= 0 && bb2 < 60 && bb3 >= 0 && bb3 < 60){B = bb1 * 3600 + bb2 * 60 + bb3;}else{MessageBox.Show("纬度输入格式不正确！", "警告");return;}if (ll1 >= 0 && ll1<360 && ll2 >= 0 && ll2 < 60 && ll3 >= 0 && ll3 < 60){L = ll1 * 3600 + ll2 * 60 + ll3;}else{MessageBox.Show("经度输入格式不正确！", "警告");return;}double b = B/p;//获取用户选项值int tyd=radioButton1.Checked?3:6;int tq=radioButton3.Checked?1:2;//按带求带号以及中央子午线经度if (tyd == 6){N =(int)( L / (6*3600.0)) + 1;L0 = (6 * N - 3)*3600.0;l1 = (L-L0)/p;}if(tyd==3){N =(int)(L / (3*3600.0)+0.5);L0 = 3 * N * 3600.0;l1 =( L-L0)/p;}//采用克氏椭球参数的计算公式if (tq == 1){double c = Math.Pow( Math.Cos(b),2);double c1= Math.Sin(b) * Math.Cos(b);double c2=Math.Cos(b);double l2 = Math.Pow(l1,2);double n = 6399698.902 - (21562.267 - (108.973 - 0.612 * c) * c) * c;double a0 = 32140.404 - (135.3302 - (0.7092 - 0.0040 * c) * c) * c;double a4 = (0.25 + 0.00252 * c) * c - 0.04166;double a6 = (0.166 * c - 0.084) * c;double a3 = (0.3333333 + 0.001123 * c) * c - 0.1666667;double a5 = 0.0083 - (0.1667 - (0.1968 + 0.0040 * c) * c) * c;x = 6367558.4969 * b - (a0 - (0.5 + (a4 + a6 * l2) * l2) * l2 * n) *c1; y = (1 + (a3 + a5 * l2) * l2) * l1 * n * c2;X = x;double y1 = y + 500000.0;for (int i = 1; y1 / i > 1; i = i * 10){Y = N * i*10 + y1;}string tuoqiu = "采用克氏椭球参数,";this.richTextBox1.Text = "按经差" + tyd + "°进行投影分带，" + tuoqiu + "其计算结果为：" + "\r\n" + "x=" + x + "\r\n" + "y=" + y + "\r\n" +"国家统一坐标为：" + "\r\n" + "X=" + X + "\r\n" + "Y=" + Y;}//采用1975国际椭球参数的计算公式if (tq == 2){double c = Math.Cos(b) * Math.Cos(b);double c1 = Math.Sin(b) * Math.Cos(b);double c2 = Math.Cos(b);double l2 = Math.Pow(l1, 2);double n = 6399596.652 - (21565.045 - (108.996 - 0.603 * c) * c) * c;double a0 = 32144.5189 - (135.3646 - (0.7034 - 0.0041 * c) * c) * c;double a4 = (0.25 + 0.00253 * c) * c - 0.04167;double a6 = (0.167 * c - 0.083) * c;double a3 = (0.3333333 + 0.001123 * c) * c - 0.1666667;double a5 = 0.00878 - (0.1702 - 0.20382 * c) * c;x = 6367452.1328 * b-(a0 - (0.5 + (a4 + a6 * l2) * l2) * l2 * n) * c1;y = (1 + (a3 + a5 * l2) * l2) * l1 * n * c2;X = x;double y1 = y + 500000.0;for (int i = 1; y1 / i > 1; i = i * 10)//Y坐标加代号{Y = N * i * 10 + y1;}string tuoqiu = "采用1975国际椭球参数,";this.richTextBox1.Text = "按经差" + tyd + "°进行投影分带，" + tuoqiu + "其计算结果为：" + "\r\n" + "x=" + x + "\r\n" + "y=" + y+"\r\n"+"国家统一坐标为：" + "\r\n" + "X=" + X + "\r\n"+"Y="+Y;}}private void button2_Click(object sender, EventArgs e){X = Convert.ToDouble(this.textBox7.Text);//获取输入数据Y = Convert.ToDouble(this.textBox8.Text);for (int i = 1; Y/i >=10; i = i * 10)//对Y坐标处理并求出中央子午线经度{y = Y - (int)(Y / i) * i-500000;L0 = 6 * (int)(Y / i) - 3;}//按6°带克氏椭球反算double bt = x / 6367558.4969*p;double BT = x / 6367558.4969;double c3=Math.Cos(BT)*Math.Cos(BT);double c4=Math.Sin(BT)*Math.Cos(BT);double Bf=(bt+(50221746+(293622+(2350+22*c3)*c3)*c3)*c4*Math.Pow(10,-10)*p)/p;double c5=Math.Pow(Math.Cos(Bf),2);double c6=Math.Sin(Bf)*Math.Cos(Bf);double Nf=6399698.902-(21562.267-(108.973-0.612*c5)*c5)*c5;double Z=y/(Nf*Math.Cos(Bf));double b2 = (0.5 + 0.003369 * c5) * c6;double b3 = 0.333333 - (0.166667 - 0.001123 * c5) * c5;double b4 = 0.25 + (0.16161 + 0.00562 * c5) * c5;double b5=0.2-(0.1667-0.0088*c5)*c5;double z2=Math.Pow(Z,2);B = (Bf*p - (1 - (b4 - 0.12 *z2) * z2) * z2 * b2 * p)/3600.0;L = L0+((1 - (b3 - b5 * z2) * z2) * Z * p)/3600.0;this.richTextBox2.Text = "按6°带克氏椭球反算后，结果为："+ "\r\n"+"B="+ B + "\r\n" + "L=" + L;}}。
高斯投影正算
高斯投影正、反算代码//高斯投影正、反算//////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;}。
高斯投影正反算
高斯投影正、反算及换带程序执行条件※数组投影选择T、换算点个数“Z=0 F≠0”、=0正算0、≠0反算※坐标系选择“54 ≠54”、=54换算为1954年北京坐标系输入54、≠54换算为1988年西安坐标系M、中央子午线经度（°′″）输入※大地坐标I、序列号B、L：大地纬度和经度（地理坐标）（°′″）※高斯平面坐标轴子午线I、序列号X、Y：高斯平面坐标（m） Z、轴子午线（°）输出※大地坐标子午收敛角N、序列号B、L：大地纬度和经度（地理坐标）（°′″） R、子午收敛角（°′″）※高斯平面坐标子午收敛角N、序列号X、Y：高斯平面坐标（m） R、子午收敛角（°′″）注：1、程序执行前必须进行数组定位。
如：Defm 10 T×2=5×2=102、Y坐标值要去掉带号及避免出现负值的500公里；4、本程序运算时，各已知数据、观测变量不会随之变化，可非常方便地进行各数据的核对;5、本程序在进行换带计算时采用的是间接换带计算法。
Prog GSXYDefm 10:TA“Z=0 F≠0”G“54 ≠54”Z:Fixm:I=0:「b」0:I=I+1◢J=2I-1:M=Z[J:L=Z[J+1:A=0=>Prog“3”:B=M:M=L+Z:Prog“3”:L=M:{BL}:M=B:Prog“2”: B=M:M=L:Prog“2”:L=M-Z:≠>X=M:Y=L:{XY}:B=X:L=Y⊿Z[J]=B:Z[J+1]=L:I<T=>Goto 0⊿G=54=>C=6399698.90178271:E=.006738525414684:≠>C=6399596.65198801:E=.006 739501819473⊿I=0:「b」0:I“N”=I+1◢J=2I-1:B=Z[J:L=Z[J+1:A≠0=>X=B:Y=L:Goto 2⊿S=sin B:G=54=>F=111134.8611B-(32 005.7799S+133.9238S∧3+.6973S∧5+.0039S∧7)cos B:≠>F=111133.0047B-(32009.857 S+133.9602S∧3+.6976S∧5+.0039S∧7)cos B⊿U=√Ecos B:V=√(1+U2:N=C÷V:W=tan B: M=cos B(Lπ÷180:X=F+NW(.5M2+1┛24(5-W2+9U2+4U∧4)M∧4+1┛720(61-58W2+W∧4)M∧6◢Y=N(M+1┛6(1-W 2+U 2)M ∧3+1┛120(5-18W 2+W ∧4+14U 2-58U 2W 2)M ∧5◢M=W ┛π(180M+60(1+3U 2+2U ∧4)M ∧3+12(2-W 2)M ∧5:Goto 3:「b 」2:W=E ﹣6X-3:G=54=>F=27.11115372595+9.024********W-.00579740442W 2-4.3532572E ﹣4W ∧3+4.857285E ﹣5W ∧4+2.15727E ﹣6W ∧5-1.9399E ﹣7W ∧6:≠>F=27.11162289465+9.024********W-.00579850656W2-4.3540029E ﹣4W ∧3+4.858357E ﹣5W ∧4+2.15769E ﹣6W ∧5-1.9404E ﹣7W ∧6⊿U=√Ecos F:V=√(1+U 2:Q=YV ÷C:W=tan F:M=F-(1+U 2)W ┛π(90Q 2-7.5(5+3W 2+U 2-9U 2W 2)Q ∧4+.25(61+90W 2+45W ∧4)Q ∧6:Prog “3”:B=M ◢M=Z+1┛(πcos F)(180Q-30(1+2W 2+U 2)Q ∧3+1.5(5+28W 2+24W ∧4)Q ∧5:Prog “3”:L=M ◢M=W ┛π(180Q-60(1+W 2-U 2)Q ∧3+12(2+5W 2+3W ∧4)Q ∧5:「b 」3:Prog “3”:R=M ◢ I<T=>Goto 1⊿“END ”概要说明：我国的经度范围西边自73°起，东边至135°，可分成6°带共11带或3°共22带。
高斯投影正反算-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秒左右。
高斯投影正反算程序
printf("高斯投影反算\nunder CGCS2000,六度带\n测试数据：x=2562083.2708,y=19512837.2851\n");
double a=6378137;
double f=1/298.257222101;
double e2=-f*f+2*f;
double x=2562083.2708;
double y=19512837.2851;
int Y=int(y);
int count=0;
while (Y>=100)
{
Y=Y/10;
count++;
}
y=y-Y*pow(10.0,count)-500000;
{
int model;
printf("高斯投影正算请输入1，反算请输入2\n");
scanf("%d",&model);
if(model==1)
{
printf("高斯投影正算\nunder WGS-84,六度带\n测试数据：B=22.155898294，L=111.285215387\n");
double W=sqrt(1-e2*sin(B)*sin(B));
double N=a/W;
double ep2=e2/(1-e2);
double ita2=ep2*cos(B)*cos(B);
double t=tan(B);
double temp1=5-t*t+9*ita2+4*ita2*ita2;
fBf=-a2*sin(2*Bf)/2+a4*sin(4*Bf)/4-a6*sin(6*Bf)/6+a8*sin(8*Bf)/8;
高斯投影正反算 c#代码
高斯投影正反算程序设计一．程序设计流程本程序的设计思路如下：（1），程序采用VS08版本作为开发平台，并采用C#语言作为开发语言，设计为WindowsForm窗体程序形式。
（2），程序主要的算法来自于教材。
但是本程序为了更加实用，添加了更多的解算基准，包括：WGS-84，国际椭球（3{{//{p ublicdoubledLongitude;publicdoubledLatitude;publicdoubledHeight;}//笛卡尔坐标//CartesianCoordinatepublicstructCRDCARTESIAN{publicdoublex;publicdoubley;publicdoublez;}publicForm1(){InitializeComponent();}privatevoidbutton1_Click(objectsender,EventArgse) {try{tt=}{}{ee=}{ee=}if(pareTo("1975国际椭球")==0){a=6378140.00;ee=}if(pareTo("2000国家大地坐标系")==0){a=6378137.0;}constdoublepai=3.1415926;doubleb=Math.Sqrt(a*a*(1-ee*ee));doublec=a*a/b;doubleepp=Math.Sqrt((a*a-b*b)/b/b);CRDGEODETICpcrdGeo;CRDCARTESIANpcrdCar;doublemidlong;//{}ai;//{}pai;//intnum;//带号midlong=0;//默认值，需要制定分带try{tt=}catch{MessageBox.Show("Choose3/6error!");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);}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;doubleB=pcrdGeo.dLatitude;doublesb=Math.Sin(B);doublecb=Math.Cos(B);doubles2b=sb*cb*2;doubles4b=s2b*(1-2*sb*sb)*2;doubles6b=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(l p,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 );}{try{tt=}{}{ee=}if(pareTo("WGS-84")==0){a=6378137.00;ee=}if(pareTo("1975国际椭球")==0){a=6378140.00;ee=}if(pareTo("2000国家大地坐标系")==0) {a=6378137.0;ee}constdoublepai=doubleb=Math.Sqrt(a*a*(1-ee*ee));//求Xtry{tt=}{}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;doublem0,m2,m4,m6,m8;doublea0,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;{}tf=Math.Tan(Bf);Vf=Math.Sqrt(1+epp*epp*Math.Cos(Bf)*Math.Cos(Bf));Nf=c/Vf;doubleynf=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:"+Conver t.ToString(pcrdGeo.dLongitude);}privatevoidlabel13_Click(objectsender,EventArgse){}}}米之间，。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
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
#define PI 3.141592653589793voiBiblioteka GaosZ_fun(){
printf("高斯投影的正算\n");
y=(1+(a3+a5*l*l)*l*l)*l*N*cos(B);
printf("x=%f\ny=%f\n",x,y);
}
void GaosF_fun()
{
printf("高斯投影的反算\n");
double B,Bf,Nf,b,b2,b3,b4,b5,Z,x,y,L0,l;
// printf("输入x y的值\n");
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.167-0.0088*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);
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;
double a0,a4,a6,a3,a5,cB2;
e2=0.006738525414683;
c=6399698.901782711;
// printf("输入x y的值\n");
// cout<<"x=";
// cin>>x;
// cout<<"y=";
// cin>>y;
// scanf("x=%f",&x); ///////////// ?????
B=Bf-(1-(b4-0.12*Z*Z)*Z*Z)*Z*Z*b2;
l=(1-(b3-b5*Z*Z)*Z*Z)*Z;
printf("B=%f\nL=%f\n",B*180/PI,(L0+l*180/PI));
}
void GaosLT_fun()
{
double B,Bf,Nf,b,b2,b3,b4,b5,Z,x,y,L01,L02,n2,l,V,L,N,c,t,e2;
// scanf("y=%f",&y);
// x=3380330.875;
// y=320089.976;
x=1944359.607;
y=240455.4563;
L01=117*3600/P;
L02=120*3600/P;
b=x/6367558.4969;
Bf=b+(50221746+ ( 293622+ (2350+22*cos(b)*cos(b))*cos(b)*cos(b) )
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.167-0.0088*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);
double l,L,B,n2,x,y,N,t,V,c,e2;
double i,j,k,n,h,a0,a4,a6,a3,a5,cB2;
int m;
e2=0.006738525414683;
c=6399698.901782711;
B=17.33557339*3600/P;
L=119.15521159*3600/P;
y=(1+(a3+a5*l*l)*l*l)*l*N*cos(B);
printf("输出高斯邻带坐标换算\n");
printf("x=%f\ny=%f\n",x,y);
}
void main()
{
int i;
printf("输入1 2 3选择投影的计算方法\n");
scanf("%d",&i);
switch (i)
// y=N*cos(B)*l+N*pow(cos(B),3)*(1-t*t+n2)*pow(l,5)/6+N*pow(cos(B),5)*(5-18*t*t
+pow(t,4)+14*n2-58*n2*t*t)*pow(l,5)/120;
x=6367558.4969*B-(a0-(0.5+(a4+a6*l*l)*l*l)*l*l*N)*sin(B)*cos(B);
// y=N*cos(B)*l+N*pow(cos(B),3)*(1-t*t+n2)*pow(l,5)/6+N*pow(cos(B),5)*(5-18*t*t
+pow(t,4)+14*n2-58*n2*t*t)*pow(l,5)/120;
x=6367558.4969*B-(a0-(0.5+(a4+a6*l*l)*l*l)*l*l*N)*sin(B)*cos(B);
{
case 1: GaosZ_fun(); break;
case 2: GaosF_fun(); break;
case 3: GaosLT_fun(); break;
default : break;
}
}
*cos(b)*cos(b)) *sin(b)*cos(b)*1e-10;
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);
*cos(b)*cos(b)) *sin(b)*cos(b)*1e-10;
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);
a5=0.0083-(0.1667-(0.1968+0.0040*cB2)*cB2)*cB2;
// x=X+N*sin(B)*cos(B)*l*l/2+N*sin(B)*pow(cos(B),3)*(5-t*t+9*n2+4*n2*n2)*pow(l,
4)/24+N*sin(B)*pow(cos(B),5)*(61-58*t*t+pow(t,4))*pow(l,6)/720;
a5=0.0083-(0.1667-(0.1968+0.0040*cB2)*cB2)*cB2;
// x=X+N*sin(B)*cos(B)*l*l/2+N*sin(B)*pow(cos(B),3)*(5-t*t+9*n2+4*n2*n2)*pow(l,
4)/24+N*sin(B)*pow(cos(B),5)*(61-58*t*t+pow(t,4))*pow(l,6)/720;
// scanf("x=%f",&x);
// scanf("y=%f",&y);