function varargout = runaways_analytic(varargin)
% returns runaway current jre*B0

coder.extrinsic('warning');
persistent f

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Input processing
if all(nargin ~= [0,6])
    error('must call with 0 or 6 inputs');
elseif nargin == 0,
    %% DEFAULT PARAMETERS
    runaways_params_default = struct(...
        'active',true, ...
        'doplot',false,...
        'dt',1e-3,...
        'ks',10e3,...
        'tau',0.1,...
        'Ecrit',1 ...
        );
    
    varargout{1} = runaways_params_default;
    return %empty call, probably to get default structures
end

% call with 6 inputs, during time loop
stap = varargin{1};
geop = varargin{2};
u = varargin{3};
it = varargin{4};
model  = varargin{5};
runaways_params = varargin{6}; % distribute inputs

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Do the actual work here
if isempty(f) || it==1
  f = zeros(model.rgrid.nrhogauss,1);
end

% init to right size
if ~runaways_params.active % if not active, return zeros
  nrhogauss = model.rgrid.nrhogauss;
  jreB = zeros(nrhogauss,1);
  return
end

%% If active, do the actual calculations

%% example of differential equation (with memory)
upl = stap.upl;
Ediff = (upl/(2*pi*model.equi.R0) - runaways_params.Ecrit);
df_dt = 1/runaways_params.tau*(-f + runaways_params.ks*(Ediff.*(Ediff>0)));
f = f + df_dt*runaways_params.dt;
jre = f;

if runaways_params.doplot
  figure(1);
  subplot(211);
  plot(df_dt)
  subplot(212)
  plot(f);
  drawnow
end
%% example of analytical equation
%jre = -exp(-beta*(runaways_params.Ecrit-upl/5);

jreB = jre*model.equi.B0;


varargout{1} = jreB;
return
