高斯反算
输入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("nt===========================计算结果===========================n");
printf("nttB=%-20.10lfttL=%-20.10lfn",FSB,DH*6-3+FSL);
printf("nt=============================================================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));
}