Home > atmlab > geodetic > geomtanpoint2d.m

geomtanpoint2d

PURPOSE ^

GEOMTANPOINT2D 2D Geometrical tangent point

SYNOPSIS ^

function [rt,latt] = geomtanpoint2d( ellipsoid, r, lat, za )

DESCRIPTION ^

 GEOMTANPOINT2D   2D Geometrical tangent point

    Calculates the 2D geometrical tangent point for arbitrary reference
    ellipsiod. That is, a spherical planet is not assumed. The tangent
    point is thus defined as the point closest to the ellipsoid (not as the
    ppoint with za=90).
   
    Geocentric coordinates are used for both sensor and tangent point
    positions. Following ARTS, 2D latitudes are not limited to [-90,90] and
    this function keeps track on the actual 2D latitude (in contrast to
    *geomtanpoint*).
 
    The algorithm used for non-spherical cases is derived by Nick Lloyd at
    University of Saskatchewan, Canada (nick.lloyd@usask.ca), and is part of
    the operational code for both OSIRIS and SMR on-board- the Odin
    satellite. His 3D code is here adopted to 2D.
 
 FORMAT  [rt,latt,lont] = geomtanpoint( ellipsoid, r, lat, lon, za, aa )

 OUT   rt         Radius of tangent point.
       latt       Latitude of tangent point.
 IN    ellipsoid  Model ellipsoid as returned by *ellipsoidmodels*.
       r          Radius of sensor position.
       lat        Latitude of sensor position.
       za         Zenith angle of sensor line-of-sight.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

geomtanpoint2d.m

SOURCE CODE ^

0001 % GEOMTANPOINT2D   2D Geometrical tangent point
0002 %
0003 %    Calculates the 2D geometrical tangent point for arbitrary reference
0004 %    ellipsiod. That is, a spherical planet is not assumed. The tangent
0005 %    point is thus defined as the point closest to the ellipsoid (not as the
0006 %    ppoint with za=90).
0007 %
0008 %    Geocentric coordinates are used for both sensor and tangent point
0009 %    positions. Following ARTS, 2D latitudes are not limited to [-90,90] and
0010 %    this function keeps track on the actual 2D latitude (in contrast to
0011 %    *geomtanpoint*).
0012 %
0013 %    The algorithm used for non-spherical cases is derived by Nick Lloyd at
0014 %    University of Saskatchewan, Canada (nick.lloyd@usask.ca), and is part of
0015 %    the operational code for both OSIRIS and SMR on-board- the Odin
0016 %    satellite. His 3D code is here adopted to 2D.
0017 %
0018 % FORMAT  [rt,latt,lont] = geomtanpoint( ellipsoid, r, lat, lon, za, aa )
0019 %
0020 % OUT   rt         Radius of tangent point.
0021 %       latt       Latitude of tangent point.
0022 % IN    ellipsoid  Model ellipsoid as returned by *ellipsoidmodels*.
0023 %       r          Radius of sensor position.
0024 %       lat        Latitude of sensor position.
0025 %       za         Zenith angle of sensor line-of-sight.
0026 
0027 % 2012-02-13 Created by Patrick Eriksson
0028 
0029 function [rt,latt] = geomtanpoint2d( ellipsoid, r, lat, za )
0030 
0031   
0032 % Check input                                                          %&%
0033 rqre_datatype( ellipsoid, @isvector );                                 %&%
0034 if length(ellipsoid) ~= 2                                              %&%
0035  error( 'The argument *ellipsoid* must be a vector of length 2' );     %&%
0036 end                                                                    %&%
0037 rqre_in_range( ellipsoid(2), 0, 1, 'element 2 of ellipsoid' );         %&%
0038 rqre_datatype( r, @istensor0 );                                        %&%
0039 rqre_datatype( lat, @istensor0 );                                      %&%
0040 rqre_datatype( za, @istensor0 );                                       %&%
0041 rqre_in_range( za, -180, 180 )
0042 
0043   
0044 %- Spherical planet
0045 %
0046 if ellipsoid(2) < 1e-7
0047 
0048   rt = r * sind( abs(za) );
0049   if za > 0
0050     latt = lat + za - 90;
0051   else
0052     latt = lat + za + 90;  
0053   end
0054   
0055   
0056 %- Ellipsoidal planet
0057 %
0058 else
0059   
0060   % Equatorial and polar radii squared:
0061   a2 = ellipsoid(1)^2;
0062   b2 = a2 * ( 1 - ellipsoid(2)^2 ); 
0063 
0064   % Code from geocentricposlos2cart adopted to 2D
0065   %
0066   coslat = cosd( lat );
0067   sinlat = sind( lat );
0068   cosza  = cosd( za );
0069   sinza  = sind( za );
0070   %
0071   x    = r * coslat; 
0072   z    = r * sinlat;  
0073   %
0074   dr   = cosza;
0075   dlat = sinza;
0076   %
0077   dx = coslat*dr - sinlat*dlat;
0078   dz = sinlat*dr + coslat.*dlat;
0079    
0080   % For testing:
0081   %[x2,y2,z2,dx2,dy2,dz2] = geocentricposlos2cart( r, lat, 0, za, 0 );
0082   %  keyboard
0083 
0084   X = [x,0,z]';
0085   
0086   xunit = [dx,0,dz]';
0087   
0088   zunit = cross( xunit, X );
0089   zunit = zunit / norm(zunit);
0090      
0091   yunit = cross( zunit, xunit );
0092   yunit = yunit / norm(yunit);
0093   
0094   w11 = xunit(1);
0095   w12 = yunit(1);
0096   w31 = xunit(3);
0097   w32 = yunit(3);
0098   yr  = dot( X, yunit );
0099   xr  = dot( X, xunit );
0100 
0101   A = w11*w11/a2 + w31*w31/b2;
0102   B = 2.0*(w11*w12/a2 + w31*w32/b2);
0103   C = w12*w12/a2 + w32*w32/b2;
0104      
0105   if B == 0.0
0106     xx = 0.0;
0107   else 
0108     K = -2.0*A/B;
0109     factor = 1.0/(A+(B+C*K)*K);
0110     xx = sqrt(factor);
0111     yy = K*xx;
0112   end 
0113      
0114   dist1 = (xr-xx)*(xr-xx) + yr*yr;
0115   dist2 = (xr+xx)*(xr+xx) + yr*yr;
0116      
0117   if dist1 > dist2
0118     xx = -xx;
0119   end
0120   
0121   xn = w11*xx + w12*yr;
0122   zn = w31*xx + w32*yr;
0123   
0124   rt   = sqrt( xn^2 + zn^2 );
0125   l    = sqrt( (xn-x)^2 + (zn-z)^2 );
0126   dlat = asind( l*sinza/rt );
0127   
0128   latt = lat + dlat;
0129 end
0130 
0131 
0132

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005