matlab大地测量高斯投影正反算程序设计实验

高斯投影正反算一、实验目的编写高斯投影正反算的程序，并对格式化文件数据进行计算，验证程序。二、实验内容：1、高斯投影正算公式已知大地坐标（B,L）及中央子午线经度L0，计算高斯平面坐标（x，y），公式如下：6222425442232)3302705861)((cos)sin(720)495)((cos)sin(24)cos()sin(2ltttBBNltBBNlBBNXx52224253223)5814185)((cos120)1)((cos6)cos(ltttBNltBNlBNy其中：B为纬度，0LLl，单位为弧度，BeaN22sin1，为卯酉圈曲率半径，Bttan,Be22'2cos，bbae22'为第二偏心率，a为旋转椭球长半轴，b为短半轴，X为子午线弧长。根据上式创建以GaussianMapDirect命名的函数，函数输入输出格式为[x,y,L0]=GaussianMapDirect(B,L,RefEllipsoid)或者[x,y,L0]=GaussianMapDirect(B,L,a,f)RefEllipsoid为椭球参数RefEllipsoid=[a,b,c,f,e2,e2_];WGS84椭球参数：长半轴a=6378137扁率f=1/298.257223563b=a(1-f)c=a*a/b;e2=(a*a-b*b)/(a*a);-1-e2_=(a*a-b*b)/(b*b);GaussianMapDirect函数如下：function[x,y,L0]=GaussianMapDirect(B,L,X)%WGS84椭球参数：a=6378137;%长半轴f=1/298.257223563;%扁率b=a*(1-f);%短半轴e=(sqrt(a^2-b^2))/a;%第一偏心率e_=(sqrt(a^2-b^2))/b;%第二偏心率%当地中央子午线决定于当地的直角坐标系统，首先确定您的直角坐标系统是3度带还是6度带投影，然后再根据如下公式推算。Q=input('请选择投影带：\n');ifQ==6%6度带：M=round((L+3)./6);%带号M=round[(L+3)/6]，即对(L+3)/6的值四舍五入取整数，L为当地经度；L0=6.*M-3;%则中央子午线经度L0=6×M-3else%3度带：M=round(L./3);%带号M=round(L/3)，即对(L/3)的值四舍五入取整数，L为当地经度；L0=3.*M;%则中央子午线经度L0=3×Mendl=(L-L0).*pi/180;N=a./(sqrt(1-(e^2).*(sin(B).^2)));t=tan(B);p=e_.*cos(B);%p表示yita-2-%计算高斯平面坐标（x，y）x=X+(N.*(sin(B)).*(cos(B)).*(l.^2))./2+(N.*(sin(B)).*((cos(B)).^3).*(5-((t).^2)+9.*(p.^2)+4.*(p.^4)).*(l.^4))./24...+(N.*sin(B).*(cos(B)).^5.*(61-58.*(t.^2)+(t.^4)+270.*(p.^2)-330.*(p.^2).*(t.^2)).*(l.^6))./720;y=N.*cos(B).*l+(N.*(cos(B)).^3.*(1-t.^2+p.^2).*(l.^3))./6+(N.*((cos(B)).^5).*(5-18.*(t.^2)+...t.^4+14.*(p.^2)-58.*(p.^2).*(t.^2)).*(l.^5))./120;L0=degree3dms(L0);endlatitude2meridian函数如下：functionx=latitude2meridian(B,a,e)%由纬度B求对应的子午线弧长x，计算公式m0=a*(1-e^2);m2=(3*(e^2)*m0)/2;m4=(5*(e^2)*m2)/4;m6=(7*(e^2)*m4)/6;m8=(9*(e^2)*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;x=a0*B-(a2*sin(2*B))/2+(a4*sin(4*B))/4-(a6*sin(6*B))/6+(a8*sin(8*B))/8;end度分秒转度函数如下：functiondegree=dms2degree(jiaodu)%度分秒(dd.mmss)－度-3-degree=fix(jiaodu);mimute=fix((jiaodu-degree)*100);second=(jiaodu-degree-mimute/100)*10000;degree=degree+mimute/60+second/3600;end度转度分秒函数如下：functiondms=degree3dms(du)%度－度分秒(dd.mmss)degree=fix(du);min=fix((du-degree)*60);second=(((du-degree)*60-min)*60);dms=degree+min/100+second/10000;end度分秒转弧度dms2rad函数如下：functionrad=dms2rad(n)deg=fix(n);%度取所给数n的整数部分min_tem=(n-deg)*100;%去掉整数部分后小数点后移两位min=fix(min_tem);%分取整sec=(min_tem-min)*100;%秒是小数点再向后移两位的数字rad=(deg+min/60.00+sec/3600)*pi/180.0;end2、高斯反算公式已知高斯平面坐标（x，y）及指定中央子午线经度L0，计算大地坐标（B,L）：64254222232)459061(720)935(242yttNMtyttNMtyNMtBBffffffffffffffff522242532230)8624285(cos1201)21(cos61cos1ytttBNytBNyBNLLfffffffffffff-4-式中，ffBeaN22sin1，3222)sin1()1(ffBeeaM，ffBe22'2cos，ffBttan，fB为根据子午线弧长x反算的底点纬度。创建以GaussianMapInverse命名的函数，函数输入输出格式为[B,L]=GaussianMapInverse(x,y,RefEllipsoid)或者[B,L]=GaussianMapInverse(x,y,a,f)GaussianMapInverse函数如下：function[L,B]=GaussianMapInverse(x,y,L0)%WGS84椭球参数：a=6378137;%长半轴f=1/298.257223563;%扁率b=a*(1-f);%短半轴e=(sqrt(a^2-b^2))/a;%第一偏心率e_=(sqrt(a^2-b^2))/b;%第二偏心率%根据中央子午线弧长x反算底点纬度BfBf=meridian2latitude(x,a,e);Nf=a./sqrt(1-(e.^2).*(sin(Bf).^2));Mf=a.*(1-e.^2)./sqrt((1-(e.^2).*(sin(Bf).^2)).^3);pf=e_.*cos(Bf);tf=tan(Bf);%已知高斯平面坐标（x，y）及指定中央子午线经度L0，计算大地坐标（B,L）B=Bf-tf.*(y.^2)./(2.*Mf.*Nf)+tf.*(5+3.*(tf.^2)+pf.^2-9.*(tf.^2).*(pf.^2)).*(y.^4)./(24.*Mf.*(Nf.^3))...+tf.*(61+90.*(tf.^2)+45.*(tf.^4)).*(y.^6)./(720.*Mf.*(Nf.^5));L=L0+y./(Nf.*cos(Bf))-(1+2.*(tf.^2)+pf.^2).*y.^3./(6.*(Nf.^3).*cos(Bf))...+(5+28.*tf.^2+24.*tf.^4+6.*pf.^2+8.*(tf.^2).*pf.^2).*y.^5./(120.*Nf.^5.*cos(Bf));B=rad2dms(B);L=rad2dms(L);-5-end子午线弧长反算函数meridian2latitude如下：functionB=meridian2latitude(x,a,e)m0=a*(1-e^2);m2=(3*(e^2)*m0)/2;m4=(5*(e^2)*m2)/4;m6=(7*(e^2)*m4)/6;m8=(9*(e^2)*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;%%纬度B的计算B0=x./a0;%B的初始值while1F=-(a2.*sin(2.*B0))./2+(a4.*sin(4.*B0))./4-(a6.*sin(6.*B0))./6+(a8.*sin(8.*B0))./8;B=(x-F)./a0;ifabs(B0-B)10.^-6break;endabs(B0-B)B0=B;endend弧度转度分秒rad2dms函数如下：functiondms=rad2dms(azimuth)%弧度转度分秒-6-dgree_tem=azimuth*180/pi;dgree=fix(dgree_tem);min_tem=((dgree_tem-dgree)*60);min=fix(min_tem);second=((min_tem-min)*60);dms=dgree+min/100+second/10000;end度分秒转弧度dms2rad函数如下：functionrad=dms2rad(n)deg=fix(n);%度取所给数n的整数部分min_tem=(n-deg)*100;%去掉整数部分后小数点后移两位min=fix(min_tem);%分取整sec=(min_tem-min)*100;%秒是小数点再向后移两位的数字rad=(deg+min/60.00+sec/3600)*pi/180.0;end3、实例计算验证首先读取文件data1.txt中的数据，计算其在相应六度带高斯投影带内的高斯平面直角坐标，并存贮在文件data2.txt中，data1.txt格式为：经度（dd.mmss）纬度（dd.mmss）大地高data2.txt格式为：x(m)y(m)中央子午线经度（dd.mmss）test61data文件程序如下:clear;clc;[filename,pathname]=uigetfile('*.*');%文件查找窗口file=fullfile(pathname,filename);%合并路径文件名-7-A=importdata(file);%将生成的文件导入工作空间，变量名为A%RefEllipsoid为椭球参数%RefEllipsoid=[a,b,c,f,e2,e2_];a=6378137;%WGS84椭球参数：长半轴f=1/298.257223563;%扁率b=a*(1-f);%WGS84椭球参数：短半轴e=(sqrt(a^2-b^2))/a;X=latitude2meridian(dms2rad(A.data(:,2)),a,e);%X为子午线弧长，有纬度B算出[x,y,L0]=GaussianMapDirect(dms2rad(A.data(:,2)),dms2degree(A.data(:,1)),X);B=[x,