function varargout = nbhcd(varargin) %#codegen

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

if nargin == 0
  %% DEFAULT PARAMETERS
  varargout{1} = RAPTORmodule('hcd_gaussian');
  return %empty call, probably to get default structures
elseif nargin == 1
  error('this use is deprecated, use RAPTORmodule() instead')
elseif nargin==6
  stap = varargin{1};
  geop = varargin{2};
  u = varargin{3};
  it = varargin{4};
  model  = varargin{5};
  nbhcd_params = varargin{6}; % distribute inputs
else
  error('must call with 0 or 6 inputs');
end

%% Compute output
Pnbe = zeros(model.rgrid.nrhogauss,1);
Pnbi = zeros(model.rgrid.nrhogauss,1);
jnb = zeros(model.rgrid.nrhogauss,1);
pfast = zeros(model.rgrid.nrhogauss,1);
dPnbe_dx = zeros(model.rgrid.nrhogauss,model.dims.nx);
dPnbi_dx = zeros(model.rgrid.nrhogauss,model.dims.nx);
dpfast_dx = zeros(model.rgrid.nrhogauss,model.dims.nx);
djnb_dx = zeros(model.rgrid.nrhogauss,model.dims.nx);
dPnbe_du = zeros(model.rgrid.nrhogauss,model.dims.nu);
dPnbi_du = zeros(model.rgrid.nrhogauss,model.dims.nu);
djnb_du = zeros(model.rgrid.nrhogauss,model.dims.nu);
dpfast_du = zeros(model.rgrid.nrhogauss,model.dims.nu);

if nbhcd_params.active
  switch model.nbhcd.name
    case 'hcd_manual'
      [Pnbe,Pnbi,jnb,pfast,...
        dPnbe_dx,dPnbi_dx,djnb_dx,dpfast_dx,...
        dPnbe_du,dPnbi_du,djnb_du,dpfast_du] = hcd_manual(stap,geop,u,it,model,model.nbhcd,nbhcd_params);
    case 'hcd_external'
      [Pnbe,Pnbi,jnb,pfast,...
        dPnbe_dx,dPnbi_dx,djnb_dx,dpfast_dx,...
        dPnbe_du,dPnbi_du,djnb_du,dpfast_du] = hcd_external(stap,geop,u,it,model,model.nbhcd,nbhcd_params);
    case 'hcd_gaussian'
      [Pnbe,Pnbi,jnb,pfast,...
        dPnbe_dx,dPnbi_dx,djnb_dx,dpfast_dx,...
        dPnbe_du,dPnbi_du,djnb_du,dpfast_du] = hcd_gaussian(stap,geop,u,it,model,model.nbhcd,nbhcd_params);
    otherwise
      error('invalid setting for nbhcd model');
  end
end

varargout{1} = Pnbe;
varargout{2} = Pnbi;
varargout{3} = jnb;
varargout{4} = pfast;
varargout{5} = dPnbe_dx;
varargout{6} = dPnbi_dx;
varargout{7} = djnb_dx;
varargout{8} = dpfast_dx;
varargout{9} = dPnbe_du;
varargout{10} = dPnbi_du;
varargout{11} = djnb_du;
varargout{12} = dpfast_du;

return
