高斯投影正算C++源代码

// PriPoint.cpp: implementation of the PriPoint class. // ////////////////////////////////////////////////////////////////////// #include<cmath> #include<iostream.h> #include "PrjPoint.h" #include "function.h" ////////////////////////////////////////////////////////////////////// // Construction/Destruction ////////////////////////////////////////////////////////////////////// CPrjPoint::CPrjPoint() { } CPrjPoint::~CPrjPoint() { } bool CPrjPoint::xy2BL() { double sinB, cosB, t, t2, N ,ng2, V, yN; double preB0, B0; double eta; y -= 500000; B0 = x / A1; do { preB0 = B0; B0 = B0 * PI / 180.0; //?????? B0 = (x - (A2 * sin(2 * B0) + A3 * sin(4 * B0) + A4 * sin(6 * B0))) / A1; eta = fabs(B0 - preB0); } while(eta > 0.000000001); B0 = B0 * PI / 180.0; B = Rad2Dms(B0); sinB = sin(B0); cosB = cos(B0); t = tan(B0); t2 = t * t; N = a / sqrt(1 - e2 * sinB * sinB); ng2 = cosB * cosB * e2 / (1 - e2); V = sqrt(1 + ng2); yN = y / N; B = B0 - (yN * yN - (5 + 3 * t2 + ng2 - 9 * ng2 * t2) * yN * yN * yN * yN / 12.0 + (61 + 90 * t2 + 45 * t2 * t2) * yN * yN * yN * yN * yN * yN / 360.0) * V * V * t / 2; L = L0 + (yN - (1 + 2 * t2 + ng2) * yN * yN * yN / 6.0 + (5 + 28 * t2 + 24 * t2 * t2 + 6 * ng2 + 8 * ng2 * t2) * yN * yN * yN * yN * yN / 120.0) / cosB; return true; } bool CPrjPoint::SetL0(double dL0) { L0 = Dms2Rad(dL0); //L0为弧度 return true; } bool CPrjPoint::SetBL(double dB, double dL) { B = Dms2Rad(dB); L = Dms2Rad(dL); // B = dB; //why??? // L = dL; return true; } bool CPrjPoint::GetBL(double *dB, double *dL) //notice the '*' { *dB = Rad2Dms(B); *dL = Rad2Dms(L); return true; } bool CPrjPoint::Setxy(double dx, double dy) { x = dx; y = dy; xy2BL(); return true; } bool CPrjPoint::Getxy(double *dx, double *dy) //notice the '*' { *dx = x; *dy = y; return true; } ///////////////////////////////////////////////////// /**************************************************************************** ** 函数名：LLtoUTM ** 输入： a 椭球体长半轴 ** f 椭球体扁率 f=(a-b)/a 其中b代表椭球体的短半轴 ** Lat 经过UTM投影之前的纬度 ** Long 经过UTM投影之前的经度 ** LongOrigin 中央经度线 ** FN 纬度起始点，北半球为0，南半球为10000000.0m ** UTMNorthing 经过UTM投影后的纬度方向的坐标 ** UTMEasting 经过UTM投影后的经度方向的坐标 ** 功能描述：UTM投影正转 ** 作者： CiCi ** 单位： 中国地质大学（北京）地球科学与资源学院 ** 创建日期：2007年3月28日 ** 版本：1.0 ***************************************************************************/ //LLtoUTM(double a,double f,const double Lat, const double Long, //double LongOrigin,double FN,double &UTMNorthing,double &UTMEasting) bool CPrjPoint::LLtoUTM() //正算 //实验证明这函数是对的 { //e表示WGS84第一偏心率,eSquare表示e的平方, double eSquare =(2*f-f*f) ; double k0 = 0.9996; double e2Square; double V, T, C, A, M; //确保longtitude位于-180.00----179.9之间 //double LongTemp = (Long+180)-int((Long+180)/360)*360-180; double LatRad = this->B;//Lat*PI/180; double LongRad =this->L;// LongTemp*PI/180; double LongOriginRad; LongOriginRad = this->L0;//LongOrigin * PI/180; e2Square = (eSquare)/(1-eSquare); V = a/sqrt(1-eSquare*sin(LatRad)*sin(LatRad)); T = tan(LatRad)*tan(LatRad); C = e2Square*cos(LatRad)*cos(LatRad); A = cos(LatRad)*(LongRad-LongOriginRad); M = a*((1-eSquare/4-3*eSquare*eSquare/64-5*eSquare*eSquare*eSquare/256)*LatRad -(3*eSquare/8+3*eSquare*eSquare/32+45*eSquare*eSquare*eSquare/1024)*sin(2*LatRad) +(15*eSquare*eSquare/256+45*eSquare*eSquare*eSquare/1024)*sin(4*LatRad) -(35*eSquare*eSquare*eSquare/3072)*sin(6*LatRad)); UTMEasting = (double)(k0*V*(A+(1-T+C)*A*A*A/6 + (5-18*T+T*T+72*C-58*e2Square)*A*A*A*A*A/120)+ 500000.0); UTMNorthing = (double)(k0*(M+V*tan(LatRad)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24 +(61-58*T+T*T+600*C-330*e2Square)*A*A*A*A*A*A/720))); //南半球纬度起点为10000000.0m //UTMNorthing=UTMNorthing+FN; return true; } bool CPrjPoint::UTMtoLL() { //UTMtoLL(int ReferenceEllipsoid, const char* UTMZone, //double&amt; Lat, double&amt; Long ) //converts UTM coords to lat/long. Equations from USGS Bulletin 1532 //East Longitudes are positive, West longitudes are negative. //North latitudes are positive, South latitudes are negative //Lat and Long are in decimal degrees. //Written by Chuck Gantz- chuck.gantz@globalstar.com double k0 = 0.9996; //double a = ellipsoid[ReferenceEllipsoid].EquatorialRadius; //a已经有了 // double eccSquared = ellipsoid[ReferenceEllipsoid].eccentricitySquared; //eccentricitySquared double eccSquared = (2*f-f*f); //right // cout<<"eccSquared "<<eccSquared<<endl; double eccPrimeSquared; double e1 = (1-sqrt(1-eccSquared))/(1+sqrt(1-eccSquared)); double N1, T1, C1, R1, D, M; double LongOrigin ; double mu, phi1, phi1Rad; double xx, yy; int ZoneNumber = 49; char* ZoneLetter; int NorthernHemisphere; //1 for northern hemispher, 0 for southern double Lat,Long; double temp;//use for some variability. int tDec; int tMin; /* x = UTMEasting ; ZoneNumber=0; for(int iLoop=0;iLoop<60;iLoop++) if(x>=833978.55644) { x-=833978.55644; ZoneNumber++; } */ xx = UTMEasting - 500000.0; //remove 500,000 meter offset for longitude yy = UTMNorthing; // ZoneNumber = strtoul(UTMZone, &amt;ZoneLetter, 10); // if((*ZoneLetter - 'N') >= 0) NorthernHemisphere = 1;//point is in northern hemisphere//只能计算北半球的 // else // { // NorthernHemisphere = 0;//point is in southern hemisphere // y -= 10000000.0;//remove 10,000,000 meter offset used for southern hemisphere // } LongOrigin = (ZoneNumber - 1)*6 - 180 + 3; //+3 puts origin in middle of zone cout<<"LongOrigin: "<<LongOrigin<<endl; eccPrimeSquared = (eccSquared)/(1-eccSquared); M = yy / k0; mu = M/(a*(1-eccSquared/4-3*eccSquared*eccSquared/64-5*eccSquared*eccSquared*eccSquared/256)); phi1Rad = mu + (3*e1/2-27*e1*e1*e1/32)*sin(2*mu) + (21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(4*mu) +(151*e1*e1*e1/96)*sin(6*mu); // phi1 = phi1Rad*rad2deg; origin N1 = a/sqrt(1-eccSquared*sin(phi1Rad)*sin(phi1Rad)); T1 = tan(phi1Rad)*tan(phi1Rad); C1 = eccPrimeSquared*cos(phi1Rad)*cos(phi1Rad); R1 = a*(1-eccSquared)/pow(1-eccSquared*sin(phi1Rad)*sin(phi1Rad), 1.5); D = xx/(N1*k0); Lat = phi1Rad - (N1*tan(phi1Rad)/R1)*(D*D/2-(5+3*T1+10*C1-4*C1*C1-9*eccPrimeSquared)*D*D*D*D/24 +(61+90*T1+298*C1+45*T1*T1-252*eccPrimeSquared-3*C1*C1)*D*D*D*D*D*D/720); // this->B = this->B * rad2deg; origin this->Bdms = Rad2Dms(Lat); Long = (D-(1+2*T1+C1)*D*D*D/6+(5-2*C1+28*T1-3*C1*C1+8*eccPrimeSquared+24*T1*T1) *D*D*D*D*D/120)/cos(phi1Rad); // Long = LongOrigin + Long * rad2deg; origin temp = Rad2Dms(Long);//myself // this->Ldms = LongOrigin*10000 + Rad2Dms(Long); this->Ldms = DmsAddDms(LongOrigin*10000.0,temp); cout<<temp<<" "<<LongOrigin*10000<<endl; return true; } ///////////////////////////////////////////////////// /////////////////////////////////////////////////// // Definition of PrjPoint_IUGG1975 /////////////////////////////////////////////////// CPrjPoint_WGS84::CPrjPoint_WGS84() { a = 6378137; b = 6356752.3142; f = 1.0/298.257223563; e2 = 1 - ((f - 1) / f) * ((f - 1) / f); e12 = (f / (f - 1)) * (f / (f - 1)) - 1; n=(a-b)/(a+b); A1 = a*(1-n+(pow(n,2)-pow(n,3))*5/4+(pow(n,4)-pow(n,5))*81/64)*PI/180; A2 = (n-pow(n,2)+(pow(n,3)-pow(n,4))*7/8+pow(n,5)*55/64)*3*a/2; A3 = (pow(n,2)-pow(n,3)+(pow(n,4)-pow(n,5))*3/4)*15*a/16; A4 = (pow(n,3)-pow(n,4)+pow(n,5)*11/16)*35*a/48; } CPrjPoint_WGS84::~CPrjPoint_WGS84() { } bool CPrjPoint::BL2XY_Gauss() //已通过测试，WGS84坐标系统 { double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval; //double a,f