%% Tutorial for time-varying geometry simulation
% Simulation with time-varying geometry is demonstrated in case of plasma
% elongation increase from 1.1 to 1.7. 
% Safety factor q95 increase is expected. 
% Two CHEASE data files are used to get equilibrium data for low and high
% plasma elongation.
%
% add RAPTOR path
close all hidden;
run(fullfile(pwd,'..','RAPTOR_path.m'));

%% Model construction and initialization
% Load default parameters
config = RAPTOR_config('ITER'); 
% Specify input data file(s)
% One can specify several input data file(s) corresponding to plasma
% equilibrium at different time steps. In case of several data
% files config.equi.filenames contains a cell array.
%% get equilibrium from IMAS
config.equi.source = 'IMAS';
config.equi.filenames = {'/equils/IMAS/METIS_4321_run2'};

%% Define time grid: time interval from 0.0 to 450.0 ms
config.grid.tgrid = [0:0.005:0.45]; 

% Build model and get params, init structures, geometry g, kinetic profiles
% v, actuator inputs U

try
[model,params,init,g,v,U] = build_RAPTOR_model(config);
catch ME
  if strcmp(ME.identifier,'RAPTOR:FileDoesNotExist')
    fprintf('Could not complete tutorial, IMAS data file does not exist.\n')
    warning(getReport(ME));
    return
  end
end


%%
% Define initial conditions for plasma state
x0 = RAPTOR_initial_conditions(model,init,g(:,1),v(:,1)); 

%% Run RAPTOR simulation   
simres = RAPTOR_predictive(x0,g,v,U,model,params);
out = RAPTOR_out(simres,model,params);
%% Plots
subplot(211)
plot(out.time,[out.q95;out.q(end,:);out.Ip(end,:)/1e6]);
xlabel('time, sec'); legend({'q_{95}','q_{edge}','Ip [MA]'},'location','northwest');
subplot(212)
plot(out.time,[out.te(1,:);out.Wth]);
xlabel('time, sec'); legend({'Te0','W_{e}'},'location','northwest');

