function varargout = vpdn_e(varargin)
% Electron particle transport
%
%   params = vpdn_e; % get default parameters
%   function [Dne,Vpe,dDne_dx,dDne_dxdot,dDne_du,dVpe_dx,dVpe_dxdot,dVpe_du] ...
%           = vpdn_e(model,geop,model,trap)
%
% returns particle diffusion coefficient Dn, pinch velocity Vn and derivatives. all outputs are on grid [rhogauss]
%
% Several options are available, set by 'cp.modeltype' options:
%   *'vpdn_manual': prescribed (constant/fixed/time-varying profile)
%   *'vpdn_chiescal': calculated from chi_e (default)
%   *'vpdn_benchmark': From Y.S. Na ITPA-IOS (April, 2017) benchmark effort
%   *'vpdn_MS': Use analytical transport model from Kim, et al PPCF 2016(58) 055002
%   *'vpdn_QLKNN10D': NN Based particle transport model from qualikiz
%   *'vpdn_QLKNN4Dkin': NN Based particle transport model from qualikiz

coder.extrinsic('vpdn_MS');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin == 0
  %% DEFAULT PARAMETERS
  varargout{1} = RAPTORmodule('vpdn_manual');
  return %empty call, to get default structures
elseif nargin == 1
  error('deprecated use, use RAPTORmodules instead')
elseif nargin==6
  stap   = varargin{1};
  geop   = varargin{2};
  trap   = varargin{3};
  it     = varargin{4};
  model  = varargin{5};
  params = varargin{6};
else
  error('must call with 0 or 6 inputs');
end

% Initialize
dzz_dx = zeros(model.rgrid.nrhogauss,model.dims.nx);
dzz_du = zeros(model.rgrid.nrhogauss,model.dims.nu);
% Set to zero matrices
if nargout>2
  dDne_dx = dzz_dx;
  dDne_dxdot = dzz_dx;
  dDne_du = dzz_du;
  dVpe_dx = dzz_dx;
  dVpe_dxdot = dzz_dx;
  dVpe_du = dzz_du;
end

switch model.vpdn_e.name
  case 'vpdn_manual'
    % Dne is prescribed
    % Vpe is prescribed
    [Dne, Vpe] = vpdn_manual(stap,geop,trap,it,model,params.vpdn_e);
  case 'vpdn_chiescal'
    % Dne is lineraly proportional to chi_e
    % Vpe is lineraly proportional to Dne
    [Dne, Vpe] = vpdn_chiescal(stap,geop,trap,it,model,params.vpdn_e);
  case 'benchmark'
    % From Y.S. Na ITPA-IOS benchmark effort
    rhogauss = model.rgrid.rhogauss;
    R0 = 6.20;
    B0 = 5.3;
    PhiN = rhogauss.^2;
    PhiNped = 0.88;
    Phib = geop.Phib;
    icore = (PhiN<=PhiNped);
    
    D0 = 0.5;
    D1 = 1.0;
    D2 = 0.11;
    V0 = 1.385;
    
    Dne(icore) = D0 + D1*(PhiN(icore)).^2;
    Dne(~icore) = D2;
    
    v0 = V0*rhogauss.*Dne/R0;
    
    rhob = sqrt(Phib/(B0*pi));
    Vpe = -v0.*(geop.g1./(geop.g0.*geop.Vp/rhob)); % = v0*<gradrho^2>/<gradrho> = <gradV^2>/<gradV>/dVdrho
    
  case 'vpdn_MS'
    if nargout<=2
      [Dne,Vpe] = vpdn_MS(stap,geop,trap,it,model,params.vpdn_e);
    else
      [Dne,Vpe,dDne_dx,dDne_dxdot,dDne_du,dVpe_dx,dVpe_dxdot,dVpe_du] = vpdn_MS(stap,geop,trap,it,model,params.vpdn_e);
    end
  case 'vpdn_QLKNN4Dkin'
    [tcoeff,dtcoeff_dx] = vpdn_QLKNN4Dkin(stap,geop,trap,it,model,params.vpdn_e);
    Dne = tcoeff.de;  Vpe = tcoeff.ve;
    dDne_dx = dtcoeff_dx.dde_dx; dVpe_dx = dtcoeff_dx.dve_dx;
  case 'vpdn_QLKNN10D'
    [tcoeff,dtcoeff_dx] = vpdn_QLKNN10D(stap,geop,trap,it,model,params.vpdn_e);
    Dne = tcoeff.de;  Vpe = tcoeff.ve;
    dDne_dx = dtcoeff_dx.dde_dx; dVpe_dx = dtcoeff_dx.dve_dx;
  otherwise
    error('not supported')
end

% outputs
varargout{1} = Dne;
varargout{2} = Vpe;

if nargout>2
  varargout{3} = dDne_dx;
  varargout{4} = dVpe_dx;
  varargout{5} = dDne_dxdot;
  varargout{6} = dVpe_dxdot;
  varargout{7} = dDne_du;
  varargout{8} = dVpe_du;
end
end
