Home > atmlab > arts > various > arts_rslope_crossing.m

arts_rslope_crossing

PURPOSE ^

ARTS_RSLOPE_CROSSING Testing of arts' rslope_crossing function

SYNOPSIS ^

function [dlat,dlat2] = arts_rslope_crossing(rp,za,r0,c)

DESCRIPTION ^

 ARTS_RSLOPE_CROSSING   Testing of arts' rslope_crossing function

    The function is a Matlab implementation of the rslope_crossing function
    found in ppath.cc of arts. The output *dlat* is as calculated inside
    arts, while *dlat2* is a comparison value that should be accurate with
    a precision of 0.0001 deg. The search for *dlat2* is limited to [-1,1].
    The search for dlat2 can return the crossing on the other side of a
    tangent point.

 FORMAT   [dlat,dlat2] = arts_rslope_crossing(rp,za,r0,c)
        
 OUT   dlat   Latitude difference calculated by Taylor expansion.
       dlat2  Same as *dlat*, but calculated by trail and error.
 IN    rp     Ppath radius.
       za     Ppath zenith angle.
       r0     Radius of level at latitude matching *rp*
       c      Slope of level [m/deg]

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

arts_rslope_crossing.m

SOURCE CODE ^

0001 % ARTS_RSLOPE_CROSSING   Testing of arts' rslope_crossing function
0002 %
0003 %    The function is a Matlab implementation of the rslope_crossing function
0004 %    found in ppath.cc of arts. The output *dlat* is as calculated inside
0005 %    arts, while *dlat2* is a comparison value that should be accurate with
0006 %    a precision of 0.0001 deg. The search for *dlat2* is limited to [-1,1].
0007 %    The search for dlat2 can return the crossing on the other side of a
0008 %    tangent point.
0009 %
0010 % FORMAT   [dlat,dlat2] = arts_rslope_crossing(rp,za,r0,c)
0011 %
0012 % OUT   dlat   Latitude difference calculated by Taylor expansion.
0013 %       dlat2  Same as *dlat*, but calculated by trail and error.
0014 % IN    rp     Ppath radius.
0015 %       za     Ppath zenith angle.
0016 %       r0     Radius of level at latitude matching *rp*
0017 %       c      Slope of level [m/deg]
0018 
0019 % 2012-02-21   Created by Patrick Eriksson.
0020 
0021 function [dlat,dlat2] = arts_rslope_crossing(rp,za,r0,c)
0022 
0023   DEG2RAD = constants( 'DEG2RAD' );
0024   RAD2DEG = constants( 'RAD2DEG' );
0025   
0026   % If r0=rp, numerical inaccuracy can give a false solution, very close
0027   % to 0, that we must throw away.
0028   dmin = 0;
0029   if( r0 == rp )
0030     dmin = 1e-12;
0031   end
0032   
0033   % The nadir angle in radians, and cosine and sine of that angle
0034   beta = DEG2RAD * ( 180 - abs(za) );
0035   cv = cos( beta );
0036   sv = sin( beta );
0037 
0038   % Convert slope to m/radian and consider viewing direction
0039   c = RAD2DEG*c;
0040   if( za < 0 )
0041     c = -c;
0042   end
0043   
0044   % The vector of polynomial coefficients
0045   p = zeros(5,1);
0046   %
0047   p(5) = ( r0 - rp ) * sv;
0048   p(4) = r0 * cv + c * sv;
0049   p(3) = -r0 * sv / 2 + c * cv;
0050   p(2) = -r0 * cv / 6 - c * sv / 2;
0051   p(1) = -c * cv / 6;
0052 
0053   % Calculate roots of the polynomial
0054   rs = roots(p);
0055   
0056   % Find the smallest root with imaginary part = 0, and real part > 0.
0057   %
0058   dlat = 1.571;
0059   %
0060   for i = 1 : length(rs)
0061     if isreal(rs(i))  &  rs(i) < dlat  &  rs(i) > dmin
0062       dlat = rs(i); 
0063     end
0064   end  
0065   
0066   % Convert back to degrees
0067   % Change also sign if zenith angle is negative
0068   if( dlat < 1.57 )
0069     
0070     dlat = RAD2DEG * dlat;
0071   
0072     if za < 0
0073       dlat = -dlat; 
0074     end
0075   else
0076     dlat = 99e99;
0077   end
0078   
0079   
0080   % Second solution
0081   if nargout > 1
0082     
0083     c = DEG2RAD * c;
0084     dr = Inf;
0085   
0086     for latt = -1 : 0.0001 : 1
0087 
0088       r1  = r0 + c*latt;
0089       r2  = (rp*sind(abs(za))) / sind(abs(za)-latt);
0090       drt = abs( r1 - r2 );
0091       %
0092       if drt < dr
0093         dr    = drt;
0094         dlat2 = latt;
0095       end
0096     end
0097     
0098     if za < 0
0099       dlat2 = -dlat2;
0100     end
0101 
0102     if dr > 1
0103       dlat2 = 99e99;
0104       fprintf( 'Closest distance for dlat2 = %.2f\n', dr );
0105     end
0106   end

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