📄 gauss.cpp
字号:
//高斯投影族
#include "stdafx.h"
#include "math.h"
#include "iostream.h"
#include "project.h"
#include "gauss.h"
//==================================================
CBGauss::CBGauss(double long_radius,double short_radius,double L0,double W)
{
/* PRJPARAM PrjParam;
PrjParam.m_bRadius_a=TRUE;
PrjParam.m_bRadius_b=TRUE;
PrjParam.m_bLatitude1=FALSE;
PrjParam.m_bLatitude2=FALSE;
PrjParam.m_bCentralLatitude=FALSE;
PrjParam.m_bCentralLongitude=TRUE;
PrjParam.m_bCentralAngle=FALSE;
PrjParam.m_bSouthLatitude=FALSE;
PrjParam.m_bZonewidth=TRUE;
PrjParam.m_bMapScale=FALSE;
*/
a = long_radius;
b = short_radius;
central_longitude = L0;
zone_width = W;
}
//-------------------------------------------------------------------------------------
//高斯投影（正解）
PRJXY CGauss::Project_xy(double latitude,double longitude)
{
double x1,x2,x3,x4,x5,x6,x7,y1;
double A1,A2,A3,B1,B2,B3,B4;
double t,t2,t4,t6,g2,g4,r,l2,l3,l5;
double e2,N,S,L,DL,B,L0,cosB2,cosB4;
PRJXY Gauss;
E_value E;
E = Get_e();
e2 = E.e2;
N = N_latitude(latitude);
S = Isodistance_S(latitude);
// int zone_number=(central_longitude+zone_width/2.0)/zone_width;
int zone_number = central_longitude/zone_width + 1;
L = Dms_to_arc(longitude);
B = Dms_to_arc(latitude);
L0 = Dms_to_arc(central_longitude);
// because GPS data is arc, so need not Dms_to_arc
/*
L=longitude;
B=latitude;
L0=central_longitude;
*/
DL = L-L0;
l2 = DL*DL;
l3 = l2*DL;
l5 = l2*l3;
cosB2 = cos(B)*cos(B);
cosB4 = cosB2*cosB2;
t = tan(B);
t2 = t*t;
t4 = t2*t2;
t6 = t2*t4;
g2 = e2*cos(B)*cos(B);
g4 = g2*g2;
A1 = 1.0/2.0*N*t;
A2 = 1.0/24.0*N*t*(5.0-t2+9.0*g2+4.0*g4);
A3 = 1.0/720.0*N*t*(61.0-58.0*t2+t4+270.0*g2-330.0*t2*g2);
B1 = N;
B2 = 1.0/6.0*N*(1.0-t2+g2);
B3 = 1.0/120.0*N*(5.0-18.0*t2+t4+14.0*g2-58.0*t2*g2);
B4 = 1.0/5040.0*N*(61.0-479.0*t2+179.0*t4-t6);
x1 = cos(B)*DL;
x2 = x1*x1;
x3 = x1*x2;
x4 = x3*x1;
x5 = x2*x3;
x6 = x2*x4;
x7 = x6*x1;
y1 = B1*x1+B2*x3+B3*x5+B4*x7+500000.0;
r = DL*sin(B)+1.0/3.0*l3*sin(B)*cosB2*(1.0+3.0*g2+4.0*g4)+1.0/15.0*sin(B)*cosB4*(2.0-t2);
Gauss.x = (S+A1*x2+A2*x4+A3*x6)/1000.0;
Gauss.y = (zone_number*1000000.0+y1)/1000.0;
Gauss.m = 1.0+1.0/2.0*(1.0+g2)*x2+1.0/24.0*(5.0-4.0*cos(B)*cos(B));
Gauss.n = Arc_to_dms(r);
Gauss.p = Gauss.m*Gauss.m;
Gauss.w = Arc_to_dms(0.0);
return(Gauss);
}
//-------------------------------------------------------------------------------------------------
//高斯投影（反解）
PRJBL CGauss::Project_bl(double x_value,double y_value)
{
double B0,B1,B2,B3,A1,A2,A3;
double N,N2,N3,N4,N5,N6,L0,DL;
double t,t2,t4,g2,g4,x,y,y2,y3,y4,y5,y6,e2;
double zone_number;
PRJBL Gauss;
x = x_value*1000.0;
y = modf(y_value/1000.0,&zone_number)*1000000.0-500000.0;
L0 = zone_number*zone_width-zone_width/2.0;
L0 = Dms_to_arc(L0);
B0 = Isodistance_B(x);
N = N_latitude(B0);
N2 = N*N;
N3 = N*N2;
N4 = N2*N2;
N5 = N2*N3;
N6 = N2*N4;
E_value E;
E = Get_e();
e2 = E.e2;
B0 = Dms_to_arc(B0);
t = tan(B0);
t2 = t*t;
t4 = t2*t2;
g2 = e2*cos(B0)*cos(B0);
g4 = g2*g2;
A1 = 1.0/2.0/N2*t*(-1.0-g2);
A2 = 1.0/24.0/N4*t*(5.0+3.0*t2+6.0*g2-6.0*t2*g2-3.0*g4-9.0*t2*g4);
A3 = 1.0/720.0/N6*t*(-61.0-90.0*t2-45.0*t4-107.0*g2+162.0*t2*g2+45.0*t4*g2);
B1 = 1.0/N/cos(B0);
B2 = -1.0/6.0/N3/cos(B0)*(1.0+2.0*t2+g2);
B3 = 1.0/120.0/N5/cos(B0)*(5.0+28.0*t2+24.0*t4+6.0*g2+8.0*t2*g2);
y2 = y*y;
y3 = y2*y;
y4 = y*y3;
y5 = y2*y3;
y6 = y3*y3;
DL = B1*y+B2*y3+B3*y5;
// no need to change to dms
Gauss.B = Arc_to_dms(B0+A1*y2+A2*y4+A3*y6);
Gauss.L = Arc_to_dms(DL+L0);
return(Gauss);
}
//==============================================
//双标准经线等角割圆锥投影（正解）
PRJXY CDLC::Project_xy(double latitude,double longitude)
{
double x0,l1,l2,l3,l4,l5,l6;
double A1,A2,A3,B1,B2,B3;
double t,t2,t4,g2,g4,r,y1;
double cosB2,cosB3,cosB4,cosB5,cosB6;
double A0,B0,C0,D0,E0;
double e2,e4,e6,e8;
double N,S,L,DL,B,L0,L2,L4;
PRJXY dlc;
E_value E;
E = Get_e();
e2 = E.e2;
e4 = e2*e2;
e6 = e4*e2;
e8 = e6*e2;
N = N_latitude(latitude);
S = Isodistance_S(latitude);
int zone_number=(central_longitude+zone_width/2.0)/zone_width;
L = Dms_to_arc(longitude);
B = Dms_to_arc(latitude);
L0 = Dms_to_arc(central_longitude);
DL = L-L0;
t = tan(B);
t2 = t*t;
t4 = t2*t2;
g2 = e2*cos(B)*cos(B);
g4 = g2*g2;
cosB2 = cos(B)*cos(B);
cosB3 = cos(B)*cosB2;
cosB4 = cos(B)*cosB3;
cosB5 = cos(B)*cosB4;
cosB6 = cos(B)*cosB5;
l1 = Dms_to_arc(zone_width/2.0);
L2 = l1*l1;
L4 = L2*L2;
A0 = 1.0/2.0+1.0/16.0*e2+3.0/128.0*e4+25.0/2048.0*e6+245.0/32768.0*e8;
B0 = 1.0/2.0-3.0/256.0*e4-5.0/512.0*e6-245.0/32768.0*e8;
C0 = -1.0/16.0*e2-3.0/128.0*e4-5.0/512.0*e6-35.0/8192.0*e8;
D0 = 3.0/256.0*e4+5.0/512.0*e6+455.0/65536.0*e8;
E0 = -5.0/2048.0*e6-105.0/32768.0*e8;
x0 = S-1.0/2.0*a*L2*(A0*B+B0/2.0*sin(2.0*B)+C0/4.0*sin(4.0*B)+D0/6.0*sin(6.0*B)+E0/8.0*sin(8.0*B));
A1 = 1.0/2.0*N*t*cosB2*(1.0-L2/2.0*cosB2*(3.0+7.0*g2));
A2 = 1.0/24.0*N*t*cosB4*((5.0-t2+9.0*g2+4.0*g4)-1.0/2.0*L2*cosB2*(33.0-27.0*t2+182.0*g2));
A3 = 1.0/720.0*N*t*cosB6*(61.0-58.0*t2+t4);
B1 = N*cos(B)*(1.0-1.0/2.0*L2*cosB2*(1.0+g2));
B2 = 1.0/6.0*N*cosB3*((1.0-t2+g2)-1.0/2.0*L2*cosB2*(3.0-9.0*t2+10.0*g2-41.0*t2*g2));
B3 = 1.0/120.0*N*cosB5*((5.0-18.0*t2+t4+14.0*g2-58.0*t2*g2)-1.0/2.0*L2*cosB2*(33.0-246.0*t2));
l2 = DL*DL;
l3 = DL*l2;
l4 = l2*l2;
l5 = l2*l3;
l6 = l2*l4;
r = sin(B)*(1.0-L2*cosB2*(1.0+3.0*g2)-1.0/2.0*L4*cosB4*(1.0+4.0*g2))*DL;
r += 1.0/3.0*sin(B)*cosB2*(1.0+3.0*g2+2.0*g4-1.0/4.0*L2*cosB2*(16.0-8.0*t2+120.0*g2+103*g2*t2))*l3;
r += 1.0/60.0*sin(B)*cosB4*(8.0-4.0*t2-12.0*t4-75.0*g2+105.0*g2*t2)*l5;
y1 = (B1*DL+B2*l3+B3*l5)+500000.0;
dlc.x = (x0+A1*l2+A2*l4+A3*l6)/1000.0;
dlc.y = (zone_number*1000000.0+y1)/1000.0;
dlc.m = 1.0+1.0/2.0*(l2-L2)*cosB2*(1.0+g2)-1.0/4.0*L2*l2*cosB4*(3.0-4.0*t2)+1.0/24.0*l4*cosB4*(5.0-4.0*t2);
dlc.n = Arc_to_dms(r);
dlc.p = dlc.m*dlc.m;
dlc.w = 0.0;
return(dlc);
}
//-------------------------------------------------------------------------------------------------
//双标准经线等角割圆锥投影（反解）
PRJBL CDLC::Project_bl(double x_value,double y_value)
{
double e2;
double B0,B1,B2,B3,A1,A2,A3;
double M,N,N2,N3,N4,N5,L0,DL,l1,L2,cosB2;
double t,t2,t4,g2,g4,x,y,y2,y3,y4,y5,y6;
double zone_number;
PRJBL dlc;
E_value E;
E = Get_e();
e2 = E.e2;
x = x_value*1000.0;
y = modf(y_value/1000.0,&zone_number)*1000000.0-500000.0;
L0 = zone_number*zone_width-zone_width/2.0;
L0 = Dms_to_arc(L0);
l1 = Dms_to_arc(zone_width/2.0);
L2 = l1*l1;
B0 = Isodistance_B(x);
M = M_longitude(B0);
N = N_latitude(B0);
N2 = N*N;
N3 = N*N2;
N4 = N2*N2;
N5 = N2*N3;
B0 = Dms_to_arc(B0);
cosB2 = cos(B0)*cos(B0);
t = tan(B0);
t2 = t*t;
t4 = t2*t2;
g2 = e2*cosB2;
g4 = g2*g2;
A1 = 1.0/2.0/M/N*(t-2.0*L2*t*cosB2*g2);
A2 = 1.0/24.0/M/N3*(t*(5.0+3.0*t2+g2-9.0*g2*t2)+4.0*L2*t*cosB2*(1.0+t2));
A3 = 1.0/720.0/M/N5*t*(61.0+90.0*t2+45.0*t4);
B1 = 1.0/N/cos(B0)*(1.0-1.0/2.0*L2*cosB2*(1.0+g2));
B2 = 1.0/6.0/N3/cos(B0)*(1.0+2.0*t2+g2+1.0/2.0*L2*cosB2*(1.0+2.0*t2-2.0*g2+26.0*g2*t2));
B3 = 1.0/120.0/N5/cos(B0)*(5.0+28.0*t2+24.0*t4+6.0*g2+8.0*t2*g2+1.0/2.0*L2*cosB2*(13.0+52.0*t2+40.0*t4));
y2 = y*y;
y3 = y2*y;
y4 = y*y3;
y5 = y2*y3;
y6 = y3*y3;
DL = B1*y-B2*y3+B3*y5;
dlc.B = Arc_to_dms(B0-A1*y2+A2*y4-A3*y6);
dlc.L = Arc_to_dms(DL+L0);
return(dlc);
}
//=======================================================
//UTM投影（正解）
PRJXY CUTM::Project_xy(double latitude,double longitude)
{
double x1,x2,x3,x4,x5,y1;
double A1,A2,B1,B2,B3;
double t,t2,t4,t6,g2,g4;
double e2,N,S,L,DL,B,L0;
PRJXY utm;
E_value E;
E = Get_e();
e2 = E.e2;
N = N_latitude(latitude);
S = Isodistance_S(latitude);
int zone_number=(central_longitude+zone_width/2.0)/zone_width;
L = Dms_to_arc(longitude);
B = Dms_to_arc(latitude);
L0 = Dms_to_arc(central_longitude);
DL = L-L0;
t = tan(B);
t2 = t*t;
t4 = t2*t2;
t6 = t2*t4;
g2 = e2*cos(B)*cos(B);
g4 = g2*g2;
A1 = 1.0/2.0*N*t;
A2 = 1.0/24.0*N*t*(5.0-t2+9.0*g2+4.0*g4);
B1 = N;
B2 = 1.0/6.0*N*(1.0-t2+g2);
B3 = 1.0/120.0*N*(5.0-18.0*t2+t4);
x1 = cos(B)*DL;
x2 = x1*x1;
x3 = x1*x2;
x4 = x3*x1;
x5 = x2*x3;
y1 = 0.9996*(B1*x1+B2*x3+B3*x5)/1000.0+500.0;
utm.x = 0.9996*(S+A1*x2+A2*x4)/1000.0;
utm.y = (zone_number*1000000.0+y1)/1000.0;
utm.m = 0.9996*(1.0+1.0/2.0*(1.0+g2)*x2+1.0/6.0*(2.0-t2)*x4-1.0/8.0*x4);
utm.n = 0.9996*(1.0+1.0/2.0*(1.0+g2)*x2+1.0/6.0*(2.0-t2)*x4-1.0/8.0*x4);
utm.p = utm.m*utm.n;
utm.w = Arc_to_dms(0.0);
return(utm);
}
//-------------------------------------------------------------------------------------------------
//UTM投影（反解）
PRJBL CUTM::Project_bl(double x_value,double y_value)
{
double B0,B1,B2,B3,A1,A2,A3;
double N,N2,N3,N4,N5,N6,L0,DL;
double t,t2,t4,g2,g4,x,y,y2,y3,y4,y5,y6,e2;
double zone_number;
PRJBL utm;
x = x_value*1000.0;
y = modf(y_value/1000.0,&zone_number)*1000000.0-500000.0;
L0 = zone_number*zone_width-zone_width/2.0;
L0 = Dms_to_arc(L0);
B0 = Isodistance_B(x);
N = N_latitude(B0);
N2 = N*N;
N3 = N*N2;
N4 = N2*N2;
N5 = N2*N3;
N6 = N2*N4;
E_value E;
E = Get_e();
e2 = E.e2;
B0 = Dms_to_arc(B0);
t = tan(B0);
t2 = t*t;
t4 = t2*t2;
g2 = e2*cos(B0)*cos(B0);
g4 = g2*g2;
A1 = 1.0/2.0/N2*t*(-1.0-g2);
A2 = 1.0/24.0/N4*t*(5.0+3.0*t2+6.0*g2-6.0*t2*g2-3.0*g4-9.0*t2*g4);
A3 = 1.0/720.0/N6*t*(-61.0-90.0*t2-45.0*t4-107.0*g2+162.0*t2*g2+45.0*t4*g2);
B1 = 1.0/N/cos(B0);
B2 = -1.0/6.0/N3/cos(B0)*(1.0+2.0*t2+g2);
B3 = 1.0/120.0/N5/cos(B0)*(5.0+28.0*t2+24.0*t4+6.0*g2+8.0*t2*g2);
y2 = y*y;
y3 = y2*y;
y4 = y*y3;
y5 = y2*y3;
y6 = y3*y3;
DL = B1*y+B2*y3+B3*y5;
utm.B = Arc_to_dms(B0+A1*y2+A2*y4+A3*y6);
utm.L = Arc_to_dms(DL+L0);
return(utm);
}
//=========================================================
//next is the main program for console
/* void main(void)
{
double B,L;
PRJXY P;
PRJBL G;
CProject *point ;
CGauss Gauss(6378245.0,6356863.0,111.0,6.0);
point = &Gauss;
cout<<"Please input coodinators for projection:"<<endl;
cout<<"Latitude=";
cin>>B;
cout<<"Longitude=";
cin>>L;
cout<<endl;
P=point->Project_xy(B,L);
cout.precision(12);
cout<<"x="<<P.x<<"(km)"<<endl;
cout<<"y="<<P.y<<"(km)"<<endl<<endl;
cout<<"m="<<P.m<<endl;
cout<<"n="<<P.n<<endl;
cout<<"p="<<(P.p-1.0)*100.0<<"(%)"<<endl;
cout<<"w="<<P.w<<"(dms)"<<endl<<endl;
G=point->Project_bl(P.x,P.y);
cout<<"B="<<G.B<<"(dms)"<<endl;
cout<<"L="<<G.L<<"(dms)"<<endl;
}
*/⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -