%% RAPTOR Exercise 1: Ramp-up current diffusion with varying heating
% In this exercise you will test the effect of differently timed heating
% during a plasma current ramp. 

%% ------------------ ITER simulation setup ------------------
% First we set up the RAPTOR simulation.
[~, current_folder,~] = fileparts(pwd); 
assert(strcmp(current_folder,'exercises'), 'You are not in the correct folder. Please open the exercices folder'); 
clear
run ../RAPTOR_path.m

%% Configuration
% load the default parameters for ITER
config = RAPTOR_config('ITER'); % load default params
config.grid.tgrid = [0:0.1:1,3:2:140];

% tranport model: choose empirical model with 
% transport reduction at low shear.
config.chi_e = RAPTORmodule('chi_FF');

% Heating actuators: EC and NBI
config.nbhcd = RAPTORmodule('hcd_gaussian');
config.echcd = RAPTORmodule('echcd_gaussian');
config.echcd.Configuration.n_units = 3; % nb of ECHCD launchers 
config.echcd.Configuration.n_rdep_external = 0; % put 0 if you want to fill rdep directly in params.echcd
config.nbhcd.Configuration.n_units = 1; % nb of NBI launchers 
config.nbhcd.Configuration.n_rdep_external = 0; % put 0 if you want to fill rdep directly in params.nbhcd

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

%% Parametrization

% Energy sources and sinks
params.prad.active = false; % impurity radiation
params.pei.active = true; % electron-ion equipartition power 
params.palpha.active = false; % alpha heating power

% Heat diffusion
params.chi_e.cano = 4; % anomalous diffusion factor
params.chi_e.cneo = 0.1; % neoclassical diffusion
params.chi_e.te_constant = 1.2; % includes a factor Te^1.2 which mimics the impact of stiffness

% NBI
% (be careful that array sizes are consistent with nb_units)
params.nbhcd.active = true;
params.nbhcd.rdep = [0.26]; %
params.nbhcd.wdep = [0.75]; % gaussian width
params.nbhcd.wdep_out = [0.75]; % gaussian width
params.nbhcd.cd_eff = [5]; % current drive efficiency factor
% the index for the power trace in the input vector is stored in the model
% structure
disp(model.nbhcd.uind(model.nbhcd.ind_power))

% ECHCD
% (be careful that array sizes are consistent with nb_units)
params.echcd.active = true;
params.echcd.rdep = [0. 0.3 0.4]; % three actuators, set rdep=0 (on-axis) for first, rdep=0.3 for second, rdep=0.4 for third
params.echcd.wdep = [.05 .05 .05];  % wdep = 0.05 for each
params.echcd.cd_eff = [0 20 20];  % current drive efficiency factor: the on-axis EC source provides pure heating, the off-axis sources provide co-ECCD
% the indices for the power trace in the input vector is stored in the model
% structure
disp(model.echcd.uind(model.echcd.ind_power))

params.debug.iterdisp = 10; % display every 10

% Plasma current ramp (Ip) & Heating powers
rampfun = @(t,tmin,ymin,tmax,ymax) max(ymin,min((ymax-ymin)/(tmax-tmin)*(t-tmin),ymax-ymin)+ymin); % anonymous function for ramps

U(1,:) = rampfun(params.tgrid,0,4e6,100,12e6); % input Ip trace
% the default Ip trace ramps linearly to 12MA at 100s, in ex 1, you will
% impose different Ip ramp rates, reaching full current at an earlier time
U(2,:) = zeros(size(U(1,:))); % placeholder, will replace this with heating later.
U(3,:) = zeros(size(U(1,:))); % placeholder, will replace this with heating later.
U(4,:) = zeros(size(U(1,:))); % placeholder, will replace this with heating later.
U(5,:) = zeros(size(U(1,:))); % placeholder, will replace this with heating later.

figure(1); 
plot(params.tgrid, U); 
title('Ohmic current & Heatings inputs'); 
xlabel('t (s)'); 
legend('U1/Ip [A]', 'U2/Pec1 [W]', 'U3/Pec2 [W]', 'U4/Pec3 [W]', 'U5/Pnb [W]');

% initial conditions
init.Ip0 = U(1,1);
init.te0 = 2e3;
init.ne0 = 1e19;

x0 = RAPTOR_initial_conditions(model,init,g,v);  % Define the initial condition

% initialize V
V = repmat(v, 1, numel(params.tgrid));

% specify density ramp: multiply initial profile by ne0(t)
negauss = model.ne.Lamgauss*V(model.ne.vind,1);
ne0 = rampfun(params.tgrid,0,1,100,5); 
negauss_varying = negauss*ne0; 
V(model.ne.vind,:) = model.ne.Lamgauss \ negauss_varying;
V(model.ni.vind,:) = 0.9*V(model.ne.vind,:);

figure(2); hold on;  
plot(params.tgrid, V(model.ne.vind(1),:), 'b'); 
plot(params.tgrid, V(model.ni.vind(1),:), 'r'); 
title('Electron & Ion density inputs'); xlabel('t (s)');
legend('n_e(\rho=0) [m-3]', 'n_i(\rho=0) [m-3]');
hold off; 

%% Run the simulation
simres = RAPTOR_predictive(x0,g,V,U,model,params);
out = RAPTOR_out(simres,model,params);
RAPTOR_plot_GUI(out);

%% ------------------ EXERCISES ------------------
% 
% # 1. By varying the Ip trace in U(1,:), investigate the effect of different plasma current ramp rates on the
% speed of penetration of inductive current, the evolution of the loop voltage profile upl and the q
% profile. 
% Note: adjust correspondingly the linear density ramp, such that the maximum density is reached when Ip reaches 12MA.
% Plot and interpret the time evolution of the edge loop voltage
% out.upl(end,:) and the internal inductance out.li3.
% Compare the different simulations (with different Ip ramp-up rates) in a plot with Ip on the x-axis and the internal inductance li3
% on the y-axis. How can we interpret the difference in li3 during ramp-up (i.e. what is the cause)  
% and what are the physics consequences?
% Compare the radial profiles upl, q and jpar at the start of the flat-top
% phase for the different cases, what can you say about the current
% distribution and stationarity?

% # 2. Returning to the original Ip time trace (reaching 12MA at 100s), now use U(5,:) to add 16.5MW 
% of NBI power starting at different times during the ramp-up.
% Examine the effect on the T_e, q and jpar profiles and explain the results.
% How does the onset of heating impact the time trace out.jpar(1,:)
% (current density on-axis)? 
% Explore how transiently negative magnetic shear during the ramp-up can
% cause the creation of an internal transport barrier (the effect can be
% highlighted by comparing simulations with params.chi_e.aitb respectively to 0 and 1).

% # 3. Use U(5,:) to add 33MW from 50s onwards and the U(2:4,:) from 75s to add at 20MW of EC
% By turning on the various entries U(2:4,:), consider respectively:
% 1. as on-axis pure heating (U(2,:))
% 2. as off-axis co-ECCD at rhodep = 0.3 (U(3,:))
% 4. as off-axis co-ECCD at rhodep = 0.4 (U(4,:))
% Investigate how off-axis ECCD leads to a transport barrier at the start of flat-top (100s) by causing a
% reversed-shear q profile.
% How does off-axis ECCD impact the minimum of the q profile? Can you
% adjust off-axis deposition radius params.echcd.rdep(3) such that
% qmin>1.5?

% # 4. Find a combination of timing for heating and current drive to keep qmin
% above 1 for as long as possible.
