%% Tutorial for simulation with ASTRA data
% Simulation with equilibrium data produced by ASTRA code is demonstrated.
% Value of plasma current Ipl, radial profiles for ion temperature Ti,
% electron density ne, ion density ni, effective charge Zeff, electron heat
% conductivity chi_e and EC heating echcd are specified according to ASTRA
% data to get comparison with ASTRA code simulation results.

% add RAPTOR path
close all hidden;
clear; run(fullfile(pwd,'..','RAPTOR_path.m'));

%% Define config
% In this section default RAPTOR config is generated. Source of equilibrium
% data is changed to 'astra'. Sources of chi_e and echcd parameters are changed to
% 'manual' to get later values from ASTRA data file.

% Load default parameters
config = RAPTOR_config;
% Specify data source
config.equi.source = 'astra';
% Specify input data file
config.equi.filenames = {'ASTRA_data_fixed_geom'};
% For fixed case specify input time for ASTRA data (sec)
% negative value - first time step
% inf - last time step
config.equi.times = 1.5;
% Specify chi_e source 
config.chi_e = RAPTORmodule('chi_manual');
% Specify echcd source
config.echcd = RAPTORmodule('hcd_manual');
% Specify time grid
config.grid.tgrid = [0.0:0.05:2.0];
% Specify way of ti calculation: independent of te
config.ti.method = 'direct';

%% Build model
% RAPTOR model is constructed with ASTRA equilibrium data. Parameter g
% contains ASTRA geometric profiles. Additional parameters are defined 
% according to ASTRA data.

% Generate model, params, init, geometry g, plasma characteristics v,
% actuator inputs U.
[model,params,init,g,v,U] = build_RAPTOR_model(config);
% Get ASTRA parameters: astra_profiles contains profiles for ne, ni, ti,
% ze; chiegauss and pechgauss - electron heat conductivity and ECH power
% density from ASTRA data recalculated to RAPTOR gauss grid; boundary
% condition for electron temperatrure Te.
[astra_profiles,chiegauss,pechgauss,~,~,~,~,te_BC_astra] = get_ASTRA_profiles(config,model);
% Boundary conditions for Te
model.te.BC.defaultValue = te_BC_astra;
% Construct kinetic profiles v with ASTRA data
v = build_kinetic_profiles(model,config,astra_profiles);
% chi_e
params.chi_e.chi_e = chiegauss;
params.chi_e.implicit = 1;
% ECH&CD
params.echcd.Paux = pechgauss;
params.echcd.jaux = zeros(size(params.echcd.Paux));
params.echcd.active = true;
params.echcd.rdep = 1;

%% Define actuator inputs
% Plasma current
U(1,:) = astra_profiles.Ip;

%% Initial conditions for state
% Plasma current
init.Ip0 = U(1,1);
x0 = RAPTOR_initial_conditions(model,init,g(:,1),v(:,1));

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

%% Plot the results
% RAPTOR and ASTRA results comparison

% Full path to ASTRA data file
astrafilepath = fullfile(config.envopts.RAPTOR_path,'equils', 'ASTRA',config.equi.filenames{1});
fprintf('loading ASTRA equilibrium from %s\n',astrafilepath);
astra = load(astrafilepath);
% Normalized radial grid (rho toroidal normalized)
astra_grid = bsxfun(@rdivide,astra.out.RHO,astra.out.RHO(end,:));
astra_grid(1,:) = 0.0;
figure;
subplot(121);
plot(astra_grid(:,end),astra.out.FP(:,end)-astra.out.FP(1,end),out.rho,out.psi(:,end)-out.psi(1,end),'ro--','LineWidth',4,'MarkerSize',5);
ylim([0.0,0.3]);
set(gca,'FontSize', 23);xlabel('\rho_{norm}');ylabel('\psi');
legend('ASTRA','RAPTOR');
subplot(122);
plot(astra_grid(:,end),astra.out.TE(:,end),out.rho,out.te(:,end)*1e-03,'ro--','LineWidth',4,'MarkerSize',5);
set(gca,'FontSize', 23);xlabel('\rho_{norm}');ylabel('T_{e}, keV');