高斯投影坐标正反算编程报告1. 编程思想进行高斯投影坐标正反算的编程需要牵涉到大量的公式,为了使程序条理更清楚,各块的数据复用性更强,这里采取了结构化的编程思想。
程序由四大块组成。
GeodesyHomework.cpp 文件用于存放main()函数,是整个程序的入口。
通过结构化的编程尽力使main()函数变得简单。
MyFunction.h 和MyFunction.cpp 用于存放计算过程中进行角度弧度换算时所要用到的一些自定的转换函数。
Zhengsuan.h 和Zhengsuan.cpp 用于存放Zhengsuan 类,在Zhengsuan 类中声明了高斯投影坐标正算所要用到的所有变量,在类的构造函数中进行成员变量的初始化及正算计算。
通过get 函数获得相应的正算结果。
Fansuan.h 和Fansuan.cpp 用于存放Fansuan 类,类似于Zhengsuan 类,Fansuan 类中声明了高斯投影坐标反算所要用到的所有变量,在类的构造函数中进行成员变量的初始化及反算计算。
通过get 函数获得相应的反算结果。
2. 计算模型高斯投影正算公式6425644223422)5861(cos sin 720)495(cos 24cos sin 2lt t B B N lt B simB N l B B N X x ''+-''+''++-''+''⋅''+=ρηηρρ5222425532233)5814185(cos 120)1(cos 6cos lt t t B N lt B N l B N y ''-++-''+''+-''+''⋅''=ηηρηρρ高斯投影反算公式()()()()22242552233642542222328624285cos 12021cos 6cos 459061720935242f f f f f ff f f ff f f ff ff ff f f f ff f ff f f t t t B N y t B N y B N y l y t t y NM t yt tNM t y N M t B B ηηηηη+++++++-=++--+++-=3. 程序框图4. 计算结果开始输入B ,L求定带号N ,中央纬度L 0,纬度差l按照实用公式计算x,y换算为国家统一坐标X ,Y输出X,Y输入国家统一坐标X,Y由Y 取定带号N ,并换算出x ,y求出中央经线L 0按照实用公式计算B ,lL=L0+l 求出大地经度L输出B ,L结束正算反算5.附录:程序代码/////主函数入口GeodesyHomework.cpp#include "MyFunction.h"#include "Zhengsuan.h"#include "Fansuan.h"#include <iostream>using namespace std;void fansuan();void zhengsuan();void main(){zhengsuan();fansuan();printf("/n over!");}void zhengsuan(){double myB,myL;cout<<"【正算】"<<endl;cout<<"请输入大地纬度B"<<endl;myB=angleToDegree();cout<<"请输入大地经度L"<<endl;myL=angleToDegree();Zhengsuan myZhengsuan1(myB,myL);printf("Radian B=%f L=%f \n",myZhengsuan1.getrB(),myZhengsuan1.getrL());myZhengsuan1.printLocation();}void fansuan(){double myX,myY;cout<<"【反算】"<<endl;cout<<"请输入国家统一坐标X Y。
例如3378627.1819 20243953.4517"<<endl; cin>>myX>>myY;Fansuan myFansuan1(myX,myY);myFansuan1.printLocation();}///自定功能函数库MyFunction.h#define PI 3.1415926#include <iostream>using namespace std;double angleToDegree(int du,int fen,float miao);double angleToDegree();//将度分秒换算为度double degreeToRadian(double degree);double degreeToRadian();//将角度换算为弧度MyFunction.cpp#include "MyFunction.h"double angleToDegree(int du,int fen,float miao){double result=0;result=miao/3600.0+fen/60.0+du;return result;}double angleToDegree(){int du,fen;float miao;double result;cout<<"请输入度分秒。
例如:30 20 00"<<endl;cin>>du>>fen>>miao;result=angleToDegree(du,fen,miao);return result;}double degreeToRadian(double degree){double result=0;result=degree/57.295779513082321;return result;}double degreeToRadian(){double result,degree;degree=angleToDegree();result=degreeToRadian(degree);return result;}///正算类Zhengsuan.h// Zhengsuan.h: interface for the Zhengsuan class.////////////////////////////////////////////////////////////////////////#if !defined(AFX_ZHENGSUAN_H__2655EA28_E810_44A3_8F14_56421A7B4466__INCL UDED_)#defineAFX_ZHENGSUAN_H__2655EA28_E810_44A3_8F14_56421A7B4466__INCLUDED_#if _MSC_VER > 1000#pragma once#endif // _MSC_VER > 1000#define rouSecond 206264.806247096355#include "MyFunction.h"#include <iostream>#include <math.h>using namespace std;class Zhengsuan{public:Zhengsuan();Zhengsuan(double fB,double fL);double getX();double getY();double getrB();double getrL();void printLocation();virtual ~Zhengsuan();private:double x;double y;//大地坐标double X;double Y;//国家统一坐标double B;double rB;int Bsecond;double L;double rL;//输入的大地纬度B,大地经度L,rB,rL为对应弧度表示值,Bsecond为换算成秒数值int n;//带号ndouble L0;//中央经线纬度L0double LDot;//纬度差L-L0int LDotSecond;//换算成秒的纬度差double l;double N;double a0;double a3;double a4;double a5;double a6;//七个计算参数};#endif// !defined(AFX_ZHENGSUAN_H__2655EA28_E810_44A3_8F14_56421A7B4466__INCLUDED _)Zhengsuan.cpp// Zhengsuan.cpp: implementation of the Zhengsuan class.////////////////////////////////////////////////////////////////////////#include "Zhengsuan.h"//////////////////////////////////////////////////////////////////////// Construction/Destruction//////////////////////////////////////////////////////////////////////Zhengsuan::Zhengsuan(){}Zhengsuan::Zhengsuan(double fB,double fL){B=fB;rB=degreeToRadian(fB);L=fL;rL=degreeToRadian(fL);Bsecond=B*3600;//初始化大地经度L,大地纬度B,Bsecond,按弧度的大地纬度rBn=(int)(L/6+1);//初始化带号nL0=6*n-3;//中央经线经度,角度单位LDot=L-L0;//经度差LDotSecond=LDot*3600;l=(LDot)*3600/rouSecond;//计算参数lN=6399698.902-(21562.267-(108.973-0.612*cos(rB)*cos(rB))*cos(rB)*cos(rB))*cos(rB) *cos(rB);//计算参数Na0=32140.404-(135.3302-(0.7092-0.004*cos(rB)*cos(rB))*cos(rB)*cos(rB))*cos(rB)*cos (rB);//计算参数a0a4=(0.25+0.00252*cos(rB)*cos(rB))*cos(rB)*cos(rB)-0.04166;//计算参数a4a6=(0.166*cos(rB)*cos(rB)-0.084)*cos(rB)*cos(rB);//计算参数a6a3=(0.3333333+0.001123*cos(rB)*cos(rB))*cos(rB)*cos(rB)-0.1666667;//计算参数a3 a5=0.0083-(0.1667-(0.1968+0.004*cos(rB)*cos(rB))*cos(rB)*cos(rB))*cos(rB)*cos(rB);//计算参数a5x=6367558.4969*Bsecond/rouSecond-(a0-(0.5+(a4+a6*l*l)*l*l)*l*l*N)*sin(rB)*cos(rB); //正算xy=(1+(a3+a5*l*l)*l*l)*l*N*cos(rB);//正算yX=x;Y=n*1000000+y+500000;//国家统一坐标}Zhengsuan::~Zhengsuan(){}double Zhengsuan::getX(){return X;}double Zhengsuan::getY(){return Y;}void Zhengsuan::printLocation(){printf("正算得国家统一坐标为:X= %8.8f Y=%8.8f \n",X,Y);}double Zhengsuan::getrB(){return rB;}double Zhengsuan::getrL(){return rL;}///反算类Fansuan.h// Fansuan.h: interface for the Fansuan class.////////////////////////////////////////////////////////////////////////#if !defined(AFX_FANSUAN_H__5B5E4668_EE81_463F_8D65_FFE2AEACE745__INCLUD ED_)#define AFX_FANSUAN_H__5B5E4668_EE81_463F_8D65_FFE2AEACE745__INCLUDED_#if _MSC_VER > 1000#pragma once#endif // _MSC_VER > 1000#define rouSecond 206264.806247096355#include <math.h>#include "MyFunction.h"#include <iostream>using namespace std;class Fansuan{public:Fansuan();Fansuan(double X,double Y);double getB();double getL();void printLocation();virtual ~Fansuan();private:double x;double y;//高斯投影坐标double X;double Y;int N;//国家统一坐标,N为带号double B,Bsecond;double L;//最后反算得到B、Ldouble L0;//中央经线经度double l,lsecond;//L=L0+l,L0=6*N-3double Bf,BfSecond,BfDegree;double beta,betaSecond,betaDegree;double Z;double Nf;double b2;double b3;double b4;double b5;//计算的8个参数};#endif// !defined(AFX_FANSUAN_H__5B5E4668_EE81_463F_8D65_FFE2AEACE745__INCLUDED_)Fansuan.cpp// Fansuan.cpp: implementation of the Fansuan class.////////////////////////////////////////////////////////////////////////#include "Fansuan.h"//////////////////////////////////////////////////////////////////////// Construction/Destruction//////////////////////////////////////////////////////////////////////Fansuan::Fansuan(){}Fansuan::Fansuan(double X,double Y){this->X=X;this->Y=Y;//初始化x,yN=(int)(Y/1000000);//取出带号L0=6*N-3;//初始化该带号的中央经线经度x=X;y=Y-1000000*N-500000;beta=x/6367558.4969;//初始化beta,弧度单位betaSecond=beta*rouSecond;//初始化beta,秒单位betaDegree=betaSecond/3600;//初始化beta,整度数单位Bf=beta+(50221746+(293622+(2350+22*cos(beta)*cos(beta))*cos(beta)*cos(beta))*co s(beta)*cos(beta))*(1e-10)*sin(beta)*cos(beta);//初始化Bf,弧度单位BfSecond=Bf*rouSecond;//初始化Bf,秒单位BfDegree=BfSecond/3600;//初始化Bf,整度数单位Nf=6399698.902-(21562.267-(108.973-0.612*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf))*cos(Bf) *cos(Bf);Z=y/(Nf*cos(Bf));b2=(0.5+0.003369*cos(Bf)*cos(Bf))*sin(Bf)*cos(Bf);b4=0.25+(0.16161+0.00562*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b3=0.333333-(0.166667-0.001123*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b5=0.2-(0.166667-0.0088*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);Bsecond=BfSecond-(1-(b4-0.12*Z*Z)*Z*Z)*Z*Z*b2*rouSecond;//计算大地经度B,单位为秒B=Bsecond/3600;//用整度数表示Blsecond=(1-(b3-b5*Z*Z)*Z*Z)*Z*rouSecond;//计算经度差l,单位为秒l=lsecond/3600;//用整度数表示lL=L0+l;//计算大地经度L}double Fansuan::getB(){return B;}double Fansuan::getL(){return L;}void Fansuan::printLocation(){printf("反算得大地坐标为:大地纬度B= %f°大地经度L=%f°\n",B,L);}Fansuan::~Fansuan(){}11。