function [alat,alon]=mapxy(x,y,hem);
%MAPXY - compute latitude and longitude from x and y
%
%   hemisphere must be 1 for north and 0 for south
%
%   Usage:
%      [latitude,longitude]=mapxy(x,y,hemisphere)
%
%   See also MAPLL, LL2XY

re = 6378137;
e2 = 0.00669437999015;

e=sqrt(e2);

%Standard parallel = latitude with no distortion
slat=71;
sn=-1.0;
xlam=0;

if hem ==1,
   slat=70;
	sn=1;
	xlam=45;
end

slat=slat/180*pi;

rho=sqrt(x.^2+y.^2);
cm=cos(slat)./sqrt(1.0-e2*(sin(slat).^2));
t=tan((pi/4)-(slat/2.))./((1.0-e*sin(slat))./(1.0+e*sin(slat))).^(e/2.);

t=rho.*t./(re*cm);

chi=(pi/2.)-2.*atan(t);

alat=chi+((e2/2.)+(5.0*e2^2/24.)+(e2^3/12.))*sin(2*chi)+...
((7.0*e2^2/48.)+(29.*e2^3/240.))*sin(4.0*chi)+...
(7.0*e2^3/120.)*sin(6.0*chi);

alat=(sn*alat/pi*180);
xpr=sn*x;
ypr=sn*y;
alon=atan2(xpr,-ypr)/pi*180-sn*xlam;
alon=sn*alon;

indice=find(alon<0);
if length(indice)>0,
   alon(indice)=alon(indice)+360;
end

indice=find(alon>360);
if length(indice)>0,
   alon(indice)=alon(indice)-360;
end
