Home > atmlab > geodetic > cartposlos2geocentric.m

cartposlos2geocentric

PURPOSE ^

CARTPOSLOS2GEOCENTRIC Converts cartesian POS/LOS to spherical coordinates

SYNOPSIS ^

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

DESCRIPTION ^

 CARTPOSLOS2GEOCENTRIC   Converts cartesian POS/LOS to spherical coordinates 

   Position is given as (x,y,z), while line-of-sight is given as
   (dx,dy,dz). The corresponding quantities in polar coordinates are
   (r,lat,lon) and (za,aa), respectively.

   See *Contents* for defintion of coordinate systems.

   If the optional arguments are given, it is ensured that latitude and
   longitude are kept constant for zenith or nadir cases, and the longitude
   and azimuth angle for N-S cases. The optional input shall be interpreted
   as the [x,y,z] is obtained by moving from [r0,lat0,lon0] in the direction
   of [za0,aa0].

 FORMAT   [r,lat,lon,za,aa] = cartposlos2geocentric(x,y,z,dx,dy,dz,...
                                                [ppc,lat0,lon0,za0,aa0] )
        
 OUT   r      Radius.
       lat    Latitude.
       lon    Longitude.
       za     Zenith angle.
       aa     Azimuth angle.
 IN    x      Coordinate in x dimension.
       y      Coordinate in y dimension.
       z      Coordinate in z dimension.
       dx     LOS component in x dimension.
       dy     LOS component in y dimension.
       dz     LOS component in z dimension.
       ppc    Propagation path constant = r0*sin(za0)
       lat0   Original latitude
       lon0   Original longitude
       za0    Orignal zenith angle
       aa0    Orignal azimuth angle

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

cartposlos2geocentric.m

SOURCE CODE ^

0001 % CARTPOSLOS2GEOCENTRIC   Converts cartesian POS/LOS to spherical coordinates
0002 %
0003 %   Position is given as (x,y,z), while line-of-sight is given as
0004 %   (dx,dy,dz). The corresponding quantities in polar coordinates are
0005 %   (r,lat,lon) and (za,aa), respectively.
0006 %
0007 %   See *Contents* for defintion of coordinate systems.
0008 %
0009 %   If the optional arguments are given, it is ensured that latitude and
0010 %   longitude are kept constant for zenith or nadir cases, and the longitude
0011 %   and azimuth angle for N-S cases. The optional input shall be interpreted
0012 %   as the [x,y,z] is obtained by moving from [r0,lat0,lon0] in the direction
0013 %   of [za0,aa0].
0014 %
0015 % FORMAT   [r,lat,lon,za,aa] = cartposlos2geocentric(x,y,z,dx,dy,dz,...
0016 %                                                [ppc,lat0,lon0,za0,aa0] )
0017 %
0018 % OUT   r      Radius.
0019 %       lat    Latitude.
0020 %       lon    Longitude.
0021 %       za     Zenith angle.
0022 %       aa     Azimuth angle.
0023 % IN    x      Coordinate in x dimension.
0024 %       y      Coordinate in y dimension.
0025 %       z      Coordinate in z dimension.
0026 %       dx     LOS component in x dimension.
0027 %       dy     LOS component in y dimension.
0028 %       dz     LOS component in z dimension.
0029 %       ppc    Propagation path constant = r0*sin(za0)
0030 %       lat0   Original latitude
0031 %       lon0   Original longitude
0032 %       za0    Orignal zenith angle
0033 %       aa0    Orignal azimuth angle
0034 
0035 % 2011-11-15   Created by Bengt Rydberg.
0036 
0037 
0038 function [r,lat,lon,za,aa] = cartposlos2geocentric(x,y,z,dx,dy,dz,...
0039                                                    ppc,lat0,lon0,za0,aa0)
0040 
0041 deg2rad = constants( 'DEG2RAD' );
0042 rad2deg = constants( 'RAD2DEG' );
0043 
0044 
0045 if nargin <= 6
0046   [r,lat,lon] = cart2geocentric( x, y, z );
0047 else
0048   [r,lat,lon] = cart2geocentric( x, y, z, lat0, lon0, za0, aa0 );    
0049 end
0050 
0051 
0052 coslat = cos( deg2rad * lat );
0053 sinlat = sin( deg2rad * lat );
0054 coslon = cos( deg2rad * lon );
0055 sinlon = sin( deg2rad * lon );
0056 dr     = coslat.*coslon.*dx + sinlat.*dz + coslat.*sinlon.*dy;
0057 
0058 
0059 % LOS angles
0060 if nargin <= 6
0061   za = rad2deg * acos( dr );
0062 else
0063   za = rad2deg * asin( ppc ./ r );
0064 end
0065 
0066 aa = zeros(size(za));
0067 
0068 % FIXME: surely at least some (if not all) cases here can be vectorised
0069 
0070 for i = 1 : length(za)
0071 
0072   % Fixes for za with optional input
0073   if nargin > 6
0074     if za0(i) < 1e-6  ||  za0(i) > 180-1e-6    % Nadir or zenith
0075       za(i) = za0(i);
0076     elseif isnan(za(i))
0077       za(i) = 90;
0078     elseif dr(i) < 0              % Expression above gives only za < 90
0079       za(i) = 180 - za(i);
0080     end
0081   end
0082 
0083   % Azimuth angle:
0084     
0085   % Take original value if N-S
0086   if nargin > 6  &&  ( abs(aa0(i)) < 1e-6  ||  abs(aa0(i)-180) < 1e-6 )
0087     % Check that not lon changed with 180 deg
0088     if lon(i) == lon0(i)
0089       aa(i) = aa0(i);
0090     else
0091       if abs(aa0(i)) < 1e-6
0092         aa(i) = 180;
0093       else
0094         aa(i) = 0;
0095       end
0096     end
0097   
0098   % Set to zero for nadir and zenith
0099   elseif za(i) < 1e-6  ||  za(i) > 180-1e-6 
0100     aa(i) = 0;
0101     
0102   
0103   % For lat = +- 90 the azimuth angle gives the longitude along which
0104   % the LOS goes
0105   elseif abs( lat(i) ) > 90-1e-8 
0106     aa(i) = rad2deg * atan2( dy(i), dx(i) );
0107 
0108   %General case
0109   else
0110     %
0111     dlat  = -sinlat(i)*coslon(i)/r(i)*dx(i) + coslat(i)/r(i)*dz(i) - ...
0112              sinlat(i)*sinlon(i)/r(i)*dy(i);
0113     dlon  = -sinlon(i)/coslat(i)/r(i)*dx(i) + coslon(i)/coslat(i)/r(i)*dy(i);
0114     aa(i) = rad2deg * acos( r(i) * dlat / sin( deg2rad * za(i) ) );
0115     %
0116     if isnan( aa(i) ) ||  ~isreal( aa(i) )
0117       if dlat >= 0
0118         aa(i) = 0;
0119       else
0120         aa(i) = 180;
0121       end
0122     elseif dlon < 0
0123       aa(i) = -aa(i);
0124     end
0125   end
0126 end
0127 
0128 end

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