Home > atmlab > geodetic > lat_crossing.m

lat_crossing

PURPOSE ^

LAT_CROSSING Calculates the point where a latitude is passed

SYNOPSIS ^

function [r,lon,l,za,aa] = lat_crossing(r0,lat0,lon0,za0,aa0,lat,x,y,z,dx,dy,dz)

DESCRIPTION ^

 LAT_CROSSING   Calculates the point where a latitude is passed

    Calculates where a 3D geometrical propagation path crosses a specified
    latitude (lat). The propagation path is specified by giving a starting
    point and the line-of-sight at that point. If these data are also at hand
    as Cartesian coordinates (x,y,z,dx,dy,dz), input these as well and the
    calculations will be quicker.

    Allowed solutions have a distance above zero (ie. l>0, note that the
    solution of l=0 is rejected).

 FORMAT [r,lon,l,za,aa] = lat_crossing(r0,lat0,lon0,za0,aa0,lat
                                                          [,x,y,z,dx,dy,dz])
        
 OUT   r     Radius of crossing point.
       lon   Longitude of crossing point.
       l     Distance to crossing point.
       za    Zenith angle at crossing point
       aa    Zenith angle at crossing point.
 IN    r0    Radius of starting point.
       lat   Latitude of starting point.
       lon   Longitude of starting point.
       za    Zenith angle at starting point.
       aa    Azimuth angle at starting point.
       lat   Latitude, for which crossing point shall be determined.
 OPT   x     x-cartesian coordinate of starting point.
       y     y-cartesian coordinate of starting point.
       z     z-cartesian coordinate of starting point.
       dx    x-component of line-of-sight.
       dy    y-component of line-of-sight.
       dz    z-component of line-of-sight.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

lat_crossing.m

SOURCE CODE ^

0001 % LAT_CROSSING   Calculates the point where a latitude is passed
0002 %
0003 %    Calculates where a 3D geometrical propagation path crosses a specified
0004 %    latitude (lat). The propagation path is specified by giving a starting
0005 %    point and the line-of-sight at that point. If these data are also at hand
0006 %    as Cartesian coordinates (x,y,z,dx,dy,dz), input these as well and the
0007 %    calculations will be quicker.
0008 %
0009 %    Allowed solutions have a distance above zero (ie. l>0, note that the
0010 %    solution of l=0 is rejected).
0011 %
0012 % FORMAT [r,lon,l,za,aa] = lat_crossing(r0,lat0,lon0,za0,aa0,lat
0013 %                                                          [,x,y,z,dx,dy,dz])
0014 %
0015 % OUT   r     Radius of crossing point.
0016 %       lon   Longitude of crossing point.
0017 %       l     Distance to crossing point.
0018 %       za    Zenith angle at crossing point
0019 %       aa    Zenith angle at crossing point.
0020 % IN    r0    Radius of starting point.
0021 %       lat   Latitude of starting point.
0022 %       lon   Longitude of starting point.
0023 %       za    Zenith angle at starting point.
0024 %       aa    Azimuth angle at starting point.
0025 %       lat   Latitude, for which crossing point shall be determined.
0026 % OPT   x     x-cartesian coordinate of starting point.
0027 %       y     y-cartesian coordinate of starting point.
0028 %       z     z-cartesian coordinate of starting point.
0029 %       dx    x-component of line-of-sight.
0030 %       dy    y-component of line-of-sight.
0031 %       dz    z-component of line-of-sight.
0032 
0033 % 2012-03-01   Created by Patrick Eriksson.
0034 
0035 function [r,lon,l,za,aa] = lat_crossing(r0,lat0,lon0,za0,aa0,lat,x,y,z,dx,dy,dz)
0036 
0037   
0038 if nargin < 7
0039   [x,y,z,dx,dy,dz]=geocentricposlos2cart(r0,lat0,lon0,za0,aa0);
0040 end
0041   
0042 
0043 % For za=0/180 and lat0+-90 there is no solution
0044 if za0 == 0  |  za0 == 180  |  abs(lat) == 90
0045   l = -1;
0046 
0047 % The expressions below can not be used for lat=0
0048 elseif abs(lat) < 1e-7    % lat=1e-8 is tested not work
0049   l = -z / dz;
0050   
0051 else
0052   
0053   t2     = tand(lat)^2;
0054   a      = t2 * ( dx*dx + dy*dy ) - dz*dz;
0055   b      = 2 * ( t2 * ( x*dx + y*dy ) - z*dz );
0056   c      = t2 * ( x*x + y*y ) - z*z;
0057   bb     = b*b;
0058   ac4    = 4*a*c;
0059   
0060   % Check if a real solution is possible
0061   if ac4 > bb
0062     l = -1;
0063   else
0064     d    = -0.5*b/a;
0065     e    = -0.5*sqrt(b*b-4*a*c)/a;      
0066     l1   = d + e;
0067     l2   = d - e;
0068 
0069     % If both l1 and l2 are > 0, we want theoretically the smallest
0070     % value. However, with lat=lat0 the "zero solution" can deviate
0071     % slightly from zero due to numerical issues, and it is not just to
0072     % pick the smallest positive value. As a solution, don't except a
0073     % final l below 1e-6 if not both l1 and l2 are inside [0,1e-6].
0074     lmin = min( l1, l2 );
0075     lmax = max( l1, l2 );
0076     if( lmin >= 0  && lmax < 1e-6 )
0077       l = lmax;
0078     else
0079       if( lmin > 1e-6 ) 
0080         l = lmin;
0081       elseif( lmax > 1e-6 )
0082         l = lmax;
0083       else
0084         l = -1;
0085       end
0086     end
0087   end
0088 end
0089 
0090 
0091 if l <= 0
0092   r   = NaN;
0093   lon = NaN;
0094   l   = NaN;
0095   za  = NaN;
0096   aa  = NaN; 
0097 
0098 else
0099  
0100   if nargout <= 3
0101     RAD2DEG = constants( 'RAD2DEG' );
0102     xp      = x+dx*l;
0103     yp      = y+dy*l;
0104     r       = sqrt( xp^2 + yp^2 + (z+dz*l)^2 );
0105     lon     = RAD2DEG * atan2( yp, xp );
0106   else
0107     [r,lat2,lon,za,aa] = cartposlos2geocentric(x+l*dx,y+l*dy,z+l*dz,dx,dy,dz);
0108     assert( abs(lat2-lat) < 1e-6 );
0109   end
0110   
0111 end

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