高斯投影正反算编程一高斯投影正反算基本公式.docx
《高斯投影正反算编程一高斯投影正反算基本公式.docx》由会员分享，可在线阅读，更多相关《高斯投影正反算编程一高斯投影正反算基本公式.docx（11页珍藏版）》请在冰豆网上搜索。
高斯投影正反算编程一高斯投影正反算基本公式
高斯投影正反算编程一．高斯投影正反算基本公式
（1）高斯正算基本公式
（2）高斯反算基本公式
以上主要通过大地测量学基础课程得到，这不进行详细的推导，只是列出基本公式指导编程的进行。
二．编程的基本方法和流程图
（1）编程的基本方法
高斯投影正反算基本上运用了所有的编程基本语句，本文中是利用C++语言进行基本的设计。
高斯正算中对椭球参数和带宽的选择主要运用了选择语句。
而高斯反算中除了选择语句的应用，在利用迭代算法求底点纬度还应用了循环语句。
编程中还应特别注意相关的度分秒和弧度之间的相互转换，这是极其重要的。
（2）相关流程图
1）正算
输入大地坐标B，L
和经差L0
选择带宽3/6度带
计算带号
计算弧长
计算平面坐标x,y
打印x,y
计算带号
计算弧长
计算平面坐标x,y
打印x,y
开始
6度带
3度带
选择椭球参数
2）反算
选择椭球参数
开始
输入自然值坐标x,y
和经差L0
利用迭代算法
求解底点纬度
利用公式计算B和L
打印B和L
三．编程的相关代码
（1）正算
#include"stdio.h"
#include"stdlib.h"
#include"math.h"
#include"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"）;
assert（r!
=NULL）;
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;
printf（"\n请输入带宽：
3度带或6度带Zonewide="）;
scanf（"%d",&Zonewide）;
getchar（）;
if（Zonewide==3）
{
Zonenumber=（int）（（g[i].L-Zonewide/2）/Zonewide+1）;
}
elseif（Zonewide==6）
{
Zonenumber=（int）g[i].L/Zonewide+1;
}
else
{
printf（"错误"）;
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;
printf（"\n所选坐标系的转换结果：
x=%lfy=%lf\n",x,y）;
fprintf（w,"%lf%lf\n",x,y）;//输出结果到文本文件//
}
fclose（r）;
fclose（w）;
system（"pause"）;
return0;
}
（2）反算
#include"stdio.h"
#include"stdlib.h"
#include"math.h"
#include"assert.h"
#definepi（4*atan（1.0））
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）;
structjin
{
doublex;
doubley;
doubleL0;
};
structjing[100];//x,y,L0为输入量：
x,y坐标和中央子午线经度//
main（intargc,double*argv[]）
{
FILE*r=fopen（"c.txt","r"）;
assert（r!
=NULL）;
FILE*w=fopen（"d.txt","w"）;
assert（r!
=NULL）;
inti=0;
while（fscanf（r,"%lf%lf%lf",&g[i].x,&g[i].y,&g[i].L0）!
=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;
doublea,b;
printf（"\n选择参考椭球：
1=克拉索夫斯基椭球，2=1975国际椭球，3=WGS-84系椭球:
"）;
scanf（"%d",&zuobiao）;
getchar（）;
if（zuobiao==1）
{
a=a1;
b=b1;
}
if（zuobiao==2）
{
a=a75;
b=b75;
}
if（zuobiao==3）
{
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;
//将弧度转化为角度//
printf（"\n所选坐标系的转换结果：
%d度%d分%lf秒%d度%d分%lf秒\n",Bdu,Bfen,Bmiao,Ldu,Lfen,Lmiao）;
fprintf（w,"%d°%d’%lf”%d°%d’%lf”\n",Bdu,Bfen,Bmiao,Ldu,Lfen,Lmiao）;//将结果输出到文本文件//
}
fclose（r）;
fclose（w）;
system（"pause"）;
return0;
}
doublesort（doublea,doubleb）
{
doublee,e1;
e=sqrt（1-（b/a）*（b/a））;
e1=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;
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;
}
升级会员 