reaction.m

function reaction = reaction(model, X, Y, U)
%function reaction = reaction(model, X, Y, U)
% 
% function computing the reaction term of a advection-diffusion-reaction
% problem, e.g. in case of a given geometry_transformation.


if isequal(model.geometry_transformation, 'spline')
  p_mu    = spline_select(model);
  [ breaks, coeffs, pieces, order ] = unmkpp(p_mu);
  p_mu_d  = mkpp(breaks, coeffs(:,1:order-1) .* repmat((order-1:-1:1),pieces,1));
  [ breaks, coeffs, pieces, order ] = unmkpp(p_mu_d);
  if order == 1
    p_mu_dd = mkpp(breaks, zeros(size(coeffs)));
  else
    p_mu_dd = mkpp(breaks, coeffs(:,1:order-1) .* repmat((order-1:-1:1),pieces,1));
  end

  div = 1 + ppval(p_mu, X);
  rec = 2 .* (ppval(p_mu_d,X) ./ div).^2 - ppval(p_mu_dd, X) ./ div;

  reaction = model.k .* rec .* U;
end

if isequal(model.name_flux, 'transport')
  reaction = model.transport_source .* U;
end;

%| \docupdate