function varargout = chi_FF(varargin)
% Ad-hoc electron heat conductivity module for RAPTOR
% Used in F.Felici PhD. Thesis (EPFL 2011)
%
%   params = chi_FF; % get default parameters
%   [chie,dchie_dx,dchie_dp] = chi_FF(x,g,v,model,module_params) % run module
%
% returns heat conductivity and derivatives. all outputs are on grid [rhogauss]
%
% chie = cp.cneo + cp.cano*Teterm*psiterm*itbterm
%

%# codegen

% F. Felici CRPP, Jan 2011, TUe 2012

% note that q*rho/s = rho/(iota*s) = -rho/(iota * diota/drho*(rho/iota))
% = 1/(diota/drho)

% 19/05/2013 PG: added a term to the chie model which deals with sawteeth. When sawteeth
% occur (q<1) in the center of the plasma, transport is increased. This is modeled by using a heavyside
% function.

% 27/05/2013 PG: added a term to scale chie with (Te/100)^constant.

% Mar 2019: S. v.Mulders fixed itb expression to avoid Inf/Inf for large q

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin == 0,
  varargout{1} = 'chi_e'; % module type
  varargout{2} = [];
  % no config for this module
  return
elseif nargin==2
  %% pass params and model
  %% model
  mm.name = mfilename;
  
  varargout{1} = mm;
  pp=struct(...
    'cneo' , 0.5 , ... % neoclassical transport
    'cano' , 7 , ... % anomalous transport
    'witb' , 3 , ... % 1/width of itb effect
    'aitb' , 1 , ... % amplitude of itb (not always used)
    'ditb' , 0, ... % shear value below which itb effect is present
    'chi_csawtooth', 0, ... % amplitude of increased transport in core due to sawtooth,  0 == off
    'chi_wsawtooth', 10,... % 1/width of sawtooth effect
    'te_constant', 0 ,... % this constant determines the power of (Te/1000)^te_constant scaling chie, 0 == off
    'aTeterm', 0, ... % importance of Te term [0: none, 1: full] (only for modeltype = qrhodlnTe)
    'arhoterm', 1,... % weight of rho term [0:none, 1:full] (only for modeltype = qrhodlnTe)
    'apsiterm', 1, ... % importance of psi term [0: none, 1: full] (only for modeltype = qrhodlnTe)
    'chiecentral', 10, ... % Central extra chie to help flat central Te
    'deltacentral',0.15, ... % width of central chie
    'chiiscal',1,... % scaling factor for chii/chie
    'implicit',true... % choose whether to treat as implicit in solver
    );

  %% params
  varargout{2} = pp;
  return %empty call, probably to get default structures
elseif nargin==6; % change this depending on number of inputs to module
  stap = varargin{1};
  %geop = varargin{2};
  %trap = varargin{3};
  %it = varargin{4};
  model = varargin{5};
  pp = varargin{6}; % distribute inputs
else
  error('must call with 0 or 6 inputs');
end

%%%% CHIE MODULE

%%% Compute chie
iota     = stap.iota;
diota_dx = stap.diota_dx;
q =   stap.q; q_squared = q.^2;
te = stap.te;
dte_dx = stap.dte_dx;
shear = stap.shear;
dshear_dx = stap.dshear_dx;

% limit effect of q on transport below q=1
iotasaw = ones(size(iota));
iotasaw(iota<1) = iota(iota<1);
diotasaw_dx = zeros(size(diota_dx)); 
diotasaw_dx(iota<1, :) = diota_dx(iota<1, :); 

% no Te term
Teterm  = ones(model.rgrid.nrhogauss,1);
rhoterm = ones(model.rgrid.nrhogauss,1);
% psiterm = r*q
psiterm = (model.rgrid.rhogauss./iotasaw); % RLW like r*q

%%
itbexp = (exp(pp.witb*(pp.ditb-shear))); % exp term in denominator
itbterm  = pp.aitb*1./(1+itbexp) + (1-pp.aitb); %F(s)

% Construction of heaviside function F(q) which increases transport in the
% core due to sawteeth activity
sawteeth_exp  = exp( pp.chi_wsawtooth *(q - 0.95) ); % q - 0.95 is done instead of q - 1, to influence the region q>1 less
sawtooth_term = 1./(1+sawteeth_exp); %F(q)

chie_central_exp = exp(-model.rgrid.rhogauss.^2/pp.deltacentral^2);
chie_central = pp.chiecentral*chie_central_exp; % to play a bit

% Te scaling chie
Te_scaling      = (abs(te(1))/1000).^(pp.te_constant);

% Generic chie expression
rho_times_q = Teterm .* rhoterm .* psiterm;
chie = pp.cneo + pp.chi_csawtooth * sawtooth_term  + pp.cano * rho_times_q .* itbterm .* Te_scaling  + chie_central;

varargout{1} = chie;
varargout{2} = pp.chiiscal*chie;
%%% GRADIENTS w.r.t. Tehat and Psihat
if nargout>1 % then compute derivatives

  dchie_dx = zeros(model.rgrid.nrhogauss,model.dims.nx); % init
  
  % psiterm = r*q
  dpsiterm_dPsihat = -bsxfun(@times,model.rgrid.rhogauss .*q_squared,...
    diotasaw_dx(:,model.psi.xind));
  
  % Te_scaling
  te_scaled = abs(te(1)/1000);
  dTeterm_dTehat = (1/1000) * pp.te_constant * Te_scaling * te_scaled.^(-1) * dte_dx(1,model.te.xind);
  
  %%
  dshear_dPsihat = dshear_dx(:,model.psi.xind);
  %clf; plot(dshear_dPsihat2,'b'); hold on; plot(dshear_dPsihat,'r--');
  %%
  % F(s)
  ditbterm_dPsihat = bsxfun(@times,pp.aitb*pp.witb.* ( 1./( 1./itbexp + 2 + itbexp ) ),dshear_dPsihat);
  
  % sawtooth_term
  dsawtooth_term_dPsihat =  bsxfun(@times,-pp.chi_csawtooth * ( 1./( 1./sawteeth_exp + 2 + sawteeth_exp ) ) .* (-pp.chi_wsawtooth./(iotasaw.^2)),diotasaw_dx(:,model.psi.xind));
  
  %%
  
  % influence of this term on whole chie
  %%dchie_dPsihat_itbterm = diag(k*(gradlnTe+2*gradlnne).*psiterm)*itbterm_dPsihat;
  % chie ~ cano*psiterm(psi)*itbterm(psi)*teterm(Te)*Te_scaling(Te)
  dchie_dTehat  = bsxfun(@times,pp.cano * rho_times_q .* itbterm,dTeterm_dTehat);
  dchie_dPsihat = bsxfun(@times,pp.cano*rho_times_q .* Te_scaling,ditbterm_dPsihat) + ...
    bsxfun(@times,pp.cano*Teterm.*rhoterm.*itbterm .* Te_scaling,dpsiterm_dPsihat) + dsawtooth_term_dPsihat;
  
  dchie_dx(:,[model.psi.xind,model.te.xind])=[dchie_dPsihat,dchie_dTehat];
  
  varargout{3} = dchie_dx;
  varargout{4} = pp.chiiscal*dchie_dx;
  % Computation of chi derivative w.r.t model parameters
end

return
