%% RAPTOR tutorial 3: Customizing heating & current drive
% In this tutorial we will add and customize auxiliary heating and current
% drive systems. We will have to specify the deposition location, as well
% as the power traces for each actuator.
%
% This tutorial uses simple EC and NBI modules where the deposition is
% modeled by gaussian distributions. The EC current drive efficiency is
% scaled by Te/ne, while the NBI current drive efficiency is constant.

%% Example 1: Add 2 EC sources at fixed location
%%
% *Preparation*
% Start by repeating the first two steps of the previous tutorial
clear; close all hidden;
run ../RAPTOR_path.m

[config] = RAPTOR_config; % load default config

%% 
% *heating & current drive parameters*
% To specify a heating and current drive actuator, we must choose a module.
% This is done through the RAPTOR_module class.

config.echcd = RAPTORmodule('echcd_gaussian');

disp(config.echcd);

disp(config.echcd.Configuration);

%%
% Inspecting the result, we see that the config.echcd contains
% * the name of the module
% * the type of the module
% * the default Configuration structure for the module
%
% Now we can change the default module configuration, e.g. to change the
% number of EC units of the module.
config.echcd.Configuration.n_units = 3;

%% 
% *Other modifications*
% make some further modifications to the configuration structure, look in RAPTOR_config.m for exact parameter definitions
config.grid.tgrid = [0:0.001:0.1]; % time grid
config.grid.rhogrid = linspace(0,1,21); % spatial grid

config.debug.verbosity = -1; % override default behaviour 
config.debug.iterplot = 0; % plot profiles every N iterations
config.debug.iterdisp = 10; % display progress every N iterations

% Create RAPTOR model,params and init structure with new parameters
[model,params,init,g,v,U] = build_RAPTOR_model(config);

%% Now we inspect the model and params structure for our module
% Note that anything in the model can not be changed without rerunning
% build_RAPTOR_model()
disp(model.echcd)

disp(params.echcd)

%% We modify the echcd parameters to suit our needs
params.echcd.active = true; % activate the module, otherwise it does nothing.
params.echcd.rdep = [0 0.4 0.6]; % two actuators, set rdep=0 for first, rdep=0.4 for second
params.echcd.wdep = [0.35, 0.35 0.2];  % gaussian deposition width
params.echcd.cd_eff = [0 1 0];  % current drive efficiency factor: pure ECH for first, co-ECCD for second

%% 
% *Define inputs*
% We must now define more inputs than the plasma current alone

% plasma current
U(1,:)= 200e3*ones(size(params.tgrid)); % input Ip trace: ramp from 40kA to 200kA
U(1,params.tgrid<0.02) = linspace(80e3,200e3,sum(params.tgrid<0.02));
% first EC actautor: on-axis heating 1MW starting at 60ms
U(2,:) = zeros(size(params.tgrid)); U(2,params.tgrid>0.06) = 1e6;
% second EC actuator power: off-axis current drive, 1MW start at 40ms
U(3,:) = zeros(size(params.tgrid)); U(3,params.tgrid>0.04) = 1e6;
% third EC actuator: off-axis heating
U(4,:) = zeros(size(params.tgrid)); U(4,params.tgrid>0.04) = 1e6;

%%
% *Initial conditions and run model*
% initial conditions
init.Ip0 = U(1,1);
% Define the initial condition for state
x0 = RAPTOR_initial_conditions(model,init,g,v); 
% run RAPTOR
simres = RAPTOR_predictive(x0,g,v,U,model,params);
%%
% *plot outputs*
out=RAPTOR_out(simres,model,params);
clf;
subplot(211); surf(out.time,out.rho,out.pauxe); 
xlabel('t [s]');ylabel('\rho'); title('P_{aux,electrons} [W/m^3]'); shading flat;
subplot(212); surf(out.time,out.rho,out.jaux); 
xlabel('t [s]');ylabel('\rho'); title('j_{aux} [A/m^2]'); shading flat;

%%
% The attentive reader will notice that the driven current increases then
% EC heating is added. This is due to the the ECCD current drive efficiency
% which increases with the higher temperature.

%% Example 2: Time-varying rdep and add NBI
% We will now let the deposition locations be specified through the input
% vector. We also add NBI.
clear; close all;
addpath(fullfile(pwd,'..','code'),'-end'); % add RAPTOR path
[config] = RAPTOR_config; % load default config
config.grid.tgrid = [0:0.001:0.1]; % time grid
config.grid.rhogrid = linspace(0,1,21); % spatial grid
config.debug.verbosity = -1;
config.debug.iterplot = 0; %10; % plot profiles every N iterations
config.debug.iterdisp = 10; % display progress every N iterations

% heating & current drive parameters

%% 
% *EC with time-varying deposition location*
% To define the deposition location as time-varying, set rdep to -1. Then
% the actuator index assigning the deposition locations are appended to
% uindex.

% set ECH/ECCD module configuration
config.echcd.Configuration.n_units = 2;
config.echcd.Configuration.n_rdep_external = 2;

%%
% *NBI*
% Configure the NBI module: use a generic hcd_gaussian module.
config.nbhcd = RAPTORmodule('hcd_gaussian');

% set NBI parameters
config.nbhcd.Configuration.n_units = 2;
config.nbhcd.Configuration.n_rdep_external = 0;

%%
% *Create model*
% rerun the config file with the new parameter set
[model,params,init,g,v] = build_RAPTOR_model(config);


%% Parameters

params.echcd.wdep = [0.3, 0.35]; % deposition width
params.echcd.cd_eff = [0 1]; % CD efficiency

