function [x,y]=gaosizs(Data)
%高斯投影坐标正算，将大地坐标（B,L）转换成高斯坐标（x,y）;
%Data为原始数据矩阵，第二列为L(经度)，第三列为B（纬度）；
disp('您选择的是高斯正算');
m=size(Data,1);
B=Data(:,3);
L=Data(:,2);
X=[];V=[];I=[];
L0=input('输入所用中央子午线L0=');
L0=DHH(L0);
for i=1:m
B(i,1)=DHH(B(i,1));
L(i,1)=DHH(L(i,1));
end
disp('1:克拉索夫斯基椭球 2:1975年国际椭球 3：WGS-84椭球');
T=0;
while(T<1||T>3)
T=input('请根据上列选择椭球类型T=');
switch T
case 1
a=6378245.0000000000;
b=6356863.0187730473;
for i=1:m
X(i,1)=111134.861*(B(i,1)*180/pi)-16036.480*sin(2*B(i,1))+16.828*sin(4*B(i,1))-0.022*sin(6*B(i,1));
end
case 2
a=6378140.0000000000;
b=6356755.2881575287;
for i=1:m
X(i,1)=111133.005*(B(i,1)*180/pi)-16038.528*sin(2*B(i,1))+16.833*sin(4*B(i,1))-0.022*sin(6*B(i,1));
end
case 3
a=6378137;
b=6356752.3142;
e=sqrt(1-b^2/a^2);
C0=a*(1-e^2/4+3*e^4/64-5*e^6/256);
C1=a*(3*e^2/8+3*e^4/32+45*e^6/1024);
C2=a*(15*e^4/256+45*e^6/1024);
C3=35*e^6/3072;
for i=1:m
X(i,1)=C0*B(i,1)-C1*sin(2*B(i,1))+C2*sin(4*B(i,1))-C3*sin(6*B(i,1));
end
otherwise
disp('T值无效 （1-3）');
end
end
e=(sqrt(a^2-b^2))/a;
e1=(sqrt(a^2-b^2))/b;
for i=1:m
V(i,1)=sqrt(1+(e1^2)*(cos(B(i,1)))^2);
end
c=(a^2)/b;
M=c./(V.^3);
N=c./V;
t=tan(B);
n=sqrt((e1^2).*(cos(B)).^2);
for i=1:m
I(i,1)=L(i,1)-L0;
end
xp1=X;
xp2=(N.*sin(B).*cos(B).*I.^2)./2;
xp3=(N.*sin(B).*((cos(B)).^3).*(5-t.^2+9*n.^2+4*n.^4).*I.^4)/24;
xp4=(N.*sin(B).*((cos(B)).^5).*(61-58*t.^2+t.^4).*I.^6)/720;
x=xp1+xp2+xp3+xp4;
yp1=N.*cos(B).*I;
yp2=N.*(cos(B)).^3.*(1-t.^2+n.^2).*I.^3/6;
yp3=N.*(cos(B)).^5.*(5-18*t.^2+t.^4+14*n.^2-58*(n.^2).*(t.^2)).*I.^5/120;
y=yp1+yp2+yp3;
format long g
plot(x(1:25),y(1:25),'bo',x(26:27),y(26:27),'rp','markersize',4);
end
%% 子程序DHH
function HD=DHH(BD)
%将度转换为弧度
D=fix(BD);
F=fix((BD-D)*100);
M=((BD-D)*100-F)*100;
HD=((D+F/60+M/3600)/180)*pi;
end
评论0