function [Tb,theta] = scat_layer(height,f,T,ka,ks,Tbo,z)

% function [Tb,theta] = scat_layer(height,f,T,ka,ks,Tbo,z)
% This function computes the brightness temperature 
% of a volume with a uniform temperature 
% whose characteristics
% are defined by absorption and scattering coefficients
% supplied by the user.  
% Isotropic scattering is assumed.  
% The numerical solution used by this function
% requires the user to set boundary conditions that define: 
% the upwelling brightness temperature
% at three angles at the bottom of the scattering layer;
% and the downwelling brightness temperature
% at three angles at the top of the scattering layer.
% The angles correspond to points used in Gaussian quadrature.
% These angles are approximately 
% 21, 49, 76 degrees (upwelling)
% and 104, 131, and 159 degrees (downwelling)
% with 0 degrees defined as straight up (+ z_hat direction).
% ----- inputs:
% height = thickness of scattering layer, m.  
% f = frequency, Hz.
% T = temperature of the scattering layer, K.  
% ka = absorption coefficient, Np m^-1.
% ks = scattering coefficient, Np m^-1.
% Tbo = boundary condition brightness temperatures, K.
%   Must be a six-element column vector. 
%   The first three elements 
%   impose the upwelling brightness temperatures at z = 0
%   at angles from zenith of 76, 49, and 21 degrees (in that order).
%   The second three elements 
%   impose the downwelling brightness temperatures at z = height
%   at angles from nadir of 14, 41, and 69 degrees (in that order).
% z = height at which Tb is desired, m.  Must be either 0 or height.
% ------ outputs:  
% Tb = brightness temperatures at z, K.    
%   Three of the elements of this six-element vector
%   will be the boundary conditions as determined in Tb for z
%   and the other three will be the brightness temperatures
%   computed given the boundary conditions.  
% theta = zenith angles of brightness temperatures, deg.
%   Defined such that theta = 0 is the +z_hat direction, deg.
% -------
% Modified for uniform temperature 
% and user-specified absorption and scattering coefficients from rain.m.
% Modified April 7, 2010, to accomodate Tbo row vector instead of column vector.

c = 2.9979e8;				% free-space velocity
ko = 2.*pi.*f./c;			% free-space wavenumber

ke = ka + ks;               % extinction coefficient, Np/m
omega = ks./ke

mu1 = 0.238619186083197;		% abscissas and weight factors
w1 = 0.467913934572691;			% for Gaussian integration n = 6
mu2 = 0.661209386466265;
w2 = 0.360761573048139;
mu3 = 0.932469514203152;
w3 = 0.171324492379170;
mu4 = -0.238619186083197;
w4 = 0.467913934572691;
mu5 = -0.661209386466265;
w5 = 0.360761573048139;
mu6 = -0.932469514203152;
w6 = 0.171324492379170;

theta = (acos([mu1 mu2 mu3 mu4 mu5 mu6])).*180./pi;

%----------------- system of 6 first order linear equations -----------------
%-------------------------- homogeneous system ------------------------------
%------------------------------ Tb' = P Tb ----------------------------------

P1 = ke./mu1.*[omega/2*w1-1, omega/2*w2, omega/2*w3, omega/2*w4, omega/2*w5, omega/2*w6];
P2 = ke./mu2.*[omega/2*w1, omega/2*w2-1, omega/2*w3, omega/2*w4, omega/2*w5, omega/2*w6];
P3 = ke./mu3.*[omega/2*w1, omega/2*w2, omega/2*w3-1, omega/2*w4, omega/2*w5, omega/2*w6];
P4 = ke./mu4.*[omega/2*w1, omega/2*w2, omega/2*w3, omega/2*w4-1, omega/2*w5, omega/2*w6];
P5 = ke./mu5.*[omega/2*w1, omega/2*w2, omega/2*w3, omega/2*w4, omega/2*w5-1, omega/2*w6];
P6 = ke./mu6.*[omega/2*w1, omega/2*w2, omega/2*w3, omega/2*w4, omega/2*w5, omega/2*w6-1];
P = [P1;P2;P3;P4;P5;P6];

[V,D] = eig(P);				% columns of V (6x6) contain eigenvectors
					% diagonal of D (6x6) contains eigenvalues

%-------------------------- particular solution ------------------------------

a = ke.*(1 - omega).*T;				               % g = ae
e = [1./mu1 1./mu2 1./mu3 1./mu4 1./mu5 1./mu6]';

r = -a.*inv(P)*e;  % Tbparticular = r

%---------------------- finding the final solution -------------------------

for h = 1:6
    for k = 1:6
        if (h == 1) | (h == 2) | (h == 3)
            F(h,k) = V(h,k);
        else					% constructing the matrix F
            F(h,k) = V(h,k).*exp(height.*D(k,k));
        end
    end
end

% zo = [0 0 0 height height height]';		% initial conditions
% Tbo = [290 290 290 60 60 60]';

d = size(Tbo);          % find dimensions of Tbo

if d(2) > d(1)          % if Tbo is a row vector instead of column vector
    Tbo = Tbo';         % change to column vector
end;

c = inv(F)*(Tbo - r);		% vector c for these boundary conditions

%-------------------------- the final solution -------------------------

for s = 1:6
    q(s,1) = c(s,1).*exp(D(s,s).*z);			% the vector q(z)
end

Tb = V*q + r;