function varargout = echcd(varargin) %#codegen

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

if nargin == 0,
  varargout{1} = RAPTORmodule('echcd_gaussian'); % load default params
  return %empty call, probably to get default structures
elseif nargin == 1
  error('this use is deprecated, use RAPTORmodule(..) instead')
  % get default parameters for each case
  def_params.modeltype = varargin{1};
  switch def_params.modeltype
    case 'echcd_gaussian'
      def_params.params = echcd_gaussian;
    case 'echcd_manual'
      def_params.params = echcd_manual;
    case 'none'
      def_params.params.active = false;
    otherwise
      error('invalid echcd model type ''%s''',def_params.modeltype)
  end
  varargout{1} = def_params; % return appropriate parameter structure
  return
elseif nargin==6
  stap = varargin{1};
  geop = varargin{2};
  u = varargin{3};
  it = varargin{4};
  model  = varargin{5};
  echcd_params = varargin{6}; % distribute inputs
else
  error('must call with 0 or 6 inputs');
end

%% Compute output
jecperunit = zeros(model.rgrid.nrhogauss,model.echcd.n_units);

if echcd_params.active
  switch model.echcd.name
    case 'echcd_gaussian'
      [Pec,jec,dPec_dx,djec_dx,dPec_du,djec_du,jecperunit] = ...
        echcd_gaussian(stap,geop,u,it,model,model.echcd,echcd_params);
    case 'hcd_manual'
      [Pec,~,jec,~,...
        dPec_dx,~,djec_dx,~,...
        dPec_du,~,djec_du,~] = ...
        hcd_manual(stap,geop,u,it,model,model.echcd,echcd_params);
    case 'hcd_external'
      [Pec,~,jec,~,...
        dPec_dx,~,djec_dx,~,...
        dPec_du,~,djec_du,~] = ...
        hcd_external(stap,geop,u,it,model,model.echcd,echcd_params);
    case 'hcd_gaussian'
      [Pec,~,jec,~,dPec_dx,~,djec_dx,~,dPec_du,~,djec_du,~,jecperunit] = ...
        hcd_gaussian(stap,geop,u,it,model,model.echcd,echcd_params);
    otherwise
      error('invalid setting %s for echcd model',model.echcd.name);
  end
else
  Pec = zeros(model.rgrid.nrhogauss,1);
  jec = zeros(model.rgrid.nrhogauss,1);
  dPec_dx = zeros(model.rgrid.nrhogauss,model.dims.nx);
  djec_dx = zeros(model.rgrid.nrhogauss,model.dims.nx);
  dPec_du = zeros(model.rgrid.nrhogauss,model.dims.nu);
  djec_du = zeros(model.rgrid.nrhogauss,model.dims.nu);
end

varargout{1} = Pec;
varargout{2} = jec;
varargout{3} = dPec_dx;
varargout{4} = djec_dx;
varargout{5} = dPec_du;
varargout{6} = djec_du;
varargout{7} = jecperunit;
return
