function  [varargout] = hcd_manual(varargin)
% RAPTOR aux heating/current drive module for manual input
% No physics, Paux and jaux just set from outside.
% Useful for benchmarking and testing.
% Factor can be provided through u if --uindices-- field provided in config

coder.extrinsic('warning');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Input processing

if nargin == 0
  %% default config
  varargout{1} = 'hcd'; % type
  varargout{2} = []; % no config
  return
elseif nargin==2
  moduleconfig = varargin{1};
  RAPTORmodel = varargin{2};
  
  mm.name = mfilename;
  mm.type = 'any';
  mm.n_units = 1;
  pp.active = true;
  pp.Pauxe = zeros(RAPTORmodel.rgrid.nrhogauss,1);
  pp.Pauxi = zeros(RAPTORmodel.rgrid.nrhogauss,1);
  pp.jaux  = zeros(RAPTORmodel.rgrid.nrhogauss,1);
  pp.pfast = zeros(RAPTORmodel.rgrid.nrhogauss,1);
  if isfield(moduleconfig,'uindices')
    mm.uind = moduleconfig.uindices;
    pp.uindices = moduleconfig.uindices;
  else 
    mm.uind = [];
  end
  
  varargout{1} = mm;
  varargout{2} = pp;
  return
end

% call with 5 inputs, during time loop
% check_params(nbhcd_params_default,nbhcd_params)
%stap = varargin{1};
%geop = varargin{2};
u = varargin{3};
it = varargin{4};
RAPTORmodel = varargin{5};
module_model  = varargin{6};
module_params = varargin{7}; % distribute inputs

if module_params.Pauxe == 0 % fill defaults
  module_params.Pauxe = zeros(RAPTORmodel.rgrid.nrhogauss,1);
end
if module_params.Pauxi == 0 % fill defaults
  module_params.Pauxi = zeros(size(module_params.Pauxe));
end
if module_params.jaux == 0 % fill defaults
  module_params.jaux = zeros(size(module_params.Pauxe));
end
if module_params.pfast == 0 % fill defaults
  module_params.pfast = zeros(size(module_params.Pauxe));
end


if it==1 % stuff to do on first call only
  % other checks here
  assert(all(size(module_params.Pauxe) == size(module_params.jaux)),...
    'Pauxe and jaux must have the same size');
  assert(all(size(module_params.Pauxi) == size(module_params.jaux)),...
    'Pauxi and jaux must have the same size');
  assert(all(size(module_params.pfast) == size(module_params.jaux)),...
    'pfast and jaux must have the same size');
  assert(size(module_params.Pauxe,1) == RAPTORmodel.rgrid.nrhogauss,...
    'Number of rows in Pauxe must correspond to nrhogauss');
  assert(size(module_params.Pauxi,1) == RAPTORmodel.rgrid.nrhogauss,...
    'Number of rows in Pauxe must correspond to nrhogauss');
  assert(size(module_params.jaux,1) == RAPTORmodel.rgrid.nrhogauss,...
    'Number of rows in jaux must correspond to nrhogauss');
  assert(size(module_params.pfast,1) == RAPTORmodel.rgrid.nrhogauss,...
    'Number of rows in pfast must correspond to nrhogauss');
end
%for all other iterations pass directly


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

% init to right size
dPauxe_dx = zeros(RAPTORmodel.rgrid.nrhogauss,RAPTORmodel.dims.nx);
dPauxi_dx = zeros(RAPTORmodel.rgrid.nrhogauss,RAPTORmodel.dims.nx);
djaux_dx = zeros(RAPTORmodel.rgrid.nrhogauss,RAPTORmodel.dims.nx);
dpfast_dx = zeros(RAPTORmodel.rgrid.nrhogauss,RAPTORmodel.dims.nx);

dPauxe_du = zeros(RAPTORmodel.rgrid.nrhogauss,RAPTORmodel.dims.nu);
dPauxi_du = zeros(RAPTORmodel.rgrid.nrhogauss,RAPTORmodel.dims.nu);
djaux_du = zeros(RAPTORmodel.rgrid.nrhogauss,RAPTORmodel.dims.nu);
dpfast_du = zeros(RAPTORmodel.rgrid.nrhogauss,RAPTORmodel.dims.nu);

nt = size(module_params.jaux,2);
flag_ufact = ~isempty(module_model.uind);
if flag_ufact 
	ufact = u(module_params.uindices,:);
else
	ufact = ones(1,nt);
end
if nt == 1
  Pauxe = ufact*module_params.Pauxe;
  Pauxi = ufact*module_params.Pauxi;
  jaux = ufact*module_params.jaux;
  pfast = ufact*module_params.pfast;
  if flag_ufact
    dPauxe_du(:,module_params.uindices) = module_params.Pauxe;
    dPauxi_du(:,module_params.uindices) = module_params.Pauxi;
    djaux_du(:,module_params.uindices) = module_params.jaux;
    dpfast_du(:,module_params.uindices) = module_params.pfast;
  end
elseif it<=nt
  Pauxe = ufact(it)*module_params.Pauxe(:,it);
  Pauxi = ufact(it)*module_params.Pauxi(:,it);
  jaux = ufact(it)*module_params.jaux(:,it);
  pfast = ufact(it)*module_params.pfast(:,it);
  if flag_ufact
    dPauxe_du(:,module_params.uindices) = module_params.Pauxe(:,it);
    dPauxi_du(:,module_params.uindices) = module_params.Pauxi(:,it);
    djaux_du(:,module_params.uindices) = module_params.jaux(:,it);
    dpfast_du(:,module_params.uindices) = module_params.pfast(:,it);
  end
else
  error('Index in Paux and jaux exceeds matrix dimensions i=%d',it)
end

% final checks
if any(Pauxe<0) || any(Pauxi<0)
  warning('negative power requested');
  % return Nans in this case, these will propagate and stop the iteration
  Pauxe = NaN*Pauxe;
  Pauxi = NaN*Pauxi;
end

varargout = {Pauxe,Pauxi,jaux,pfast,...
  dPauxe_dx,dPauxi_dx,djaux_dx,pfast,...
  dPauxe_du,dPauxi_du,djaux_du,pfast};

return

