Home > atmlab > geodetic > geomtanpoint.m

geomtanpoint

PURPOSE ^

GEOMTANPOINT Geometrical tangent point

SYNOPSIS ^

function [rt,latt,lont] = geomtanpoint( ellipsoid, r, lat, lon, za, aa )

DESCRIPTION ^

 GEOMTANPOINT   Geometrical tangent point

    Calculates the 3D 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.
 
    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.
 
 FORMAT  [rt,latt,lont] = geomtanpoint( ellipsoid, r, lat, lon, za, aa )

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

geomtanpoint.m

SOURCE CODE ^

0001 % GEOMTANPOINT   Geometrical tangent point
0002 %
0003 %    Calculates the 3D 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.
0010 %
0011 %    The algorithm used for non-spherical cases is derived by Nick Lloyd at
0012 %    University of Saskatchewan, Canada (nick.lloyd@usask.ca), and is part of
0013 %    the operational code for both OSIRIS and SMR on-board- the Odin
0014 %    satellite.
0015 %
0016 % FORMAT  [rt,latt,lont] = geomtanpoint( ellipsoid, r, lat, lon, za, aa )
0017 %
0018 % OUT   rt         Radius of tangent point.
0019 %       latt       Latitude of tangent point.
0020 %       lont       Longitude of tangent point.
0021 % IN    ellipsoid  Model ellipsoid as returned by *ellipsoidmodels*.
0022 %       r          Radius of sensor position.
0023 %       lat        Latitude of sensor position.
0024 %       lon        Longitude of sensor position.
0025 %       za         Zenith angle of sensor line-of-sight.
0026 %       aa         Azimuth angle of sensor line-of-sight.
0027 
0028 % 2012-02-10 Created by Patrick Eriksson
0029 
0030 function [rt,latt,lont] = geomtanpoint( ellipsoid, r, lat, lon, za, aa )
0031 
0032   
0033 % Check input                                                          %&%
0034 rqre_datatype( ellipsoid, @isvector );                                 %&%
0035 if length(ellipsoid) ~= 2                                              %&%
0036  error( 'The argument *ellipsoid* must be a vector of length 2' );     %&%
0037 end                                                                    %&%
0038 rqre_in_range( ellipsoid(2), 0, 1, 'element 2 of ellipsoid' );         %&%
0039 rqre_datatype( r, @istensor0 );                                        %&%
0040 rqre_datatype( lat, @istensor0 );                                      %&%
0041 rqre_datatype( lon, @istensor0 );                                      %&%
0042 rqre_datatype( za, @istensor0 );                                       %&%
0043 rqre_datatype( aa, @istensor0 );                                       %&%
0044 % Values of r, lat, lon, za and aa checked by geocentricposlos2cart    %&%
0045 
0046 [x,y,z,dx,dy,dz] = geocentricposlos2cart( r, lat, lon, za, aa );
0047 
0048 
0049 %- Spherical planet
0050 %
0051 if ellipsoid(2) < 1e-7
0052    
0053   ppc = r * sind( abs(za) );
0054   
0055   l_tan = sqrt( r*r - ppc*ppc );
0056    
0057   [rt,latt,lont] = cart2geocentric( x+dx*l_tan, y+dy*l_tan, z+dz*l_tan );
0058  
0059   
0060 %- Ellipsoidal planet
0061 %
0062 else
0063   
0064   % Equatorial and polar radii squared:
0065   a2 = ellipsoid(1)^2;
0066   b2 = a2 * ( 1 - ellipsoid(2)^2 ); 
0067   
0068   X = [x,y,z]';
0069   
0070   xunit = [dx,dy,dz]';
0071   
0072   zunit = cross( xunit, X );
0073   zunit = zunit / norm(zunit);
0074      
0075   yunit = cross( zunit, xunit );
0076   yunit = yunit / norm(yunit);
0077   
0078   w11 = xunit(1);
0079   w12 = yunit(1);
0080   w21 = xunit(2);
0081   w22 = yunit(2);
0082   w31 = xunit(3);
0083   w32 = yunit(3);
0084   yr  = dot( X, yunit );
0085   xr  = dot( X, xunit );
0086         
0087   A = (w11*w11 + w21*w21)/a2 + w31*w31/b2;
0088   B = 2.0*((w11*w12 + w21*w22)/a2 + (w31*w32)/b2);
0089   C = (w12*w12 + w22*w22)/a2 + w32*w32/b2;
0090      
0091   if B == 0.0
0092     xx = 0.0;
0093   else 
0094     K = -2.0*A/B;
0095     factor = 1.0/(A+(B+C*K)*K);
0096     xx = sqrt(factor);
0097     yy = K*x;
0098   end 
0099      
0100   dist1 = (xr-xx)*(xr-xx) + (yr-yy)*(yr-yy);
0101   dist2 = (xr+xx)*(xr+xx) + (yr+yy)*(yr+yy);
0102      
0103   if dist1 > dist2
0104     xx = -xx;
0105   end
0106 
0107   [rt,latt,lont] = cart2geocentric( w11*xx + w12*yr, ...
0108                                     w21*xx + w22*yr, ...
0109                                     w31*xx + w32*yr );
0110 end
0111 
0112 
0113

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