MATLAB高斯投影正反算-白衣卿相

MATLAB高斯投影正反算-白衣卿相

正算.m文件代码*************************************************************** function GaoSiTouYingZS

%UNTITLED Summary of this function goes here

% Detailed explanation goes here

%******* 6度带投影正算******

DB1=[17,33,55.7339];

DL1=[119,15,52.1159];

L0 = 117;

B=DB1(1)+DB1(2)/60+DB1(3)/3600;%把B，l矩阵度分秒转换成度的小数形式

B=B/(180/3.141592654); %角度转弧度

L=DL1(1)+DL1(2)/60+DL1(3)/3600;

l=L-L0;

%克拉索夫斯基椭球参数

a=6378245.0; %长半轴

b=6356863.0187730473; %短半轴

c=6399698.901782711; %极点处曲率半径

ro=206264.806247096355; %角度转弧度常数

E1=(a^2-b^2)/a^2; %E2=e^2

E2=(a^2-b^2)/b^2; %E2=(e')^2

t=tan(B);

YeTa2=E2*cos(B)*cos(B);

v=sqrt(1+E2*cos(B)*cos(B));

m0=a*(1-E1);

m2=3/2*E1*m0;

m4=5/4*E1*m2;

m6=7/6*E1*m4;

m8=9/8*E1*m6;

a0=m0+m2/2+3/8*m4+5/16*m6+35/128*m8;

a2=m2/2+m4/2+15/32*m6+7/16*m8;

a4=m4/8+3/16*m6+7/32*m8;

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)+a8*sin(8*B)/8;

N=c/v;

x=X+N*sin(B)*cos(B)*(l*3600)^2/(2*ro^2)+N*sin(B)*cos(B)^3*(5-t^2+9*YeTa2+4*YeTa2^2)*(l*36 00)^4/(24*ro^4)+N*sin(B)*cos(B)^5*(61-58*t^2+t^4)*(l*3600)^6/(720*ro^6);

y=N*cos(B)*l*3600/ro+N*cos(B)^3*(1-t^2+YeTa2)*(l*3600)^3/(6*ro^3)+N*cos(B)^5*(5-18*t^2+t

^4+14*YeTa2-58*YeTa2*t^2)*(l*3600)^5/(120*ro^5);

end

**************************************************************************

反算.m文件代码******************************************************** function GaoSiTouYingFs

%UNTITLED3 Summary of this function goes here

% Detailed explanation goes here

x=1944359.607;

y=240455.4563;

L0 = 117;

%克拉索夫斯基椭球参数

a=6378245.0; %长半轴

b=6356863.0187730473; %短半轴

c=6399698.901782711; %极点处曲率半径

ro=206264.806247096355; %角度转弧度常数

E1=(a^2-b^2)/a^2; %E2=e^2

E2=(a^2-b^2)/b^2; %E2=(e')^2

e=0.0001; %计算底点纬度-->Bf

Bf1=2;

Bf2=3;

while Bf2-Bf1>e

Bf1=x/111134.8611;

Bf1=Bf1/(180/3.141592654); %角度转弧度

F_Bf2=-16036.4803*sin(2*Bf1)+16.8281*sin(4*Bf1)-0.022*sin(6*Bf1);

Bf2=(x-F_Bf2)/111134.8611;

end

t=tan(Bf2);

YeTa2=E2*cos(Bf2)*cos(Bf2);

v=sqrt(1+E2*cos(Bf2)^2);

N=c/v;

M=c/(v^3);

B=Bf2-t*y^2/(2*M*N)+t*(5+3*t^2+YeTa2-9*YeTa2*t^2)*y^4/(24*M*N^3)-t*(61+90*t^2+45*t^4) *y^6/(720*M*N^5);

l=y/(N*cos(Bf2))-(1+2*t^2+YeTa2)*y^3+(5+28*t^2+24*t^4+6*YeTa2+8*YeTa2*t^2)*y^5/(120*N^ 5*cos(Bf2));

end

*****************************************************************************