经纬度坐标系转换为UTM坐标系（matlab）

如题所说，直接上程序。验证自己做一下，结果应该是对的。诚不我欺。
（注意：程序名与函数名保持一致！！！）

function[shuchu]=lat_lon2utm(lat_shuru,lon_shuru)
%地理经纬度坐标转换为UTM坐标
size_shuzu=size(lat_shuru);
for i=1:size_shuzu(2)
%输入经纬度
% lat=29.819206;
% lon=116.133243;
lat=lat_shuru(i);
lon=lon_shuru(i);
%%     %%unit：km
a=6378.137; 
e=0.0818192;
k0=0.9996;
E0=500;
N0=0;
Zonenum=fix(lon/6)+31;
lamda0=(Zonenum-1)*6-180+3;%degree
lamda0=lamda0*pi/180;%radian
phi=lat*pi/180;
lamda=lon*pi/180;%radian

v=1/sqrt(1-e^2*(sin(phi)^2));
A=(lamda-lamda0)*cos(phi);
T=tan(phi)*tan(phi);
C=e^2*cos(phi)*cos(phi)/(1-e^2);
s=(1-e^2/4-3*e^4/64-5*e^6/256)*phi-(3*e^2/8+3*e^4/32+45*e^6/1024)*sin(2*phi)+(15*e^4/256+45*e^6/1024)*sin(4*phi)-35*e^6/3072*sin(6*phi);
UTME=E0+k0*a*v*(A+(1-T+C)*A^3/6+(5-18*T+T^2)*A^5/120);
UTMN=N0+k0*a*(s+v*tan(phi)*(A^2/2+(5-T+9*C+4*C^2)*A^4/24+(61-58*T+T^2)*A^6/720));
shuchu(1,i)=UTME;
shuchu(2,i)=UTMN;
end