function [x,y]=gaosizs(L,B,L0)
% 功能：将大地坐标解算到高斯平面上。L为大地经度，B为大地纬度，L0为中央子午线。
% B0=24.271152150;L0=110.174395360;
% [x,y]=gaosizs1(L0,B0,111)
% 椭球选择：1 克拉索夫斯基椭球 2 1975国际椭球 3 WGS84椭球 4 CGCS2000椭球
% 请输入所用椭球编号：3
% x =
% 2.705667588948869e+06
% y =
% 4.285722813950625e+05
disp('椭球选择：1 克拉索夫斯基椭球 2 1975国际椭球 3 WGS84椭球 4 CGCS2000椭球');
T=input('请输入所用椭球编号：');
switch(T)
case 1
a=6378245;f=1/298.3;
case 2
a=6378140;f=1/298.257;
case 3
a=6378137;f=0.003352810664;
case 4
a =6378137;f = 1/298.257222101;
otherwise
disp('输入有误，请输入1到4之间的数！');
end
B2=DHH(B);
l=DHH(L)-DHH(L0);
X=HC(B,a,f);
b=a*(1-f); %半短轴
e1=sqrt(a^2-b^2)/a; %第一偏心率
e2=sqrt(a^2-b^2)/b; %第二偏心率
V=sqrt(1+e2^2*(cos(B2))^2);
W=sqrt(1-e1^2*(sin(B2))^2);
n=sqrt(e2^2*(cos(B2))^2);
t=tan(B2);
c=a^2/b;
M=c/V^3; %子午圈曲率半径
N=c/V; %卯酉圈曲率半径
x1=X;
x2=N/2*t*(cos(B2))^2*l^2;
x3=N/24*t*(5-t^2+9*n^2+4*n^4)*(cos(B2))^4*l^4;
x4=N/720*t*(61-58*t^2+t^4)*(cos(B2))^6*l^6;
x=x1+x2+x3+x4;
y1=N*cos(B2)*l;
y2=N/6*(1-t^2+n^2)*(cos(B2))^3*l^3;
y3=N/120*(5-18*t^2+t^4+14*n^2-58*n^2*t^2)*(cos(B2))^5*l^5;
y=y1+y2+y3;
y=y+500000;
end %高斯正算，将国家大地坐标解算为高斯平面坐标
function X=HC(B,a,f) %计算子午线的弧长，B为维度（弧度），a为半长轴，f为扁率
b=a*(1-f) ;%半短轴
e2=sqrt(a^2-b^2)/b; %第一偏心率
c=a^2/b;
B1=DHH(B);
a0=c*(1-3/4*e2^2+45/64*e2^4-175/256*e2^6+11025/16384*e2^8);
a1=c*(-3/8*e2^2+15/32*e2^4-525/1024*e2^6+2205/4096*e2^8);
a2=c*(15/256*e2^4-105/1024*e2^6+2205/16384*e2^8);
a3=c*(-35/3072*e2^6+105/4096*e2^8);
a4=c*(315/131072*e2^8);
X=a0*B1+a1*sin(2*B1)+a2*sin(4*B1)+a3*sin(6*B1)+a4*sin(8*B1);
end
function HD=DHH(B)
%DHH是将度化算为弧度的函数
%度的表达形式形如（三十度四十五分二十三秒：30.4523）
D=fix(B);
F=fix((B-D)*100);
M=((B-D)*100-F)*100;
HD=((D+F/60+M/3600)/180)*pi;
end %弧度转换，将度分秒换算成相对应的弧度值
评论0