function varargout = hcd_external(varargin)

switch nargin
    case 0
        % call to get default config
        config = get_default_config;
        varargout{1} = 'hcd'; % type
        varargout{2} = config;
        return
    case 2
        % call to get model and params for given config
        moduleconfig = varargin{1};
        model = varargin{2};
        
        [module,params] = setup_module(moduleconfig,model);
        varargout{1} = module;
        varargout{2} = params;
        return
end
nexp = 7; % expected input arguments
assert(nargin==nexp,'RAPTOR:WrongModuleNargin',...
    'Wrong number of input arguments for this function, must be 0,1 or %d',nexp)

% stap = varargin{1};
% geop = varargin{2};
u = varargin{3};
% it = varargin{4};
model  = varargin{5};
hcd_model = varargin{6};
hcd_params = varargin{7}; % distribute inputs

% indices
uind = hcd_model.uind;
uind_pi = uind(hcd_model.ind_pi);
uind_pe = uind(hcd_model.ind_pe);
uind_j = uind(hcd_model.ind_j);
uind_pfast = uind(hcd_model.ind_pfast);

% init matrices
dpe_dx = zeros(model.rgrid.nrhogauss,model.dims.nx);
dpi_dx = zeros(model.rgrid.nrhogauss,model.dims.nx);
dj_dx = zeros(model.rgrid.nrhogauss,model.dims.nx);
dpfast_dx = zeros(model.rgrid.nrhogauss,model.dims.nx);

dpe_du = zeros(model.rgrid.nrhogauss,model.dims.nu);
dpi_du = zeros(model.rgrid.nrhogauss,model.dims.nu);
dj_du = zeros(model.rgrid.nrhogauss,model.dims.nu);
dpfast_du = zeros(model.rgrid.nrhogauss,model.dims.nu);

% powers and jacobians
dpe_du(:,uind_pe) = hcd_model.Gam;
dpi_du(:,uind_pi) = hcd_model.Gam;
dj_du(:,uind_j) = hcd_model.Gam;
dpfast_du(:,uind_pfast) = hcd_model.Gam;

pe = hcd_model.Gam*u(uind_pe);
pi = hcd_model.Gam*u(uind_pi);
j = hcd_model.Gam*u(uind_j);
pfast = hcd_model.Gam*u(uind_pfast);

varargout{1} = pe;
varargout{2} = pi;
varargout{3} = j;
varargout{4} = pfast;

varargout{5} = dpe_dx;
varargout{6} = dpi_dx;
varargout{7} = dj_dx;
varargout{8} = dpfast_dx;

varargout{9} = dpe_du;
varargout{10} = dpi_du;
varargout{11} = dj_du;
varargout{12} = dpfast_du; % one unit
end

function config = get_default_config
config = struct('rhoknots',linspace(0,1,11)',... % rho grid for spline knots
    'rho_in',linspace(0,1,21)',... % rho grid for specifying input
    'order',3);
end

function [mm,pp] = setup_module(myconfig,model)

%% settings to model
mm.name = mfilename();
mm.type = 'echcd';
mm.n_units = 1;
% spline matrices
mm.rho_in = myconfig.rho_in; % basis for input
mm.Gam_in = splineval(myconfig.rho_in,myconfig.rhoknots,myconfig.order)';
mm.Gam = splineval(model.rgrid.rhogauss,myconfig.rhoknots,myconfig.order)';

assert(diff(size(mm.Gam_in))<=0,...
  'insufficient rho_in points or too many rhoknots specified. Mapping matrix size: [%d,%d]\n',size(mm.Gam_in));
% indices
nindices     = uint32( size(mm.Gam,2) );
mm.uind      = uint32( 1:(4*nindices) );
mm.ind_pe    = uint32(   1:nindices   );
mm.ind_pi    = uint32( nindices   + (1:nindices) );
mm.ind_j     = uint32( 2*nindices + (1:nindices) );
mm.ind_pfast = uint32( 3*nindices + (1:nindices) );

%% settings to params
pp.active = true;
end