%% RAPTOR Tutorial 2: Change input, display control and output analysis
%
% In this second tutorial we will customize how the simulation results are
% presented during the run, and will make some figures after the run.
% This will help underestand how profile data is stored and extracted in
% RAPTOR. We will also modify the input signal to have something slightly
% more interesting.

%% Preparation
% Start by repeating the first two steps of the previous tutorial
clear; close all hidden;
run(fullfile(pwd,'..','RAPTOR_path.m')); % add RAPTOR path
[config] = RAPTOR_config; % load default config

% build model for these params.
[model,params,init,g,v,U] = build_RAPTOR_model(config);

%% Manually customize display options
% Now we will customize the display options in the params structure to suit
% our needs
params.debug.verbosity = -1; % -1 indicates we manually set the output detail
params.debug.iterdisp = 5; % display every 5 iterations only (default=1)
params.debug.iterplot = 1; % plot profile evolution every iteration (default=0, no plot)

%% Run RAPTOR with new settings

% Choose a more interesting input: modify U
nt = numel(params.tgrid);
% the first row of U is the plasma current as a function of time. Assign
% it here:
U(1,params.tgrid <0.02) = linspace(80e3,250e3,sum(params.tgrid<0.02));
U(1,params.tgrid>=0.02) = 250e3;

% Update init to have correct plasma current at first time step
init.Ip0 = U(1,1);
%%
% Define the initial condition
x0 = RAPTOR_initial_conditions(model,init,g,v); 
%%
% Run the simulation
simres = RAPTOR_predictive(x0,g,v,U,model,params);

%% 
% The out function does the post-processing to compute all the outputs of
% interest.
out = RAPTOR_out(simres,model,params);
disp(out);

%% Customized display of results

% surf plots
subplot(231); surf(out.time,out.rho,out.psi); xlabel('t [s]'); ylabel('\rho'); title('\psi[Wb]');
subplot(232); surf(out.time,out.rho,out.q); xlabel('t [s]'); ylabel('\rho'); title('q');
subplot(233); surf(out.time,out.rho,out.Ip); xlabel('t [s]'); ylabel('\rho'); title('I_p');
subplot(234); surf(out.time,out.rho,out.upl); xlabel('t [s]'); ylabel('\rho'); title('U_{pl} [V]');
subplot(235); surf(out.time,out.rho,out.te); xlabel('t [s]'); ylabel('\rho'); title('T_e [eV]');
subplot(236); surf(out.time,out.rho,out.tep); xlabel('t [s]'); ylabel('\rho'); title('dT_e/d\rho [eV]');

%%
% note that some quantities in this output are also given on the grid of gauss
% points: out.rhogauss, though the normal grid version is also available.
% Beware, the normal grid point values are just a linear interpolation/extrapolation of
% the gauss point values.
clf;
hp=plot(out.rhogauss,out.jbsgauss(:,end),'x',out.rho,out.jbs(:,end),'o'); 
grid on; ylabel('[A/m^2]'); xlabel('\rho'); title('j_{BS}');
legend(hp,{'gauss grid points','main grid points'},'location','northeast')

%% Interactive plotting tool
% There exists an interactive plotting tool to plot several RAPTOR results 
% at the same time, allowing one to adjust the time of the plotted profiles
% with the mouse

outs{1} = out; % store old out structure

% Do a new RAPTOR run with different input parameters
params.debug.iterplot = 0; % do not plot this time
U = bsxfun(@times,U,linspace(1,2,numel(params.tgrid))); % scale U
simres = RAPTOR_predictive(x0,g,v,U,model,params);
outs{2} = RAPTOR_out(simres,model,params); % different structure

%%
% Plot 'overview'
hf = RAPTOR_plot_GUI(outs,'overview');
% drag your mouse over the time traces to see the time evolution
figure(hf); % make visible to export

%%
