坐标正算c语言,[转载]C语言版——高斯坐标正算程序

高斯正算

已知克拉索夫斯基椭球上一点大地坐标为

(B=22.365435915，L=113.022979765)，按六度带投

影求其高斯平面坐标(用函数实现)。

#include

#include

#define PI 3.1415926535897932

void main()

{

long double AngleToRadian(long double alfa);

long double MYQQLBJN(long double a,long double b,long double

B);

long double HCX(long double a,long double b,long double B);

void GSZS(long double a,long double b,long double l,long double

B);

long double a,b,B,L;

a=6378245.0000000000;b=6356863.0187730473;

B=22.365435915; L=113.022979765;

GSZS(a,b,L,B);

}

void GSZS(long double a,long double b,long double L,long double

B)

{

int N0; long double

e1,t,u2,L1,l,x,y;

long double

X;

long double N;

N0=(int)((L+3)/6.+0.5);

B=AngleToRadian(B); L=AngleToRadian(L); L1=AngleToRadian(N0*6-3);

l=L-L1;

e1=(sqrt(a*a-b*b))/b; //求椭圆的第二偏心率

t=tan(B); u2=e1*e1*cos(B)*cos(B);

X=HCX(a,b,B);

N=MYQQLBJN(a,b,B);

x=X+N/2.0*t*cos(B)*cos(B)*l*l+N/24.0*t*(5-t*t+9*u2+4*pow(u2,2))*pow(cos(B),4)*pow(l,4)+N/720.0*t*(61-58*t*t+pow(t,4))*pow(cos(B),6)*pow(l,6);

y=N*cos(B)*l+N/6.0*(1-t*t+u2)*pow(cos(B),3)*pow(l,3)+N/120.0*(5-18*t*t+pow(t,4)+14*u2-58*u2*t*t)*pow(cos(B),5)*pow(l,5);

printf("nt===============计算结果===============n");

printf("ntx=%-20.4lfy=%-20.4lfn",x,(y+500000+N0*1000000));

printf("nt=====================================n");

}

long double AngleToRadian(long double

alfa)

{

long double alfa1,alfa2;

long double HS;

alfa = alfa + 0.00000000000001;

alfa1=floor(alfa)+floor((alfa-floor(alfa))*100.)/60;

alfa2=(alfa*100.-floor(alfa*100.0))/36;

alfa1+=alfa2;

HS=alfa1/180.*PI;

return(HS);

}

long double MYQQLBJN(long double a,long double b,long double

B)

{

long double

e,n0,n2,n4,n6,n8,sin2B,sin4B,sin6B,sin8B,N;

e=(sqrt(a*a-b*b))/a;

n0=a;

n2=1.0/2.0*e*e*n0; n4=3.0/4.0*e*e*n2;

n6=5.0/6.0*e*e*n4; n8=7.0/8.0*e*e*n6;

sin2B=pow(sin(B),2); sin4B=pow(sin(B),4);

sin6B=pow(sin(B),6); sin8B=pow(sin(B),8);

N=n0+n2*sin2B+n4*sin4B+n6*sin6B+n8*sin8B;

return

N;

}

long double HCX(long double a,long double b,long double B)

{

long double e,m0,m2,m4,m6,m8,a0,a2,a4,a6,a8,X;

e=(sqrt(a*a-b*b))/a; //求椭圆的第一偏心率

m0=a*(1-e*e);

m2=3.0/2.0*e*e*m0; m4=5.0/4.0*e*e*m2;

m6=7.0/6.0*e*e*m4; m8=9.0/8.0*e*e*m6;

a0=m0+m2/2.0+3.0/8.0*m4+5.0/16.0*m6+35.0/128.0*m8;

a2=m2/2.0+m4/2.0+15.0/32.0*m6+7.0/16.0*m8;

a4=m4/8.0+3.0/16.0*m6+7.0/32.0*m8;

a6=m6/32.0+m8/16.0;

a8=m8/128.0;

X=a0*B-a2/2.0*sin(2*B)+a4/4.0*sin(4*B)-a6/6.0*sin(6*B)+a8/8.0*sin(8*B);

return X;