% SATVIEWDEPTH   Satellite viewing depth
%
%    Calculates how far down into the atmosphere "can see". The viewing
%    depth is based on the optical thickness, calculated from the top of
%    the atmosphere and downwards. The atmosphere top is defined by the 
%    highest altitude for which absorption is provided.
%
%    Refraction is not considered. In this case, the propagation path is
%    best described by the tangent altitude. Down-looking cases correspond
%    to negative tangent altitudes.
%
% FORMAT   zs = satviewdepth(z,Abs,ztan,tau,Re)
%        
% OUT   zs    Viewing altitude, for each combination of frequency and tangent
%             altitudes.
% IN    z     Altitudes for absorption.
%       Abs   Absorption, where each column corresponds to an altitude.
%       ztan  Tangent altitude(s).
% OPT   tau   Definition of view limit, in optical thickness. Default is 1.
%       Re    Earth radius. Defualt is taken from *constants*.

% 2006-04-26   Created by Parick Eriksson.


function zs = satviewdepth(z,Abs,ztan,tau,Re)


tau_DEFAULT = 1;
Re_DEFAULT  = constants( 'EARTH_RADIUS' );
%
set_defaults;


if ~issorted(z)
  error('Argument *z* must be a sorted increasing vector.');
end

if length(z) ~= size(Abs,2)
  error('Size mismatch between *z* and *Abs*.');
end


%= Return argument
%
zs = zeros( length(ztan), size(Abs,1) );


%= Loop tangent altitudes
%
for j = 1:length(ztan)


  %= Matrix of cumulative optical depths
  %
  Tau = zeros( size(Abs) );
  %
  ind = 1:size(Abs,2);

  %= Loop downwards
  %
  for i = size(Abs,2)-1:-1:1

     %- Above tangent altitude
     %
     if z(i) > ztan(j)  
       l        = sqrt( (Re+z(i+1)).^2 - (Re+ztan(j))^2 ) - ...
                                         sqrt( (Re+z(i)).^2 - (Re+ztan(j))^2 );
       Tau(:,i) = Tau(:,i+1) +   l * (Abs(:,i+1)+Abs(:,i))/2; 

   
     %- Below tangent altitude
     %
     else
       l            = sqrt( (Re+z(i+1)).^2 - (Re+ztan(j))^2 );
       Tau(:,i)     = Tau(:,i+1) + l * (Abs(:,i+1)+interp1(z,Abs',ztan(j))')/2;
       ind          = i:ind(end);
       break
     end
  end


  %= Interpolate in tau
  %
  for i = 1:size(Abs,1)
    zs(j,i) = interp1( Tau(i,ind), z(ind), tau );
  end

end