function [C,z0]=gaussianPlume(Q, u_ref, h, varargin)
%gaussianPlume  Steady-state gaussian plume distribution model
%
%   gaussianPlume models the dispersion of a continuous point source, i.e.
%   plume, in various conditions and terrains. The output of gaussianPlume
%   is a 3-dimensional matrix containing the concentrations of the emitted
%   substance over a field with the first dimension (y) representing the
%   cross-wind axis, the second dimension (x) the downwind distance, and
%   the third dimension (z) the vertical axis. The origin is set at the
%   base of the stack unless Briggs plume rise model is used in which case
%   the horizontal (x axis) coordinate is offset such that x=0 marks the 
%   maximum rise of the plume downwind. All units are in <mgs> (meters,
%   grams, seconds) except where noted.
%
%   Most of the equations were taken from the ISC3 User's Manual, Volume
%   II, available online at
%   http://www.epa.gov/scram001/userg/regmod/isc3v2.pdf
%
%   The EPA now has newer software for modeling atmospheric dispersion
%   including AERMOD and CALPUFF
%   http://www.epa.gov/scram001/dispersion_prefrec.htm#aermod as these are
%   their preferred/recommended models.
%
%   C = gaussianPlume(Q) returns the steady-state Gaussian distribution
%   model of a single, continuous point source emitting at a rate of Q
%   grams per second for a 50m physical stack height with no calculations
%   for plume rise, in rural terrain with stability class "F" in the 
%   Guifford-Pasquiill scale. Wind speed is assumed to be 1m/s at the stack
%   tip (50m).
%
%   C = gaussianPlume(Q, u_ref) sets the wind at 1m to be u_ref
%
%   C = gaussianPlume(Q, u_ref, h) sets the stack height to h
%
%   [C, z0] = gaussianPlume(...) returns the effective stack height in
%       meters as z0.
%
%   C = gaussianPlume(Q, u_ref, h, ...) allows you to set certain options.
%   Options are set by 'option', <option value> pairs. The list of
%   permissible options are listed below:
%
%       {'h_ref', h_ref} Sets the reference height for the wind speed u_ref
%           in meters. Default is 1m.
%       {'stability', 'class'}    Sets Guifford-Pasquill stability class to
%           class 'class'. 'class' must be one of 'A', 'B', 'C', 'D', 'E', 
%           or 'F'. Default is 'F'.
%       {'terrain', 't_type'}   Sets the terrain to one of 'rural' or
%           'urban'. Default is 'urban'
%       {'p', p}    Sets the scaling factor for change in wind as a
%           function of altitude to p. Default is 0.4.
%       {'plum_rise_model', 'model'}  Sets the model to use for calculation
%           of plume rise. 'model' is one of 'none', 'CONCAWE',
%           'CarlsonMoses', 'Holland', or 'Briggs'. Default is 'none'.
%       {'mw', mw}  Sets the molecular weight of the exhaust for plume rise
%           calculations. Specify in atomic mass units. Default is average
%           of air, 30g/mole.
%       {'amb_temp', Ta} Sets the ambient temperature in degrees Celsius.
%           Default is 25C. Alternatively, you may specify {'Ta',Ta}.
%       {'stack_temp', Ts} Sets the stack temperature in degrees Celsius.
%           Default is 200C. Alternatively, you may specify {'Ts',Ts}.
%       {'stack_diameter', ds} Sets the stack diameter (in meters). Default
%           is 10m. Alternatively, you may specify {'ds',ds}.
%       {'specific_heat', Cp} Sets the specific heat of the exhaust gas in
%           J/degree Celsius/g. Default is 1.020 (constant pressure Cp for
%           dry air). Alternatively, you may specify {'Cp',Cp}.
%       {'amb_pres', Pa} Sets ambient pressure in millibars. Default is
%           1010mb. Alternatively, you may specify {'Pa',Pa}.
%       {'stack_pres', Ps} Sets stack tip pressure in millibars. Default is
%           1010mb. Alternatively, you may specify {'Ps',Ps}.
%       {'stack_velocity', vs} Sets the stack exit velocity in m/s. Default
%           is 1m/s. Alternativelt, you may specify {'vs',vs}.
%       {'lapse', eta} Sets the lapse rate for stable conditions (class E or
%           F). Alternatively, you may specify {'eta',eta}.
%       {'reflection', true} enables ground reflection whereas
%       {'reflection', false} disables ground reflection
%       {'deposition', true} models dry deposition whereas
%       {'deposition', false} disables dry deposition modeling
%       {'term_velocity', vt} Sets terminal/settling velocity in m/s.
%           Default is settling velocity for PM2.5, 0.5cm/s. Synonyms
%           include {'settling_velocity', vt} or simply {'vt', vt}.
%       {'X', x} Sets sampling points at the downwind distances. x is a
%           vector specifying the downwind distances to evaluate for C in
%           meters. Default is 1 to 5km in 100m resolution.
%       {'Y', y} Same as above except for the y-axis (cross-wind)
%       {'Z', z} Same as above except for the vertical axis
%
%   Example
%   ------------------------------
%   C=gaussianPlume(500, 5, 50, 'h_ref', 10, 'reflection', false, ...
%       'deposition', false, 'amb_pres', 1103, 'amb_temp', 22);
%       Computes the steady-state Gaussian-distribution of an emission at a
%       rate of 500g/s with average wind speed of 5m/s at 10m from a
%       physical stack height of 50m. No plume rise, reflection, nor
%       deposition will be modeled. The data will be sampled in 100m
%       resolution across the downwind and crosswind axes from 1 to 5km and
%       1 to 1km in the vertical axis.


