function [ delay, Rcompl ] = genrefl( epsilon, mu, d, f )

% Function GENREFL ( Generalized reflection coefficient )
%
% Abstract
%
%    Computes the frequency-dependent complex electromagnetic reflectivity
%    of a plane parallel stratigraphy at normal incidence.
%
% Copyright
%
%    Copyright (2010), Istituto Nazionale di Astrofisica.
%
% Required Reading
%
%    None.
%
% Keywords
%
%    PROPAGATION
%
% Brief I/O
%
%    Variable  I/O  Description
%    --------  ---  --------------------------------------------------
%    epsilon    I   Relative complex permittivity.
%    mu         I   Relative complex permeability.
%    d          I   Depth of layer interfaces.
%    f          I   Frequencies for which the coefficient is to be computed.
%    delay      O   Time delay of the echo from layer interfaces.
%    Rcompl     O   Generalized reflection coefficient of the layers.
%
% Detailed Input
%
%    epsilon     is an array containing the relative complex permittivity
%                of the plane parallel layers, at every frequency for which
%                the reflection coefficient needs to be computed. Because
%                of the convention in the formulas used to compute the
%                reflection coefficient, the imaginary part of the complex
%                permittivity must be non-negative.
%
%    mu          is an array containing the relative complex permeability
%                of the plane parallel layers, at every frequency for which
%                the reflection coefficient needs to be computed. Because
%                of the convention in the formulas used to compute the
%                reflection coefficient, the imaginary part of the complex
%                permittivity must be non-negative.
%
%    d           is a vector containing the depth of interfaces between the
%                plane parallel layers. Depths are positive although they
%                are measured in the minus z direction.
%
%    f           is a vector containingthe set of frequencies for which the
%                reflection coefficient is to be computed.
%
% Detailed Output
%
%    delay       is an array containing the time delay of echoes from
%                interfaces within the planarly layered medium, for all
%                frequencies for which the reflection coefficient is to be
%                computed.
%
%    Rcompl      is a vector containing the generalized reflection
%                coefficient of the layers as a function of frequency.
%                Expressed as a phasor, it is the exact solution of
%                Maxwell's equations for a planarly layered medium in which
%                the electromagnetic wave arrives at normal incidence.
%
% Parameters
%
%    c           is the velocity of light in vacuo, set here to be equal to
%                299792458 m/s.
%
% Exceptions
%
%    1) The two input arrays epsilon and mu must have the same size, as
%       they contain the dielectric and magnetic properties of the layers.
%       The number of rows and columns of the two array are compared, and
%       if they differ, the error message 'Genrefl:ArraySizeMismatch' is
%       issued.
%
%    2) The two input arrays must have the same number of lines as there
%       are layers. The number of lines of the two arrays is compared with
%       the length of the vector d, and if they do not match the error
%       message 'Genrefl:WrongNumberOfArrayLines' is issued.
%
%    3) The two input arrays must have the same number of columns as there
%       are frequencies. The number of columns of the two arrays is
%       compared with the length of the vector f, and if they do not match
%       the error message 'Genrefl:WrongNumberOfArrayColumns' is issued.
%
%    4) The formula used to compute the generalized reflection coefficient
%       does not work for negative frequencies. If negative values are
%       present in vector f, the error message
%       'Genrefl:NegativeFrequencies' is issued.
%
%    5) Because of the convention in the formulas used to compute the
%       reflection coefficient, the imaginary part of both complex
%       permittivity and complex permeability must be non-positive. The
%       function checks the sign of the imaginary part of the permittivity
%       and permeability vectors, and issues the error message
%       'Genrefl:WrongConvention' if this is not true.
%
% Files
%
%    None.
%
% Details
%
%    The planarly layered medium for which the reflection coefficient is
%    computed is made of a number of regions of uniform permittivity and
%    permeability, the first and the last of which are semi-infinite. Such
%    regions and their interfaces are numbered according to the following
%    scheme:
%
%       Region 1         epsilon(1),   mu(1)
%       ----------------------------------------   -d(1),   R(1)
%       Region 2         epsilon(2),   mu(2)
%       ----------------------------------------   -d(2),   R(2)
%       Region 3         epsilon(3),   mu(3)
%       ----------------------------------------   -d(3),   R(3)
%            .
%            .
%            .
%       ----------------------------------------   -d(N-2), R(N-2)
%       Region N-1       epsilon(N-1), mu(N-1)
%       ----------------------------------------   -d(N-1), R(N-1)
%       Region N         epsilon(N),   mu(N)
%
%    where epsilon is the complex permittivity, mu the complex
%    permeability, d is the depth of the interface and R the Fresnel
%    reflection coefficient of the interface.
%
%    Following Chew (1990), a complex generalized reflection coefficient
%    for the planarly layered medium can be computed by means of the
%    following recursive relation:
%
%                       R(i) + Rcompl(i+1) * exp( 2 * j * k(i+1) * ( d(i+1) - d(i) ) )
%       Rcompl(i) = ------------------------------------------------------------------
%                   1 + R(i) * Rcompl(i+1) * exp( 2 * j * k(i+1) * ( d(i+1) - d(i) ) )
%
%    where j is the imaginary unit, k(i) is the wavenumber within the i-th
%    region, and Rcompl(i) is the complex reflection coefficient for all
%    layers from the N-th to the i-th. 
%
%    Rcompl(1) is then the generalized reflection coefficient of the layers
%    at the frequency corresponding to the wavenumber k. Expressed as a
%    phasor, it is the exact solution of Maxwell's equations for a planarly
%    layered medium in which the electromagnetic wave arrives at normal
%    incidence.
%
% Examples
%
%    1) The following script demonstrates the use of the function to
%       produce a simulated radar echo for a ground penetrating radar. The
%       script makes use of the functions FFTVARS, RICKER and GETLAYERS, to
%       which the reader is referred for explanations on their use:
%
%     clc
%     clear all
%     close all
% 
%     %
%     % The following physical constants are used in the computations:
%     %
%     %    c           Velocity of light in vacuo.
% 
%     c = 299792458;
% 
%     %
%     % The following parameters describe the radar setup.
%     %
%     %    deltat      Duration of the receiver window.
%     %    Ns          Number of samples in the echo.
%     %    fc          Central frequency in the transmitted signal.
% 
%     deltat = 2e-7;
%     Ns     = 512;
%     fc     = 1e8;
% 
%     %
%     % Several numerical parameters are computed:
%     %
%     %    dt          Sampling interval of the time window.
%     %    fs          Sampling frequency of the time window.
%     %    t           Times of acquisition of values in FFT computations.
%     %    f           Frequency spectrum for FFT computations.
% 
%     dt       = deltat / Ns;
%     fs       = 1 / dt;
%     [ t, f ] = fftvars( fs, Ns );
% 
%     %
%     % A discretely sampled Ricker wavelet and its spectrum are computed
%     % according to the above parameters:
% 
%     V = ricker ( fc, t );
%     S = fft ( V );
% 
%     %
%     % The electromagnetic parameters of the layered medium are defined:
% 
%     epsilon = [ 1, 3.15-0.003i, 9 ];
%     mu      = ones( size( epsilon ) );
%     d       = [ 2, 4.5, Inf ];
%     [ fdump, epsilon ] = meshgrid( f, epsilon );
%     [ fdump, mu      ] = meshgrid( f, mu      );
%     clear fdump
% 
%     %
%     % The complex reflectivity of the layered medium is computed:
% 
%     [ delay, Rcompl ] = genrefl( epsilon, mu, d, abs( f ) );
%     Rcompl( f < 0 ) = conj( Rcompl( f < 0 ) );
% 
%     %
%     % The simulated echo is obtained as the convolution of the transmitted
%     % signal with the reflectivity of the layered medium:
% 
%     signal = ifft( Rcompl .* S );
% 
%     %
%     % Results are displayed:
% 
%     plot( t ( 1 : Ns ), real( signal ( 1 : Ns  ) ) )
%     hold on
%     stem( delay,  max( abs( real( signal ) ) ) .* ones( size( delay ) ), 'r.-' )
%     stem( delay, -max( abs( real( signal ) ) ) .* ones( size( delay ) ), 'r.-' )
%     hold off
%     title( 'Received signal' )
%     legend( 'Signal', 'Interface positions' )
%     xlabel( 'Time (s)' )
%     ylabel( 'Voltage (V)')
%
% Restrictions
%
%    It is important to remember that the arrays describing the
%    electromagnetic properties of the layers must have a line for every
%    layer and a column for every frequency, in case electromagnetic
%    properties are frequency-dependent. Furthermore, lines in the array
%    must be ordered in order of increasing depth.
%
% Literature References
%
%    Chew, Weng Cho. Waves and fields in inhomogeneous media. Van Nostrand
%    Reinhold, New York, 1990.
%
% Author and Institution
%
%    R. Orosei   (INAF)
%
% Version
%
%    Version 1.0.0, 26-OCT-2010 (RO)
%
%       First release.
%
% Index Entries
%
%    None at this time.
%
% Revisions
%
%    None.
%

