高斯正反算计算函数

1、./servey.h/ #ifndef SERVEY_H/ #define SERVEY_H#include #include #include const double PI = 3.149323846;const double epsilon = 0.;/角度（度、分、秒）化弧度(带符号）double angle_to_radian (double alfa) double alfa1,alfa2,fsign,fbeta;if( fabs(alfa) =0)alfa_sign = 1;elsealfa_sign = -1;alfa = fabs(alfa); double alfa1,al

2、fa2;double A = floor(alfa+epsilon);double B = floor(alfa-A)*100+epsilon);alfa1 = A+B/60;alfa2=(alfa*100.-floor(alfa*100.+epsilon)/36.;alfa1+=alfa2;return (alfa_sign*alfa1);/弧度化度分秒double radian_to_angle(double alfa)double alfa_sign; /alfa的正负号if(alfa=0)alfa_sign = 1;elsealfa_sign = -1;alfa = fabs(alfa

3、);double alfa1=alfa*180./PI+epsilon;double alfa2=floor(alfa1+epsilon)+floor(alfa1-floor(alfa1+epsilon)*60+epsilon)/100; double alfa3=(alfa1-angle_to_degree(alfa2)*0.36;alfa2+=alfa3;return (alfa_sign*alfa2);/高斯正算函数/Ellipse为枚举类型/Cent_Meridian为中央子午线经度/x,y为返回的高斯平面坐标int BL_xy(int ellipse, double Cent_Mer

4、idian, double B, double L, double& x, double& y)double a, b, e1, e2; /e1为第一偏心率，e2为第二偏心率 switch (ellipse) /选椭球 case 0 : /54椭球a=6378245.; /长半轴b=6356863.0187730473; /短半轴e1=sqrt(pow(a,2.)-pow(b,2.)/a;e2=sqrt(pow(a,2.)-pow(b,2.)/b;break; case 1 : /80椭球a=6378140.;b=6356755.2881575287;e1=sqrt(pow(a,2)-pow(

5、b,2)/a;e2=sqrt(pow(a,2)-pow(b,2)/b;break;case 2 : /84椭球a=6378137.;b=6356752.3142;e1=sqrt(pow(a,2)-pow(b,2)/a;e2=sqrt(pow(a,2)-pow(b,2)/b;break;default : break; double l=angle_to_degree(L)-angle_to_degree(Cent_Meridian); /l为经差 /角度化弧度B=angle_to_radian(B);L=angle_to_radian(L);l=l*PI/180.;double X0; /X0

6、为当l=0时，从赤道起算的子午线弧长 /计算子午线弧长X的系数double A0=1.0+3.0/4*pow(e1,2.0)+45.0/64*pow(e1,4.0)+350.0/512*pow(e1,6.0)+11025.0/16384*pow(e1,8.0);double A2=-(3.0/4*pow(e1,2.0)+60.0/64*pow(e1,4.0)+525.0/512*pow(e1,6.0)+17640.0/16384*pow(e1,8.0)/2.0;double A4=(15.0/64*pow(e1,4.0)+210.0/512*pow(e1,6.0)+8820.0/16384*p

7、ow(e1,8.0)/4;double A6=-(35.0/512*pow(e1,6.0)+2520.0/16384*pow(e1,8.0)/6;double A8=(315.0/16384.0*pow(e1,8.0)/8; /计算子午线弧长XX0=a*(1.0-pow(e1,2.0)*(A0*B+A2*sin(2*B)+A4*sin(4*B)+A6*sin(6*B)+A8*sin(8*B);double t=tan(B);double anke=e2*cos(B);double N=a/sqrt(1.0-pow(e1,2.0)*pow(sin(B),2.0); /N为投影点的卯酉圈曲率半径/

8、坐标计算x=X0+1.0/2*N*t*pow(cos(B),2.0)*pow(l,2.0)+1.0/24*N*t*(5.0-pow(t,2.0)+9.0*pow(anke,2.0)+4.0*pow(anke,4.0)*pow(cos(B),4.0)*pow(l,4.0)+1.0/720*N*t*(61.0-58.0*pow(t,2.0)+pow(t,4.0)+270*pow(anke,2.0)-330.0*pow(anke,2.0)*pow(t,2.0)*pow(cos(B),6.0)*pow(l,6.0);y=N*cos(B)*l+1.0/6*N*(1-pow(t,2.0)+pow(anke

9、,2.0)*pow(cos(B),3.0)*pow(l,3.0) +1.0/120*N*(5.0-18.0*pow(t,2.0)+pow(t,4.0)+14.0*pow(anke,2.0)-58*pow(anke,2)*pow(t,2)*pow(cos(B),5.0)*pow(l,5.0);y+=500000.0; /add 500kmreturn 0;/高斯反算函数 int xy_BL(int ellipse, double Cent_Meridian, double x, double y, double& B, double& L)double a, b, e1, e2; /e1为第一偏

10、心率，e2为第二偏心率 switch (ellipse) /选椭球 case 0 : /54椭球a=6378245.; /长半轴b=6356863.0187730473; /短半轴e1=sqrt(pow(a,2.)-pow(b,2.)/a;e2=sqrt(pow(a,2.)-pow(b,2.)/b;break; case 1 : /80椭球a=6378140.;b=6356755.2881575287;e1=sqrt(pow(a,2)-pow(b,2)/a;e2=sqrt(pow(a,2)-pow(b,2)/b;break;case 2 : /84椭球a=6378137.;b=6356752.

11、3142;e1=sqrt(pow(a,2)-pow(b,2)/a;e2=sqrt(pow(a,2)-pow(b,2)/b;break;default : break;double Bf; /底点纬度/以下计算Bf的系数B0,K0,K2,K4,K6double E=e1*e1;double A0=1.0+3./4*E+45./64*E*E+350./512*E*E*E+11025./16384*E*E*E*E +43659./65536*E*E*E*E*E;double B0=x/(a*(1-E)*A0);double K0=(3./4*E+45./64*E*E+350./512*E*E*E+1

12、1025./16384*pow(e1,8.)/2;double K2=-(63./64*E*E+1108./512*pow(e1,6.)+58239./16384*pow(e1,8.)/3;double K4=(604./512*E*E*E+68484./16384*E*E*E*E)/3;double K6=-(26328./16384*pow(e1,8)/3; double si=sin(B0)*sin(B0);Bf=B0+sin(2.*B0)*(K0+si*(K2+si*(K4+K6*si);double tf=tan(Bf);double ankef=e2*cos(Bf);/double

13、 Nf=a/sqrt(1.0-E*pow(sin(Bf),2.0);double c=a/sqrt(1-E);double Vf=sqrt(1+ankef*ankef); double Nf=c/Vf;double Mf=a*(1.-E)/pow(1.-E*sin(Bf),1.5);y=y-500000.; double T=tf*tf, an=ankef*ankef, y2=y*y, MN=Mf*Nf; /B=Bf-tf*y2/(2*MN)+tf*(5.+3.*T+an-9.*an*T)*y2*y2/(24.*MN*Nf*Nf)/-tf*(61.+90*T+45.*T*T)*pow(y2,3.)/(720.*Mf*pow(Nf,5.);/另一种计算B的算法double q=y/Nf;B=Bf-pow(Vf,2.0)*tf*(1.0+(-(5+3.0*T+an*(1-9.0*T)+(61+45.0*T*(2+T)*pow(q,2.0)/30.0)*pow(q,2.0)/12.0)*pow(q,2.0)/2.0; double