%BFPSP  Parameters for B-spline basis functions
%
% Constructs necessary data for evaluation of B-splines of order 5 with 0
% value at the edge and natural or not-a-knot boundary conditions for the
% first derivative.
% Assumes same knots distribution for p' and TT'.
%
% P = BFPSP(T[,EC]) where T is the knot distribution corresponding to values
% of FN (values between 0 and 1) and EC sets the boundary conditions for P'
% and TT' at FN=0,1. It can be either 'a' for not-a-knot (continuous f''')
% or 'n' for natural (f''=0), the default value is 'a'.
%
% See also: BFSP, CSDEC, BSPSUM, BSPBASE
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

function [bfp] = bfpsp(t,ec)

% arguments
if nargin < 2
	ec = 'a';
end
nt = numel(t);
if nt < 2, error('Give more than one knot.'); end

d = 4; % Degree of the polynomials for the integral bf

% augment knot sequence
dt = mean(diff(t));
taug = [(-d:-1)'*dt+t(1);t(:);t(nt)+(1:d)'*dt];

% compute the first derivative for the augmented base function set
b = bspsum(taug, eye(nt+d-1), t, 1);

% set up constraint matrix for specified edge conditions: T*c = 0
if nt < 4, error('Cases with fewer than 4 knots not implemented yet');
elseif ec == 'n' % f'''(t(1)) = f'''(t(nt)) = 0
 T = bspsum(taug, eye(nt+d-1), t([1,nt]), 3) * dt^3;
elseif ec == 'a' % f'''' equal at some mid-points
 T = [1 -1 0 0; 0 0 1 -1] * bspsum(taug, eye(nt+d-1), ...
     [t(1)+t(2) t(2)+t(3) t(nt-2)+t(nt-1) t(nt-1)+t(nt)]/2, 4) * dt^4;
elseif ec == 'z' % f'(t(nt)) = 0
 T = bspsum(taug, eye(nt+d-1), t([1 nt]), 3) * dt^3;
else
 error('Invalid edge condition.')
end
% Add condition f(t(end)) = 0
T = [bspbase(taug,d+1,t(nt)).';T];

% constraint regression
[q,~] = qr(T'); q = q(:,size(T,1)+1:end);
b = b*q; M = q/b;

bfp.taug = taug;
bfp.M    = M;

