大地坐标向笛卡尔坐标转换
高斯正反算
采用不同椭球实现高斯克里格投影，将经纬度坐标转换为高斯平面坐标：正算
高斯平面坐标转换为不同椭球下的经纬度坐标：反算
http://robekeane.iteye.com/blog/1441566
http://blog.163.com/superliuwhu@126/blog/static/68557075200992654039167/
1 void GaussProjectDirect(double a,double efang,double B,double L,double L0,double& x,double &y,double& R)//高斯投影正算克氏 2 { 3 4 double b=aefangtob(a,efang); 5 double e2=seconde(a,b); 6 double W=sqrt(1-efang*sin(B)*sin(B));printf("W=%f",W); 7 double N=a/W;printf("N=%f",N); 8 double M=a*(1-efang)/pow(W,3);printf("M=%f",M); 9 double t=tee(B); 10 double eitef=eitefang(a,b,B); 11 double l=L-L0; 12 //主曲率半径计算 13 double m0,m2,m4,m6,m8,n0,n2,n4,n6,n8; 14 m0=a*(1-efang); n0=a; 15 m2=3.0/2.0*efang*m0; n2=1.0/2.0*efang*n0; 16 m4=5.0/4.0*efang*m2; n4=3.0/4.0*efang*n2; 17 m6=7.0/6.0*efang*m4; n6=5.0/6.0*efang*n4; 18 m8=9.0/8.0*efang*m6; n8=7.0/8.0*efang*n6; 19 //子午线曲率半径 20 double a0,a2,a4,a6,a8; 21 a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128; 22 a2=m2/2+m4/2+15*m6/32+7.0/16.0*m8; 23 a4=m4/8.0+3.0/16.0*m6+7.0/32.0*m8; 24 a6=m6/32+m8/16; 25 a8=m8/128; 26 27 double X=a0*B-a2/2*sin(2*B)+a4/4*sin(4*B)-a6/6*sin(6*B)+a8/8*sin(8*B); 28 x=X+N/2*t*cos(B)*cos(B)*l*l+N/24*t*(5-t*t+9*eitef+4*pow(eitef,2))*pow(cos(B),4)*pow(l,4)+N/720*t*(61-58*t*t+pow(t,4))*pow(cos(B),6)*pow(l,6); 29 y=N*cos(B)*l+N/6*(1-t*t+eitef)*pow(cos(B),3)*pow(l,3)+N/120*(5-18*t*t+pow(t,4)+14*eitef-58*eitef*t*t)*pow(cos(B),5)*pow(l,5); 30 R=sqrt(M*N); 31 } 32 33 34 //高斯投影反算 35 36 37 void GaussProjectInvert(double a,double efang,double x,double y,double L0,double &B,double& L,double& R) 38 { 39 double b=aefangtob(a,efang); 40 41 42 double m0,m2,m4,m6,m8,n0,n2,n4,n6,n8; 43 m0=a*(1-efang); n0=a; 44 m2=3.0/2.0*efang*m0; n2=1.0/2.0*efang*n0; 45 m4=5.0/4.0*efang*m2; n4=3.0/4.0*efang*n2; 46 m6=7.0/6.0*efang*m4; n6=5.0/6.0*efang*n4; 47 m8=9.0/8.0*efang*m6; n8=7.0/8.0*efang*n6; 48 49 50 //子午线曲率半径 51 double a0,a2,a4,a6,a8; 52 a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128; 53 a2=m2/2+m4/2+15*m6/32+7.0/16.0*m8; 54 a4=m4/8.0+3.0/16.0*m6+7.0/32.0*m8; 55 a6=m6/32+m8/16; 56 a8=m8/128; 57 58 59 double X=x; 60 double FBf=0; 61 double Bf0=X/a0,Bf1=0; 62 while((Bf0-Bf1)>=0.0001) 63 { Bf1=Bf0; 64 FBf=a0*Bf0-a2/2*sin(2*Bf0)+a4/4*sin(4*Bf0)-a6/6*sin(6*Bf0)+a8/8*sin(8*Bf0); 65 Bf0=(X-FBf)/a0; 66 } 67 double Bf=Bf0; 68 double Vf=bigv(a,b,Bf); 69 double tf=tee(Bf); 70 double Nf=bign(a,b,Bf); 71 double eiteffang=eitefang(a,b,Bf); 72 double Bdu=rad_deg(Bf)-1/2.0*Vf*Vf*tf*(pow((y/Nf),2)-1.0/12*(5+3*tf*tf+eiteffang-9*eiteffang*tf*tf)*pow((y/Nf),4)+1.0/360.0*(61+90*tf*tf+45*tf*tf)*pow((y/Nf),6))*180/PI; 73 double ldu=1.0/cos(Bf)*(y/Nf+1.0/6.0*(1+2*tf*tf+eiteffang)*pow((y/Nf),3)+1.0/120.0*(5+28*tf*tf+24*tf*tf+6*eiteffang+8*eiteffang*tf*tf)*pow((y/Nf),5))*180.0/PI; 74 75 76 B=deg_int(Bdu); 77 L=L0+deg_int(ldu); 78 double W=sqrt(1-efang*sin(B)*sin(B));printf("W=%f ",W); 79 double N=a/W;printf("N=%f ",N); 80 double M=a*(1-efang)/pow(W,3);printf("M=%f ",M); 81 R=sqrt(M*N); 82 83 84 }
