function varargout = lhhcd(varargin) %#codegen

%% DEFAULT PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Input processing

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

%% Compute output
if lhhcd_params.active
    switch model.lhhcd.name
            case 'hcd_manual'
            [Plh,~,jlh,~,...
              dPlh_dx,~,djlh_dx,~,dPlh_du,~,djlh_du] = hcd_manual(stap,geop,u,it,model,model.lhhcd,lhhcd_params);
            case 'hcd_external'
            [Plh,~,jlh,~,...
              dPlh_dx,~,djlh_dx,~,...
              dPlh_du,~,djlh_du,~] = hcd_external(stap,geop,u,it,model,model.lhhcd,lhhcd_params);
            case 'hcd_gaussian'
            [Plh,~,jlh,~,...
              dPlh_dx,~,djlh_dx,~,...
              dPlh_du,~,djlh_du,~] = hcd_gaussian(stap,geop,u,it,model,model.lhhcd,lhhcd_params);
          otherwise
            error('invalid setting for lhhcd model');
    end
    
else
    Plh = zeros(model.rgrid.nrhogauss,1);
    jlh = zeros(model.rgrid.nrhogauss,1);
    dPlh_dx = zeros(model.rgrid.nrhogauss,model.dims.nx);
    djlh_dx = zeros(model.rgrid.nrhogauss,model.dims.nx);
    dPlh_du = zeros(model.rgrid.nrhogauss,model.dims.nu);
    djlh_du = zeros(model.rgrid.nrhogauss,model.dims.nu);
end

varargout{1} = Plh;
varargout{2} = jlh;
varargout{3} = dPlh_dx;
varargout{4} = djlh_dx;
varargout{5} = dPlh_du;
varargout{6} = djlh_du;
return