c白塞尔大地主题解算

#include

#include

double hudu(double,double,double); /*度分秒转换为弧度*/

double du(double); /*弧度转换为度*/

double fen(double); /*弧度转换为分*/

double miao(double); /*弧度转换为秒*/

#define PI 3.1415926

void main (void)

{

int k;

printf("请选择大地主题算法，若执行正算，请输入1；若执行反算，请输入2。\n");

scanf("%d",&k);

/*大地主题正算*/

if(k==1)

{

double ax,ay,az,bx,by,bz,cx,cy,cz,S,dz,ez,fz,B1,B2,L1,L2,A1,A2;

int dx,dy,ex,ey,fx,fy;

double e2,W1,sinu1,cosu1,sinA0,coto1,sin2o1,cos2o1,sin2o,cos2o,A,B,C,r,t,o0,o,g,sinu2,q;

/*输入度分秒数据*/

printf("请输入大地线起点纬度度分秒\n");

scanf("%lf%lf%lf",&ax,&ay,&az);

printf("请输入大地线起点经度度分秒\n");

scanf("%lf%lf%lf",&bx,&by,&bz);

printf("请输入大地方位角度分秒\n");

scanf("%lf%lf%lf",&cx,&cy,&cz);

printf("请输入大地线长度\n");

scanf("%lf",&S);

/*调用函数*/

B1=hudu(ax,ay,az);

L1=hudu(bx,by,bz);

A1=hudu(cx,cy,cz);

/*白塞尔大地主题解算*/

e2=0.006693421622966;

W1=sqrt(1-e2*sin(B1)*sin(B1));

sinu1=sin(B1)*(sqrt(1-e2))/W1;

cosu1=cos(B1)/W1;

sinA0=cosu1*sin(A1);

coto1=cosu1*cos(A1)/sinu1;

sin2o1=2*coto1/(coto1*coto1+1);

cos2o1=(coto1*coto1-1)/(coto1*coto1+1);

A=6356863.020+(10718.949-13.474*(1-sinA0*sinA0))*(1-sinA0*sinA0);

B=(5354.469-8.978*(1-sinA0*sinA0))*(1-sinA0*sinA0);

C=(2.238*(1-sinA0*sinA0))*(1-sinA0*sinA0)+0.006;

r=691.46768-(0.58143-0.00144*(1-sinA0*sinA0))*(1-sinA0*sinA0);

t=(0.2907-0.0010*(1-sinA0*sinA0))*(1-sinA0*sinA0);

o0=(S-(B+C*cos2o1)*sin2o1)/A;

sin2o=sin2o1*cos(2*o0)+cos2o1*sin(2*o0);

cos2o=cos2o1*cos(2*o0)-sin2o1*sin(2*o0);

o=o0+(B+5*C*cos2o)*sin2o/A;

g=(r*o+t*(sin2o-sin2o1))*sinA0;

/*求B2*/

sinu2=sinu1*cos(o)+cosu1*cos(A1)*sin(o);

B2=atan(sinu2/(sqrt(1-e2)*sqrt(1-sinu2*sinu2)));

/*求L2*/

q=atan(sin(A1)*sin(o)/(cosu1*cos(o)-sinu1*sin(o)*cos(A1)));

/*判断q*/

if(sin(A1)>0 && tan(q)>0)

q=fabs(q);

else if(sin(A1)>0 && tan(q)<0)

q=PI-fabs(q);

else if(sin(A1)<0 && tan(q)<0)

q=-fabs(q);

else

q=fabs(q)-PI;

L2=L1+q-g/3600/180*PI;

/*求A2*/

A2=atan(c

osu1*sin(A1)/(cosu1*cos(o)*cos(A1)-sinu1*sin(o)));

/*判断A2*/

if(sin(A1)<0 && tan(A2)>0)

A2=fabs(A2);

else if(sin(A1)<0 && tan(A2)<0)

A2=PI-fabs(A2);

else if(sin(A1)>0 && tan(A2)>0)

A2=PI+fabs(A2);

else

A2=2*PI-fabs(A2);

/*调用函数*/

dx=(int)(du(B2));

dy=(int)(fen(B2));

dz=miao(B2);

ex=(int)(du(L2));

ey=(int)(fen(L2));

ez=miao(L2);

fx=(int)(du(A2));

fy=(int)(fen(A2));

fz=miao(A2);

printf("大地线终点纬度度分秒分别为：\n%d\n%d\n%lf\n",dx,dy,dz);

printf("大地线终点经度度分秒分别为：\n%d\n%d\n%lf\n",ex,ey,ez);

printf("终点大地方位角度分秒分别为：\n%d\n%d\n%lf\n",fx,fy,fz);

}

/*大地主题反算*/

else if(k==2)

{

double fx,fy,fz,gx,gy,gz,hx,hy,hz,ix,iy,iz,jz,kz,B1,B2,L1,L2,S,A1,A2;

int jx,jy,kx,ky;

double e2,W1,W2,sinu1,sinu2,cosu1,cosu2,L,a1,a2,b1,b2,g,g2,g0=0.0,r,p,q,sino,coso,o,sinA0,x,t1,t2,A,B,C,y;

/*输入度分秒数据*/

printf("请输入大地线起点纬度度分秒\n");

scanf("%lf%lf%lf",&fx,&fy,&fz);

printf("请输入大地线起点经度度分秒\n");

scanf("%lf%lf%lf",&gx,&gy,&gz);

printf("请输入大地线终点纬度度分秒\n");

scanf("%lf%lf%lf",&hx,&hy,&hz);

printf("请输入大地线终点经度度分秒\n");

scanf("%lf%lf%lf",&ix,&iy,&iz);

/*调用函数*/

B1=hudu(fx,fy,fz);

L1=hudu(gx,gy,gz);

B2=hudu(hx,hy,hz);

L2=hudu(ix,iy,iz);

/*白塞尔大地主题解算*/

e2=0.006693421622966;

W1=sqrt(1-e2*sin(B1)*sin(B1));

W2=sqrt(1-e2*sin(B2)*sin(B2));

sinu1=sin(B1)*sqrt(1-e2)/W1;

sinu2=sin(B2)*sqrt(1-e2)/W2;

cosu1=cos(B1)/W1;

cosu2=cos(B2)/W2;

L=L2-L1;

a1=sinu1*sinu2;

a2=cosu1*cosu2;

b1=cosu1*sinu2;

b2=sinu1*cosu2;

/*逐次趋近法求解A1*/

g=0;

r=L;

while(1)

{

p=cosu2*sin(r);

q=b1-b2*cos(r);

A1=atan(p/q);

/*判断A1*/

if(p>0 && q>0)

A1=fabs(A1);

else if(p>0 && q<0)

A1=PI-fabs(A1);

else if(p<0 && q<0)

A1=PI+fabs(A1);

else

A1=2*PI-fabs(A1);

sino=p*sin(A1)+q*cos(

A1);

coso=a1+a2*cos(r);

o=atan(sino/coso);

/*判断o*/

if(coso>0)

o=fabs(o);

else

o=PI-fabs(o);

sinA0=cosu1*sin(A1);

x=2*a1-(1-sinA0*sinA0)*coso;

t1=(33523299-(28189-70*(1-sinA0*sinA0))*(1-sinA0*sinA0))*1e-10;

t2=(28189-94*(1-sinA0*sinA0))*1e-10;

g2=(t1*o-t2*x*sino)*sinA0;

if(fabs(g2-g0)<=1e-100)

break;

else

{

r=L+g2;

g0=g2;

}

}

/*求解S*/

A=6356863.020+(10708.949-13.474*(1-sinA0*sinA0))*(1-sinA0*sinA0);

B=10708.938-17.956*(1-sinA0*sinA0);

C=4.487;

y=((1-sinA0*sinA0)*(1-sinA0*sinA0)-2*x*x)*coso;

S=A*o+(B*x+C*y)*sino;

/*求解A2*/

A2=atan(cosu1*sin(r)/(b1*cos(r)-b2));

/*判断A2*/

if(p<0 && q<0)

A2=fabs(A2);

else if(p<0 && q>0)

A2=PI-fabs(A2);

else if(p>0 && q>0)

A2=PI+fabs(A2);

else

A2=2*PI-fabs(A2);

/*调用函数*/

jx=(int)(du(A1));

jy=(int)(fen(A1));

jz=miao(A1);

kx=(int)(du(A2));

ky=(int)(fen(A2));

kz=miao(A2);

printf("起点大地方位角度分秒分别为：\n%d\n%d\n%lf\n",jx,jy,jz);

printf("终点大地方位角度分秒分别为：\n%d\n%d\n%lf\n",kx,ky,kz);

printf("大地线长度为：%lf\n",S);

}

/*数据错误*/

else

printf("数据错误，请重新输入\n");

}

/*度分秒转换为弧度*/

double hudu(double x,double y,double z)

{

double A0;

A0=(x+y/60+z/3600)*PI/180;

return A0;

}

/*弧度转换为度*/++++++++++++++

double du(double B0)

{

double x0;

x0=(int)(B0*180/PI);

return x0;

}

/*弧度转换为分*/

double fen(double C0)

{

double _y,y0;

_y=(int)(C0*180/PI);

y0=(fabs)((int)((C0*180/PI-_y)*60));

return y0;

}

/*弧度转换为秒*/

double miao(double D0)

{

double _z1,_z2,z0;

_z1=(int)(D0*180/PI);

_z2=(int)((D0*180/PI-_z1)*60);

z0=(fabs)((double)(((D0*180/PI-_z1)*60-_z2)*60));

return z0;

}