% BTEX10 Flow with axial dispersion through a one-region pipe of uniform
% cross-section (ported from Physiome JSim model 079 by DH  9/1/12).
function BTEX10
  close all

  % Parameters
  Fp    = 1/60;   % ml/g/sec
  hVp   = 0.05;   % ml/g
  Gp    = 0.2/60; % ml/g/sec
  Dp    = 1e-5;   % cm^2/sec

  tmin  = 0;      % sec
  tmax  = 30;     % sec
  tdelt = 0.1;    % sec

  xmin  = 0;      % cm
  L     = 0.1;    % cm
  xdelt = L/60;   % cm

  % Spacing in x, t
  xspan = xmin:xdelt:L;
  tspan = tmin:tdelt:tmax;

  % Solver Parameters
  options = ddeset('RelTol', 1e-7, 'AbsTol', 1e-7);

  % Inflowing Concentration
  Cin = LagNormal(5, .3, 1.3, tspan);

  sol = pdepe(0, @pdefun, @icfun, @bcfun, xspan, tspan, options);
  
  figure('name', 'Conc_Spatial', 'NumberTitle', 'off', 'Units', 'normalize', ...
         'Position', [0 0 .5 1])
  subplot(2, 1, 1)
  plot(tspan, Cin, 'k', tspan, sol(:, end), 'r')
  title('Concentration')
  xlabel('Time, sec')
  ylabel('Concentration, mM')
  axis([tmin tmax 0 .4])
  legend('Cin mM', 'Cout mM')

  subplot(2, 1, 2)
  plot(xspan, sol(find(tspan == 3), :), 'k', ...
       xspan, sol(find(tspan == 4), :), 'r', ...
       xspan, sol(find(tspan == 5), :), 'm', ...
       xspan, sol(find(tspan == 6), :), 'y', ...
       xspan, sol(find(tspan == 7), :), 'g', ...
       xspan, sol(find(tspan == 8), :), 'b')
  title('Cp(t,x) Spatial Profiles at time t')
  xlabel('Distance along axis, cm')
  ylabel('Concentration mM at position x')
  legend('Cp(3,x) mM', 'Cp(4,x) mM', 'Cp(5,x) mM', ...
         'Cp(6,x) mM', 'Cp(7,x) mM', 'Cp(8,x) mM')
  
  function [c, f, s] = pdefun(x, t, u, DuDx)
    c = 1;
    f = Dp * DuDx;
    s = -Fp * L / hVp * DuDx - Gp / hVp * u;
  end % pdefun

  function u = icfun(x)
    u = 0;
  end % icfun

  function [pl, ql, pr, qr] = bcfun(xl, ul, xr, ur, t)
    cint = interp1(tspan, Cin, t);
    pl = -Fp * L / hVp * (ul - cint);
    ql = 1;
    pr = 0;
    qr = 1;
  end % bcfun
end % BTEX10

% LagNormal Lagged Normal Density Function.
% pdf = LagNormal(tSpan, tMean, RD, skewn) returns a vector of the same
% size as tSpan with the relative likelihoods that a random variable will
% take the values of a Lagged Normal Probability Density Function. The PDF
% is characterized by the parameters tMean (mean transit time), RD
% (relative dispersion) and skewn(skewness) as described in the reference.
% tSpan should be a monotonously increasing or decreasing vector. For a
% sufficiently long interval that includes the peak, the area under the
% density function approaches 1.
%
% Reference:
% Bassingthwaighte, J.B., F.H. Ackerman and E.H. Wood: Applications of the
% lagged normal density curve as a model for arterial dilution curves.
% Circ. Res. 18:398-415, 1966.
function pdf = LagNormal(tMean, RD, skewn, tspan)
  tau = RD * tMean * (skewn / 2)^(1 / 3);
  sigma = sqrt((RD * tMean) ^ 2 - tau ^ 2);
  mean = tMean - tau;
  
  options = ddeset('RelTol', 1e-7, 'AbsTol', 1e-7);
  sol =  ode23(@dCdt, tspan, 0, options);
  pdf = interp1(sol.x, sol.y, tspan);
    
  function Ret = dCdt(t, C)
    Cin = exp(-.5 * ((t - mean) / sigma) ^ 2) / (sigma * sqrt(2 * pi));
    Ret = (Cin - C) / tau;
  end % dCdt
end % LagNormal