matlab高斯正反算程序6,高斯正反算MATLAB程序

以下是一个基于 MATLAB 的高斯正反算程序，其中包括高斯投影坐标的正算和反算。

高斯正算：

function [x,y]=gaussfor(l,b,L0,E0,a,f)
%Bessel椭球体参数a=6377397.155m,f=1/299.1528128
%中央经度L0,中央经差E0
%WGS-84椭球体参数a=6378137m,f=1/298.2572236
%调用程序
%[x,y]=gaussfor(l,b,L0,E0,a,f)
%输入：l，b，L0，E0，a，f
%输出：x，y
%  弧度制转换为度分秒制
d2mf=@(d)mfix(d*180/pi,0);%度转度分
mf2s=@(mf)(mf-fix(mf))*60;%度分转度秒
dms2mf=@(d,m,s)d+sign(d)*(abs(m)+abs(s)/60)/60;%度分秒转度分

% 常数
L0 = d2mf(L0);
B= d2mf(b);
L= d2mf(l);
E= E0;

% 椭球参数
e2 = f*(2-f);
ee2 = e2/(1-e2);
n = f/(2-f);
A = a/(1+n)*(1+n^2/4+n^4/64);
alpha = [1/2*n-2/3*n^2+5/16*n^3+41/180*n^4;-13/48*n^2+3/5*n^3+557/1440*n^4;61/240*n^3-103/140*n^4;49561/161280*n^4];

% 投影坐标
t = tan(B*pi/180);
eta2 = ee2*cos(B*pi/180)^2;
x0 = A*(1-e2)*((1+1/4*ee2+1/64*ee2^2)*B-rad(alpha,2)*sin(2*B)+1/8*rad(alpha,4)*sin(4*B)-1/256*rad(alpha,6)*sin(6*B));
x = x0 + A*t.*(rad(alpha,1)+rad(alpha,3)*cos(2*B)+rad(alpha,5)*cos(4*B)+rad(alpha,7)*cos(6*B));

y = A*(1-e2)*(rad(alpha,0)*L-rad(alpha,2)/2*sin(2*L)+rad(alpha,4)/4*sin(4*L)-rad(alpha,6)/6*sin(6*L)) + A*t.^2/2.*(cos(B*pi/180).^2.*rad(alpha,2)+1/3*rad(alpha,4)*cos(B*pi/180).^4+1/5*rad(alpha,6)*cos(B*pi/180).^6);
% 转换为大地坐标系下的坐标
x = x + E;
end

function rad = rad(alpha,i)
rad = 1;
for j=1:i
    rad = rad*alpha(j);
end
end

function x = mfix(x,n)
if x<0
    x = fix(x)-1/n;
else
    x = fix(x);
end
end

高斯反算：

function [l,b]=gaussinv(x,y,L0,E0,a,f)
%Bessel椭球体参数a=6377397.155m,f=1/299.1528128
%中央经度L0,中央经差E0
%WGS-84椭球体参数a=6378137m,f=1/298.2572236
%调用程序
%[l,b]=gaussinv(x,y,L0,E0,a,f)
%输入：x，y，L0，E0，a，f
%输出：l，b

% 弧度制转换为度分秒制
d2mf=@(d)mfix(d*180/pi,0);%度转度分
mf2s=@(mf)(mf-fix(mf))*60;%度分转度秒
dms2mf=@(d,m,s)d+sign(d)*(abs(m)+abs(s)/60)/60;%度分秒转度分

% 常数
L0 = d2mf(L0);
E= E0;

% 椭球参数
e2 = f*(2-f);
ee2 = e2/(1-e2);
n = f/(2-f);
A = a/(1+n)*(1+n^2/4+n^4/64);
alpha = [1/2*n-2/3*n^2+5/16*n^3+41/180*n^4;-13/48*n^2+3/5*n^3+557/1440*n^4;61/240*n^3-103/140*n^4;49561/161280*n^4];

% 投影坐标
y = y - E;
Bf = y/A;
Mf = A*((1-e2/4-3/64*e2^2-5/256*e2^3)*Bf-rad(alpha,2)/2*(sin(2*Bf)+1/2*rad(alpha,1)*sin(4*Bf)+1/4*rad(alpha,3)*sin(6*Bf)+1/6*rad(alpha,4)*sin(8*Bf)));
Nf = A./sqrt(1-e2*sin(Bf*pi/180)^2);
t = tan(Bf*pi/180);
eta2 = ee2*cos(Bf*pi/180)^2;
r = x/A;
D = r./Nf;
C = eta2*cos(Bf*pi/180)^2;
lmd = L0+dms2mf(D-D.^3/6*(1+t^2+C+14/15*t^2*(1+C)-9/40*C^2).*cos(Bf*pi/180));
B = Bf-dms2mf((D.^2/2-t^2/2-2*t^4/3-37/96*t^6)*sin(Bf*pi/180).*cos(Bf*pi/180).*rad(alpha,2)+D.^4/24*(5-4*t^2+14*C-28/3*t^2.*C+6*C.^2-3*t^2.*C.^2-9/5*t^4).*sin(Bf*pi/180).*cos(Bf*pi/180).*rad(alpha,4)-D.^6/720*(61+662/45*t^2+1320/45*t^4+720/45*t^6).*sin(Bf*pi/180).*cos(Bf*pi/180).*rad(alpha,6));
% 转换为度分秒制
l = dms2mf(lmd);
b = dms2mf(B);

end

function rad = rad(alpha,i)
rad = 1;
for j=1:i
    rad = rad*alpha(j);
end
end

function x = mfix(x,n)
if x<0
    x = fix(x)-1/n;
else
    x = fix(x);
end
end