收藏
下载资源下载资源 加入VIP,免费下载

大地坐标与空间直角坐标的转换程序代码.docx

上传人：b****6 文档编号：5736072 上传时间：2022-12-31 格式：DOCX 页数：13 大小：17.83KB
下载 相关 举报
大地坐标与空间直角坐标的转换程序代码.docx_第1页
第1页 / 共13页
大地坐标与空间直角坐标的转换程序代码.docx_第2页
第2页 / 共13页
大地坐标与空间直角坐标的转换程序代码.docx_第3页
第3页 / 共13页
大地坐标与空间直角坐标的转换程序代码.docx_第4页
第4页 / 共13页
大地坐标与空间直角坐标的转换程序代码.docx_第5页
第5页 / 共13页
点击查看更多>>
下载资源
资源描述

大地坐标与空间直角坐标的转换程序代码.docx

《大地坐标与空间直角坐标的转换程序代码.docx》由会员分享，可在线阅读，更多相关《大地坐标与空间直角坐标的转换程序代码.docx（13页珍藏版）》请在冰豆网上搜索。

大地坐标与空间直角坐标的转换程序代码.docx

大地坐标与空间直角坐标的转换程序代码

#include"stdio.h"

#include"math.h"

#include"stdlib.h"

#include"iostream"

#definePI3.1415926535897323

doublea,b,c,e2,ep2;

intmain（）

{

intm,n,t;

doubleRAD（doubled,doublef,doublem）;

voidRBD（doublehd）;

voidBLH_XYZ（）;

voidXYZ_BLH（）;

voidB_ZS（）;

voidB_FS（）;

voidGUS_ZS（）;

voidGUS_FS（）;

printf（"大地测量学\n"）;

sp1:

printf（"请选择功能:

\n"）;

printf（"1.大地坐标系到大地空间直角坐标的转换\n"）;

printf（"2.大地空间直角坐标到大地坐标系的转换\n"）;

printf（"3.贝塞尔大地问题正算\n"）;

printf（"4.贝塞尔大地问题反算\n"）;

printf（"5.高斯投影正算\n"）;

printf（"6.高斯投影反算\n"）;

printf（"0.退出程序\n"）;

scanf（"%d",&m）;

if（m==0）exit（0）;

sp2:

printf（"请选择椭球参数（输入椭球序号）:

\n"）;

printf（"1.克拉索夫斯基椭球参数\n"）;

printf（"2.IUGG_1975椭球参数\n"）;

printf（"3.CGCS_2000椭球参数\n"）;

printf（"0.其他椭球参数（自行输入）\n"）;

scanf（"%d",&n）;

switch（n）

{

case1:

a=6378245.0;b=6356863.0188;c=6399698.9018;e2=0.00669342162297;ep2=0.00673852541468;break;

case2:

a=6378140.0;b=6356755.2882;c=6399596.6520;e2=0.00669438499959;ep2=0.00673950181947;break;

case3:

a=6378137.0;b=6356752.3141;c=6399593.6259;e2=0.00669438002290;ep2=0.00673949677547;break;

case0:

{

printf（"请输入椭球参数：

\n"）;

printf（"长半径a="）;scanf（"%lf",&a）;

printf（"短半径b="）;scanf（"%lf",&b）;

c=a*a/b;

ep2=（a*a-b*b）/（b*b）;

e2=（a*a-b*b）/（a*a）;

break;

}

default:

printf（"\n\n输入错误!

\n请重新输入!

\n\n"）;gotosp2;

}

while

（1）

{

switch（m）

{

case1:

BLH_XYZ（）;break;

case2:

XYZ_BLH（）;break;

case3:

B_ZS（）;break;

case4:

B_FS（）;break;

case5:

GUS_ZS（）;break;

case6:

GUS_FS（）;break;

default:

printf（"\n\n输入错误!

\n请重新输入!

\n\n"）;gotosp1;

}

printf（"是否继续进行此功能计算?

\n\n"）;

printf（"（若继续进行此功能计算，则输入1；\n若选择其他功能进行计算，则输入2；\n若退出，则输入0.）\n"）;

scanf（"%d",&t）;

switch（t）

{

case1:

break;

case2:

gotosp1;

case0:

exit（0）;

}

}

}

doubleRAD（doubled,doublef,doublem）

