Changes between Version 6 and Version 7 of shelf_thickness

Timestamp:

10/19/10 10:21:23 (14 years ago)

Author:

schodlok

Comment:

--

shelf_thickness

@@ -87 +87 @@
 md.bed=md.surface-md.thickness;
 }}}
+
+ 
+== stereographic conversion part 1 ==
+
+{{{
+#!m
+  function [x,y] = map2ll(lat,lon,sgn)  
+%     
+% NAME:            map2ll
+%===============================================================================    
+%   
+% USAGE:           [x,y] = map2ll(lat,lon,sgn)
+%
+% DESCRIPTION:
+%                 Converts from geodetic latitude and longitude to Polar 
+%                 Stereographic (X,Y) coordinates for the polar regions.
+% 
+% INPUT:    
+%                 Lat         = Geodetic Latitude  (degrees, +90 to -90)  
+%                 Lon         = Geodetic Longitude (degrees, 0 to 360)
+%                 SGN         = Sign of latitude
+%                                +1 : north latitude
+%                                -1 : south latitude
+%
+%
+% OUTPUT:         X  =  Polar Stereographic X Coordinate (km)
+%                 Y  =  Polar Stereographic Y Coordinate (km)
+%            
+%
+% AUTHOR:         The original equations are from Snyder, J. P.,
+%                 1982,  Map Projections Used by the U.S. Geological Survey,
+%                 Geological Survey Bulletin 1532, U.S. Government Printing 
+%                 Office.  See JPL Technical Memorandum 3349-85-101 for 
+%                 further details.
+%
+%                 Original FORTRAN program written by C. S. Morris,
+%                 April 1985, Jet Propulsion Laboratory, California
+%                 Institute of Technology
+%
+%                 IDL conversion by Helmut Schottmueller, September 1995,
+%                 Institute for environmental physics, University of Bremen
+%
+%                 Matlab conversion by Michael P. Schodlok, December 2003,
+%                 Alfred Wegener Institute for Polar and Marine Research, Bremerhaven
+%
+%                 Error correction for RHO: 1997-11-03 by LML
+%====================================================================================
+
+%-------------
+% CHECK INPUTS
+%-------------
+if nargin ~= 3
+   error('map2ll.m: must pass 3 parameters lat, lon, hemisphere')
+end %if
+
+% CHECK lat,lon,hems dimensions and verify consistent
+[ma,na] = size(lat);
+[mb,nb] = size(lon);
+[mc,nc] = size(sgn);
+
+% CHECK THAT lat & lon HAVE SAME SHAPE
+if (ma~=mb) | (na~=nb)
+   error('check dims: lat & lon must have same dimensions')
+end %if
+
+% CHECK THAT sgn IS A SCALAR
+if     (mc~=1  & np~=1 ) 
+   error('check sgn: sgn is either 1 or -1')
+end %if   
+
+
+%======
+% BEGIN
+%======
+
+
+% Conversion constant from degrees to radians
+  cde  = 57.29577951;
+% Standard latitude for the SSM/I grid with no distorsion
+  slat = 70.;
+% Radius of the earth in kilometers
+  re   = 6378.273;
+% Eccentricity of the Hughes ellipsoid squared
+  ex2   = .006693883;
+% Eccentricity of the Hughes ellipsoid
+  ex    =  sqrt(ex2);
+ 
+  if sgn == 1 
+    delta = 45.;
+  else
+    delta = 0.0;
+  end %if
+
+  latitude  = abs(lat) * pi/180.;
+  longitude = (lon + delta) * pi/180.;
+ 
+% compute X and Y in grid coordinates.
+  T = tan(pi/4-latitude/2) ./ ((1-ex*sin(latitude))./(1+ex*sin(latitude))).^(ex/2);
+
+   if (90 - slat) <  1.e-5 
+    rho = 2.*re*T/sqrt((1.+ex)^(1.+ex)*(1.-ex)^(1.-ex));
+   else
+    sl  = slat*pi/180.;
+    tc  = tan(pi/4.-sl/2.)/((1.-ex*sin(sl))/(1.+ex*sin(sl)))^(ex/2.);
+    mc  = cos(sl)/sqrt(1.0-ex2*(sin(sl)^2));
+    rho = re*mc*T/tc;
+   end % if
+  
+  
+  y = -rho .* sgn .* cos(sgn.*longitude);
+  x =  rho .* sgn .* sin(sgn.*longitude);
+
+  [cnt1,cnt2] = find(latitude >= pi / 2.);
+  
+ if cnt1
+   x(cnt1,1) = 0.0
+   y(cnt1,1) = 0.0
+ end %if
+
+ return
+%=========================================================================
+
+}}}
+
+== stereographic conversion part 2 ==
+
+   function [lat,lon] = map2xy(x,y,sgn)
+%     
+% NAME:            map2xy
+%===============================================================================    
+%   
+% USAGE:           [lat,lon] = map2xy(x,y,sgn)
+%   
+% DESCRIPTION:
+%                 Converts from geodetic latitude and longitude to Polar 
+%                 Stereographic (X,Y) coordinates for the polar regions.
+% 
+% INPUT:
+%                 X  =  Polar Stereographic X Coordinate (km)
+%                 Y  =  Polar Stereographic Y Coordinate (km)   
+%                 SGN         = Sign of latitude
+%                                +1 : north latitude
+%                                -1 : south latitude
+%
+%
+% OUTPUT:         
+%                 Lat         = Geodetic Latitude  (degrees, +90 to -90)  
+%                 Lon         = Geodetic Longitude (degrees, 0 to 360)
+%            
+%
+% AUTHOR:         The original equations are from Snyder, J. P.,
+%                 1982,  Map Projections Used by the U.S. Geological Survey,
+%                 Geological Survey Bulletin 1532, U.S. Government Printing 
+%                 Office.  See JPL Technical Memorandum 3349-85-101 for 
+%                 further details.
+%
+%                 Original FORTRAN program written by C. S. Morris,
+%                 April 1985, Jet Propulsion Laboratory, California
+%                 Institute of Technology
+%
+%                 IDL conversion by Helmut Schottmueller, September 1995,
+%                 Institute for environmental physics, University of Bremen
+%
+%                 Matlab conversion by Michael P. Schodlok, December 2003,
+%                 Alfred Wegener Institute for Polar and Marine Research, Bremerhaven
+%
+%                 Error correction for RHO: 1997-11-03 by LML
+%====================================================================================
+
+%-------------
+% CHECK INPUTS
+%-------------
+if nargin ~= 3
+   error('map2ll.m: must pass 3 parameters lat, lon, hemisphere')
+end %if
+
+% CHECK lat,lon,hems dimensions and verify consistent
+[ma,na] = size(x);
+[mb,nb] = size(y);
+[mc,nc] = size(sgn);
+
+% CHECK THAT lat & lon HAVE SAME SHAPE
+if (ma~=mb) | (na~=nb)
+   error('check dims: lat & lon must have same dimensions')
+end %if
+
+% CHECK THAT sgn IS A SCALAR
+if     (mc~=1  & np~=1 ) 
+   error('check sgn: sgn is either 1 or -1')
+end %if   
+
+
+%======
+% BEGIN
+%======
+
+
+% Conversion constant from degrees to radians
+  cde  = 57.29577951;
+% Standard latitude for the SSM/I grid with no distorsion
+  slat = 70.;
+% Radius of the earth in kilometers
+  re   = 6378.273;
+% Eccentricity of the Hughes ellipsoid squared
+  ex2   = .006693883;
+% Eccentricity of the Hughes ellipsoid
+  ex    =  sqrt(ex2);
+
+  if sgn == 1 
+    delta = 45.;
+  else
+    delta = 0.0;
+  end %if
+  
+  
+  sl  = slat*pi/180.;
+  rho = sqrt(x.^2 + y.^2);
+  cm = cos(sl) / sqrt(1.0 - ex2 * (sin(sl)^2));
+  T = tan((pi / 4.0) - (sl / 2.0)) / ((1.0 - ex * sin(sl)) / (1.0 + ex * sin(sl)))^(ex / 2.0);
+  
+  if  abs(slat-90.) < 1.e-5
+    T = rho * sqrt((1. + ex)^(1. + ex) * (1. - ex)^(1. - ex)) / 2. / re;
+  else
+    T = rho * T / (re * cm);
+  end %if
+ 
+  chi = (pi / 2.0) - 2.0 * atan(T);
+  lat = chi + ((ex2 / 2.0) + (5.0 * ex2^2.0 / 24.0) + (ex2^3.0 / 12.0)) * ...
+         sin(2 * chi) + ((7.0 * ex2^2.0 / 48.0) + (29.0 * ex2^3 / 240.0)) * ...
+         sin(4.0 * chi) + (7.0 * ex2^3.0 / 120.0) * sin(6.0 * chi) ;
+
+  lat = sgn * lat;
+  lon = atan2(sgn * x,-sgn * y);
+  lon = sgn * lon;
+  
+  [res1,res2] = find(rho <= 0.1);
+  if res1
+    lat(res1,1) = 90. * sgn;
+    lon(res1,1) = 0.0;
+  end %if
+
+  lon = lon * 180. / pi;
+  lat = lat * 180. / pi;
+  lon = lon - delta; 
+ 
+ return
+%=========================================================================
+ 
+}}}
+
+