%% RAPTOR tutorial 6: Simulating other species
% The default setting in RAPTOR is to simulate only flux diffusion and
% electron temperature evolution. Here we explain how to add also an
% equation for evolving ion temperature and/or electron density.
%%
run(fullfile(pwd,'..','RAPTOR_path.m')); % add RAPTOR path
close all hidden; clear;

%% Adding Ti equation

% Load some defaults
[config] = RAPTOR_config; % load default params
config.debug.iterdisp = 10;
config.grid.tgrid = [0:2e-3:0.4];
config.echcd.Configuration.n_units = 1;
config.nbhcd.Configuration.n_units = 1;

%% 
% Use BgB transport model since that outputs ion and electron transport
% coefficients separately
config.chi_e = RAPTORmodule('chi_BgB');

%% 
% Ti defined as independent state to be solved for
config.ti.method = 'state';
[model,params,init,g,v,U] = build_RAPTOR_model(config); % generate model structure for these params.
params.pei.active = true; % source electron-ion equipartitions, switch it on.

%%
% Inputs
params.echcd.active = true;
params.nbhcd.active = true;
params.nbhcd.frac_to_electrons = 0.5; 

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

U(1,:) = rampfun(params.tgrid,0,80e3,0.2,150e3); % input Ip trace: constant 80kA
U(model.echcd.uind,params.tgrid>=0.1 & params.tgrid<=0.15) = 500e3; % briefly add some ECH (electron only)
U(model.nbhcd.uind,:) = rampfun(params.tgrid,0.2,0e3,0.25,500e3); % input NBI: heat those ions
init.Ip0 = U(1); % initial condition for current
x0 = RAPTOR_initial_conditions(model,init,g(:,1),v(:,1)); % Define the initial condition

%%
% run RAPTOR and get outputs
simres = RAPTOR_predictive(x0,g,v,U,model,params);
out = RAPTOR_out(simres,model,params);

%%
% Plots
subplot(211)
plot(out.time,out.Ip(end,:)/1e6,'b',...
out.time,out.Pec/1e6,'k--',...
out.time,out.Pnbe/1e6,'k',...
out.time,out.ne(1,:)/1e19,'m'); legend('I_p [MA]','P_{EC} [MW]','P_{NB,el} [MW]','n_{e0,19}');
subplot(212)
plot(out.time,out.te(1,:),'b',out.time,out.ti(1,:),'r');legend('T_{e0}','T_{i0}');

%% Particle density evolution

%% Solve for ne with manual transport coefficients, fix ni/ne fraction
% Load the default configuration structure.
[config] = RAPTOR_config; % load default params

% To solve ne equation, set these options:
config.ne.method = 'state';

config.ni.method = 'nescal'; % also specify ion density method: scale w.r.t. ne

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

% inspect transport model used for ne transport coefficieents;
disp(model.vpdn_e); % manual model
disp(params.vpdn_e);
%%
init.ni0scal = 0.9; init.niescal = 1; % scaling factor (core,edge)
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);
out = RAPTOR_out(simres,model,params);
% Plots
hf=figure(1); set(hf,'position',[50 50 700 300]); clf(hf);
subplot(131); plot(out.time,out.ne/1e19); xlabel('time [s]'); title('n_e at various rho [10^{19}m^{-3}]');
subplot(132); plot(out.rho,out.ne(:,1:5:end)/1e19); xlabel('\rho_{tor}'); title('n_e profiles [10^{19}m^{-3}]');
subplot(133); plot(out.time,out.ne(1,:)/1e19,'b',out.time,out.ni(1,:)/1e19,'r--'); xlabel('time [s]'); title('ne(-),ni(--)');

%% Solve for ne with manual transport coefficients, compute ni from Zeff and quasineutrality
% To solve for ni via Zeff and quasineutrality, we need to introduce one
% impurity species (since we have two constraints).

% Change the method for n to 'qnze' (quasi-neutrality and Z-effective).
config.ni.method = 'qnze';
config.n1.method = 'qnze'; % solve for one impurity species as well.
config.atom.Z1 = 6; config.atom.A1 = 12; % choose carbon impurity

[model,params,init,g,v,U] = build_RAPTOR_model(config);
init.ze0 = 2.0; init.zee = 3.0; % set Zeff at axis and edge.
v = build_kinetic_profiles(model,init); % update Zeff in v
x0 = RAPTOR_initial_conditions(model,init,g(:,1),v(:,1)); 

% Run RAPTOR
simres = RAPTOR_predictive(x0,g,v,U,model,params);
% To compute some physical outputs, we run the RAPTOR_out function.
out = RAPTOR_out(simres,model,params);

hf=figure(1); set(hf,'position',[50 50 700 300]); clf(hf);
subplot(131); plot(out.rho,[out.ne(:,end),out.ni(:,end)]/1e19); xlabel('\rho_{tor}'); title('n_e,n_i [10^{19}m^{-3}]');
subplot(132); plot(out.rho,out.ze,'r',out.rho,out.ne(:,end)./out.ni(:,end)); xlabel('\rho_{tor}'); title('Z_{eff}, n_i/n_e');
subplot(133); plot(out.time,out.ne(1,:)/1e19,'b',out.time,out.ni(1,:)/1e19,'r--'); xlabel('time [s]'); title('ne(-),ni(--)');
