#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#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;
do
{
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);
}
}
加载中，请稍候......
加载中…
喜欢