%% RAPTOR Tutorial 4: Time varying kinetic profiles
% Example of how to input time-varying kinetic profiles ne,ni,ti,zeff
% These are provided externally to the simulation.
clear;close all hidden;
run ../RAPTOR_path.m

% load the default parameters,
[config] = RAPTOR_config; % load default config

% HCD
config.echcd = RAPTORmodule('echcd_gaussian');

% generate the following structures for these params: model, params, init, geometry g, kinetic profiles v
[model,params,init,g,v,U] = build_RAPTOR_model(config);

params.echcd.active = true; % activate the module, otherwise it does nothing.
params.echcd.rdep = [0 0.4]; % two actuators, set rdep=0 for first, rdep=0.4 for second
params.echcd.wdep = [0.35, 0.35];  % wdep =0.35 for each
params.echcd.cd_eff = [0 1];  % current drive efficiency factor: pure ECH for first, co-ECCD for second
params.echcd.uindices = uint32([2 3]); % index of power for each actuator in input vector

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

% simple inputs: constant current plasma
U(1,:) = 200e3*ones(size(params.tgrid)); % input Ip trace
U(2:3,:) = zeros(2,numel(params.tgrid));
% initial conditions
init.Ip0 = U(1,1);
% Define the initial condition
x0 = RAPTOR_initial_conditions(model,init,g(:,1),v(:,1)); 

% run basic RAPTOR simulation
simres = RAPTOR_predictive(x0,g,v,U,model,params);
out0 = RAPTOR_out(simres,model,params);

%% time varying ne: scale ne by time-dependent factor
necoeffs = v(model.ne.vind); %NB these are the spline coefficients, not the profile itself!
% initialize V0
nt = numel(params.tgrid);
V0 = repmat(v,1,nt);

% factor
nefactor = linspace(1,3,nt);
% clip to 2
nefactor(round(numel(nefactor)/2):end) = nefactor(round(numel(nefactor)/2)); 

% substitute new time-varying ne0
V0(model.ne.vind,:) = bsxfun(@times,necoeffs,nefactor);

% run RAPTOR simulation
simres = RAPTOR_predictive(x0,g,V0,U,model,params);
out1 = RAPTOR_out(simres,model,params);

%% time varying Zeff
ze0 = v(model.ze.vind); % get original Zeff profile coefficients
V0 = repmat(v,1,nt); % reset to old values

zefactor = 1+linspace(1,0,numel(ze0))'*linspace(0,2,nt); 
% Zeff ramp by factor 3 in core only, keep edge fixed
V0(model.ze.vind,:) = bsxfun(@times,ze0,zefactor);
simres = RAPTOR_predictive(x0,g,V0,U,model,params);

out2 = RAPTOR_out(simres,model,params);

%% plot all three cases
RAPTOR_plot_GUI({out0,out1,out2})
% end of tutorial