function [alat,alon]=mapxy(x,y,hem,varargin);
%MAPXY - compute latitude and longitude from x and y
%
%   hemisphere must be 'n' for north and 's' for south
%   One can specify the standard parallel and center meridien.
%   Default values are (70,45) for northern hemisphere and (71,0) for southerm hemisphere.
%
%   Usage:
%      [latitude,longitude]=mapxy(x,y,hemisphere)
%      [latitude,longitude]=mapxy(x,y,hemisphere,parallel,meridien)
%
%   Example:
%      [lat,long]=mapxy(x,y,'s')
%      [lat,long]=mapxy(x,y,'n',71,39)
%
%   See also MAPLL, LL2XY

%check number of arguments.
if ~((nargin==3) | (nargin==5)),
	mapxyerrorusage();
end

%hem: must be either 's' or 'n'
if ~ischar(hem),
	error('mapxy error message: hemisphere argument should be ''n'' or ''s''');
end

%check hem is either 'n' or 's': 
if ~(strcmpi(hem,'n') | strcmpi(hem,'s')),
	error('mapxy error message: hem should be ''s'' or ''n''');
end

%set sn: 
if strcmpi(hem,'s'),
	slat=71;
	sn=-1.0;
	xlam=0;
else
	slat=70;
	sn=1.0;
	xlam=45;
end

%set defaults for standard parallels and centre meridians.
if nargin==5,
	slat=varargin{1};
	xlam=varargin{2};
end

re = 6378137;
e2 = 0.00669437999015;

e=sqrt(e2);

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


end
	

function mapxyerrorusage()
	help mapxy
end