%% RAPTOR tutorial: ad-hoc transport model tuning
% In this tutorial we will vary the parameters of the 
% electron transport model to see their effect on the plasma.
%
% We use the heating and current drive setup from tutorial 3 again, but
% modify some parmeters: in particular we choose a smaller time grid to
% handle the rapid rise of the pedestal.
%%

close all hidden;
clear; run('../RAPTOR_path.m');

[config] = RAPTOR_config; % load default config
config.grid.tgrid = 0:0.001:0.2;
config.grid.rhogrid = linspace(0,1,21); % spatial grid
config.debug.iterdisp = 20; % display progress every N iterations

config.echcd.Configuration.n_units = 1;
config.nbhcd.Configuration.n_units = 1;

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

% set ECH/ECCD parameters
params.echcd.active = true;
params.echcd.rdep = [0.4]; %
params.echcd.wdep = [0.35]; % deposition width
params.echcd.cd_eff = [1]; % CD efficiency

% set NBI parameters
params.nbhcd.active = true;
params.nbhcd.rdep = [0]; %
params.nbhcd.wdep = [0.8]; % broad heating
params.nbhcd.cd_eff = [0]; % no current drive

% 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));

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

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

% Run RAPTOR
[simres] = RAPTOR_predictive(x0,g,v,U,model,params);
out = RAPTOR_out(simres,model,params);