%
% The following physical constant is used in the computations:
%
%    c           Velocity of light in vacuo.

c = 299792458;

%
% Several auxiliary quantities are computed:
%
%    Nl          Number of layers.
%    Nf          Number of frequencies.

[ Nl,    Nf    ] = size( epsilon );

%
% Input arrays containing the dielectric and magnetic properties of the
% plane parallel layers must have the same size.

[ Nlmu, Nfmu ] = size( mu );

if Nl ~= Nlmu || Nf ~= Nfmu
    error( 'Genrefl:ArraySizeMismatch', ...
         [ 'The two input arrays epsilon and mu do not have the same size.\n', ...
           'epsilon is a %d X %d array, while mu is a %d X %d array.' ], ...
            Nl, Nf, Nlmu, Nfmu )
end

%
% Both epsilon and mu must have as many lines as there are layers, and as
% many columns as there are frequencies for which the reflection
% coefficient is to be computed.

if Nl ~= length( d ) || Nlmu ~= length( d )
    error( 'Genrefl:WrongNumberOfArrayLines', ...
         [ 'The two input arrays epsilon and mu must have as many lines as there are layers.\n', ...
           'epsilon has %d lines, mu has %d lines, and there are %d layers.' ], ...
            Nl, Nlmu, length( d ) )
