高斯投影正反算编程一高斯投影正反算基本公式文档格式.docx
《高斯投影正反算编程一高斯投影正反算基本公式文档格式.docx》由会员分享,可在线阅读,更多相关《高斯投影正反算编程一高斯投影正反算基本公式文档格式.docx(11页珍藏版)》请在冰豆网上搜索。
3度带
选择椭球参数
2)反算
输入自然值坐标x,y
利用迭代算法
求解底点纬度
利用公式计算B和L
打印B和L
三.编程的相关代码
(1)正算
#include"
stdio.h"
stdlib.h"
math.h"
assert.h"
#definepi(4*atan(1.0))
inti;
structjin
{
doubleB;
doubleL;
doubleL0;
};
structjing[100];
main(intargc,double*argv[])
FILE*r=fopen("
a.txt"
"
r"
);
assert(r!
=NULL);
FILE*w=fopen("
b.txt"
w"
inti=0;
while(fscanf(r,"
%lf%lf%lf"
&
g[i].B,&
g[i].L,&
g[i].L0)!
=EOF)
{
doublea,b;
intzuobiao;
printf("
\n请输入坐标系:
北京54=1,西安80=2,WGS84=3:
"
scanf("
%d"
zuobiao);
getchar();
if(zuobiao==1)
a=6378245;
b=6356863.0187730473;
}
if(zuobiao==2)
a=6378140;
b=6356755.2881575287;
if(zuobiao==3)
a=6378137;
b=6356752.3142;
}//选择坐标系//
doublef=(a-b)/a;
doublee,e2;
e=sqrt(2*f-f*f);
e2=sqrt((a/b)*(a/b)-1);
//求椭球的第一,第二曲率//
doublem0,m2,m4,m6,m8;
doublea0,a2,a4,a6,a8;
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;
a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;
a2=m2/2+m4/2+15*m6/32+7*m8/16;
a4=m4/8+3*m6/16+7*m8/32;
a6=m6/32+m8/16;
a8=m8/128;
doubleBmiao,Lmiao,L0miao;
Bmiao=(int)(g[i].B)*3600.0+(int)((g[i].B-(int)(g[i].B))*100.0)*60.0+(g[i].B*100-(int)(g[i].B*100))*100.0;
Lmiao=(int)(g[i].L)*3600.0+(int)((g[i].L-(int)(g[i].L))*100.0)*60.0+(g[i].L*100-(int)(g[i].L*100))*100.0;
L0miao=(int)(g[i].L0)*3600.0+(int)((g[i].L0-(int)(g[i].L0))*100.0)*60.0+(g[i].L0*100-(int)(g[i].L0*100))*100.0;
doubledb;
db=pi/180.0/3600.0;
doubleB1,L1,l;
B1=Bmiao*db;
L1=Lmiao*db;
l=L1-L0miao*db;
//角度转化为弧度//
doubleT=tan(B1)*tan(B1);
doublen=e2*e2*cos(B1)*cos(B1);
doubleA=l*cos(B1);
doubleX,x,y;
X=a0*(B1)-a2*sin(2*B1)/2+a4*sin(4*B1)/4-a6*sin(6*B1)/6+a8*sin(8*B1)/8;
//求弧长//
doubleN=a/sqrt(1-e*e*sin(B1)*sin(B1));
intZonewide;
intZonenumber;
\n请输入带宽:
3度带或6度带Zonewide="
Zonewide);
if(Zonewide==3)
Zonenumber=(int)((g[i].L-Zonewide/2)/Zonewide+1);
elseif(Zonewide==6)
Zonenumber=(int)g[i].L/Zonewide+1;
else
错误"
exit(0);
}//选择带宽//
double
FE=Zonenumber*1000000+500000;
//改写为国家通用坐标//
y=FE+N*A+A*A*A*N*(1-T*T+n*n)/6+A*A*A*A*A*N*(5-18*T*T+T*T*T*T+14*n*n-58*n*n*T*T)/120;
x=X+tan(B1)*N*A*A/2+tan(B1)*N*A*A*A*A*(5-T*T+9*n*n+4*n*n*n*n)/24+tan(B1)*N*A*A*A*A*A*A*(61-58*T*T+T*T*T*T)/720;
\n所选坐标系的转换结果:
x=%lfy=%lf\n"
x,y);
fprintf(w,"
%lf%lf\n"
//输出结果到文本文件//
fclose(r);
fclose(w);
system("
pause"
return0;
}
(2)反算
doubleX,Y,B1,B2,B3,F,t;
doublem0,m2,m4,m6,m8;
doublea0,a2,a4,a6,a8,a1,b1;
doubleBB,LL,Bf;
doublee,e1;
intd,m,s,i,zuobiao;
doublesort(double,double);
doublex;
doubley;
doubleL0;
//x,y,L0为输入量:
x,y坐标和中央子午线经度//
c.txt"
d.txt"
g[i].x,&
g[i].y,&
=EOF)//文件为空,无法打开//
doublea1=6378245.0000000000;
//克拉索夫斯基椭球参数//
doubleb1=6356863.0187730473;
doublea75=6378140.0000000000;
//1975国际椭球参数//
doubleb75=6356755.2881575287;
doublea84=6378137.0000000000;
//WGS-84系椭球参数//
doubleb84=6356752.3142000000;
doubleM,N;
//mouyou圈曲率半径,子午圈曲率半径//
doublet,n;
doubleA,B,C;
doubleBB,LL,Bf,LL0,BB0;
\n选择参考椭球:
1=克拉索夫斯基椭球,2=1975国际椭球,3=WGS-84系椭球:
a=a1;
b=b1;
a=a75;
b=b75;
a=a84;
b=b84;
}//选择参考椭球,求解第一偏心率e,第二偏心率e1//
Bf=sort(a,b);
//调用求解底点纬度的函数//
doubleq=sqrt(1-e*e*sin(Bf)*sin(Bf));
doubleG=cos(Bf);
M=a*(1-e*e)/(q*q*q);
N=a/q;
doubleH,I;
A=g[i].y/N;
H=A*A*A;
I=A*A*A*A*A;
t=tan(Bf);
n=e1*cos(Bf);
B=t*t;
C=n*n;
BB0=Bf-g[i].y*t*A/(2*M)+g[i].y*t*H/(24*M)*(5+3*B+C-9*B*C)-g[i].y*t*I/(720*M)*(61+90*B+45*B*B);
LL0=g[i].L0*pi/180.0+A/G-H/(6*G)*(1.0+2*B+C)+I/(120*G)*(5.0+28*B+24*B*B+6*C+8*B*C);
//利用公式求解经纬度//
intBdu,Bfen,Ldu,Lfen;
doubleBmiao,Lmiao;
Ldu=int(LL0/pi*180);
Lfen=int((LL0/pi*180)*60-Ldu*60);
Lmiao=LL0/pi*180*3600-Ldu*3600-Lfen*60;
Bdu=int(BB0/pi*180);
Bfen=int((BB0/pi*180)*60-Bdu*60);
Bmiao=BB0/pi*180*3600-Bdu*3600-Bfen*60;
//将弧度转化为角度//
%d度%d分%lf秒%d度%d分%lf秒\n"
Bdu,Bfen,Bmiao,Ldu,Lfen,Lmiao);
fprintf(w,"
%d°
%d’%lf”%d°
%d’%lf”\n"
//将结果输出到文本文件//
doublesort(doublea,doubleb)
e=sqrt(1-(b/a)*(b/a));
e1=sqrt((a/b)*(a/b)-1);
B1=g[i].x/a0;
do
{
F=-a2*sin(2*B1)/2+a4*sin(4*B1)/4-a6*sin(6*B1)/6+a8*sin(8*B1)/8;
B2=(g[i].x-F)/a0;
B3=B1;
B1=B2;
}while(fabs(B3-B2)>
10e-10);
//利用迭代算法求解底点纬度//
returnB2;