matlab大地測(cè)量高斯投影正反算程序設(shè)計(jì)實(shí)驗(yàn)

1、 高斯投影正反算 高斯投影正反算 一、實(shí)驗(yàn)?zāi)康?一、實(shí)驗(yàn)?zāi)康?編寫(xiě)高斯投影正反算的程序，并對(duì)格式化文件數(shù)據(jù)進(jìn)行計(jì)算，驗(yàn)證程序。 二、實(shí)驗(yàn)內(nèi)容： 二、實(shí)驗(yàn)內(nèi)容： 1、高斯投影正算公式 高斯投影正算公式 已知大地坐標(biāo)（B,L）及中央子午線經(jīng)度 L0，計(jì)算高斯平面坐標(biāo)（x，y） ，公式如下： 6 2 2 2 4 2 54 4 2 2 3 2) 330 270 58 61 )( ( cos ) sin( 720) 4 9 5 )( ( cos

2、) sin( 24 ) cos( ) sin( 2l t t t B B Nl t B B N l B B N X x? ?? ?? ? ? ? ?? ? ? ? ? ?5 2 2 2 4 2 5 3 2 2 3 ) 58 14 18 5 )( ( cos 120 ) 1 )( ( cos 6 ) cos( l t t t B N l t B N l B N y ? ? ? ? ? ? ? ? ? ? ? ?其中： B 為緯度， 0 L

3、 L l ? ? ，單位為弧度，B ea N 2 2 sin 1?? ，為卯酉圈曲率半徑， B t tan ? , B e 2 2 ' 2 cos ? ? ， bb a e2 2' ? ? 為第二偏心率，a 為旋轉(zhuǎn)橢球長(zhǎng)半軸，b 為短半軸，X 為子午線弧長(zhǎng)。 根據(jù)上式創(chuàng)建以 GaussianMapDirect 命名的函數(shù)，函數(shù)輸入輸出格式為 [x,y,L0] = GaussianMapDirect(B, L, RefEl

4、lipsoid) 或者[x,y,L0] = GaussianMapDirect(B, L, a, f) RefEllipsoid 為橢球參數(shù) RefEllipsoid = [a, b, c, f, e2, e2_]; WGS84 橢球參數(shù)：長(zhǎng)半軸 a = 6378137 扁率 f = 1/298.257223563 b=a(1-f) c = a*a/b; e2 = (a*a-b*b)/(a*a); - 2 - %計(jì)算高斯平面坐標(biāo)

5、（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;

6、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); end latitude2meridian 函數(shù)如下： 函數(shù)如下： function x = latitud

7、e2meridian(B,a,e) %由緯度 B 求對(duì)應(yīng)的子午線弧長(zhǎng) x，計(jì)算公式 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