{

doublee;

doublesign=（d<0.0）?

-1.0:

1.0;

if（d==0）

{

sign=（f<0.0）?

-1.0:

1.0;

if（f==0）

{

sign=（m<0.0）?

-1.0:

1.0;

}

}

if（d<0）

d=d*（-1.0）;

if（f<0）

f=f*（-1.0）;

if（m<0）

m=m*（-1.0）;

e=sign*（d*3600+f*60+m）*PI/（3600*180）;

returne;

}

voidRBD（doublehd）

{

intt;

intd,f;

doublem;

doublesign=（hd<0.0）?

-1.0:

1.0;

if（hd<0）

hd=fabs（hd）;

hd=hd*3600*180/PI;

t=int（hd/3600）;

d=sign*t;

hd=hd-t*3600;

f=int（hd/60）;

m=hd-f*60;

printf（"%d'%d'%lf'\n",d,f,m）;

}

voidBLH_XYZ（）

{

doubleB,L,H,N,W;

doubled,f,m;

doubleX,Y,Z;

printf（"请输入大地坐标（输入格式为角度（例如：

30'40'50'））：

\n"）;

printf（"大地经度L="）;

scanf（"%lf'%lf'%lf'",&d,&f,&m）;

L=RAD（d,f,m）;

printf（"大地纬度B="）;

scanf（"%lf'%lf'%lf'",&d,&f,&m）;

B=RAD（d,f,m）;

printf（"大地高H="）;

scanf（"%lf",&H）;

W=sqrt（1-e2*sin（B）*sin（B））;

N=a/W;

X=（N+H）*cos（B）*cos（L）;

Y=（N+H）*cos（B）*sin（L）;

Z=（N*（1-e2）+H）*sin（B）;

printf（"\n\n转换后得到大地空间直角坐标为：

\n\n"）;

printf（"X=%lf\nY=%lf\nZ=%lf\n\n",X,Y,Z）;

}

voidXYZ_BLH（）

{

doubleB,L,H,N,W;

doubleX,Y,Z;

doubletgB0,tgB1;

printf（"请输入大地空间直角坐标：

\n"）;

printf（"X="）;

scanf（"%lf",&X）;

printf（"Y="）;

scanf（"%lf",&Y）;

printf（"Z="）;

scanf（"%lf",&Z）;

printf（"\n\n转换后得到大地坐标为：

\n\n"）;

L=atan（Y/X）;

printf（"大地经度为：

L="）;

RBD（L）;

printf（"\n"）;

tgB0=Z/sqrt（X*X+Y*Y）;

tgB1=（1/sqrt（X*X+Y*Y））*（Z+a*e2*tgB0/sqrt（1+tgB0*tgB0-e2*tgB0*tgB0））;

while（fabs（tgB0-tgB1）>5*pow（10,-10））

{

tgB0=tgB1;

tgB1=（1/sqrt（X*X+Y*Y））*（Z+a*e2*tgB0/sqrt（1+tgB0*tgB0-e2*tgB0*tgB0））;

}

B=atan（tgB1）;

printf（"大地纬度为：

B="）;

RBD（B）;

printf（"\n"）;

W=sqrt（1-e2*sin（B）*sin（B））;

N=a/W;

H=sqrt（X*X+Y*Y）/cos（B）-N;

printf（"大地高为：

H=%lf\n\n",H）;

}

voidB_ZS（）

