%% RAPTOR NTM simulation tutorial
% NTMs can also be simulated, either by prescribing the time-evolution of
% the modes directly (affecting transport) or by solving MRE
% self-consistently.

%% Prescribed NTM time-evolution
% This simulation also includes sawteeth, you will see the effect!

disp('this tutorial does not work and is temporarily disabled')
return

close all hidden; clear all;
run('../RAPTOR_path.m');
[config] = RAPTOR_config; % load default params

config.saw = saw('porcelli');

% load an ntm model.
config.ntm = ntm('prescribed'); 
config.debug.iterdisp = 50;

% build RAPTOR structure
config.init.te0=100;
config.grid.tgrid = [0.:0.001:1];
% generate model structure for these params.
[model,params,init,g,v,U] = build_RAPTOR_model(config);

% define input
rampfun = @(t,tmin,ymin,tmax,ymax) ...
    max(ymin,min((ymax-ymin)/(tmax-tmin)*(t-tmin),ymax-ymin)+ymin); % anonymous function for ramps
init.Ip0 = 80e3;
U(1,:) = rampfun(config.grid.tgrid,0,init.Ip0,0.05,160e3);

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

% Prescribe NTM evolution
% make v0 vector time-varying
nt = numel(params.tgrid); vtvar = v*ones(1,nt);

% parameters
params.saw.s1crit = 0.3;
params.ntm.chiefact_A(:) = 2; % factor increase in chie

% define NTM appearance times and width
w43 = zeros(1,nt);
w43(params.tgrid > 0.2 & params.tgrid<0.3) = 0.05; % 4/3 NTM
w32 = zeros(1,nt);
w32(params.tgrid > 0.4 & params.tgrid<0.5) = 0.05; % 3/2 NTM
w21 = zeros(1,nt);
w21(params.tgrid > 0.6 & params.tgrid<0.7) = 0.05; % 2/1 NTM
w31 = zeros(1,nt);
w31(params.tgrid > 0.8 & params.tgrid<0.9) = 0.05; % 3/1 NTM

% assign to V
% Make sure that the indices match with the m and n parameters
vtvar(model.ntm.vind(1),:) = w21;
vtvar(model.ntm.vind(2),:) = w32;
vtvar(model.ntm.vind(3),:) = w43;
vtvar(model.ntm.vind(4),:) = w31;

% must turn these off for technical reasons
params.neos.implicit = false;
params.chi_e.implicit = false;

simres = RAPTOR_predictive(x0,g,vtvar,U,model,params);
out = RAPTOR_out(simres,model,params);

% inspect result
RAPTOR_plot_GUI(out,'overview');

plot(out.time,[out.w_32; out.w_43; out.w_31;out.w_21]);

%% NTM simulation using Modified Rutherford equation
% Includes ECRH power to stabilise the mode

close all hidden;
run('../RAPTOR_path.m');
[config] = RAPTOR_config; % load default params

% Specify the ECRH launcher, make sure that the last NTM launcher can be
% moved
config.echcd = RAPTORmodule('echcd_gaussian');
config.echcd.Configuration.n_units = 3;
config.echcd.Configuration.n_rdep_external = 1;

% load the Modified Rutherford Equation ntm model.
config.ntm = ntm('mre'); 
config.ntm.params.m=[2 3];
config.ntm.params.n=[1 2];
config.ntm.params.chiefact_A=[2 1]; % factor increase in chie
config.ntm.params.chiefact_w=2*ones(1,config.ntm.n_ntm);
config.ntm.params.w_d_e=0.01*ones(1,config.ntm.n_ntm);
config.ntm.params.w_d_i=0.01*ones(1,config.ntm.n_ntm);
config.ntm.params.a2_bs=0.1*ones(1,config.ntm.n_ntm);
config.ntm.params.a_ggj=ones(1,config.ntm.n_ntm);
config.ntm.params.conf_deg=zeros(1,config.ntm.n_ntm);
config.ntm.params.w_cd_ntm=0.01*ones(1,config.ntm.n_ntm);
config.ntm.params.rhodeltap_cl0=-config.ntm.params.m;


% build RAPTOR structure
config.debug.iterdisp = 50;
config.init.te0=100;
config.grid.tgrid = 0.:0.001:1;
% generate model structure for these params.
[model,params,init,g,v,U] = build_RAPTOR_model(config);

%  ECRH
params.echcd.active=1;
params.echcd.rdep=[0      0.4100 -1];
params.echcd.wdep=[0.3000 0.3500 0.35];
params.echcd.cd_eff=[0 1 1];
params.echcd.eta_eccd0=7.0000e+14;

%%
% define input
nt = numel(params.tgrid);

rampfun = @(t,tmin,ymin,tmax,ymax) ...
    max(ymin,min((ymax-ymin)/(tmax-tmin)*(t-tmin),ymax-ymin)+ymin); % anonymous function for ramps
init.Ip0 = 80e3;

P_central1 = zeros(1,size(config.grid.tgrid,2));
P_central1(:,params.tgrid >= 0.05 & params.tgrid <= 1.1) = 1.0e6;


NTM_stabpow2=zeros(1,size(config.grid.tgrid,2));
NTM_stabpow2(params.tgrid >= 0.5  & params.tgrid <= 0.6) = 0.1e6;

NTM_stabpow3=zeros(1,size(config.grid.tgrid,2));
rhodep3 = 0.45*ones(1,size(config.grid.tgrid,2));

U(1:5,:) = ...
    [rampfun(config.grid.tgrid,0,init.Ip0,0.05,160e3); ...
     P_central1;...
     NTM_stabpow2;...
     NTM_stabpow3;...
     rhodep3]; % rho dep

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

% ntm
% make v0 vector time-varying
vtvar = v*ones(1,nt);

%% 
% define NTM appearance times and width

% Seed the 3/2 to 5 cm which provides a lower bound between 0.4 s and 0.5 s
w32 = zeros(1,nt);
w32(params.tgrid >= 0.4 & params.tgrid <= 0.45) = 0.05; % 3/2 NTM, island held for some time
% Seed the 2/1 at 0.1 s to 15 cm
w21 = zeros(1,nt);
w21(params.tgrid == 0.1) = 0.15; % 2/1 NTM

% assign to V
% Make sure that the indices match with the m and n parameters
vtvar(model.ntm.vind(1),:) = w21;
vtvar(model.ntm.vind(2),:) = w32;

% run RAPTOR
params.neos.implicit = false;
params.chi_e.implicit = false;

params.ntm.doplot = false;
simres = RAPTOR_predictive(x0,g,vtvar,U,model,params);
% generate output structure

out = RAPTOR_out(simres,model,params);
%%
% Inspect result
RAPTOR_plot_GUI(out,'overview');

figure
hold on
plot(out.time,out.w_21,'x')
plot(out.time,out.w_32,'or')
xlabel('Time t [s]')
ylabel('Island width  w [m]')
legend('2/1','3/2')