%% Changing parameters for the default ad-hoc transport model
% RAPTOR includes an ad-hoc expression for \chi_e
% 
% $\chi_e = \chi_{neo} + c_{ano}\rho q F(s)H +
% \chi_{central}\exp(-\rho^2/\delta_{central})$
%
% where 
%
% $F(s) = \frac{a_{ic}}{1+\exp(w_{ic}(d_{ic}-s))} + 1 - a_{ic}$
%
% and $H$ is an extra factor to reduce the transport at the edge to create
% an H-mode pedestal.
%
% This model is written as
% 
% * $\chi_{prof} = \exp(-(\rho-\rho_{\chi,edge})^2/w_{\chi,edge}.^2)$
% * $H = a_{LH}[1-1/(1+\exp(w_{LH}(d_{LH}-LH_{crit}))]\chi_{prof}(\rho)$
% * $LH_{crit} = -\frac{\partial p_{e}}{\partial\rho}$
%
% The relevant parameters can be set from the params structure
params.chi_e

%%
% For example, to remove the confinement enhancement at low shear we can
% set $a_{ic}=0$
params.chi_e.aitb = 0;

% simulate
[simres_noenhance] = RAPTOR_predictive(x0,g,v,U,model,params);
out_noenhance = RAPTOR_out(simres_noenhance,model,params);
RAPTOR_plot_GUI({out,out_noenhance},'overview')

%% Bohm Gyrobohm
% We can also use the Bohm-gyroBohm model
% first initialize new default parameters
config.chi_e = RAPTORmodule('chi_BgB');

% make new model and parameters
[model,params,init,g,v,U] = build_RAPTOR_model(config);
%
% we now have a new set of parameters to change, but leave them fixed for now
params.chi_e

% run simulation with new model
params.debug.iterdisp = 10; % plot progress to see pedestal growing
[simres_BgB] = RAPTOR_predictive(x0,g,v,U,model,params);
out_BgB = RAPTOR_out(simres_BgB,model,params);
RAPTOR_plot_GUI({out_noenhance,out_BgB},'overview')

%% Manual chie
% We can also choose a fixed chi_e profile manually, by setting
% |'modeltype'| to |'manual'| and assigning the profile to |'cneo'|.

config.chi_e = RAPTORmodule('chi_manual');

[model,params,init,g,v,U] = build_RAPTOR_model(config);
disp(params.chi_e);

params.chi_e.chi_e = 10; % constant value
params.debug.iterplot = 0; % no plots

[simres_manual1] = RAPTOR_predictive(x0,g,v,U,model,params);
out_manual1 = RAPTOR_out(simres_manual1,model,params);

params.chi_e.chi_e = (5+10*model.rgrid.rhogauss.^2); % spatially varying profile
params.debug.iterplot = 0; % no plots

[simres_manual2] = RAPTOR_predictive(x0,g,v,U,model,params);
out_manual2 = RAPTOR_out(simres_manual2,model,params);

%% plot some comparisons
figure; clf;
subplot(211)
plot(out_manual1.rho,out_manual1.te(:,end)/1e3,'r',...
     out_manual2.rho,out_manual2.te(:,end)/1e3,'k'...
     ); ylabel('T_e [keV]'), xlabel('\rho')

subplot(212)
plot(out_manual1.rho,out_manual1.chie(:,end),'r',...
     out_manual2.rho,out_manual2.chie(:,end),'k'...
     ); ylabel('\chi_e'), xlabel('\rho')
shg;
%% Time-varying Manual Chie
%
% It is also possible to assign to |'cneo'| a time-varying cie_e profile. 
% This profile should be provided as a matrix of dimension:$\rho_{gauss}$
% by t. Here is an example of a profile which is based on a parabolic
% function with time dependent gradient. A default value for tgrid is
% used of dimension 2 by 1 so that the previous examples hold because
% these have dimension # of rhogauss by 1.
%

% Create a matrix of the right size with a parabol which gradient increases
% in time:

[X,Y] = meshgrid(25*params.tgrid,model.rgrid.rhogauss);
params.chi_e.chi_e = 5 + 2.^(X).*Y.^2;

[simres_manual_time] = RAPTOR_predictive(x0,g,v,U,model,params);
out_manual_time = RAPTOR_out(simres_manual_time,model,params);

% Plot the result
subplot(211)
plot(out_manual_time.rho,out_manual_time.te(:,end)/1e3,'r');
    ylabel('T_e [keV]'), xlabel('\rho')
    
subplot(212) % Chi at the beginning and the end
plot(out_manual_time.rho,out_manual_time.chie(:,end),'r',...
     out_manual_time.rho,out_manual_time.chie(:,1),'b'...
    ); ylabel('\chi_e'), xlabel('\rho')

return


%% Analytical formula for chi_e with caracteristic scalelength
% L-mode simulation
config.chi_e = chi_e('MS');
% A new set of parameters for chi_e module
config.chi_e.params

% Time grid
config.grid.tgrid = [0:0.001:0.15];
% Define reference He factor(==tauEe/tauE_H98) on RAPTOR time grid
config.chi_e.params.hh = 0.2*ones(size(config.grid.tgrid));
% Make new model and parameters
[model,params,init,g,v,U] = build_RAPTOR_model(config);

% Run simulation with new model
simres_MS = RAPTOR_predictive(x0,g,v,U,model,params);
out_MS = RAPTOR_out(simres_MS,model,params);

% Profile of electron temperature
plot(out_MS.rho,out_MS.te(:,end),'m','LineWidth',2); hold on
% Define inded for rho_inv on gauss grid
[~,ind_rinv] = min((out_MS.rhogauss - params.chi_e.rhoinv).^2);
% rho_ped index
[~,ind_rped] = min((out_MS.rhogauss - params.chi_e.rhoped(1)).^2);
% Te on gauss grid
tegauss = eval_te(simres_MS.X(:,end),simres_MS.G(:,1),simres_MS.V(:,1),model,'gauss');
% dlnTe/drho
dlnTe_drhogauss = bsxfun(@rdivide,diff(log(tegauss(:,end))),diff(out_MS.rhogauss(:,end)));
dlnTe_drhogauss = [dlnTe_drhogauss' dlnTe_drhogauss(end)]';
% Get lambdaTe value from simulation
slambdaTe = -dlnTe_drhogauss(ind_rinv:ind_rped,end);
% Get muTe value from simulation
smuTe = -tegauss(ind_rped:end,end).*dlnTe_drhogauss(ind_rped:end,end);
muTe = mean(smuTe);

% Add core part
coeffl = tegauss(ind_rped)/exp(-params.chi_e.lambdaTe(1)*out_MS.rhogauss(ind_rped));
te_core = coeffl*exp(-params.chi_e.lambdaTe(1)*out_MS.rhogauss(ind_rinv:ind_rped));
plot(out_MS.rhogauss(ind_rinv:ind_rped),te_core,'b--','LineWidth',2);

% Add pedestal part
coeffm = tegauss(end) + muTe*out_MS.rhogauss(end);
te_edge = coeffm - muTe*out_MS.rhogauss(ind_rped:end);
plot(out_MS.rhogauss(ind_rped:end),te_edge,'r--','LineWidth',2);    

xlabel('rho');legend('Te\_prof','labmdaTe','muTe','Location','SouthWest');
title('Te profile');