c语言坐标反算,C语言版——高斯坐标反算程序

高斯反算

输入1975年国际椭球体下的6度带坐标x和y，例

如(2506658.3680，19711330.6495)求其大地经纬度L，

B。(现用函数处理，且用指针类型的数据作函数参数)

#include

#include

#define PI 3.1415926535897932

#define P 206264.806247096355

void main()

{

long double RadianToAngle(long double

alfa);

long double Bf(long double a,long double b,long double

x);

void GSFS(long double a,long double b,long double Bf,long double

y,long double *B,long double *l);

long double a,b,x,y,L0;

long double FSB,FSL;

long double

*pointer_B,*pointer_L;

long double

DH;

pointer_B=&FSB; pointer_L=&FSL;

a=6378140.0000000000;b=6356755.2881575287;

printf("请输入x坐标和y坐标:(输入时x坐标和y坐标之间用逗号隔开)\n");

scanf("%lf,%lf",&x,&y);

FSB=Bf(a,b,x);

DH=floor(y/1000000); y=y-DH*1000000-500000;

GSFS(a,b,FSB,y,pointer_B,pointer_L);

FSB=RadianToAngle(FSB);

FSL=RadianToAngle(FSL);

L0=DH*6-3;

printf("\n\t===========================计算结果===========================\n");

printf("\n\t\tB=%-20.10lf\t\tL=%-20.10lf\n",FSB,DH*6-3+FSL);

printf("\n\t=============================================================\n");

}

long double RadianToAngle(long double alfa)

{

long double alfa1,alfa2;

alfa=alfa*180./PI;

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

alfa2=(alfa*60.-floor(alfa*60.))*0.006;

alfa1+=alfa2;

return(alfa1);

}

long double Bf(long double a,long double b,long double x)

{

long double

e,m0,m2,m4,m6,m8;

long double a0,a2,a4,a6,a8;

long double B,FB,Bn1,Bn;

int n=0;

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;

B=x/a0;

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

do

{

Bn=B;

B=(x+FB)/a0;

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

Bn1=B;

}

while(fabs((Bn1-Bn)*P)>=0.000000000001);

return(B);

}

void GSFS(long double a,long double b,long double Bf,long double

y,long double *B,long double *l)

{

long double

e,e1,V,t,u2;

long double

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

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

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

V=sqrt(1+e1*e1*cos(Bf)*cos(Bf));

t=tan(Bf);

u2=e1*e1*cos(Bf)*cos(Bf);

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(Bf),2); sin4B=pow(sin(Bf),4);

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

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

*B=Bf-1.0/2.0*V*V*t*((y/N)*(y/N)-1.0/12.0*(5+3*t*t+u2-9*u2*t*t)*pow((y/N),4)+1.0/360.0*(61+90*t*t+45*pow(t,4))*pow((y/N),6));

*l=1.0/cos(Bf)*((y/N)-1.0/6.0*(1+2*t*t+u2)*pow((y/N),3)+1.0/120.0*(5+28*t*t+24*pow(t,4)+6*u2+8*u2*t*t)*pow(y/N,5));

}