高斯投影正反算程序

// GKprojection.cpp : 定义控制台应用程序的入口点。
//


#include "math.h"
#define pi 3.141592653

double Torad(double);
double Todouble(double);

int _tmain(int argc, _TCHAR* argv[])
{
int model;
printf("高斯投影正算请输入1，反算请输入2\n");
scanf("%d",&model);
if(model==1)
{
printf("高斯投影正算\nunder WGS-84,六度带\n测试数据：B=22.155898294，L=111.285215387\n");
double a=6378137,f=0.00335281066474;
double B=22.155898294,L=111.285215387;
B=Torad(B);
double l=L-(int(L/6)*6+3);
int n=int(L/6)+1;
l=Torad(l);
double e2=-f*f+2*f;
double m0=a*(1-e2);
double m2=3*e2*m0/2;
double m4=5*e2*m2/4;
double m6=7*e2*m4/6;
double m8=9*e2*m6/8;
double a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;
double a2=m2/2+m4/2+15*m6/32+7*m8/16;
double a4=m4/8+3*m6/16+7*m8/32;
double a6=m6/32+m8/16;
double a8=m8/128;
double X=a0*B-a2*sin(2*B)/2+a4*sin(4*B)/4-a6*sin(6*B)/6+a8*sin(8*B)/8;
double W=sqrt(1-e2*sin(B)*sin(B));
double N=a/W;
double ep2=e2/(1-e2);
double ita2=ep2*cos(B)*cos(B);
double t=tan(B);
double temp1=5-t*t+9*ita2+4*ita2*ita2;
double temp2=61-58*t*t+pow(t,4);
double temp3=1-t*t+ita2;
double temp4=5-18*t*t+pow(t,4)+14*ita2-58*t*t*ita2;
double x=X+N*sin(B)*cos(B)*l*l/2+N*sin(B)*pow(cos(B),3)*temp1*pow(l,4)/24+N*sin(B)*pow(cos(B),5)*temp2*pow(l,6)/720;
double y=N*cos(B)*l+N*pow(cos(B),3)*temp3*pow(l,3)/6+N*pow(cos(B),5)*temp4*pow(l,5)/120;
y=500000+y;
printf("计算结果x=%f,y=%d%f\n",x,n,y);
}
else
{

printf("高斯投影反算\nunder CGCS2000,六度带\n测试数据：x=2562083.2708,y=19512837.2851\n");
double a=6378137;
double f=1/298.257222101;
double e2=-f*f+2*f;
double x=2562083.2708;
double y=19512837.2851;
int Y=int(y);
int count=0;
while (Y>=100)
{
Y=Y/10;
count++;
}
y=y-Y*pow(10.0,count)-500000;
double m0=a*(1-e2);
double m2=3*e2*m0/2;
double m4=5*e2*m2/4;
double m6=7*e2*m4/6;
double m8=9*e2*m6/8;
double a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;
double a2=m2/2+m4/2+15*m6/32+7*m8/16;
double a4=m4/8+3*m6/16+7*m8/32;
double a6=m6/32+m8/16;
double a8=m8/128;
double Bf0=x/a0;
double Bf=Bf0;
double fBf=-a2*sin(2*Bf)/2+a4*sin(4*Bf)/4-a6*sin(6*Bf)/6+a8*sin(8*Bf)/8;
Bf=(x-fBf)/a0;
while (fabs(Bf-Bf0)>Torad(0.000000000001))
{
Bf0=Bf;
fBf=-a2*sin(2*Bf)/2+a4*sin(4*Bf)/4-a6*sin(6*Bf)/6+a8*sin(8*Bf)/8;
Bf=(x-fBf)/a0;
}
double W=sqrt(1-e2*sin(Bf)*sin(Bf));
double N=a/W;
double M=a*(1-e2)/pow(W,3);
double ep2=e2/(1-e2);
double ita2=ep2*cos(Bf)*cos(Bf);
double t=tan(Bf);
double temp1=5+3*t*t+ita2-9*ita2*t*t;
double temp2=61+90*t*t+45*pow(t,4);
double temp3=1+2*t*t+ita2;
double temp4=5+28*t*t+6*ita2+24*ita2*ita2+8*t*t*ita2;
double B=Bf-t*y*y/(2*M*N)+t*temp1/(24*M*pow(N,3))-t*temp2*pow(y,6)/(720

*M*pow(N,5));
double l=y/(N*cos(Bf))-temp3*pow(y,3)/(6*pow(N,3)*cos(Bf))+temp4*pow(y,5)/(120*pow(N,5)*cos(Bf));
B=Todouble(B);
double L=Todouble(Torad(Y*6.0)+l);
printf("计算结果是B=%f,L=%f\n",B,L);

}

return 0;
}

double Torad(double s)
{
int A=int(s);
double temp1=(s-A)*100;
int B=int(temp1);
double C=(temp1-B)*100;
double result=(A+B/60.0+C/3600)*pi/180;
return result;
}

double Todouble(double s)
{
s=s*180/pi;
int A=int(s);
double temp=(s-A)*60;
int B=int(temp);
double C=(temp-B)*60;
double result=A+B*0.01+C*0.0001;
return result;
}