/ComputeControl.m

 1function [ut] = ComputeControl(ct, s_int, kd)
 2    % Integrate controls using the Trapezoidal rule and approxmiate derivatives.
 3    % kd is a function handle to a function that calculates kd - kdsquared
 4    % for delta^2 CCM, kdfourth for delta^4 CCM.
 5    
 6    ut = zeros(1/s_int + 1,1);
 7    for s = s_int:s_int:1
 8        e = round(s/s_int + 1); % end index for current s
 9        % calculate end derivatives using skewed stencils
10        derivative0 = (ct(2) - ct(1))/s_int;
11        derivativen = (ct(e) - ct(e - 1))/s_int;
12        ut(e) = 0;
13        ut(e) = ut(e) + kd(derivative0);
14        ut(e) = ut(e) + kd(derivativen);
15        for i = 2:e-1
16            % calculate derivative using centered stencil
17            derivative = (ct(i + 1) - ct(i - 1))/2/s_int;
18            ut(e) = ut(e) + 2*kd(derivative);
19        end
20        ut(e) = s_int/2 * ut(e);
21    end
22end