mapxy

PURPOSE ^

MAPXY - compute latitude and longitude from x and y

SYNOPSIS ^

function [alat,alon]=mapxy(x,y,hem);

DESCRIPTION ^

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [alat,alon]=mapxy(x,y,hem);
0002 %MAPXY - compute latitude and longitude from x and y
0003 %
0004 %   hemisphere must be 1 for north and 0 for south
0005 %
0006 %   Usage:
0007 %      [latitude,longitude]=mapxy(x,y,hemisphere)
0008 %
0009 %   See also MAPLL, LL2XY
0010 
0011 re = 6378137;
0012 e2 = 0.00669437999015;
0013 
0014 e=sqrt(e2);
0015 
0016 %Standard parallel = latitude with no distortion
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

Generated on Sun 29-Mar-2009 20:22:55 by m2html © 2003