{

doubleL1,B1,A1,s,d,f,mi;

doubleu1,u2,m,M,k2,alfa,bt,r,kp2,alfap,btp,rp;

doublesgm0,sgm1,lmd,lmd1,lmd2,A2,B2,l,L2;

printf（"请输入已知点的大地坐标（输入格式为角度（例如：

30'40'50'），下同）：

\nL1="）;

scanf（"%lf'%lf'%lf'",&d,&f,&mi）;

L1=RAD（d,f,mi）;

printf（"\nB1="）;

scanf（"%lf'%lf'%lf'",&d,&f,&mi）;

B1=RAD（d,f,mi）;

printf（"请输入大地方位角：

\nA1="）;

scanf（"%lf'%lf'%lf'",&d,&f,&mi）;

A1=RAD（d,f,mi）;

printf（"请输入该点至另一点的大地线长：

\ns="）;

scanf（"%lf",&s）;

u1=atan（sqrt（1-e2）*tan（B1））;

m=asin（cos（u1）*sin（A1））;

M=atan（tan（u1）/cos（A1））;

m=（m>0）?

m:

m+2*PI;

M=（M>0）?

M:

M+PI;

k2=ep2*cos（m）*cos（m）;

alfa=（1-k2/4+7*k2*k2/64-15*k2*k2*k2/256）/b;

bt=k2/4-k2*k2/8+37*k2*k2*k2/512;

r=k2*k2/128-k2*k2*k2/128;

sgm0=alfa*s;

sgm1=alfa*s+bt*sin（sgm0）*cos（2*M+sgm0）+r*sin（2*sgm0）*cos（4*M+2*sgm0）;

while（fabs（sgm0-sgm1）>2.8*PI/180*pow（10,-7））

{

sgm0=sgm1;

sgm1=alfa*s+bt*sin（sgm0）*cos（2*M+sgm0）+r*sin（2*sgm0）*cos（4*M+2*sgm0）;

}

sgm0=sgm1;

 

A2=atan（tan（m）/cos（M+sgm0））;

A2=（A2>0）?

A2:

A2+PI;

A2=（A1>PI）?

A2:

A2+PI;

u2=atan（-cos（A2）*tan（M+sgm0））;

lmd1=atan（sin（u1）*tan（A1））;

lmd1=（lmd1>0）?

lmd1:

lmd1+PI;

lmd1=（m

lmd1:

lmd1+PI;

lmd2=atan（sin（u2）*tan（A2））;

lmd2=（lmd2>0）?

lmd2:

lmd2+PI;

lmd2=（m

（（（M+sgm0）

lmd2:

lmd2+PI）:

（（（M+sgm0）>PI）?

lmd2:

lmd2+PI）;

lmd=lmd2-lmd1;

B2=atan（sqrt（1+ep2）*tan（u2））;

kp2=e2*cos（m）*cos（m）;

alfap=（e2/2+e2*e2/8+e2*e2*e2/16）-e2/16*（1+e2）*kp2+3*e2*kp2*kp2/128;

btp=e2*（1+e2）*kp2/16-e2*kp2*kp2/32;

rp=e2*kp2*kp2/256;

l=lmd-sin（m）*（alfap*sgm0+btp*sin（sgm0）*cos（2*M+sgm0）+rp*sin（2*sgm0）*cos（4*M+2*sgm0））;

L2=L1+l;

printf（"\n\n得到另一点的大地坐标和大地线在该点的大地方位角为：

\n\n"）;

printf（"L2="）;

RBD（L2）;printf（"\n"）;

printf（"B2="）;

RBD（B2）;printf（"\n"）;

printf（"A2="）;

RBD（A2）;

printf（"\n"）;

}

voidB_FS（）

