0001 function [alat,alon]=mapxy(x,y,hem);
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 re = 6378137;
0012 e2 = 0.00669437999015;
0013
0014 e=sqrt(e2);
0015
0016
0017 slat=71;
0018 sn=-1.0;
0019 xlam=0;
0020
0021 if hem ==1,
0022 slat=70;
0023 sn=1;
0024 xlam=45;
0025 end
0026
0027 slat=slat/180*pi;
0028
0029 rho=sqrt(x.^2+y.^2);
0030 cm=cos(slat)./sqrt(1.0-e2*(sin(slat).^2));
0031 t=tan((pi/4)-(slat/2.))./((1.0-e*sin(slat))./(1.0+e*sin(slat))).^(e/2.);
0032
0033 t=rho.*t./(re*cm);
0034
0035 chi=(pi/2.)-2.*atan(t);
0036
0037 alat=chi+((e2/2.)+(5.0*e2^2/24.)+(e2^3/12.))*sin(2*chi)+...
0038 ((7.0*e2^2/48.)+(29.*e2^3/240.))*sin(4.0*chi)+...
0039 (7.0*e2^3/120.)*sin(6.0*chi);
0040
0041 alat=(sn*alat/pi*180);
0042 xpr=sn*x;
0043 ypr=sn*y;
0044 alon=atan2(xpr,-ypr)/pi*180-sn*xlam;
0045 alon=sn*alon;
0046
0047 indice=find(alon<0);
0048 if length(indice)>0,
0049 alon(indice)=alon(indice)+360;
0050 end
0051
0052 indice=find(alon>360);
0053 if length(indice)>0,
0054 alon(indice)=alon(indice)-360;
0055 end