%% RAPTOR tutorial: ITER scenario with alpha particle heating
% In this tutorial we simulate the alpha particle heating for an
% ITER scenario.

%% ITER scenario parameters
clear; close all hidden;
run(fullfile(pwd,'..','RAPTOR_path.m')); % add RAPTOR path
% ramp function
rampfun = @(t,tmin,ymin,tmax,ymax) max(ymin,min((ymax-ymin)/(tmax-tmin)*(t-tmin),ymax-ymin)+ymin); % anonymous function for ramps
[config] = RAPTOR_config('ITER'); % load default ITER params

% Physics modules for standard ITER scenario
% chi_e
config.chi_e = RAPTORmodule('chi_BgB');

% ECHCD
config.echcd = RAPTORmodule('echcd_gaussian');
config.echcd.Configuration.n_units = 3;

% NBI
% add physics-based NBI module
config.nbhcd = RAPTORmodule('hcd_gaussian');
config.nbhcd.Configuration.n_units = 1;

% Timing   
t_0    = 0.0; t_end  = 10;
Ts = 1;    % slow sample time

% time grid
config.grid.tgrid = unique([t_0  : Ts :t_end]);
                  
% spatial grid
nrho=21;
config.grid.rhogrid = [0:0.05:0.9,0.925:0.025:1];
config.te.BC.defaultValue = 300;
config.ni.method = 'nescal';
config.ti.method = 'tescal';

% Build model 
% generate model, params, init, geometry g, kinetic profiles v
[model,params,init,g,v,U] = build_RAPTOR_model(config);

% chi_e
params.chi_e.cneo = 0.3;
params.chi_e.aegb = 1.8 * params.chi_e.aegb;
params.chi_e.aeb  = 0.3 * params.chi_e.aeb;

params.pbrem.active = true;
params.pei.active = true;
params.prad.active = false;

% ECHCD
params.echcd.active = true;
params.echcd.rdep = [0.05 0.2 0.3]; %
params.echcd.wdep = [.1 .15 .15];  %
params.echcd.cd_eff = [5 5 15];  % current drive efficiency factor
% {JL} Forse da rimuovere params.echcd.uindices = uint32([3 4 5]); % index of power for each
%actuator in input vector

% NBI
% add physics-based NBI module
params.nbhcd.uindices=uint32([2]);
params.nbhcd.rdep = [0.23]; % 
params.nbhcd.wdep = [0.15]; % broad heating
params.nbhcd.wdep_out = [1.1]; % broad heating
params.nbhcd.cd_eff = [4]; % current drive
params.nbhcd.active = true;
params.nbhcd.frac_to_electrons = 0.72;

% switch on alpha heating
params.palpha.active = true;
params.palpha.calpha = 1.0; % factor scaling alpha power (see P_alpha.m for details)

%% Actuator inputs
% plasma current
Ipflattop = 12e6;

U(1,:) = Ipflattop;

% NBI
Pnb_max = 33e6;
Pec_max = 30e6;
U(model.nbhcd.uind,:) = Pnb_max; % NBI

% ECCD
% rdep = [0.05 0.2 0.3]; %
U(3,:) = Pec_max;

%% Kinetic profile settings
% initial conditions 
init.newidth = 2;
init.ne0 = 1e20;
init.ni0scal = 0.9;
init.niescal = 0.9;
init.ti0scal = 1;
init.tiescal = 1;

init.te0 = 15e3;
init.tewidth = 0.6; 
init.ze0 = 1.6;
init.zee = 1.6;
init.jpow = 1;
init.Ip0 = U(1,1);

% Construct kinetic profiles according to the new settings.
v = build_kinetic_profiles(model,init);
% Define the initial condition
x0 = RAPTOR_initial_conditions(model,init,g(:,1),v(:,1));

% hmode
v(model.hmode.vind.activation,:) = 1;
v(model.hmode.vind.rhoped,:) = 0.9;
v(model.hmode.vind.gteped,:) = -(4.3e3-model.te.BC.defaultValue)/0.1;

%% Run RAPTOR
% H-mode simulation
simres = RAPTOR_predictive(x0,g,v,U,model,params,'verbosity',2);
% output
out = RAPTOR_out(simres,model,params);

%% Plot results
subplot(321); plot(out.time,out.U(1,:)/1e6); title('Ip [MA]');
subplot(322); plot(out.time,out.U(2:end,:)/1e6); title('P_{aux} [MW]')
subplot(323); plot(out.time,out.Palpha(end,:)/1e6); title('P_{alpha} [MW]')
subplot(324); plot(out.time,out.f_ss); title('f_{ss} (steady stateness)')
subplot(325); [cs,h] = contour(out.time,out.rho,out.q,[1 3/2 2 3 4],'color','k'); clabel(cs,h,'labelspacing',72); title('rational q locations');
xlabel('t [s]'); ylabel('\rho');
subplot(326); [cs,h] = contour(out.time,out.rho,out.te/1e3,[0.5 1 2 4 8 16 32],'color','k'); clabel(cs,h,'labelspacing',72); title('T_e contours [keV]')
xlabel('t [s]'); ylabel('\rho');

