高斯投影正反算c代码.docx
- 文档编号:9135859
- 上传时间:2023-02-03
- 格式:DOCX
- 页数:11
- 大小:16.78KB
高斯投影正反算c代码.docx
《高斯投影正反算c代码.docx》由会员分享,可在线阅读,更多相关《高斯投影正反算c代码.docx(11页珍藏版)》请在冰豆网上搜索。
高斯投影正反算c代码
文档编制序号:
[KK8UY-LL9IO69-TTO6M3-MTOL89-FTT688]
高斯投影正反算c代码
高斯投影正反算程序设计
一.程序设计流程
本程序的设计思路如下:
(1),程序采用VS08版本作为开发平台,并采用C#语言作为开发语言,设计为WindowsForm窗体程序形式。
(2),程序主要的算法来自于教材。
但是本程序为了更加实用,添加了更多的解算基准,包括:
WGS-84,国际椭球1975,克氏椭球,和2000国家大地坐标系。
(3),程序为了更方便的读取数据和输出数据,故需要自己定义了固定的数据输入格式和数据输出格式或形式,请老师注意查看。
二.代码
usingSystem;
usingSyst
usingSystem.ComponentModel;
usingSystem.Data;
usingSystem.Drawing;
usingSystem.Text;
namespaceGauss
{
publicpartialclassForm1:
Form
{
//大地坐标
//GeodeticCoordinate
publicstructCRDGEODETIC
{
publicdoubledLongitude;
publicdoubledLatitude;
publicdoubledHeight;
}
//笛卡尔坐标
//CartesianCoordinate
publicstructCRDCARTESIAN
{
publicdoublex;
publicdoubley;
publicdoublez;
}
publicForm1()
{
InitializeComponent();
}
privatevoidbutton1_Click(objectsender,EventArgse)
{
doubleee=0;
doublea=0;
stringtt;
try
{
}
catch
{
MessageBox.Show("GaussInverse:
Choosedatumerror!
");
return;
}
if(tt.CompareTo("克氏椭球")==0)
{
a=6378245.00;
}
if(tt.CompareTo("WGS-84")==0)
{
a=6378137.00;
}
if(tt.CompareTo("1975国际椭球")==0)
{
a=6378140.00;
ee=
}
if(tt.CompareTo("2000国家大地坐标系")==0)
{
a=6378137.0;
}
constdoublepai=3.1415926;
doubleb=Math.Sqrt(a*a*(1-ee*ee));
doublec=a*a/b;
doubleepp=Math.Sqrt((a*a-b*b)/b/b);
CRDGEODETICpcrdGeo;
CRDCARTESIANpcrdCar;
doublemidlong;
//求纬度
string[]temp;
double[]tempradius=newdouble[3];
for(inti=0;i<3;i++)
{
tempradius[i]=Convert.ToDouble(temp[i]);
}
pcrdGeo.dLatitude=tempradius[0]/180.0*pai+tempradius[1]/180.0/60.0*pai+tempradius[2]/180/60.0/60*pai;
//求经度
for(inti=0;i<3;i++)
{
tempradius[i]=Convert.ToDouble(temp[i]);
}
pcrdGeo.dLongitude=tempradius[0]/180.0*pai+tempradius[1]/180.0/60.0*pai+tempradius[2]/180/60.0/60*pai;
intdeglon=Convert.ToInt32(pcrdGeo.dLongitude*180/pai);
//求中央经度
intnum;//带号
midlong=0;//默认值,需要制定分带
try
{
}
catch
{
MessageBox.Show("Choose3/6error!
");
return;
}
if(tt.CompareTo("3度带")==0)
{
num=Convert.ToInt32(deglon/6+1);
midlong=(6*num-3)/180.0*pai;
}
if(tt.CompareTo("6度带")==0)
{
num=Convert.ToInt32((deglon+1.5)/3);
midlong=num*3*pai/180;
}
doublelp=pcrdGeo.dLongitude-midlong;
doubleN=c/Math.Sqrt(1+epp*epp*Math.Cos(pcrdGeo.dLatitude)*Math.Cos(pcrdGeo.dLatitude));
doubleM=c/Math.Pow(1+epp*epp*Math.Cos(pcrdGeo.dLatitude)*Math.Cos(pcrdGeo.dLatitude),1.5);
doubleita=epp*Math.Cos(pcrdGeo.dLatitude);
doublet=Math.Tan(pcrdGeo.dLatitude);
doubleNscnb=N*Math.Sin(pcrdGeo.dLatitude)*Math.Cos(pcrdGeo.dLatitude);
doubleNcosb=N*Math.Cos(pcrdGeo.dLatitude);
doublecosb=Math.Cos(pcrdGeo.dLatitude);
doubleX;
doublem0,m2,m4,m6,m8;
doublea0,a2,a4,a6,a8;
m0=a*(1-ee*ee);
m2=3.0/2.0*m0*ee*ee;
m4=5.0/4.0*ee*ee*m2;
m6=7.0/6.0*ee*ee*m4;
m8=9.0/8.0*ee*ee*m6;
a0=m0+m2/2.0+3.0/8.0*m4+5.0/16.0*m6+35.0/128.0*m8;
a2=m2/2+m4/2+15.0/32.0*m6+7.0/16.0*m8;
a4=m4/8.0+3.0/16.0*m6+7.0/32.0*m8;
a6=m6/32.0+m8/16.0;
a8=m8/128.0;
doubleB=pcrdGeo.dLatitude;
doublesb=Math.Sin(B);
doublecb=Math.Cos(B);
doubles2b=sb*cb*2;
doubles4b=s2b*(1-2*sb*sb)*2;
doubles6b=s2b*Math.Sqrt(1-s4b*s4b)+s4b*Math.Sqrt(1-s2b*s2b);
X=a0*B-a2/2.0*s2b+a4*s4b/4.0-a6/6.0*s6b;
pcrdCar.x=Nscnb*lp*lp/2.0+Nscnb*cosb*cosb*Math.Pow(lp,4)*(5-t*t+9*ita*ita+4*Math.Pow(ita,4))/24.0
+Nscnb*Math.Pow(cosb,4)*Math.Pow(lp,6)*(61-58*t*t+Math.Pow(t,4))/720.0+X;
pcrdCar.y=Ncosb*Math.Pow(lp,1)+Ncosb*cosb*cosb*(1-t*t+ita*ita)/6.0*Math.Pow(lp,3)+Ncosb*Math.Pow(lp,5)*Math.Pow(cosb,4)*(5-18*t*t
+Math.Pow(t,4)+14*ita*ita-58*ita*ita*t*t)/120.0;
if(pcrdCar.y<0)
pcrdCar.y+=500000;
richTextBox1.Text="Results:
\nX:
\t"+Convert.ToString(pcrdCar.x)+"\nY:
\t"+Convert.ToString(pcrdCar.y);
}
privatevoidbutton2_Click(objectsender,EventArgse)
{
doubleee=0;
doublea=0;
stringtt;
intnum;//带号
stringytext;//利用y值求带号和中央经线
try
{
}
catch
{
MessageBox.Show("GaussInverse:
Choosedatumerror!
");
return;
}
if(tt.CompareTo("克氏椭球")==0)
{
a=6378245.00;
}
if(tt.CompareTo("WGS-84")==0)
{
a=6378137.00;
}
if(tt.CompareTo("1975国际椭球")==0)
{
a=6378140.00;
}
if(tt.CompareTo("2000国家大地坐标系")==0)
{
a=6378137.0;
}
doubleb=Math.Sqrt(a*a*(1-ee*ee));
doublec=a*a/b;
doubleepp=Math.Sqrt((a*a-b*b)/b/b);
CRDGEODETICpcrdGeo;
CRDCARTESIANpcrdCar;
doublemidlong=0;
//求X,Y和带号
pcrdCar.x=Convert.ToDouble(textBox4.Text);
ytext=textBox5.Text;
stringtemp=ytext.Substring(0,2);
num=Convert.ToInt32(temp);
ytext=ytext.Remove(0,2);
pcrdCar.y=Convert.ToDouble(ytext)-500000;
try
{
}
catch
{
MessageBox.Show("Choose3/6error!
");
return;
}
if(tt.CompareTo("3度带")==0)
{
midlong=num*3*pai/180;
}
if(tt.CompareTo("6度带")==0)
{
midlong=(6*num-3)*pai/180;
}
b=Math.Sqrt(a*a*(1-ee*ee));
c=a*a/b;
epp=Math.Sqrt(a*a-b*b)/b;
doublem0,m2,m4,m6,m8;
doublea0,a2,a4,a6,a8;
m0=a*(1-ee*ee);
m2=3.0/2.0*m0*ee*ee;
m4=5.0/4.0*ee*ee*m2;
m6=7.0/6.0*ee*ee*m4;
m8=9.0/8.0*ee*ee*m6;
a0=m0+m2/2.0+3.0/8.0*m4+5.0/16.0*m6+35.0/128.0*m8;
a2=m2/2+m4/2+15.0/32.0*m6+7.0/16.0*m8;
a4=m4/8.0+3.0/16.0*m6+7.0/32.0*m8;
a6=m6/32.0+m8/16.0;
a8=m8/128.0;
doubleBf,B;
Bf=pcrdCar.x/a0;
B=0.0;
while(Math.Abs(Bf-B)>1E-10)
{
B=Bf;
doublesb=Math.Sin(B);
doublecb=Math.Cos(B);
doubles2b=sb*cb*2;
doubles4b=s2b*(1-2*sb*sb)*2;
doubles6b=s2b*Math.Sqrt(1-s4b*s4b)+s4b*Math.Sqrt(1-s2b*s2b);
Bf=(pcrdCar.x-(-a2/2.0*s2b+a4/4.0*s4b-a6/6.0*s6b))/a0;
}
doubleitaf,tf,Vf,Nf;
itaf=epp*Math.Cos(Bf);
tf=Math.Tan(Bf);
Vf=Math.Sqrt(1+epp*epp*Math.Cos(Bf)*Math.Cos(Bf));
Nf=c/Vf;
doubleynf=pcrdCar.y/Nf;
pcrdGeo.dLatitude=Bf-1.0/2.0*Vf*Vf*tf*(ynf*ynf-1.0/12.0*Math.Pow(ynf,4)*(5+3*tf*tf+itaf*itaf-9*Math.Pow(itaf*tf,2))+
1.0/360.0*(61+90*tf*tf+45*Math.Pow(tf,4))*Math.Pow(ynf,6));
pcrdGeo.dLongitude=(ynf/Math.Cos(Bf)-(1+2*tf*tf+itaf*itaf)*Math.Pow(ynf,3)/6.0/Math.Cos(Bf)+
(5+28*tf*tf+24*Math.Pow(tf,4)+6*itaf*itaf+8*Math.Pow(itaf*tf,2))*Math.Pow(ynf,5)/120.0/Math.Cos(Bf));
pcrdGeo.dLongitude=pcrdGeo.dLongitude+midlong;
//pcrdGeo.dLatitude=pcrdGeo.dLatitude;
richTextBox2.Text="Results:
\nLatitude:
"+Convert.ToString(pcrdGeo.dLatitude)+"\nLongtitude:
"+Convert.ToString(pcrdGeo.dLongitude);
}
privatevoidlabel13_Click(objectsender,EventArgse)
{
}
}
}
三.程序运行结果分析
通过选取书上的具体实例进行测试,本程序的精度大体满足要求,一般正算的精度在0.01米和0.001米之间,反算的精度在0.0001秒左右。
以下是程序运行的截图。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 投影 正反 代码