% Check for missing arguments. Set to empty so we don't get undefined
% errors.
if(nargin<2)
    u_ref=[];
end

if(nargin<3)
    h=[];
end

if(isempty(u_ref))
    % Default wind speed is 1m/s
    u_ref=1;
end

if(isempty(h))
    % Default stack height is 50m
    h=50;
end

% Set default values
stability='F';              % Guifford-Pasquill stability class
plume_rise_model='none';    % Plume rise models: 'none', 'CONCAWE', 'Briggs'
reflection=false;           % Model ground reflection?
deposition=false;           % Model dry deposition?
terrain='urban';            % Terrain type, 'urban' or 'rural'
p=0.4;                      % Wind speed variation as a function of altitude, from 0.07 to 0.6
mw=30;                      % Average molecular weight of the efflux, typically air (g/mol)
Ta=25;                      % Ambient temperature (degrees C)
Ts=200;                     % Stack tip temperature (degrees C)
ds=10;                      % Stack tip diameter (meters)
Cp=1.020;                   % Specific heat of efflux (J/g degrees Kelvin). 1.020 is for air (http://www.efunda.com/materials/common_matl/show_gas.cfm?MatlName=Air0C)
Pa=1010;                    % Ambient pressure (millibars = 100 N/m^2 = 100 [kg*m/s^2]/m^2 = 1e5 [g/m/s^2])
Ps=1010;                    % Stack tip pressure (millibar)
x_grid=[1 100:100:5e3];     % Downwind distance of each sampled point (m)
y_grid=[1 100:100:4e3];     % Crosswind distance of each sampled point (m)
z_grid=[1 100:100:1e3];     % Altitude above ground level of each sampled point (m)
vt=0.005;                   % Settling (terminal) velocity in m/s
vs=1.5;                     % Stack exit velocity (m/s)
h_ref=1;                    % Altitude above ground level of measured windspeed (u_ref) (in m)
eta=0.035;                  % Lapse rate (degrees C/m)

% The rest of the arguments are set by optional 'string' 'value' pairs
for argnum=1:2:length(varargin)
    switch(varargin{argnum})
        case 'stability'
            stability=varargin{argnum+1};            
        case 'plume_rise_model'
            plume_rise_model=varargin{argnum+1};            
        case 'reflection'
            reflection=varargin{argnum+1};        
        case 'deposition'
            deposition=varargin{argnum+1};
        case 'terrain'
            terrain=varargin{argnum+1};        
        case 'p'
            p=varargin{argnum+1};            
        case 'mw'
            mw=varargin{argnum+1};            
        case 'amb_temp'
            Ta=varargin{argnum+1};            
        case 'Ta'
            Ta=varargin{argnum+1};            
        case 'stack_temp'
            Ts=varargin{argnum+1};            
        case 'Ts'
            Ts=varargin{argnum+1};            
        case 'stack_diameter'
            ds=varargin{argnum+1};            
        case 'ds'
            ds=varargin{argnum+1};            
        case 'specific_heat'
            Cp=varargin{argnum+1};            
        case 'Cp'
            Cp=varargin{argnum+1};            
        case 'amb_pres'
            Pa=varargin{argnum+1};
        case 'Pa'
            Pa=varargin{argnum+1};
        case 'stack_pres'
            Ps=varargin{argnum+1};            
        case 'Ps'
            Ps=varargin{argnum+1};            
        case 'X'
            x_grid=varargin{argnum+1};            
        case 'Y'
            y_grid=varargin{argnum+1};            
        case 'Z'
            z_grid=varargin{argnum+1};            
        case 'term_velocity'
            vt=varargin{argnum+1};            
        case 'vt'
            vt=varargin{argnum+1};            
        case 'settling_velocity'
            vt=varargin{argnum+1};            
        case 'stack_velocity'
            vs=varargin{argnum+1};            
        case 'vs'
            vs=varargin{argnum+1};            
        case 'h_ref'
            h_ref=varargin{argnum+1};            
        case 'lapse'
            eta=varargin{argnum+1};            
        case 'eta'
            eta=varargin{argnum+1};            
        otherwise
            warn('gaussianPlume:args', ['Unrecognized option: ', varargin{argnum}]);
    end    
end

% Setup some physical constants
Rc=8.3144e3;  % Universal gas constant in mJ/mol/K = [g*m^2/s^2/mol/K]
gc=9.80616;     % Graviational acceleration, m/s^2

% For internal calculations, convert stack and ambient temperatures 
% into degrees Kelvin
Ta=Ta+273;
Ts=Ts+273;

% Force x_grid, y_grid, and z_grid to be row vectors
if(size(x_grid, 1)==1)
    x_grid=x_grid';
end
if(size(y_grid, 1)==1)
    y_grid=y_grid';
end
if(size(z_grid, 1)==1)
    z_grid=z_grid';
end

% Determine wind velocity (us) at tip of stack using power law for wind
% speed measured at height h_ref (m) to be u_ref (m/s)
us=((h/h_ref)^p)*u_ref; % in m/s

% ISC3 manual warns of issues with us<1m/s
if(us<1)
    warning('gaussianPlume:badus','Stack height wind speed, u_s, is %.3fm/s which is less than 1m/s.', us);
end

% Compute mass flow using ideal gas law assumption
% PV=nRT (ideal gas law) <==> n=PV/RT
% m_dot=[m^2]*[m/s]*[g/m/s^2]*[g/mol]/[(mJ/mol/K)*K]
%      =[m^3][1/s]*[g/m/s^2]*[g/mol]/[g*m^2/s^2/mol]
%      =[g/s]
m_dot=pi*(ds/2)^2*vs*(Ps*1e5)*mw/(Rc*Ts);   % g/s
% Now compute heat flux
% Qh=[g/s]*[J/(g*degK)]*degK
%   =[J/s]
Qh=m_dot*Cp*(Ts-Ta);    % J/s

% Compute plume rise
switch(plume_rise_model)
    case 'none'
        % The rise is then set to zero
        plume_rise=0;
        x_f=0;
    case 'CONCAWE'
        % TODO: Verify
        plume_rise=4.71*(Qh^0.44/us^0.694);
        x_f=0;
    case 'Holland'
        % TODO: Verify
        plume_rise=vs*ds/us*(1.5+0.01*Qh/(vs*ds));
        x_f=0;
    case 'CarlsonMoses'
        % TODO: Verify
        plume_rise=0.029*vs*ds/us+2.62*(Qh^0.5/us);
        x_f=0;
    case 'Briggs'
        % Compute buoyancy and momentum factors
        % Fb=[m/s^2]*[m/s]*[m^2]*[degK]/[degK]
        %   =[m^4]/[s^3]
        Fb=gc*vs*ds^2*(Ts-Ta)/(4*Ts);
        % Fm=[m^2/s^2]*[m^2]*[degK]/[degK]
        %   =[m^4]/[s^2]
        Fm=vs^2*ds^2*Ta/(4*Ts);
        % Compute stability parameter if stable
        if(stability=='E' || stability=='F')
            % s = [m/s^2]*[degK/m]/[degK]
            %   = [1/s^2]
            s=gc*eta/Ta; % 1/s^2
            % Can't perform unit analysis, factor 0.019... is unknown
            dTc=0.019582*Ts*vs*sqrt(s);
            % Check if buoyancy dominated (thermal) or momentum (kinetic)
            if( (Ts-Ta)>=dTc)
                % Buoyancy dominates
                x_f=2.0715*us/sqrt(s);
                plume_rise=2.6*(Fb/(us*s))^(1/3);
            else
                x_f=0;
                % Note that we also evaluate unstable/neutral and select
                % the lower value
                plume_rise=min(1.5*(Fm/(us*sqrt(s)))^(1/3), 3*ds*vs/us);
            end
        else
            % Unstable or neutral              
            if(Fb<55)
                % Check for buoyancy dominated or momentum
                dTc=0.0297*Ts*vs^(1/3)/ds^(2/3);
                if( (Ts-Ta)>=dTc )
                    % Buoyancy dominated
                    x_f=49*Fb^(5/8);
                    plume_rise=21.425*Fb^(3/4)/us;
                else
                    % Momentum dominated
                    x_f=0;
                    plume_rise=3*ds*vs/us;
                    % Check
                    if(vs/us<=4)
                        warning('gaussianPlume:Briggs', 'Momentum rise model selected but may not be accurate due to low stack exit velocity');
                    end
                end
            else
                % Brigg's equations seem to diverage at Fb==55
                dTc=0.00575*Ts*vs^(2/3)/ds^(1/3);
                if( (Ts-Ta)>=dTc )
                    % Offset horizontal downwind distance for distance to
                    % max height (rise)
                    x_f=119*Fb^(2/5);
                    plume_rise=38.71*Fb^(3/5)/us;
                else
                    x_f=0;
                    plume_rise=3*ds*vs/us;
                    if(vs/us<=4)
                        warning('gaussianPlume:Briggs', 'Momentum rise model selected but may not be accurate due to low stack exit velocity');
                    end
                end
            end
        end            
    otherwise
        warning('gaussianPlume:riseModel', ['Unrecognized plume rise model specified: ', plume_rise_model]);
        plume_rise=0;
        x_f=0;
end

% Set effective stack height
z0=h+plume_rise;

% Offset the downwind direction to take into account distance to max rise
% TODO: Substitute appropriate equations for distance less than final rise
% per ISC3v2 manual for buoyancy dominated conditions
x_grid=x_grid-x_f;

% Compute the dispersion coefficients
switch(terrain)
    case 'rural'
        % Pasquill-Gifford curves
        switch(stability)            
            case 'A'
                % [c, d] coefficients
                coeffs_y=[24.1670, 2.5334];
                % [x a b] matrix
                coeffs_z=[0.10 122.800 0.94470;
                          0.16 158.080 1.05420;
                          0.21 170.220 1.09320;
                          0.26 179.520 1.12620;
                          0.31 217.410 1.26440;
                          0.41 258.890 1.40940;
                          0.51 346.750 1.72830;
                          3.11 453.850 2.11660;
                          inf  nan     nan];
            case 'B'
                coeffs_y=[18.3330, 1.8096];
                coeffs_z=[0.20 90.673 0.93198;
                          0.40 98.483 0.98332;
                          inf  109.300 1.09710];
            case 'C'
                coeffs_y=[12.5000, 1.0857];
                coeffs_z=[inf 61.141 0.91465];
            case 'D'
                coeffs_y=[8.3330, 0.72382];
                coeffs_z=[0.31 34.459 0.86974;
                          1.01 32.093 0.81066;
                          3.01 32.093 0.64403;
                          10.01 33.504 0.60486;
                          30.00 36.650 0.56589;
                          inf   44.053 0.51179];
            case 'E'
                coeffs_y=[6.2500, 0.54287];
                coeffs_z=[0.10 24.260 0.83660;
                          0.31 23.331 0.81956;
                          1.01 21.628 0.75660;
                          2.01 21.628 0.63077;
                          4.01 22.534 0.57154;
                          10.01 24.703 0.50527;
                          20.01 26.970 0.46713;
                          40.00 35.420 0.37615;
                          inf 47.618 0.29592];                      
            case 'F'
                coeffs_y=[4.1667, 0.36191];
                coeffs_z=[0.21 15.209 0.81558;
                          0.71 14.457 0.78407;
                          1.01 13.953 0.68465;
                          2.01 13.953 0.63227;
                          3.01 14.823 0.54503;
                          7.01 16.187 0.46490;
                          15.01 17.836 0.41507;
                          30.01 22.651 0.32681;
                          60.00 27.074 0.27436;
                          inf   34.219 0.21716];
            otherwise
                error('gaussianPlume:stability', ['Unknown stability class ', stability]);                
        end
        % Construct sigma_y vector along the x-axis
        % Note that x should be in kilometers (x_grid/1.e3)
        sigma_y=465.11628.*(x_grid./1e3).*tan(0.017453293.*(coeffs_y(1)-coeffs_y(2).*log(x_grid./1e3))); % m
        % Construct sigma_z vector along the x-axis
        prev_boundary=0; % in km    
        % Pre-allocate (should be same size as x_grid since all tables end
        % with 'inf')
        sigma_z=nan(size(x_grid));
        for section=1:size(coeffs_z, 1)
            idx=find(prev_boundary<=(x_grid./1e3) & (x_grid./1e3)<coeffs_z(section, 1));
            sigma_z(idx)=coeffs_z(section, 2).*(x_grid(idx)./1e3).^coeffs_z(section, 3);
            prev_boundary=coeffs_z(section, 1);
        end
        % Move the vector into the proper dimension (cols=x)
        sigma_y=shiftdim(sigma_y, -1);
        % Clip all values>5e3 for stability classes A-C
        switch(stability)
            case {'a', 'b', 'c'}
                sigma_y(sigma_y>5e3)=5e3;
        end        
        %sigma_z=shiftdim(sigma_z, -1);        
    case 'urban'
        % Pasquill-Gifford with urban fit (McElroy-Pooler)
        switch(stability)
            case 'A'
                coeffs_y=0.32;
                coeffs_z=[0.24 1 0.001 0.5];
            case 'B'
                coeffs_y=0.32;
                coeffs_z=[0.24 1 0.001 0.5];
            case 'C'
                coeffs_y=0.22;
                coeffs_z=[0.20 1 0 0];
            case 'D'
                coeffs_y=0.16;
                coeffs_z=[0.14 1 0.0003 -0.5];
            case 'E'
                coeffs_y=0.11;
                coeffs_z=[0.08 1 0.0015 -0.5];
            case 'F'
                coeffs_y=0.11;
                coeffs_z=[0.08 1 0.0015 -0.5];
            otherwise
                error('gaussianPlume:stability', ['Unrecognized stability class ', stability]);
        end
        % Construct sigma_y along x-axis
        sigma_y=coeffs_y(1).*x_grid.*(1+0.0004.*x_grid).^(-0.5);
        sigma_y=shiftdim(sigma_y, -1);
        % Construct sigma_z along x-axis
        sigma_z=coeffs_z(1).*x_grid.*(coeffs_z(2)+coeffs_z(3).*x_grid).^coeffs_z(4);        
        sigma_z=shiftdim(sigma_z, -1);
    otherwise
        error('gaussianPlume:terrain', ['Unrecognized terrain option ', terrain]);        
end

% Now we have a choice when we compute Gaussian distribution. We
% can either save memory and use a for-loop, or speed the code up
% at the expense of more memory requirement. Since memory is cheap
% these days, we'll go ahead and optimize for speed in favor of
% size.
        
% Now replicate sigma vectors into 3D matrices
sigma_y=repmat(sigma_y, [length(y_grid) 1 length(z_grid)]);
sigma_z=repmat(sigma_z, [length(y_grid) 1 length(z_grid)]);

% Create wind matrix (scaled according to altitude)
u_matrix=((z_grid./h_ref).^p).*u_ref; % [m/s]
% Move over into z-dimension
u_matrix=shiftdim(u_matrix, -2);
u_matrix=repmat(u_matrix, [length(y_grid) length(x_grid) 1]);

% At this point, we will scale x_grid/y_grid/z_grid into full 3-D matrices
% Because we now manipulate them, first save their original lengths
sx=size(x_grid, 1);
sy=size(y_grid, 1);
sz=size(z_grid, 1);
x_grid=repmat(shiftdim(x_grid, -1), [sy 1 sz]);
y_grid=repmat(y_grid, [1 sx sz]);
z_grid=repmat(shiftdim(z_grid, -2), [sy sx 1]);

% Setup deposition
if(deposition)
    % Note here we had to massage x_grid into full 3-D matrix
    dz_dep=vt.*x_grid./u_matrix;
else
    dz_dep=zeros(size(x_grid));
end

if(reflection)
    % Reflection scales everything by 2
    r=2;
else
    r=1;
end

% Now compute steady-state Gaussian dispersion
% (Note everything inside the exp() becomes unitless)
% C=[g/s]/[m/s*m*m]=[g/m^3]
C=r.*Q./ ...
    (2.*pi.*u_matrix.*sigma_y.*sigma_z).*exp( ...
        (-y_grid.^2)./(2.*sigma_y.^2) - ((z_grid-z0-dz_dep).^2)./(2.*sigma_z.^2));
