%% RAPTOR tutorial: radiation module from ADAS data
% In this tutorial central ECRH is used with BgB transport model
% for the electron heat transport
%
% Three impurities are considered from light to heavy
% Total radiations computed with corresponding cooling factors
% which depend on Te

close all hidden;
clear;
run(fullfile(pwd,'..','RAPTOR_path.m'));

% To create a matlab structure with cooling factors and average charge for a list of impurities:
%  - Specify the list of imputities in 'chargemass.m' 
%  - Use only 'make_Lz_table_with_H0.m' which uses also neutral deuterium for computations including charge exchange

%  make_Lz_table_with_H0.m

% It will create a structure and a corresponding matlab function with writestruct
% A default matlab function 'Lz_zave_HDT0_default.m' is used in 'prad_adas.m' and 'qnze.m'

[config] = RAPTOR_config; % load default config
config.grid.tgrid = 0:0.005:0.05;
config.grid.rhogrid = linspace(0,1,21); % spatial grid
config.debug.iterdisp = 20; % display progress every N iterations
config.chi_e = RAPTORmodule('chi_BgB');

config.echcd.Configuration.n_units = 1;

config.ni.method = 'qnze'; % quasineutrality
config.n1.method = 'nescal'; 
config.n2.method = 'nescal'; 
config.n3.method = 'nescal'; 
config.ze.method = 'qnze'; % quasineutrality

config.atom.Z1=7;
config.atom.A1=14.007;
config.atom.Z2=29;
config.atom.A2=63.546;
config.atom.Z3=74;
config.atom.A3=183.84;

config.prad=RAPTORmodule('prad_adas');

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

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

% Plasma current
U(1,:)= 200e3*ones(size(params.tgrid)); % input Ip 
U(1,params.tgrid<0.02) = linspace(80e3,200e3,sum(params.tgrid<0.02));

% EC power ramp (one can track the W radiation peaking at 1.5 keV)
U(2,:) = zeros(size(params.tgrid)); U(2,params.tgrid>0.02) = linspace(0,2e6,sum(params.tgrid>0.02));

% Initial condition
init.Ip0 = U(1,1);
init.n1scal = 0.01; 
init.n2scal = 0.001; 
init.n3scal = 1e-3;  % Tungsten contamination

v = build_kinetic_profiles(model,init);
x0 = RAPTOR_initial_conditions(model,init,g(:,1),v(:,1)); 

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

% Shows the radiated power against the electron temperature
% The peak of W radiation is not exactly around 1.5 keV due
% to peaked density profiles. In contrast, the W cooling factor
% also plotted, behaves as expected from ADAS data

figure
plot(out.te,out.prad,'-','Linewidth',2)
set(gca,'FontSize',20)
xlabel('T_e (eV)')
ylabel('P_{rad} (W/m^3)')
  
figure
plot(out.te,out.prad./out.ne./out.n3,'x')
set(gca,'FontSize',20)
xlabel('T_e (eV)')
ylabel('L_Z (W m^3)')
