function varargout = sn_e(varargin)
% Electron particle sources
%
%   params = sn_e; % get default parameters
%   function [trap.Sne, trap.dSne_dx, trap.dSne_dxdot, trap.dSne_du] = sn_e(it,model,params);
%
% returns source for electrons and derivatives. all outputs are on grid [rhogauss]
%
% Several options are available, set by 'cp.modeltype' options:
%   *'manual': prescribed (fixed by default/time-varying)
%   *'benchmark': From Y.S. Na ITPA-IOS (April, 2017) benchmark effort

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin == 0,
    %% DEFAULT PARAMETERS
    module_params = struct(...
    'modeltype' , 'manual' , ... % default module
    'params', sn_manual); % load default params

    varargout{1} = module_params;
    return %empty call, probably to get default structures
elseif nargin == 1
    % get default parameters for each case
    def_params.modeltype = varargin{1}; 
    switch def_params.modeltype
        case 'manual'
            def_params.params = sn_manual;
        case 'benchmark'
            def_params.params = [];
        otherwise
            error('invalid Sne model type ''%s''',def_params.modeltype)
    end
    varargout{1} = def_params; % return appropriate parameter structure
    return
elseif nargin==3; % change this depending on number of inputs to module
    it     = varargin{1};
    model  = varargin{2};
    params = varargin{3}; % distribute inputs
else
    error('must call with 0 or 3 inputs');
end

% Initalize
zz = zeros(model.rgrid.nrhogauss,1);
dzz_dx = zeros(model.rgrid.nrhogauss,model.dims.nx);
dzz_du = zeros(model.rgrid.nrhogauss,model.dims.nu);
% Set to zero matrices
Sne = zz;
dSne_dx = dzz_dx;
dSne_dxdot = dzz_dx;
dSne_du = dzz_du;

switch model.sn_e
    case 'manual'
        % Sne is prescribed
        Sne = sn_manual(params.sn_e,it,model);
    case 'benchmark'
        % From Y.S. Na ITPA-IOS benchmark effort
        rhogauss = model.rgrid.rhogauss;
        PhiN = rhogauss.^2;
        
        S0 = 7.5e20;
        C = 0.25e24;
        d = 0.225;
        Sne_edge = S0*exp(15*(PhiN - 1));
        Sne_core = C*d^2* PhiN.^(6.5).*(1-PhiN).^(8.5)./(d^2+(PhiN-0.5).^2) ;
      
        Sne = (Sne_core+Sne_edge); 
end

% outputs
varargout{1} = Sne;
if nargout>1
    varargout{2} = dSne_dx;
    varargout{3} = dSne_dxdot;
    varargout{4} = dSne_du;
end
return

%% Manual Sne
function varargout = sn_manual(varargin)
if nargin == 0,
    %% DEFAULT PARAMETERS
    module_params=struct(...
        'sne_pr',   zeros(21,1), ... % prescribed Sne: constant (default), fixed profile (rho_gauss), time-varying
        'rho', linspace(0,1,21) ...
        );
    varargout{1} = module_params;
    return %empty call, probably to get default structures
elseif nargin==3; % change this depending on number of inputs to module
    sp = varargin{1};
    it = varargin{2};
    model = varargin{3};
else
    error('must call with 0 or 3 inputs');
end

assert(size(sp.sne_pr,1)==numel(sp.rho),'size of profile must match number of rho elements')

if (size(sp.sne_pr,2)==1)&&(size(sp.sne_pr,1)~=1)
    % Fixed profiles
    Sne = sp.sne_pr;
else
    % Time-varying Dne
    Sne = sp.sne_pr(:,it);
end

varargout{1} = interp1(sp.rho,Sne,model.rgrid.rhogauss);
return