function  [varargout] = RAPTOR_cvp(varargin)
%
% cvp_params = RAPTOR_cvp % no inputs to get default parameter structure
% [duk_dp] = RAPTOR_cvp(it) % 1 input to get instantaneous duk_dp
% [U,dU_dp,dU_dt,nnp] = RAPTOR_cvp(p0,params) % 2 inputs to evaluate U,dU_dp,dU_dt for all times

persistent module_params
persistent dU_dp
%% DEFAULT PARAMETERS

module_params = struct(...
    'pscale'  ,  1e6 * ones(3,1) , ...
    'psporder',  [1 0 0] , ...
    'tknots'  ,  3 , ...
    'uind'    ,  [], ...
    'transp_matrix',[]...
);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Input processing

if nargin == 0,
    varargout{1} = module_params;
    return %empty call, probably to get default structures
elseif nargin==1; % change this depending on number of inputs to module
    [it] = deal(varargin{:}); % distribute inputs
elseif nargin==2; % change this depending on number of inputs to module
    [p0,params] = deal(varargin{:}); % distribute inputs
    
    % check input parameters
    module_params_in = params.cvp; % module parameter structure
    module_params = local_procinput(module_params, {module_params_in} ); %convert list to structure, replace default by user-defined values.
else
    error('must call with 0, 1 or 2 inputs');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Do the actual work here

if nargin==2 
    % in this case, it is called to generate the full u time evolution 
    % (before the RAPTOR loop) and initialize dU_dp
    par = module_params;
   
    % check inputs
    if numel(par.pscale) ~= numel(par.psporder); error('cvp.pscale must match cvp.psporder'); end
    if ~isvector(p0) || size(p0,1) ~=numel(p0); error('p0 must be a column vector'); end
    
    % knot grid for splines
    if numel(par.tknots)==1 % number of knots specified a a scalar
        tknots = linspace(params.tgrid(1),params.tgrid(end),par.tknots)';
    else
        % knot location explicitly given
        tknots = par.tknots;
    end

    % control vector parametrization
    [dU_dp,dU_dt,nnp] = cvp_setup(params.tgrid,tknots,par.psporder,par.pscale);
    U = zeros(numel(par.pscale),numel(params.tgrid));
    for it = 1:numel(params.tgrid) 
        U(:,it) = dU_dp(:,:,it)*p0;
    end % calculate U for this parameter vector p0
    
    % npp is an array giving the number of parameters parametrizing each input
    np = sum(nnp); % total number of parameters

    varargout = {U,dU_dp,dU_dt,nnp}; 

elseif nargin==1 % This is the call from inside the RAPTOR time loop, only return the duk_dt for the appropriate time
    if isempty(dU_dp); error('RAPTOR_cvp called without initializing it first! Run RAPTOR_cvp(p0,params) first!'); end
    %varargout = {dU_dp_free(:,:,it)}; % duk_dp
    varargout = {dU_dp(:,:,it)}; % duk_dp
    return
end
    


