function [varargout] = hcd_gaussian(varargin) %#codegen
% RAPTOR auxiliary heating/current drive module
% Generic module for gaussians hcd, fixed current drive efficiency per watt
%
% rdep, wdep and cd_eff set through parameters
coder.extrinsic('warning');

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

if nargin == 0
  % call for default config
  config = struct('n_units',2,...
    'n_rdep_external',0);
  varargout{1} = 'hcd'; % type
  varargout{2} = config;
  return
elseif nargin==2
  myconfig = varargin{1};
  
  [mm,pp] = setup_module(myconfig);
  varargout{1} = mm;
  varargout{2} = pp;
  return %empty call, probably to get default structures
end

%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

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

% init to right size
nrhogauss = numel(model.rgrid.rhogauss);
pfast = zeros(nrhogauss,1);
dpe_dx = zeros(nrhogauss,model.dims.nx);
dpi_dx = zeros(nrhogauss,model.dims.nx);
djcd_dx = zeros(nrhogauss,model.dims.nx);
dpfast_dx = zeros(nrhogauss,model.dims.nx);
dpe_du = zeros(nrhogauss,model.dims.nu);
dpi_du = zeros(nrhogauss,model.dims.nu);
djcd_du = zeros(nrhogauss,model.dims.nu);
dpfast_du = zeros(nrhogauss,model.dims.nu);

iupow = hcd_model.uind(hcd_model.ind_power); % indices (in u) of u for POWER
iurdep = hcd_model.uind(hcd_model.ind_rhodep); % indices (in u) of u for POWER
nunits = hcd_model.n_units;

rdep = zeros(1,nunits); % init

irdep = hcd_model.rdep_ext; % indices of rdep taken externally
rdep(~irdep) = hcd_params.rdep(~irdep); % some deposition locations are directly passed as fixed parameters
if any(irdep)
  rdep(irdep) = u(iurdep); % other deposition locations are time varying, taken from u
end

wdep_in = hcd_params.wdep_in;
cd_eff = hcd_params.cd_eff;
wdep_out = hcd_params.wdep_out;

rr = model.rgrid.rhogauss;

% Heating
% gaussian profiles
Pprofile = zeros(nrhogauss,numel(rdep));

if isempty(wdep_out) || all(wdep_out == wdep_in); % symmetric profile
  Pprofile = exp(-bsxfun(@times,bsxfun(@minus,rr,rdep).^2,1./wdep_in.^2));
else % asymmetric profile
  for ii=1:numel(rdep);
    iout = (rr>rdep(ii));
    Pprofile(~iout,ii)  = Pprofile(~iout,ii) + exp(-4*bsxfun(@times,bsxfun(@minus,rr(~iout),rdep(ii)).^2,1./wdep_in(ii).^2)); % -4 removed
    Pprofile(iout,ii)   = Pprofile( iout,ii) + exp(-4*bsxfun(@times,bsxfun(@minus,rr( iout),rdep(ii)).^2,1./wdep_out(ii).^2));
  end
end

% normalize to unit total power
Ptot = int_Vtot(Pprofile,geop,model);
dp_dupow = hcd_params.frac_absorbed * bsxfun(@times,Pprofile,1./Ptot);

% Current drive
% NB efficiency model
jcdeff = hcd_params.eta_jcd0;

