%% RAPTOR tutorial - Variable timegrid
% It is possible to use a time-varying grid, typically to have 
% shorter time steps in the early (transient) part of the shot and
% fewer time points later.
%
% This is illustrated with a simple ITER-like simulation below
%% Set up the simulation

close all hidden;
clear; run('../RAPTOR_path.m');

addpath(fullfile(pwd,'..','code'),'-end'); % add RAPTOR path
[config] = RAPTOR_config('ITER'); % load default ITER params

config.nbhcd.Configuration.n_units = 1;

config.echcd = RAPTORmodule('none');

% create model, params, init, geometry g, kinetic profiles v.
[model,params,init,g,v,~] = build_RAPTOR_model(config);

% set NBI parameters
params.nbhcd.active = true;
params.nbhcd.rdep = 0.25; % 
params.nbhcd.wdep = 0.25; % broad heating
params.nbhcd.wdep_out = 0.75; % broad heating
params.nbhcd.cd_eff = 0; % current drive

% switch off alpha power and losses
params.palpha.active = false;
params.pbrem.active = false;
params.pei.active = false;
params.prad.params.active = false;
params.debug.verbosity = 2;

%% Simulation with variable time step
% define a monotonic, but not evenly spaced time grid:
tgrid0 = [0:0.5:3,4:2:80];
params.tgrid = tgrid0;

% inputs
rampfun = @(t,tmin,ymin,tmax,ymax) ...
    max(ymin,min((ymax-ymin)/(tmax-tmin)*(t-tmin),ymax-ymin)+ymin); % anonymous function for ramps
% Construct actuator inputs on the time grid tgrid0.
U0(1,:) = rampfun(tgrid0,0,4e6,80,12e6); % Ip
U0(2,:) = rampfun(tgrid0,30,0,50,40e6); % NBI power

%%
init.te0  = 3000;
init.jpow = 0;
% Construct kinetic profiles according to the new settings.
v = build_kinetic_profiles(model,init);

% initial condition
init.Ip0 = U0(1,1);
x0 = RAPTOR_initial_conditions(model,init,g(:,1),v(:,1)); 
simres0 = RAPTOR_predictive(x0,g,v,U0,model,params);
out0 = RAPTOR_out(simres0,model,params);

%% Simulation with constant time step
% Repeat the same simulation but with a fixed time step.

tgrid1 = 0:0.1:80;
params.tgrid = tgrid1;

% inputs
rampfun = @(t,tmin,ymin,tmax,ymax) max(ymin,min((ymax-ymin)/(tmax-tmin)*(t-tmin),ymax-ymin)+ymin); % anonymous function for ramps
% Construct actuator inputs on the time grid tgrid1.
U1(1,:) = rampfun(tgrid1,0,4e6,80,12e6); % Ip
U1(2,:) = rampfun(tgrid1,30,0,50,40e6); % NBI power

% initial condition
init.Ip0 = U1(1,1);
params.debug.verbosity = -1; % let user define details
params.debug.iterdisp = 10; % disp every 10
params.debug.iterplot = 00; % no plots

x0 = RAPTOR_initial_conditions(model,init,g(:,1),v(:,1)); 
simres1 = RAPTOR_predictive(x0,g,v,U1,model,params);
out1 = RAPTOR_out(simres1,model,params);

%% Compare
% Comparing the two simulations confirms their equivalence
subplot(311)
hh=plot(tgrid1,out1.iota(1:5:end,:),'b.-',tgrid0,out0.iota(1:5:end,:),'ro');
ylabel('iota []'); xlabel('t[s]');
legend(hh([1,end]),'fixed dt','variable dt');

subplot(312)
plot(tgrid1,out1.te(1:5:end,:)/1e3,'b.-',tgrid0,out0.te(1:5:end,:)/1e3,'ro');
ylabel('T_e [keV]'); xlabel('t[s]');

subplot(313)
plot(tgrid1,out1.psi(1:5:end,:),'b.-',tgrid0,out0.psi(1:5:end,:),'ro');
ylabel('\psi [Wb]'); xlabel('t[s]');