end

if Nf ~= length( f ) || Nfmu ~= length( f )
    error( 'Genrefl:WrongNumberOfArrayColumns', ...
         [ 'The two input arrays epsilon and mu must have as many columns as there are frequencies.\n', ...
           'epsilon has %d columns, mu has %d columns, and there are %d frequencies.' ], ...
            Nf, Nfmu, length( f ) )
end

%
% The formula used to compute the generalized reflection coefficient does
% not work for negative frequencies. Negative frequencies are thus not
% allowed.

if any( f < 0 )
    error( 'Genrefl:NegativeFrequencies', ...
         [ 'Some input frequencies are negative.\n', ...
           'In this function frequencies are required to be non-negative.' ] )
end

%
% Because of the convention in the formulas used to compute the reflection
% coefficient, the imaginary part of both complex permittivity and complex
% permeability must be non-negative. This has been verified through
% testing.

if any( any( imag( epsilon ) > 0 ) ) || any( any( imag( mu ) > 0 ) )
    error( 'Genrefl:WrongConvention', ...
         [ 'The imaginary part of either permittivity or permeability is positive.\n', ...
           'In this function it is required that both have non-positive imaginary part.' ] )
end

% More auxiliary quantities are computed:
%
%    fgrid       Array containing the frequency vector in each line.
%    dgrid       Array containing the depth vector in each column.

[ fgrid, dgrid ] = meshgrid( f, d );
dgrid = [ zeros( 1, Nf ); dgrid ];

%
% The time delay of the echo from each layer interface is computed for
% every frequency, in case dielectric properties are frequency-dependent.

delay = cumsum( 2 .* diff( dgrid ) .* real( sqrt( epsilon .* mu ) ) ./ c );

%
% The Fresnel reflection coefficients at normal incidence are computed for
% each interface and for each frequency.

R = ( sqrt( epsilon( 1 : Nl - 1, : ) .* mu( 2 : Nl, : ) )          ...
    - sqrt( epsilon( 2 : Nl, : )     .* mu( 1 : Nl - 1, : ) ) ) ./ ...
    ( sqrt( epsilon( 1 : Nl - 1, : ) .* mu( 2 : Nl, : ) )          ...
    + sqrt( epsilon( 2 : Nl, : )     .* mu( 1 : Nl - 1, : ) ) )    ;

%
% The wave number is computed for each layer and for each frequency.

k   = 2 .* pi ./ c .* sqrt( epsilon .* mu ) .* fgrid;

%
% The recursive formula taken from Chew (1990) is now used to compute the
% frequency-dependent complex electromagnetic reflectivity of the layers at
% normal incidence. The sign of the complex exponential has been changed
% w.r.t. the formula listed in Chew (1990) because of the implementation of
% the FFT which would cause it to produce a mirrored signal.

Rcompl = R ( end, : );

for i = Nl - 2 : -1 : 1
    Rcompl = (     R( i, : ) +  Rcompl .* exp( - 2 .* 1i .* k ( i + 1, : ) .* ( d ( i + 1 ) - d ( i ) ) ) ) ...
          ./ ( 1 + R( i, : ) .* Rcompl .* exp( - 2 .* 1i .* k ( i + 1, : ) .* ( d ( i + 1 ) - d ( i ) ) ) ) ;
end

%
% The time delay due to the depth of the first interface is introduced
% according to the shift theorem for Fourier transforms.

Rcompl = Rcompl .* exp( - 2 .* pi .* 1i .* 2. * d ( 1 ) ./ c .* f );