% take into account CD efficiency per actuator
djcd_dupow = bsxfun(@times,jcdeff*cd_eff,dp_dupow);
jcdperunit = bsxfun(@times,djcd_dupow,u(iupow)');
jcd = sum(jcdperunit,2);

% sum beams to get total profile
efrac = hcd_params.frac_to_electrons;

dpe_dupow = bsxfun(@times,dp_dupow,efrac);
dpi_dupow = bsxfun(@times,dp_dupow,(1-efrac));

pe_perunit = bsxfun(@times,dpe_dupow,u(iupow)');
pi_perunit = bsxfun(@times,dpi_dupow,u(iupow)');

pe = sum(pe_perunit,2);
pi = sum(pi_perunit,2);

if nargout>2
  
  if ~isempty(iurdep) % if some inputs are rdep
    disp('this term needs fixing following efrac introduction')
    % Pnb = Pnbprofile/Ptot : use chain rule
    % dPnb_durdep = dPnbprofile_durdep / Ptot
    %              - Pnbprofile/Ptot^2*(dPtot/durdep)
    % Ptot = int(V' Pnbprofile) dr
    %  -> dPtot/durdep = int(V'dPnbprofile_durdep) dr
    
    dPnbprofile_durdep = (Pprofile(:,irdep) .* (8*bsxfun(@times,bsxfun(@minus,rr,rdep(irdep)),1./wdep_in(irdep).^2) ));
    dPnb_durdep = bsxfun(@times,...
      bsxfun(@times,dPnbprofile_durdep,1./Ptot(irdep)) ...
      - bsxfun(@times,Pprofile(:,irdep),int_Vtot(dPnbprofile_durdep,geop,model)./Ptot(irdep).^2),...
      u(iupow(irdep))');
    
    dPe_durdep = efrac*dPnb_durdep;
    dPi_durdep = (1-efrac)*dPnb_durdep;
  else
    dPe_durdep = []; % init
    dPi_durdep = []; % init
  end
  
  dpe_du(:,[iupow,iurdep]) = [dpe_dupow,dPe_durdep]; % assign to complete dPnb_du
  dpi_du(:,[iupow,iurdep]) = [dpi_dupow,dPi_durdep]; % assign to complete dPnb_du
  
  % cd
  if ~isempty(iurdep) % if some inputs are rdep
    djnb_durdep = bsxfun(@times,jcdeff*cd_eff(irdep),dPnb_durdep);
  else
    djnb_durdep = [];
  end
  djcd_du(:,[iupow,iurdep]) = [djcd_dupow,djnb_durdep]; % assign to complete dPnb_du
  
  % sensitivity
  djcd_dx(:,model.te.xind) = zeros(nrhogauss,numel(model.te.xind));
  
  % checks
  if any(u(iupow)<0);
    warning('negative power requested');
    % return Nans in this case, these will propagate and stop the iteration
  end
  %
end
%% assign outputs
varargout{1} = pe;
varargout{2} = pi;
varargout{3} = jcd;
varargout{4} = pfast;
varargout{5} = dpe_dx;
varargout{6} = dpi_dx;
varargout{7} = djcd_dx;
varargout{8} = dpfast_dx;
varargout{9} = dpe_du;
varargout{10} = dpi_du;
varargout{11} = djcd_du;
varargout{12} = dpfast_du;
varargout{13} = jcdperunit;

end


function [mm,pp] = setup_module(myconfig)
n_units = myconfig.n_units;
n_rdep = myconfig.n_rdep_external;

%% Model
mm.name = mfilename;
mm.type = 'any';
mm.n_units = uint32(n_units);
mm.rdep_ext = false(1,n_units); mm.rdep_ext((end-n_rdep+1):end)=true;
mm.ind_power = uint32(1:n_units); % indices of uind for power
mm.ind_rhodep = uint32(n_units + (1:n_rdep)); % indices of uind for rhodep
mm.uind = uint32(1:(n_units + n_rdep)); % u indices

%% Parameters
pp = struct(...
  'active',true,... % activation switch
  'rdep', linspace(0,0.4,n_units), ...  % centers of deposition gaussian, set to -1 if passed through input
  'wdep_in', 0.1*ones(1,n_units), ... % full width of gaussian for part towards center (to have asymmetric distribution)
  'wdep_out', 0.4*ones(1,n_units), ... % full width of gaussian for part towards edge (to have asymmetric distribution)
  'cd_eff', zeros(1,n_units), ... %  current drive effect
  'eta_jcd0', 1.,... % Crrent drive efficiency scaling in A/m^2 / W
  'frac_to_electrons', 1., ... % fraction of power to electrons
  'frac_absorbed',1. ... % fraction of power absorbed
  );
end
