function epsilon = matzler( f, T )

% ATTENZIONE: C'E' ANCORA DEL LAVORO DA FARE, POICHE' IL METODO DI INTERPOLAZIONE
% DEI VALORI DELLA PARTE REALE IN FUNZIONE DELLA TEMPERATURA PRODUCE DEGLI
% SBALZI TROPPO NETTI. RICORDA CHE OGNI MODIFICA VA DOCUMENTATA. PER ORA
% STO SPERIMENTANDO CON UNA QUARTICA

% Function MATZLER ( Relative dielectric permittivity of pure water ice )
%
% Abstract
%
%    Computes the relative dielectric permittivity of pure water ice (ice
%    Ih) for any temperature and frequency according to formulas provided
%    in Matzler (1998).
%
% Copyright
%
%    Copyright (2010), Istituto Nazionale di Astrofisica.
%
% Required Reading
%
%    None.
%
% Keywords
%
%    DIELECTRIC PERMITTIVITY
%    WATER ICE
%
% Brief I/O
%
%    Variable  I/O  Description
%    --------  ---  --------------------------------------------------
%    f          I   Frequency.
%    T          I   Absolute temperature.
%    epsilon    O   Relative complex permittivity of water ice.
%
% Detailed Input
%
%    f           is an array of dimension <=2 containing frequencies for
%                which the relative complex permittivity of water ice needs
%                to be computed. Frequency is expressed in Hz. f is
%                required to have the same size as T.
%
%    T           is an array of dimension <=2 containing absolute
%                temperatures for which the relative complex permittivity
%                of water ice needs to be computed. Temperature is
%                expressed in K. T is required to have the same size as f.
%
% Detailed Output
%
%    epsilon     is an array having the same size as f and T, containing
%                the values of the relative complex permittivity of water
%                ice corresponding to values of f and T.
%
% Parameters
%
%    None.
%
% Exceptions
%
%    1) Input arrays need to have the same size for the function to work
%       properly. The size of input arrays is compared. If input arrays
%       have different sizes, the error message 'Matzler:ArraySizeMismatch'
%       is issued.
%
% Files
%
%    None.
%
% Details
%
%    See Matzler (1998) for a detailed explanation of how the formulas used
%    in this function were obtained.
%
% Examples
%
%    1) The following script demonstrates the use of the function to
%       produce a graph identical to Figure 2 in Matzler (1998):
%
%     clc
%     clear all
%     close all
% 
%     %
%     % The following physical constants are used in the computations:
%     %
%     %    c           Velocity of light in vacuo.
% 
%     c = 299792458;
% 
%     %
%     % Arrays of equal size are created to store frequencies and temperatures
% 
%     f = logspace( 6, 12, 100 );
%     T = [ 100, 150, 200, 240, 270 ];
% 
%     [ f, T ] = meshgrid( f, T );
% 
%     %
%     % The relative dielectric permittivity of water ice and the corresponding
%     % absorption coefficient is computed.
% 
%     epsilon = matzler( f, T );
% 
%     gamma_a = -2 .* pi .* imag( epsilon ) .* f ./ ( c .* sqrt( real( epsilon ) ) );
% 
%     %
%     % Data are plotted as in Figure 2 of Matzler, C. 1998. Microwave properties
%     % of ice and snow. In: B. Schmitt et al. (eds), Solar System Ices, 241-257.
% 
%     loglog( f'./1e9, gamma_a')
%     axis square
%     axis( [ 1e-3, 1e3, 1e-8, 1e2 ] )
%     xlabel( 'Frequency (GHz)' )
%     ylabel( 'Absorption coefficient (m^{-1})' )
%
% Restrictions
%
%    The formula provided in Matzler (1998) is valid for any temperature
%    between 0 K and 273 K, and for frequencies greater than 0 and below
%    1000 GHz.
%
% Literature References
%
%    Matzler, C. 1998. Microwave properties of ice and snow. In: B. Schmitt
%    et al. (eds), Solar System Ices, 241-257.
%
% Author and Institution
%
%    R. Orosei   (INAF)
%
% Version
%
%    Version 1.0.0, 25-OCT-2010 (RO)
%
%       First release.
%
% Index Entries
%
%    None at this time.
%
% Revisions
%
%    Version 1.0.0, 25-OCT-2010 (RO)
%
%       Initial release.
%

%
% Input arrays need to have the same size for the function to work
% properly.

[ rowsf, colsf ] = size( f );
[ rowsT, colsT ] = size( T );

if rowsf ~= rowsT || colsf ~= colsT
    error( 'Matzler:ArraySizeMismatch', ...
         [ 'Input arrays f and T must have the same size.\n', ...
           'f has size %d X %d, T has size %d X %d.', ], ...
            rowsf, colsf, rowsT, colsT )
end

%
% Input frequencies must be positive to be physically meaningful.

if any( f ) < 0
    error( 'Matzler:NegativeInputFrequencies', ...
           'Input arrays f must contain only non-negative elements.' )
end

%
% Compute the real part of the relative dielectric permittivity,
% interpolating for temperatures outside the range of available
% measurements.

% epsilonre = 3.1884 + 9.1e-4 .* ( T - 273 );
% epsilonre( epsilonre < 3.1 ) = 3.1;

b = 3.1;
a = ( 3.1884 - 3.1 ) / 273^4;

epsilonre = a .* T.^4 + b;

%
% The formula for the imaginary part of the relative dielectric
% permittivity given in Matzler (1998) is written for frequencies expressed
% in GHz: an auxiliary variable containing frequency in GHz is used.

fGHz = f ./ 1e9;

%
% The imaginary part of the relative dielectric permittivity is computed
% according to equations (2), (3), (A1), (A2) and (A3) in Matzler (1998).

T_0   = 300;
theta = T_0 ./ T - 1;
alpha = ( 0.00504 + 0.0062 .* theta ) .* exp( -22.1 .* theta );

B_1 = 0.0207;
B_2 = 1.16e-11;
b   = 335;

beta_M    = B_1 ./ T .* exp( b ./ T ) ./ ( exp( b ./ T ) - 1 ).^2 + B_2 .* fGHz.^2;
Deltabeta = exp( -10.02 + 0.0364 .* ( T - 273 ) );
beta      = beta_M + Deltabeta;

epsilonim = alpha ./ fGHz + beta .* fGHz;

%
% The relative dielectric permittivity is computed, in the convention that
% the imaginary part of a dielectric medium is negative.

epsilon = epsilonre - 1i .* epsilonim;