Changeset 15126

Timestamp:

05/28/13 11:55:54 (12 years ago)

Author:

adhikari

Message:

CHG: xy2lambert function is added for reverse calculation of lambert2xy

Location:

issm/trunk-jpl/src/m/lambert

Files:

• lambert2xy.m (modified) (3 diffs)
• xy2lambert.m (added)

issm/trunk-jpl/src/m/lambert/lambert2xy.m

@@ -1 @@
-function [x,y] = lambert2xy(lat,lon,projection_center_lat,projection_center_lon)  
+function [x,y] = lambert2xy(lat,lon,sgn,projection_center_lat,projection_center_lon)  
 %LAMBERT2XY - converts lat long from Lambert Azimuthal to Polar Stereographic
 %
@@ -7 +7 @@
 %
 %   Usage:
-%      [x,y] = ll2xy(lat,lon)
-%      [x,y] = ll2xy(lat,lon,projection_center_lat,projection_center_lon)
+%      [x,y] = lambert2xy(lat,lon,sgn)
+%      [x,y] = lambert2xy(lat,lon,sgn,projection_center_lat,projection_center_lon)
 %
 %      - provide lat in [-90,90] and lon in [-180,180].
-%
-%      - default: North hemisphere [projection center lat = 90 lon=0]
-%                 South hemisphere [projection center lat = -90 lon=0]
+
+%      - sgn = +1 N hemisphere [default projection center lat = 90 lon=0]
+%              -1 S hemisphere [default projection center lat = -90 lon=0]
 
 %Get projection_center_lat and projection_center_lon 
-if nargin==4,
+if nargin==5,
         latitude0  = projection_center_lat;
         longitude0 = projection_center_lon;
-elseif nargin==2,
-        if lat>0,
+elseif nargin==3,
+        if sgn==1,
                 latitude0 = 90; longitude0 = 0;
                 disp('Info: creating coordinates in polar stereographic (Projection center lat: 90N lon: 0)');
-        elseif lat<0,
+        elseif sgn==-1,
                 latitude0 = -90; longitude0 = 0;
                 disp('Info: creating coordinates in polar stereographic (Projection center lat: 90S lon: 0)');
+        else
+                error('Sign should be either +1 or -1');
         end
 else
@@ -58 +60 @@
 
 % Calculation of x and y
-if(lat==90)
-        rho=a*sqrt(qp-q);
-        x=rho*sin(lam-lam0);
-        y=-rho*cos(lam-lam0);
-elseif(lat==-90)
-        rho=a*sqrt(qp+q);
-        x=rho*sin(lam-lam0);
-        y=rho*cos(lam-lam0);
+if(abs(lat)==90)
+        x=0.0; 
+        y=0.0; 
 else
         x=(B*D)*(cos(b)*sin(lam-lam0));
         y=(B/D)*((cos(b0)*sin(b))-(sin(b0)*cos(b)*cos(lam-lam0)));
 end
+