用c语言坐标正反算,高斯正反算程序代码——C语言

最新推荐文章于 2021-05-18 16:33:11 发布

weixin_39862985

最新推荐文章于 2021-05-18 16:33:11 发布

阅读量1k

点赞数

文章标签：用c语言坐标正反算

#include

#define PI 3.141592653589793

double DMS2RAD(double dmsAngle)

{

int degAngle,minAngle,intSignOfDms;

double radAngle,secAngle;

intSignOfDms = 1;

if(dmsAngle<0) intSignOfDms=-1;

dmsAngle=fabs(dmsAngle);

degAngle=(int)(dmsAngle+0.0001);

minAngle=(int)((dmsAngle-degAngle)*100+0.0001);

secAngle=(dmsAngle-degAngle-minAngle/100.0)*10000.0;

radAngle=(degAngle+minAngle/60.0+secAngle/3600.0)*PI/180.0;

radAngle=radAngle*intSignOfDms;

return radAngle;

}

double RAD2DMS(double radAngle)

{

int degAngle,minAngle,intSignOfRad;

double secAngle,dmsAngle;

intSignOfRad = 1;

if(radAngle<0) intSignOfRad=-1;

radAngle=fabs(radAngle);

secAngle=radAngle*180.0/PI*3600.0;

degAngle=(int)(secAngle/3600+0.0001);

minAngle=(int)((secAngle-degAngle*3600.0)/60.0+0.0001);

secAngle=secAngle-degAngle*3600.0-minAngle*60.0;

dmsAngle=degAngle+minAngle/100.0+secAngle/10000.0;

dmsAngle=dmsAngle*intSignOfRad;

return dmsAngle;

}

void a0a2a4a6a8(double a,double e,double *Coeficient_a0)

{

double

m0,m2,m4,m6,m8;

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

m2=3*e*e*m0/2;

m4=5*e*e*m2/4;

m6=7*e*e*m4/6;

m8=9*e*e*m6/8;

Coeficient_a0[0]=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;

Coeficient_a0[1]=m2/2+m4/2+15*m6/32+7*m8/16;

Coeficient_a0[2]=m4/8+3*m6/16+7*m8/32;

Coeficient_a0[3]=m6/32+m8/16;

Coeficient_a0[4]=m8/128;

}

void main()

