%% RAPTOR tutorial: Runaways
% Add description here
%
%

clear; close all hidden;
run ../RAPTOR_path.m

[config] = RAPTOR_config; % load default config

%% 
% *runaway parameters*
% To specify a runaway model, we need to choose a module
config.runaways = runaways('runaways_analytic');

% turn off hcd modules
config.echcd = RAPTORmodule('none');
config.nbhcd = RAPTORmodule('none');

% Create RAPTOR model,params and init structure with new parameters
[model,params,init,g,v,U] = build_RAPTOR_model(config);

% time grid
dt = 1e-3;
params.tgrid = 0:dt:0.3;
% For this analytical model, we need to prescribe the time step externally
params.runaways.dt = dt;
params.runaways.Ecrit = 0.1;
params.runaways.ks = 100e3;
params.runaways.doplot = false;

%%  Inputs
% plasma current
U = 200e3*ones(size(params.tgrid)); % input Ip trace: ramp from 40kA to 200kA
U(1,params.tgrid<0.02) = linspace(80e3,200e3,sum(params.tgrid<0.02));

%%
% *Initial conditions and run model*
% initial conditions
init.Ip0 = U(1,1);
% Define the initial condition for state
x0 = RAPTOR_initial_conditions(model,init,g,v); 
% run RAPTOR
simres = RAPTOR_predictive(x0,g,v,U,model,params);
%%
% *plot outputs*
out=RAPTOR_out(simres,model,params);

%%
clf;
subplot(311); surf(out.time,out.rho,out.jre); 
xlabel('t [s]');ylabel('\rho'); title('j_{re} [A/m^2]'); shading flat;
subplot(312); surf(out.time,out.rho,out.upl); 
xlabel('t [s]');ylabel('\rho'); title('upl'); shading flat;
subplot(313); plot(out.time,[out.Ip(end,:);out.Ibs(end,:);out.Ire(end,:)]); 
legend('Ip','bs','re');
