高斯投影坐标正反算编程报告
高斯投影正反算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
高斯正反算实验报告(王震阳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.实验总结这次实验是实现高斯投影的正反算，方法很多，实现并不复杂，但是计算公式复杂，变量繁多，稍有不慎，就会造成计算错误。
高斯投影正反算-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秒左右。
高斯投影正反算实习
大地测量学编程实习报告姓名：鲁尼学号：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应该等于投影前从赤道量至该点的子午弧长。
高斯投影正反算程序
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;
- 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. 计算结果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;//计算参数a3a5=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(){}。