{ int h,k;

double a,Alfa;

double dmslat,dmslon,dmsl0;

double a0 ,a2 ,a4, a6, a8;

double

radlat,radlon,radl0,l;

double b,t,sb,cb,ita,e1,e;

double X,l0;

double N,c,v;

double coor_x,coor_y;

double Bf0,Bf1,dB,FBf,bf;

double itaf, tf;

double Nf,Mf;

double B,L,dietaB,dietal;

double sinBf,cosBf;

double *Coeficient_a0;

Coeficient_a0=(double

*)malloc(5*sizeof(double));

printf("正算请选择1，反算请选择2：\n");

scanf("%d",&k);

if(k==1) //正算

{ printf("请选择坐标系:\n");

printf("

选择WGS-84坐标系,请按1\n");

printf("

选择BJ-54坐标系,请按2\n");

printf("

选择GDZ-80坐标系,请按3\n");

printf("

其他坐标系,请按4\n");

scanf("%d",&h);

if(h==1)

a=6378137,Alfa=1.0/298.257223563;

if(h==2)

a=6378245,Alfa=1.0/298.3;

if(h==3)

a=6378140,Alfa=1.0/298.257;

if(h==4)

{

printf("输入椭球长轴：");

scanf("%lf",&a);

printf("输入椭球扁率分母：");

scanf("%lf",&Alfa);

Alfa=1.0/Alfa;

}

printf("请输入已知点纬度:\n");

scanf("%lf",&dmslat);

printf("请输入已知点经度:\n");

scanf("%lf",&dmslon);

printf("请输入中央子午线经度

:\n");

scanf("%lf",&dmsl0);

radlat=DMS2RAD(dmslat);

radlon=DMS2RAD(dmslon);

radl0=DMS2RAD(dmsl0);

l=radlon-radl0;

b=a*(1-Alfa);

sb=sin(radlat);

cb=cos(radlat);

t=sb/cb;

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

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

ita=e1*cb;

a0a2a4a6a8(a,e,Coeficient_a0);

a0=Coeficient_a0[0];

a2=Coeficient_a0[1];

a4=Coeficient_a0[2];

a6=Coeficient_a0[3];

a8=Coeficient_a0[4];

X=a0*radlat-sb*cb*((a2-a4+a6)+(2*a4-16*a6/3)*sb*sb+16*a6*sb*sb*sb*sb/3.0);

c=a*a/b;

v=sqrt(1+e1*e1*cb*cb);

N=c/v;

coor_x=X+N*sb*cb*l*l/2+N*sb*cb*cb*cb*(5-t*t+9*ita*ita+4*ita*ita*ita*ita)*l*l*l*l/24+N*sb*cb*cb*cb*cb*cb*(61-58*t*t+t*t*t*t)*l*l*l*l*l*l/720;

coor_y=N*cb*l+N*cb*cb*cb*(1-t*t+ita*ita)*l*l*l/6+N*cb*cb*cb*cb*cb*(5-18*t*t+t*t*t*t+14*ita*ita-58*ita*ita*t*t)*l*l*l*l*l/120;

printf("coor_x=%.4lf\n",coor_x);

printf("coor_y=%.4lf\n",coor_y);

}

if(k==2)//反算

{

printf("请选择坐标系:\n");

printf(" 选择WGS-84坐标系,请按1\n");

printf(" 选择BJ-54坐标系,请按2\n");

printf(" 选择GDZ-80坐标系,请按3\n");

printf(" 其他坐标系,请按4\n");

scanf("%d",&h);

if(h==1) a=6378137,Alfa=1.0/298.257223563;

if(h==2) a=6378245,Alfa=1.0/298.3;

if(h==3) a=6378140,Alfa=1.0/298.257;

if(h==4)

{

scanf("输入椭球长轴：%lf",&a);

scanf("输入椭球扁率分母：%lf",&Alfa);

Alfa=1.0/Alfa;

}

printf("请输入l0:\n");

scanf("%lf",&l0);

l0=l0*PI/180.0;

printf("请输入x坐标:\n");

scanf("%lf",&coor_x);

printf("请输入y坐标:\n");

scanf("%lf",&coor_y);

b=a*(1-Alfa);

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

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

a0a2a4a6a8(a,e,Coeficient_a0);

a0=Coeficient_a0[0];

a2=Coeficient_a0[1];

a4=Coeficient_a0[2];

a6=Coeficient_a0[3];

a8=Coeficient_a0[4];

X=coor_x;

Bf0=X/a0;

{

sinBf=sin(Bf0);cosBf=cos(Bf0);

FBf=-sinBf*cosBf*((a2-a4+a6)+(2.0*a4-16.0*a6/3.0)*sinBf*sinBf+

(16.0/3.0)*a6*(sinBf*sinBf*sinBf*sinBf));

Bf1=(X-FBf)/a0;

dB=Bf1-Bf0;

Bf0=Bf1;

}

while(fabs(dB*180.0/PI*3600.0)>0.00001);

bf=Bf1;

c=a*a/b;

v=sqrt(1+e1*e1*cos(bf)*cos(bf));

Nf=c/v;

Mf=c/(v*v*v);

tf=sin(bf)/cos(bf);

itaf=e1*cos(bf);

dietaB=tf*coor_y*coor_y/(2*Mf*Nf)-tf*(5+3*tf*tf+itaf*itaf-9*tf*tf*itaf*itaf)*coor_y*coor_y*coor_y*coor_y/(24*Mf*Nf*Nf*Nf)+(61+90*tf*tf+45*tf*tf*tf*tf)*coor_y*coor_y*coor_y*coor_y*coor_y*coor_y/(720*Mf*Nf*Nf*Nf*Nf*Nf);

dietal=coor_y/(Nf*cos(bf)+(1+2*tf*tf+itaf*itaf)*cos(bf)*coor_y*coor_y/(6*Nf))+(5+44*tf*tf+32*tf*tf*tf*tf-2*itaf*itaf-16*itaf*itaf*tf*tf)/(360*Nf*Nf*Nf*Mf*Mf*cos(bf));

B=bf-dietaB;

L=l0+dietal;

dmslat=RAD2DMS(B);

dmslon=RAD2DMS(L);

printf("已知点的纬度(ddd.mmss)为：.9lf\n",dmslat);

printf("已知点的经度(ddd.mmss)为：.9lf\n",dmslon);

}