{

doubleL1,B1,L2,B2,s,A1,A2,du,f,mi,m0,m,M;

doublel,u1,u2,alfa,bt,r,lmd0,dit_lmd,lmd,sgm,dit_sgm,sgm0,sgm1,alfap,btp,rp,k2,kp2;

printf（"请输入第一个点大地坐标（输入格式为角度（例如：

30'40'50'），下同）：

\n大地经度L1="）;

scanf（"%lf'%lf'%lf'",&du,&f,&mi）;

L1=RAD（du,f,mi）;

printf（"大地纬度B1="）;

scanf（"%lf'%lf'%lf'",&du,&f,&mi）;

B1=RAD（du,f,mi）;

printf（"\n请输入第二个点大地坐标：

\n大地经度：

L2="）;

scanf（"%lf'%lf'%lf'",&du,&f,&mi）;

L2=RAD（du,f,mi）;

printf（"大地纬度：

B2="）;

scanf（"%lf'%lf'%lf'",&du,&f,&mi）;

B2=RAD（du,f,mi）;

l=L2-L1;

u1=atan（sqrt（1-e2）*tan（B1））;

u2=atan（sqrt（1-e2）*tan（B2））;

sgm0=acos（sin（u1）*sin（u2）+cos（u1）*cos（u2）*cos（l））;

m0=asin（cos（u1）*cos（u2）*sin（l）/sin（sgm0））;

dit_lmd=0.003351831*sgm0*sin（m0）;

lmd0=l+dit_lmd;

dit_sgm=sin（m0）*dit_lmd;

sgm1=sgm0+dit_sgm;

m=asin（cos（u1）*cos（u2）*sin（lmd0）/sin（sgm1））;

A1=atan（sin（lmd0）/（cos（u1）*tan（u2）-sin（u1）*cos（lmd0）））;

A1=（A1>0）?

A1:

A1+PI;

A1=（m>0）?

A1:

A1+PI;

M=atan（sin（u1）*tan（A1）/sin（m））;

M=（M>0）?

M:

M+PI;

k2=ep2*cos（m）*cos（m）;

alfa=（1-k2/4+7*k2*k2/64-15*k2*k2*k2/256）/b;

bt=k2/4-k2*k2/8+37*k2*k2*k2/512;

r=k2*k2/128-k2*k2*k2/128;

kp2=e2*cos（m）*cos（m）;

alfap=（e2/2+e2*e2/8+e2*e2*e2/16）-e2/16*（1+e2）*kp2+3*e2*kp2*kp2/128;

btp=e2*（1+e2）*kp2/16-e2*kp2*kp2/32;

rp=e2*kp2*kp2/256;

sgm0=acos（sin（u1）*sin（u2）+cos（u1）*cos（u2）*cos（l））;

sgm1=acos（sin（u1）*sin（u2）+cos（u1）*cos（u2）*cos（l+sin（m）*（alfap*sgm0+btp*sin（sgm0）*cos（2*M+sgm0））））;

while（fabs（sgm0-sgm1）>1*PI/180*pow（10,-8））

{

sgm0=sgm1;

sgm1=acos（sin（u1）*sin（u2）+cos（u1）*cos（u2）*cos（l+sin（m）*（alfap*sgm0+btp*sin（sgm0）*cos（2*M+sgm0））））;

}

sgm=sgm1;

lmd=l+sin（m）*（alfap*sgm+btp*sin（sgm）*cos（2*M+sgm））;

s=（sgm-bt*sin（sgm）*cos（2*M+sgm）-r*sin（2*sgm）*cos（4*M+2*sgm））/alfa;

A1=atan（sin（lmd）/（cos（u1）*tan（u2）-sin（u1）*cos（lmd）））;

A1=（A1>0）?

A1:

A1+PI;

A1=（m>0）?

A1:

A1+PI;

A2=atan（sin（lmd）/（sin（u2）*cos（lmd）-tan（u1）*cos（u2）））;

A2=（A2>0）?

A2:

A2+PI;

A2=（m<0）?

A2:

A2+PI;

printf（"\n\n得到两点间大地线长S和大地正反方位角A1、A2如下：

\n\n"）;

printf（"s=%lf\n",s）;

printf（"A1="）;

RBD（A1）;printf（"\n"）;

printf（"A2="）;

RBD（A2）;printf（"\n"）;

}

voidGUS_ZS（）

