高斯投影正反算编程Word文件下载.docx
《高斯投影正反算编程Word文件下载.docx》由会员分享，可在线阅读，更多相关《高斯投影正反算编程Word文件下载.docx（13页珍藏版）》请在冰豆网上搜索。
三．编程的相关代码
（1）正算
#iｎｃluｄe"
stdio.ｈ"
#ｉncludｅ "
sｔdｌiｂ.h＂
#include＂math.ｈ＂
#includｅ"
asｓｅrt.h＂
#ｄeｆinepｉ（4*atan（１．0））
ｉntｉ;
sｔructjin
{
dｏuble Ｂ；
dｏuｂle L；
ｄｏuｂleL0；
};
ｓtrｕctjing[１00];
main（intarｇc,double *aｒｇｖ[]）
ＦILE*r=fopeｎ（＂a．ｔxt"
，"
r"
）;
aｓseｒt（r！
=ＮＵLL）；
ﻩＦILE *ｗ=fopen（"
b．ｔxt"
＂w"
aｓsert（r!
=ＮULL）；
ﻩｉnｔｉ=0；
ｗhilｅ（fscaｎf（r，"
%lf％lｆ％ｌf＂，&
ｇ[i］．Ｂ,＆g[i].Ｌ,&
g[ｉ].L0）!
=ＥOＦ）
ﻩ｛
double ａ,b；
intzｕoｂｉａo；
priｎtf（"
\ｎ请输入坐标系:
北京5４=1，西安80＝２，WＧS84=3:
"
ﻩｓｃaｎf（"
%d"
&
zuobiａo）；
getｃhaｒ（）；
if（zuobiao==1）
{
a＝6３78245；
b=6３５６863.018773047３;
ﻩ｝
ｉf（ｚｕobiaｏ==2）
ﻩ {
a=63７８14０;
b=6３5675５.２881５７5287;
｝
iｆ（zuoｂiaｏ==3）
ﻩ{
ﻩa＝63７8137;
ﻩ b=６３５６７52．3142;
} /／选择坐标系／/
doubleｆ=（ａ-ｂ）/ａ;
doｕｂlｅ e,ｅ2;
ﻩe＝sqｒt（2*ｆ－f*f）；
ｅ２=ｓqrt（（a/b）*（a/b）-1）;
/／求椭球的第一,第二曲率//
ﻩ
ｄouｂlem0,m2，ｍ4，m6,ｍ8;
doｕble a0,a2,a4,ａ６,a8;
m０＝ａ＊（1-ｅ*e）；
ｍ２=3*e*e*m0／２;
m4=５*ｅ*e*m2/４;
ﻩ m6=７*ｅ*ｅ＊m４/6;
m8=９*e*e＊m6／8；
ａ０＝m0+ｍ2/2＋3*m4/8+5*m6/16＋３5*m８/128;
a２=ｍ2/2+m4/２＋1５*m6／32+7*m8/16;
ﻩ a４=m4/８+３＊ｍ6/1６+７*m８／３2;
a6=m6／３2+m8／１６;
a８=ｍ8／1２8；
ﻩ
douｂle Bmiao,Lｍiａo, L０mｉao；
Bmiao＝（int）（g[i].B）＊３600.0+（int）（（ｇ[i].B-（int）（g［i]．B））*10０.0）*60.0+（g[i].Ｂ＊100-（inｔ）（g[i]．Ｂ*１00））*100．0;
Lｍiａo＝（int）（ｇ[ｉ].L）＊36０0.０+（iｎt）（（ｇ[i］.L－（int）（ｇ[i].L））＊10０.0）*60．0+（g[i]．Ｌ*100－（int）（ｇ［i].L*100））*100.0；
ﻩＬ0mｉao=（ｉnt）（g[i]．L０）*360０.0+（int）（（g[ｉ］.L０-（int）（g［i].L0））*１00.０）*60.0+（g[i]．Ｌ0*１00-（ｉｎt）（g［ｉ].L０*１０0））＊100．0;
dｏuｂledb;
ﻩ db=pi/180．0/3600．0；
ﻩdｏublｅ B1,Ｌ1,l；
B1=Bmiao*dｂ;
L1= Lｍiaｏ*db;
l＝L1－Ｌ0miao*ｄb;
//角度转化为弧度//
douｂｌe T＝ｔａn（B１）*ｔan（B1）;
ｄoublｅn＝e２*e２*cos（B１）*coｓ（Ｂ1）;
douｂleA=l*cｏs（B1）;
doubleX,x,ｙ;
X=a0＊（Ｂ1）－a2*sin（2*B1）／２＋a４*sｉn（4＊B１）/4-ａ6*sin（6*B1）／６＋ａ8＊sin（８＊Ｂ1）/８；
//求弧长／／
ｄoubleN=a/sｑrt（１-e＊ｅ*sｉｎ（Ｂ1）*ｓｉn（B1））;
inｔZonｅwide；
ﻩintZｏｎeｎumｂer;
printf（＂\ｎ请输入带宽：
3度带或６度带Zｏｎewide＝"
）；
sｃanf（"
＆Ｚoneｗidｅ）;
ｇeｔchar（）;
ﻩiｆ（Zonewｉde==3）
{
ﻩZoｎenumber=（ｉnt）（（g[i]．Ｌ-Zoｎｅwiｄe/2）/Zoｎewｉde+1）;
else if（Zｏnewiｄe==６）
ﻩ Zonｅnumｂer=（ｉnt）g［i].L／Zonewide+１;
}
eｌｓe
ﻩ｛
ﻩ pｒｉntf（"
错误＂）;
exiｔ（０）；
ﻩ}／/选择带宽//
doｕble
ﻩ FE＝Zoｎenumbｅｒ*10000００+5000００;
//改写为国家通用坐标//
ﻩｙ＝FE+N＊Ａ+A*A*Ａ*Ｎ＊（1－T*T+n＊n）/6+A＊Ａ*A*A*A*Ｎ＊（5-18＊T*T+T*T＊Ｔ*Ｔ+1４*n＊ｎ-58*ｎ*n＊T*T）/120;
ﻩx＝Ｘ＋tan（B1）*Ｎ＊A＊A/２+tan（B1）*N*A*A＊A*A＊（5-T＊T+９*n*n+4*n*n*ｎ＊ｎ）/2４+ｔan（B1）＊N*Ａ*Ａ*A＊A*A*A*（6１-58＊T*T+T*T*T*T）/72０;
ﻩｐriｎtf（＂\n所选坐标系的转换结果:
x＝％lfy=%lf\n＂,x,ｙ）;
ﻩｆprｉntｆ（ｗ,＂%lｆ%ｌｆ＼n"
x,ｙ）；
//输出结果到文本文件/／
ﻩﻩ
ﻩ}
ﻩfclｏse（r）；
fclose（ｗ）；
ﻩsystｅｍ（"
pause"
ﻩrｅtｕｒn 0;
}
（２）反算
#include"
stｄiｏ．h"
#iｎｃlｕde"
stdlib.h"
＃ iｎclude"
ｍath.h"
# ｉｎcluｄe"
asｓeｒt.h"
#dｅfｉnｅpｉ（4*atan（1．0））
doｕbｌeＸ,Y,B1，B2,Ｂ３，F，t;
dｏublｅ m０,m２,m4,m6，m8;
doｕblea０，a2,a4，a6,a8,ａ1,b1;
doubleBＢ,ＬL，Ｂf;
doｕｂle e，e1；
inｔ ｄ,m,s,ｉ,ｚｕobiao；
doubleｓorｔ（ｄoubｌe，double）;
stｒuｃt jiｎ
dｏuble x;
doubley;
dｏuｂleＬ０；
struct jｉn g［100］;
／/ｘ，ｙ，L0为输入量：
x,y坐标和中央子午线经度//
main（intａrgc, ｄoubｌe *argv[］）
ＦILE*r＝fｏpen（＂c.tｘt"
"
aｓsｅrt（r!
=ＮＵＬL）;
FＩＬE*w＝ｆopｅn（"
ｄ．txt"
assｅｒｔ（ｒ!
＝NＵLL）;
inti＝０;
ｗｈile（fｓcａnf（ｒ,"
%lｆ％lf％lｆ"
，&
g[i］．x,&
g[i]．y，&
g［i]．L0）!
=EOＦ）／/文件为空,无法打开//
ﻩdoublea1=6378245.００0000000０；
//克拉索夫斯基椭球参数//
dｏublｅｂ1=６35６8６3.01877３０4７3;
doｕｂｌea75=637８140.０0０0000000;
//１９75国际椭球参数／/
ﻩ doｕbｌｅb75=6356７５５．２8815７５2８7;
ｄoubｌｅａ8４＝６37８13７.0000000000;
/／ＷGS－84系椭球参数//
ﻩdouｂleb8４＝6３５6７5２.３1420000００;
doｕbｌe M,N;
//mouyou圈曲率半径，子午圈曲率半径//
ﻩdoｕblet,n;
ｄoｕｂlｅ A,B,C;
doubｌeＢB，LＬ,Bf,ＬL0,BB0;
ﻩ
dｏublea，b;
ﻩpriｎtf（"
\ｎ选择参考椭球：
1=克拉索夫斯基椭球,2＝1９75国际椭球,3＝WGS-8４系椭球:
scaｎｆ（＂%ｄ"
zuｏbiaｏ）;
ｇeｔchaｒ（）;
ﻩiｆ（zｕobiao==1）
ﻩa=a1；
b=b1;
ﻩ}
if（ｚuoｂiao==2）
a＝a75;
ｂ=ｂ75;
if（zｕobｉaｏ==3）
{
ﻩﻩa=a84;
ﻩﻩb=b8４;
ﻩ }//选择参考椭球,求解第一偏心率e,第二偏心率e1//
Bf=sort（a,b）;
/／调用求解底点纬度的函数／/
dｏｕbleq＝sqｒｔ（1-e*e*ｓiｎ（Ｂf）*ｓin（Bf））;
doubｌeＧ＝coｓ（Bｆ）；
ﻩ M=a*（1-e*e）/（q*q＊ｑ）;
N=a/q;
douｂlｅH,Ｉ;
A=g[i].ｙ/N;
ﻩ H=A*A＊Ａ;
ﻩI＝A*A＊A*A*A;
t=taｎ（Bf）；
ﻩn=e1*cｏs（Bf）;
B=t*t;
ﻩＣ＝n*n;
ＢＢ0＝Bf-ｇ[ｉ］.y＊t＊A／（2*M）+ｇ[i].y*ｔ＊Ｈ/（24*M）＊（5+3＊B＋C-９*B＊C）-g[i].y*t*I/（72０*Ｍ）*（6１+90*B＋45*B*B）;
ﻩLL0=g[i]．L0＊pi／１80．0＋A／G－H/（6*Ｇ）*（1.0＋2＊Ｂ+C）+Ｉ／（１20*G）*（5.0+28＊B+24*B*B+6*C+8*B*Ｃ）;
//利用公式求解经纬度//
ﻩｉntBdu,Bfeｎ,Ldu,Lｆeｎ;
ﻩdoｕblｅBmiao,Lmiaｏ;
Ｌdu=ｉnt（LＬ０/pi*18０）;
ﻩLｆen=int（（LL0/pi*180）*6０-Ｌdu*6０）;
ﻩＬmiao＝ＬL0/pi*1８0＊3600-Ldｕ*3600-Lｆen*60；
ﻩBdu＝int（BB0/pｉ＊180）;
ﻩBfen=int（（BＢ０/ｐｉ＊１80）*60-Bdu*6０）；
ﻩBmiａo=ＢB０／pi*180*3６0０-Bdu*３600-Bfeｎ
升级会员 