source: issm/trunk/src/m/utils/LatLong/mapxy.m@ 2734

Last change on this file since 2734 was 2734, checked in by seroussi, 15 years ago

Minor, improved headers for latlong/xy conversions

File size: 1.8 KB
Line 
1function [alat,alon]=mapxy(x,y,hem,varargin);
2%MAPXY - compute latitude and longitude from x and y
3%
4% hemisphere must be 'n' for north and 's' for south
5% One can specify the standard parallel and center meridien.
6% Default values are (70,45) for northern hemisphere and (71,0) for southerm hemisphere.
7%
8% Usage:
9% [latitude,longitude]=mapxy(x,y,hemisphere)
10% [latitude,longitude]=mapxy(x,y,hemisphere,parallel,meridien)
11%
12% Example:
13% [lat,long]=mapxy(x,y,'s')
14% [lat,long]=mapxy(x,y,'n',71,39)
15%
16% See also MAPLL, LL2XY
17
18%check number of arguments.
19if ~((nargin==3) | (nargin==5)),
20 mapxyerrorusage();
21end
22
23%hem: must be either 's' or 'n'
24if ~ischar(hem),
25 error('mapxy error message: hemisphere argument should be ''n'' or ''s''');
26end
27
28%check hem is either 'n' or 's':
29if ~(strcmpi(hem,'n') | strcmpi(hem,'s')),
30 error('mapxy error message: hem should be ''s'' or ''n''');
31end
32
33%set sn:
34if strcmpi(hem,'s'),
35 slat=71;
36 sn=-1.0;
37 xlam=0;
38else
39 slat=70;
40 sn=1.0;
41 xlam=45;
42end
43
44%set defaults for standard parallels and centre meridians.
45if nargin==5,
46 slat=varargin{1};
47 xlam=varargin{2};
48end
49
50re = 6378137;
51e2 = 0.00669437999015;
52
53e=sqrt(e2);
54
55slat=slat/180*pi;
56
57rho=sqrt(x.^2+y.^2);
58cm=cos(slat)./sqrt(1.0-e2*(sin(slat).^2));
59t=tan((pi/4)-(slat/2.))./((1.0-e*sin(slat))./(1.0+e*sin(slat))).^(e/2.);
60
61t=rho.*t./(re*cm);
62
63chi=(pi/2.)-2.*atan(t);
64
65alat=chi+((e2/2.)+(5.0*e2^2/24.)+(e2^3/12.))*sin(2*chi)+...
66((7.0*e2^2/48.)+(29.*e2^3/240.))*sin(4.0*chi)+...
67(7.0*e2^3/120.)*sin(6.0*chi);
68
69alat=(sn*alat/pi*180);
70xpr=sn*x;
71ypr=sn*y;
72alon=atan2(xpr,-ypr)/pi*180-sn*xlam;
73alon=sn*alon;
74
75indice=find(alon<0);
76if length(indice)>0,
77 alon(indice)=alon(indice)+360;
78end
79
80indice=find(alon>360);
81if length(indice)>0,
82 alon(indice)=alon(indice)-360;
83end
84
85
86end
87
88
89function mapxyerrorusage()
90 help mapxy
91end
Note: See TracBrowser for help on using the repository browser.