%% Tutorial running RAPTOR using IMAS IDS inputs

%error('this tutorial is not yet ready for use')

% add RAPTOR path

close all hidden;
%load ../data/raptor_ws

%% IMAS shot and run number
%shot = 110005;
shot = 100002;
run = 1;

if exist('ids_get_slice','file')
    inIMAS = true;
else
    inIMAS = false;
end

%% Initializations
% RAPTOR time grid
%timegrid = [300:310];
%timegrid = [50:50:647];
%timegrid = [5:0.5:50];

%% Read IDSs from local database
if ~exist('imas_open','file');
  warning('imas_open.m not found. Can not complete tutorial. Please check IMAS installation\n');
  return
end

fprintf('opening IMAS database\n')
idx_in = imas_open('ids',shot,run);

%% To get time vector
if ~exist('cp_tot')
  fprintf('get time vector\n')

  cp_tot = ids_get(idx_in,'core_profiles');
  ps_tot = ids_get(idx_in,'pulse_schedule');
  %timegrid = cp_tot.time(2:2:30);
end
%timegrid=[10:1:20]
%timegrid=[cp_tot.time(2):1:cp_tot.time(end)];
timegrid=[cp_tot.time(2):1:80];

%% Load default configuration structure for ITER (default settings)
[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.

%% Define time grid: time interval from 0.0 to 450.0 ms

% Metis resolution is 21 so Raptor resolution needs to be lower
config.grid.rhogrid = linspace(0,1,16);
config.grid.tgrid = 0; % just one time
config.ti.method = 'direct'; % By default Ti would be a scaling from Te => Here we directly give it explicitly instead

config.echcd = echcd('echcd_gaussian');
config.ichcd = ichcd('ichcd_manual');
config.nbhcd = nbhcd('nbhcd_manual');
config.lhhcd = lhhcd('lhhcd_manual');

% initial initialization of hcd configuration
ec = config.echcd.params;
nb = config.nbhcd.params;
ic = config.ichcd.params;
lh = config.lhhcd.params;
if inIMAS
core_sources  = ids_get_slice(idx_in,'core_sources',timegrid(1),1); % needed once outside the time loop
core_profiles = ids_get_slice(idx_in,'core_profiles',timegrid(1),1);
end
[config.echcd.params,config.nbhcd.params,config.lhhcd.params,config.ichcd.params] = get_HCD_from_IDS(core_sources,ec,nb,ic,lh);

config.init.Ip0 = core_profiles.global_quantities.ip;

%% Load default configuration structure
[model,params,init,~,~,U] = build_RAPTOR_model(config);

%% Time loop
for it=1:numel(timegrid)
    if it < numel(timegrid)
        params.tgrid = timegrid(it:it+1);
    else
        % extrapolate
        params.tgrid = timegrid(it)+[0,diff(timegrid(it-1:it))];
    end
    %% get time slices
    t_now = timegrid(it);
    if inIMAS
        core_profiles  = ids_get_slice(idx_in,'core_profiles',t_now,1);
        equilibrium    = ids_get_slice(idx_in,'equilibrium',t_now,1);
        pulse_schedule = ids_get_slice(idx_in,'pulse_schedule',t_now,1);
        core_sources   = ids_get_slice(idx_in,'core_sources',t_now,1);
    end
    %% Heating & Current Drive sources
    ec = params.echcd;
    nb = params.nbhcd;
    ic = params.ichcd;
    lh = params.lhhcd;
    [params.echcd,params.nbhcd,params.lhhcd,params.ichcd] = get_HCD_from_IDS(core_sources,ec,nb,ic,lh);
    
    %% Read Ip waveform from pulse_schedule
    U(1,1) = pulse_schedule.flux_control.i_plasma.reference.data;
    
    %% Read waveforms from input pulse_schedule IDS
    U = get_HCD_waveform_from_IDS(pulse_schedule,params,U);
    
    %% Read input kinetic profiles from input core_profiles IDS
    v = read_core_profiles_ids(model,params,core_profiles);
    %%
    g = read_geometry_ids(equilibrium,model.geom);
    
    %% Define initial conditions for x state (only for 1st time step)
    if it==1
      x = RAPTOR_initial_conditions(model,init,g(:,1),v(:,1)); 
    else
      x = xprev;
    end
   
    %% 
    RAPTOR_plot_xgv(x,g,v,model);% check inputs going into RAPTOR
    %% Run RAPTOR
    simres_t = RAPTOR_predictive(x,g,v,U*[1 1],model,params);
    xprev = simres_t.X(:,2);
    assert(~any(isnan(xprev)),'RAPTOR time step failed, exitflag=%d',simres_t.exitflags(end-1))
    
    %% store results
    simres.X(:,it) = simres_t.X(:,1);
    simres.G(:,it) = simres_t.G(:,1);
    simres.V(:,it) = simres_t.V(:,1);
    simres.Xdot(:,it) = simres_t.Xdot(:,1);
    simres.U(:,it) = simres_t.U(:,1);
    simres.staps{it} = simres_t.staps{1};
    simres.geops{it} = simres_t.geops{1};
    simres.traps{it} = simres_t.traps{1};
    simres.exitflags{it} = simres_t.exitflags(1);
    simres.newtres(it) = simres_t.newtres(1);
    simres.df_dx{it} = simres_t.df_dx{1};
    simres.df_dxdot{it} = simres_t.df_dxdot{1};
    simres.df_dxk{it} = simres_t.df_dxk{1};
    simres.df_dxk1{it} = simres_t.df_dxk1{1};
    simres.df_du{it} = simres_t.df_du{1};
    simres.df_dd{it} = simres_t.df_dd{1};

end %% End time loop

%% Close file of local database
if inIMAS
imas_close(idx_in);
end

%% Compute some Physical outputs
params.tgrid = timegrid; % restore true timegrid
out = RAPTOR_out(simres,model,params);

%% Plots
subplot(311)
plot(ps_tot.time,(ps_tot.ec.antenna{1}.power.reference.data+ps_tot.ec.antenna{2}.power.reference.data)*1.e-6)
xlabel('time, sec'); legend({'EC power waveform (MW)'},'location','northwest');
set(gca,'xlim',[5. 100.]);
set(gca,'ylim',[0 30])
nice
subplot(312)
plot(out.time,[out.q95;out.q(end,:);out.Ip(end,:)/1e6]);
xlabel('time, sec'); legend({'q_{95}','q_{edge}','Ip [MA]'},'location','northwest');
set(gca,'xlim',[5. 100.]);
nice
subplot(313)
plot(out.time,[out.te(1,:);out.Wth/1e3]);
xlabel('time, sec'); legend({'Te0','W_{e}'},'location','northwest');
set(gca,'xlim',[5. 100.]);
nice