{

doubleB,L,x3,x6,y3,y6,Y3,Y6,du,f,mi,X,N,n,t;

doubleAt,Bt,Ct,Dt,m3,m6,l3,l6,W,L03,L06;

intDH3,DH6;

printf（"请输入大地坐标（输入格式为角度（例如：

30'40'50'））：

\n大地经度L="）;

scanf（"%lf'%lf'%lf'",&du,&f,&mi）;

L=RAD（du,f,mi）;

printf（"\n大地纬度B="）;

scanf（"%lf'%lf'%lf'",&du,&f,&mi）;

B=RAD（du,f,mi）;

At=1+3*e2/4+45*e2*e2/64+175*e2*e2*e2/256;

Bt=3*e2/4+15*e2*e2/16+525*e2*e2*e2/512;

Ct=15*e2*e2/64+105*e2*e2*e2/256;

Dt=35*e2*e2*e2/512;

X=a*（1-e2）*（At*B-Bt*sin（2*B）/2+Ct*sin（4*B）/4-Dt*sin（6*B）/6）;

W=sqrt（1-e2*sin（B）*sin（B））;

N=a/W;

n=sqrt（ep2）*cos（B）;

t=tan（B）;

DH3=（L-（1.5*PI/180））/（3*PI/180）+1;

DH6=L/（6*PI/180）+1;

L03=DH3*（3*PI/180）;

L06=DH6*（6*PI/180）-（3*PI/180）;

l3=L-L03;

l6=L-L06;

m3=cos（B）*l3;

m6=cos（B）*l6;

x3=X+N*t*（m3*m3/2+（5-t*t+9*n*n+4*n*n*n*n）*m3*m3*m3*m3/24+（61-58*t*t+t*t*t*t）*m3*m3*m3*m3*m3*m3/720）;

x6=X+N*t*（m6*m6/2+（5-t*t+9*n*n+4*n*n*n*n）*m6*m6*m6*m6/24+（61-58*t*t+t*t*t*t）*m6*m6*m6*m6*m6*m6/720）;

y3=N*（m3+（1-t*t+n*n）*m3*m3*m3/6+（5-18*t*t+t*t*t*t+14*n*n-58*n*n*t*t）*m3*m3*m3*m3*m3/120）;

y6=N*（m6+（1-t*t+n*n）*m6*m6*m6/6+（5-18*t*t+t*t*t*t+14*n*n-58*n*n*t*t）*m6*m6*m6*m6*m6/120）;

Y3=DH3*1000000+500000+y3;

Y6=DH6*1000000+500000+y6;

printf（"\n\n得到的高斯平面坐标为:

\n\n"）;

printf（"对于3度带:

\n纵坐标x=%.3lf\n横坐标y=%.3lf（通用坐标Y=%.3lf）\n\n",x3,y3,Y3）;

printf（"对于6度带:

\n纵坐标x=%.3lf\n横坐标y=%.3lf（通用坐标Y=%.3lf）\n\n",x6,y6,Y6）;

}

voidGUS_FS（）

{

doublex,y,Y,B,B0,B1,Bf,Vf,tf,Nf,nf,L,At,Bt,Ct,Dt,L3,L6;

longDH;

printf（"请输入高斯平面坐标：

\n\n"）;

printf（"纵坐标X="）;

scanf（"%lf",&x）;printf（"\n"）;

printf（"自然坐标y="）;

scanf（"%lf",&y）;printf（"\n"）;

printf（"通用坐标Y="）;

scanf（"%lf",&Y）;printf（"\n"）;

At=1+3*e2/4+45*e2*e2/64+175*e2*e2*e2/256;

Bt=3*e2/4+15*e2*e2/16+525*e2*e2*e2/512;

Ct=15*e2*e2/64+105*e2*e2*e2/256;

Dt=35*e2*e2*e2/512;

B0=x/（a*（1-e2）*At）;

B1=（x-a*（1-e2）*（-Bt*sin（2*B0）/2+Ct*sin（4*B0）/4-Dt*sin（6*B0）/6））/（a*（1-e2）*At）;

while（fabs（B1-B0）>1*pow（10,-8））

{

B0=B1;

B1=（x-a*（1-e2）*（-Bt*sin（2*B0）/2+Ct*sin（4*B0）/4-Dt*sin（6*B0）/6））/（a*（1-e2）*At）;

}

Bf=B1;

nf=sqrt（ep2）*cos（Bf）;

tf=tan（Bf）;

Vf=sqrt（1+ep2*cos（Bf）*cos（Bf））;

Nf=c/Vf;

B=Bf-Vf*Vf*tf/2*（（y/Nf）*（y/Nf）-（5+3*tf*tf+nf*nf-9*nf*nf*tf*tf）*pow（（y/Nf）,4）/12+（61+90*tf*tf+45*tf*tf）*pow（（y/Nf）,6）/360）;

L=（（y/Nf）-（1+2*tf*tf+nf*nf）*（y/Nf）*（y/Nf）*（y/Nf）/6+（5+28*tf*tf+24*pow（tf,4）+6*nf*nf+8*nf*nf*tf*tf）*pow（（y/Nf）,5）/120）/cos（Bf）;

DH=Y/1000000;

L3=3*PI/180*double（DH）+L;

L6=6*PI/180*double（DH）-3*PI/180+L;

printf（"\n\n得到的大地坐标为:

\n\n"）;

printf（"

展开阅读全文
相关资源
猜你喜欢
相关搜索

当前位置：首页 > 经管营销

copyright@ 2008-2022 冰豆网网站版权所有

经营许可证编号:鄂ICP备2022015515号-1