function  [varargout] = echcd_gaussian(varargin) %#codegen
% RAPTOR auxiliary heating/current drive module for EC heating & current drive
% ECRH/ECCD modeled as gaussians, current drive efficiency dependent on
% local Te/ne
%
% rdep, wdep and cd_eff set through parameters
% for variable rdep, set rdep to -1 and pass rdep through input vector
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 % call for setup
  myconfig = varargin{1};
  model = varargin{2};
  [mm,pp] = setup_module(myconfig,model);
  varargout{1} = mm;
  varargout{2} = pp;
  return
end

%% Run-time call
stap = varargin{1};
geop = varargin{2};
u = varargin{3};
%it = varargin{4};
model  = varargin{5};
echcd_model = varargin{6};
echcd_params = varargin{7}; % distribute inputs

if ~isfield(echcd_params,'basis_fcts')
  basis_fcts = 'gaussian';
else
  basis_fcts = echcd_params.basis_fcts;
end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% init to right size
ngauss = numel(model.rgrid.rhogauss);
dPec_dx = zeros(ngauss,model.dims.nx);
djec_dx = zeros(ngauss,model.dims.nx);
dPec_du = zeros(ngauss,model.dims.nu);
djec_du = zeros(ngauss,model.dims.nu);

% Plasma profile information needed for evaluating CD
te     = stap.te;
dte_dx = stap.dte_dx;
ne = stap.ne;
dne_dx = stap.dne_dx;

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

rdep = zeros(1,nunits); % init

irdep = echcd_model.rdep_ext; % indices of rdep taken externally
rdep(~irdep) = echcd_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

assert(all(rdep>=0 & rdep<=1),'rdep must be between 0 and 1')

cd_eff = echcd_params.cd_eff; %current drive efficiency

rr = model.rgrid.rhogauss; %spatial grid in units of normalised rho_toroidal

% Heating
switch basis_fcts 
  case 'gaussian'
    wdep = echcd_params.wdep; %deposition width
    Pecprofile = exp(-4*bsxfun(@times,bsxfun(@minus,rr,rdep).^2,1./wdep.^2));
  case 'spline'
    wdeploc = echcd_params.rdep_outer/(echcd_params.nr_fcts-2); % nr of knot points should be 1 below the number of splines requested
    rdep = ( 0:wdeploc:echcd_params.rdep_outer )';
    S = splineval(rr, rdep, 3, [1 0 0])'; % 3rd order
    Pecprofile = repmat(S, 1, echcd_params.copy_basis);  
  case 'spline_located'
    % wdeploc = (echcd_params.rdep_outer-echcd_params.rdep_inner)/(echcd_params.nr_fcts-2); % nr of knot points should be 1 below the number of splines requested
    rdep = linspace(echcd_params.rdep_inner, echcd_params.rdep_outer, echcd_params.nr_fcts+1)';
    S = splineval(rr, rdep, 3, [1 0 0])'; % 3rd order
    S = S(:, 2:end-1);
    Pecprofile = repmat(S, 1, echcd_params.copy_basis); 
end 

Ptot = int_Vtot(Pecprofile,geop,model); % normalize to unit total power
dPec_dupow = Pecprofile*diag(1./Ptot);

Pecperunit = bsxfun(@times,dPec_dupow,u(iupow)'); % matrix of power per source
Pec = sum(Pecperunit,2); % total

% ECCD efficiency model
djeccdeff_dte = echcd_params.eta_eccd0 * exp(-echcd_params.cd_trap*(rr.^2)/0.5^2)./ne;
jeccdeff = djeccdeff_dte.*te; % profile of efficiency_cd(rho)

% take into account CD efficiency per actuator
djec_dupow = bsxfun(@times,jeccdeff*cd_eff,dPec_dupow); % matrix of efficiency per source
jecperunit = bsxfun(@times,djec_dupow,u(iupow)'); % jec profile per source
jec = sum(jecperunit,2);

if nargout>2
  if ~isempty(iurdep) % if some inputs are rdep
    % Pec = Pecprofile/Ptot : use chain rule
    % dPec_durdep = dPecprofile_durdep / Ptot
    %              - Pecprofile/Ptot^2*(dPtot/durdep)
    % Ptot = int(V' Pecprofile) dr
    %  -> dPtot/durdep = int(V'dPecprofile_durdep) dr
    
    dPecprofile_durdep = (Pecprofile(:,irdep) .* (8*bsxfun(@times,bsxfun(@minus,rr,rdep(irdep)),1./wdep(irdep).^2) ));
    dPec_durdep = bsxfun(@times,...
      bsxfun(@times,dPecprofile_durdep,1./Ptot(irdep)) ...
      - bsxfun(@times,Pecprofile(:,irdep),int_Vtot(dPecprofile_durdep,geop,model)./Ptot(irdep).^2),...
      u(iupow(irdep))');
    
    djec_durdep = bsxfun(@times,jeccdeff*cd_eff(irdep),dPec_durdep);
  else
    dPec_durdep = zeros(ngauss,0); % init
    djec_durdep = zeros(ngauss,0); % init
  end
  
  dPec_du(:,[iupow,iurdep]) = [dPec_dupow,dPec_durdep]; % assign to complete dPec_du
  dPec_dx = zeros(ngauss,model.dims.nx);
  
  % cd
  djec_du(:,[iupow,iurdep]) = [djec_dupow,djec_durdep]; % assign to complete dPec_du
  
  % sensitivity
  
  djeccdeff_dne = -jeccdeff./ne;
  djeccdeff_dx= bsxfun(@times,djeccdeff_dte,dte_dx) + bsxfun(@times,djeccdeff_dne,dne_dx);
  
  djec_dx = bsxfun(@times, djeccdeff_dx, dPec_du(:,iupow)*diag(echcd_params.cd_eff)*u(iupow) );
  
  % checks
  if any(u(iupow)<0)
    warning('negative EC power requested');
  end
  %
  
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;
end

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

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

%% Parameters

% Defaults per machine
switch model.tokamak
  case 'TCV'
    rdep = linspace(0,0.4,n_units);
    wdep = 0.1*ones(1,n_units);
    cd_eff = linspace(0,1,n_units);  
  case 'AUG'
    rdep = 0*ones(1,n_units); %use rhodep, not rdep
    wdep = 0.3*ones(1,n_units);
    cd_eff = 0*ones(1,n_units);
  case 'ITER'
    rdep = zeros(1,n_units); rdep(1:3) = [0.4 0.55 .2];
    wdep = .05*ones(1,n_units);
    cd_eff = 10*ones(1,n_units);  % current drive efficiency factor: pure ECH for first, co-ECCD for second  
  otherwise
    % use standard default
    rdep = linspace(0,0.6,n_units);
    wdep = 0.1*ones(1,n_units);
    cd_eff = linspace(0,1,n_units);
end
    

pp = struct(...
  'active',true,... % activation switch
  'rdep', rdep, ...  % centers of deposition gaussian, set to -1 if passed through input
  'wdep', wdep, ... % full 1/e width of gaussian
  'cd_eff', cd_eff, ... %  current drive effect
  'eta_eccd0', 7e14,... % coefficient for CD efficiency calculation
  'cd_trap', 0 ...  % Effect of trapped particles on efficiency as a function of rho
  );



end