params.nbhcd.active = true;
params.nbhcd.rdep = [0 0.4]; %
params.nbhcd.wdep = [0.8 0.4]; % broad heating
params.nbhcd.wdep_out = [0.8 0.6]; % broad heating
params.nbhcd.cd_eff = [0 0]; % no current drive
params.nbhcd.uindices = uint32([6 7]); % index in input vector

%%
% *Define inputs*:
% Must define additional inputs for rdep and for the NBI power.

% plasma current
U(1,:)= 200e3*ones(size(params.tgrid)); % input Ip trace: ramp from 80kA to 200kA
U(1,params.tgrid<0.02) = linspace(80e3,200e3,sum(params.tgrid<0.02));

% indices for EC power
iupow = model.echcd.uind(model.echcd.ind_power); % indices to put power
iurhodep = model.echcd.uind(model.echcd.ind_rhodep); % indices to put rhodep

% first EC actautor:  1MW starting at 60ms
U(iupow(1),:) = zeros(size(params.tgrid)); U(2,params.tgrid>0.06) = 1e6;
% second EC actuator power:  1MW start at 40ms
U(iupow(2),:) = zeros(size(params.tgrid)); U(3,params.tgrid>0.04) = 1e6;

% deposition location for first EC actuator: on-axis with small movement
U(iurhodep(1),:) = zeros(size(params.tgrid)); 
U(iurhodep(1),params.tgrid>0.08) = linspace(0,0.2,sum(params.tgrid>0.08));

% deposition location for first EC actuator: off-axis with sweep inward
U(iurhodep(2),:) = 0.4*ones(size(params.tgrid)); 
U(iurhodep(2),params.tgrid>0.05) = linspace(0.4,0.2,sum(params.tgrid>0.05));

% NBI power: ramp slowly
iupow = model.nbhcd.uind;
U(iupow(1),:) = 1e6*ones(size(params.tgrid)); 
U(iupow(1),params.tgrid<0.06) = linspace(0,500e3,sum(params.tgrid<0.06));

% Second NBI source: off
U(iupow(2),:) = 0e6*ones(size(params.tgrid)); 

%%
% initial condition
init.Ip0 = U(1,1);
x0 = RAPTOR_initial_conditions(model,init,g,v); 
% Run RAPTOR
simres = RAPTOR_predictive(x0,g,v,U,model,params);
%% 
% *plot results*
out = RAPTOR_out(simres,model,params);
subplot(321); plot(out.time,out.U(1,:)); title('Ip');
subplot(322); plot(out.time,out.U([2:3,6],:)); title('P_{aux}')
subplot(323); plot(out.time,out.U([4:5],:)); title('\rho_{dep}'); set(gca,'Ylim',[0,1]);
subplot(324); pcolor(out.time,out.rho,bsxfun(@times,out.Vp,out.pec)); shading flat; title('p_{EC}*dV/d\rho')
subplot(325); [cs,h] = contour(out.time,out.rho,out.q,[1 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],'color','k'); clabel(cs,h,'labelspacing',72); title('T_e contours [keV]')
xlabel('t [s]'); ylabel('\rho');

%% Example 3: Completely manual specification of heating & current drive
clear; close all;
[config] = RAPTOR_config; % load default config
config.grid.tgrid = [0:0.001:0.1]; % time grid
config.grid.rhogrid = linspace(0,1,21); % spatial grid
config.debug.iterplot = 10; % plot profiles every N iterations
config.debug.iterdisp = 10; % display progress every N iterations

% Manual heating and current drive specification
config.echcd = RAPTORmodule('none');
config.nbhcd = RAPTORmodule('hcd_manual'); % define the manual module
rho = linspace(0,1,21); % arbitrary rho grid
disp(config.nbhcd)

% *Create model*
% rerun with the new parameter set
[model,params,init,g,v,U] = build_RAPTOR_model(config);
disp(model.nbhcd)
disp(params.nbhcd)
%%
% We now set the nbh deposition parameters manually: 
% specify them on rhogauss
params.nbhcd.Pauxe  = 1e6*(1-model.rgrid.rhogauss); % heating profile to electrons
params.nbhcd.Pauxi  = 0.8*(1-model.rgrid.rhogauss); % heating profile to ions (ignored anyway since we don't solve for ions now)
params.nbhcd.jaux  = zeros(model.rgrid.nrhogauss,1);
params.nbhcd.pfast = zeros(model.rgrid.nrhogauss,1); % placeholders for future fast ion things

init.te0 = 300;

% define inputs
U(1,:)= 200e3*ones(size(params.tgrid)); % input Ip trace: ramp from 40kA to 200kA
U(1,params.tgrid<0.02) = linspace(80e3,200e3,sum(params.tgrid<0.02));

% initial conditions
init.Ip0 = U(1,1);
x0 = RAPTOR_initial_conditions(model,init,g,v);

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

% *plot outputs*
out = RAPTOR_out(simres,model,params);
figure;
subplot(221); surf(out.time,out.rho,out.pauxe); 
xlabel('t [s]');ylabel('\rho'); title('P_{aux,electrons} from NBI [W/m^3]'); shading flat;
subplot(222); surf(out.time,out.rho,out.jaux); 
xlabel('t [s]');ylabel('\rho'); title('j_{aux} NBI [A/m^2]'); shading flat;
subplot(223); plot(out.time,[out.Pnbe;out.Pnbi]/1e6); 
xlabel('t [s]');ylabel('[MW]'); title('Total power to electrons (-), ions (--)');

return