第1页 / 共13页
第2页 / 共13页
第3页 / 共13页
第4页 / 共13页
第5页 / 共13页
第6页 / 共13页
第7页 / 共13页
第8页 / 共13页
C++ 1.大地坐标与空间直角坐标的相互转化 （2）大地问题正反算 （3）高斯投影正反算.docx
1. （1）大地坐标与空间直角坐标的相互转化
（2）大地问题正反算
（3）高斯投影正反算
#include "stdio.h"
#include "math.h"
#include "stdlib.h"
#include "iostream"
#define PI 3.1415926535897323
double a,b,c,e2,ep2;
int main()
{
int m,n,t;
double RAD(double d,double f,double m);
void RBD(double hd);
void BLH_XYZ();
void XYZ_BLH();
void B_ZS();
void B_FS();
void GUS_ZS();
void GUS_FS();
printf("
sp1:printf("请选择功能:\n");
大地测量学
\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)
{
case
1:a=6378245.0;b=6356863.0188;c=6399698.9018;e2=0.00669342162297;ep2=0.0067385254146
8;break;
case
2:a=6378140.0;b=6356755.2882;c=6399596.6520;e2=0.00669438499959;ep2=0.0067395018194
7;break;
case
3:a=6378137.0;b=6356752.3141;c=6399593.6259;e2=0.00669438002290;ep2=0.0067394967754
7;break;
case 0:{
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");goto sp2 ;
}
while(1)
{
switch(m)
{
}
case 1:BLH_XYZ();break;
case 2:XYZ_BLH();break;
case 3:B_ZS();break;
case 4:B_FS();break;
case 5:GUS_ZS();break;
case 6:GUS_FS();break;
default:printf("\n\n 输入错误!\n 请重新输入!\n\n");goto sp1 ;
printf("是否继续进行此功能计算? \n\n");
printf("（ 若继续进行此功能计算，则输入 1；\n
）\n");
若退出， 则输入 0.
则输入 2；\n
若选择其他功能进行计算，
scanf("%d",&t);
switch(t)
{
case 1:break;
case 2:goto sp1;
case 0:exit(0);
}
}
}
double RAD(double d,double f,double m)
{
double e;
double sign=(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);
return e;
}
void RBD(double hd)
{
int t;
int d,f;
double m;
double sign=(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);
}
void BLH_XYZ()
{
double B,L,H,N,W;
double d,f,m;
double X,Y,Z;
printf(" 请输入大地坐标（输入格式为角度(例如：30'40'50')）：\n");
printf("
scanf("%lf'%lf'%lf'",&d,&f,&m);
L=RAD(d,f,m);
printf("
scanf("%lf'%lf'%lf'",&d,&f,&m);
B=RAD(d,f,m);
printf("
scanf("%lf",&H);
大地高 H=");
大地经度 L=");
大地纬度 B=");
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
printf("X=%lf\nY=%lf\nZ=%lf\n\n",X,Y,Z);
转换后得到大地空间直角坐标为：\n\n");
}
void XYZ_BLH()
{
double B,L,H,N,W;
double X,Y,Z;
double tgB0,tgB1;
请输入大地空间直角坐标：\n");
printf("
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("
RBD(L);
printf("\n");
大地经度为： L=");
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.0,-10))
{
tgB0=tgB1;
tgB1=(1/sqrt(X*X+Y*Y))*(Z+a*e2*tgB0/sqrt(1+tgB0*tgB0-e2*tgB0*tgB0));
大地纬度为：B=");
}
B=atan(tgB1);
printf("
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);
}
void B_ZS()
{
double L1,B1,A1,s,d,f,mi;
double u1,u2,m,M,k2,alfa,bt,r,kp2,alfap,btp,rp;
double sgm0,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.0,-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=(m0)?lmd2:lmd2+PI;
lmd2=(mPI)?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*sg
m0));
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");
}
void B_FS()
{
double L1,B1,L2,B2,s,A1,A2,du,f,mi,m0,m,M;
double l,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.0,-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");
}
void GUS_ZS()
{
double B,L,x3,x6,y3,y6,Y3,Y6,du,f,mi,X,N,n,t;
double At,Bt,Ct,Dt,m3,m6,l3,l6,W,L03,L06;
int DH3,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));
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
