经纬度转高斯平面坐标

function [ xy ] = BL2XY（ BL ）
% 此函数用于进行高斯正算
%WGS84椭球
a=6378137;
f=1/298.257223563;
b=a-a*f;
%c=a^2/b;
% k =6;  %6度带
e=sqrt（a^2-b^2）/a;
ee=e/sqrt（1-e^2）;
% 计算中央子午线经度
% if k==6
%     L0=6*（fix（L/6）+1）-3;
% elseif k==3
%     L0=3*（fix（（L-1.5）/3）+1）;
% end
L0=6*（fix（L/6）+1）-3;
LP=（L-L0）/180*pi; %单位为角度，需转换为弧度
B=B/180*pi;       

% 计算子午圈弧长
ap=1+3/4*e^2+45/64*e^4+175/256*e^6+11025/16384*e^8+43659/65536*e^10;
bp=3/4*e^2+15/16*e^4+525/512*e^6+2205/2068*e^8+72765/65536*e^10;
cp=15