1 function reaction = reaction_term(X, Y, U, params)
3 if isfield(params, 'geometry_transformation')
5 if isequal(params.geometry_transformation, 'spline')
6 % x_hill = params.geometry_transformation_spline_x;
7 % y_hill = params.geometry_transformation_spline_y;
8 % y_hill(2) = params.hill_height;
9 p_mu = spline_select(params);
10 [ breaks, coeffs, pieces, order ] = unmkpp(p_mu);
11 p_mu_d = mkpp(breaks, coeffs(:,1:order-1) .* repmat(order-1:-1:1,pieces,1));
12 [ breaks, coeffs, pieces, order ] = unmkpp(p_mu_d);
14 p_mu_dd = mkpp(breaks, zeros(size(coeffs)));
16 p_mu_dd = mkpp(breaks, coeffs(:,1:order-1) .* repmat(order-1:-1:1,pieces,1));
19 div = 1 + ppval(p_mu, X);
20 rec = 2 .* (ppval(p_mu_d,X) ./ div).^2 - ppval(p_mu_dd, X) ./ div;
22 reaction = 0.001 .* rec .* U;
27 % TO BE ADJUSTED TO NEW SYNTAX