Home > atmlab > geodetic > r_crossing.m

r_crossing

PURPOSE ^

R_CROSSING Calculates the point where a radius is passed

SYNOPSIS ^

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

DESCRIPTION ^

 R_CROSSING   Calculates the point where a radius is passed

    Calculates where a 3D geometrical propagation path crosses a specified
    radius (r). 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 [lat,lon,l,za,aa] = r_crossing(r0,lat0,lon0,za0,aa0,r[,x,y,z,dx,dy,dz])
        
 OUT   lat   Latitude 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.
       r     Radius, 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 ^

r_crossing.m

SOURCE CODE ^

0001 % R_CROSSING   Calculates the point where a radius is passed
0002 %
0003 %    Calculates where a 3D geometrical propagation path crosses a specified
0004 %    radius (r). The propagation path is specified by giving a starting point
0005 %    and the line-of-sight at that point. If these data are also at hand as
0006 %    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 [lat,lon,l,za,aa] = r_crossing(r0,lat0,lon0,za0,aa0,r[,x,y,z,dx,dy,dz])
0013 %
0014 % OUT   lat   Latitude of crossing point.
0015 %       lon   Longitude of crossing point.
0016 %       l     Distance to crossing point.
0017 %       za    Zenith angle at crossing point
0018 %       aa    Zenith angle at crossing point.
0019 % IN    r0    Radius of starting point.
0020 %       lat   Latitude of starting point.
0021 %       lon   Longitude of starting point.
0022 %       za    Zenith angle at starting point.
0023 %       aa    Azimuth angle at starting point.
0024 %       r     Radius, for which crossing point shall be determined.
0025 % OPT   x     x-cartesian coordinate of starting point.
0026 %       y     y-cartesian coordinate of starting point.
0027 %       z     z-cartesian coordinate of starting point.
0028 %       dx    x-component of line-of-sight.
0029 %       dy    y-component of line-of-sight.
0030 %       dz    z-component of line-of-sight.
0031 
0032 % 2012-03-01   Created by Patrick Eriksson.
0033 
0034 function [lat,lon,l,za,aa] = r_crossing(r0,lat0,lon0,za0,aa0,r,x,y,z,dx,dy,dz)
0035 
0036   
0037 if nargin < 7
0038   [x,y,z,dx,dy,dz]=geocentricposlos2cart(r0,lat0,lon0,za0,aa0);
0039 end
0040   
0041   
0042 if ( r0 >= r  &  za0 <= 90 )  |  r0*sind(za0) > r
0043   lat = NaN;
0044   lon = NaN;
0045   l   = NaN;
0046   za  = NaN;
0047   aa  = NaN; 
0048 
0049 else
0050 
0051   p  = x.*dx + y.*dy + z.*dz;
0052   pp = p .* p;
0053   q  = x.*x + y.*y + z.*z - r.*r;
0054   sq = sqrt( pp - q );
0055   l1 = -p + sq;
0056   l2 = -p - sq;
0057 
0058   lmin = min( l1, l2 );
0059   lmax = max( l1, l2 );
0060 
0061   % If r0 equals r solutions just above zero can appear (that
0062   % shall be rejected). So we pick the biggest solution if lmin is
0063   % negative or just above zero.
0064   l = zeros(size(lmax));
0065   l(lmin < 1e-6) = lmax(lmin < 1e-6);
0066   l(lmin >= 1e-6) = lmin(lmin >= 1e-6);
0067   assert(all(l > 0));  
0068  
0069   if nargout <= 3
0070     RAD2DEG = constants( 'RAD2DEG' );
0071     lat     = RAD2DEG .* asin( ( z+dz.*l ) ./ r );
0072     lon     = RAD2DEG .* atan2( y+dy.*l, x+dx.*l );
0073   else
0074     [r2,lat,lon,za,aa] = cartposlos2geocentric(x+l.*dx,y+l.*dy,z+l.*dz,dx,dy,dz);
0075     assert( all(abs(r2-r) < 1e-6) );
0076   end
0077 
0078 end   
0079 end

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