# RAPTOR Tutorial 2: Change input, display control and output analysis

In this second tutorial we will customize how the simulation results are presented during the run, and will make some figures after the run. This will help underestand how profile data is stored and extracted in RAPTOR. We will also modify the input signal to have something slightly more interesting.

## Contents

- Preparation
- Customize display options
- Run RAPTOR with new settings
- Customized display of results
- Interactive plotting tool

## Preparation

Start by repeating the first two steps of the previous tutorial

    clear; close all hidden;
    run(fullfile(pwd,'..','RAPTOR_path.m')); % add RAPTOR path
    [config] = RAPTOR_config; % load default config

    % build model for these params.
    [model,params,init,g,v,U] = build_RAPTOR_model(config);

## Customize display options

Now we will customize the display options in the params structure to suit our needs

    params.debug.iterdisp = 5; % display every 5 iterations only (default=1)
    params.debug.iterplot = 1; % plot profile evolution every iteration (default=0, no plot)

## Run RAPTOR with new settings

    % Choose a more interesting input: modify U
    nt = numel(params.tgrid);
    % the first row of U is the plasma current as a function of time. Assign
    % it here:
    U(1,params.tgrid <0.02) = linspace(80e3,250e3,sum(params.tgrid<0.02));
    U(1,params.tgrid>=0.02) = 250e3;

    % Update init to have correct plasma current at first time step
    init.Ip0 = U(1,1);
    % Define the initial condition
    x0 = RAPTOR_initial_conditions(model,init,g,v);

    simres = RAPTOR_predictive(x0,g,v,U,model,params);

    % For convenience, a genetic function to directly extract interesting data is available:
    out = RAPTOR_out(simres,model,params);
    disp(out);

    figure(findall(0,'name','RAPTOR step plot')); % to bring figure forward

         it  telaps newt     res   t[ms]  dt[ms]  Ip[kA] Icd[kA] Ibs[kA] Ioh[kA]      qe    qmin      q0   Vl[V] Te0[keV] Ti0[keV] ne0[e19]   f_ss
          1    0.56    4 5.7e-10       0       1      80       0    1.41    78.5    15.3    5.55    5.55 3.4e+00    0.21    0.21    1.00 1.1e+01
          6     2.2    2 9.7e-09       5       1     125       0    3.08     122    9.84    3.99    3.99 6.1e+00    0.31    0.31    1.00 2.1e+01
         11     3.8    2 1.2e-09      10       1     169       0    4.76     165    7.24    3.31    3.31 5.9e+00    0.44    0.44    1.00 2.2e+01
         16     5.3    2 2.6e-10      15       1     214       0    6.31     208    5.73    2.84    2.84 5.7e+00    0.59    0.59    1.00 2.2e+01
         21     6.9    3 6.7e-12      20       1     250       0    7.65     242    4.91    2.46     2.5 4.3e+00    0.76    0.76    1.00 1.7e+01
         26     8.8    2 9.2e-10      25       1     250       0     8.3     242    4.91     2.1    2.24 3.2e+00    0.90    0.90    1.00 1.3e+01
         31      10    2 6.6e-11      30       1     250       0    8.39     241    4.91    1.87    2.02 2.8e+00    1.01    1.01    1.00 1.1e+01
         36      12    2 9.8e-12      35       1     250       0    8.27     242    4.91    1.69    1.81 2.6e+00    1.08    1.08    1.00 9.8e+00
         41      13    2 2.4e-12      40       1     250       0     8.1     242    4.91    1.53    1.61 2.4e+00    1.14    1.14    1.00 8.8e+00
         46      15    2 8.1e-13      45       1     250       0    7.93     242    4.91    1.42    1.43 2.2e+00    1.19    1.19    1.00 8.0e+00
         51      17    2 3.3e-13      50       1     250       0    7.79     242    4.91    1.27    1.27 2.1e+00    1.24    1.24    1.00 7.3e+00
                  it: [1x51 double]
                time: [1x51 double]
               ntime: 51
                 rho: [16x1 double]
            rhogauss: [60x1 double]
                nrho: 16
           nrhogauss: 60
                   U: [4x51 double]
                 psi: [16x51 double]
                psiN: [16x51 double]
              rhopol: [16x51 double]
                psip: [16x51 double]
                 upl: [16x51 double]
            dupldrho: [16x51 double]
                f_ss: [1x51 double]
                   H: [16x51 double]
                  g0: [16x51 double]
                  g1: [16x51 double]
                  g2: [16x51 double]
                  g3: [16x51 double]
               g23or: [16x51 double]
                  Vp: [16x51 double]
                   F: [16x51 double]
                Phib: [1x51 double]
               kappa: [16x51 double]
             epsilon: [16x51 double]
                rhob: [1x51 double]
                  te: [16x51 double]
                 te0: [1x51 double]
                 tep: [16x51 double]
                  ne: [16x51 double]
                 nep: [16x51 double]
                  ti: [16x51 double]
                  ni: [16x51 double]
                  n1: [16x51 double]
                  n2: [16x51 double]
                  n3: [16x51 double]
               nitot: [16x51 double]
                  ze: [16x51 double]
                   p: [16x51 double]
                iota: [16x51 double]
                   q: [16x51 double]
               shear: [16x51 double]
              qgauss: [60x51 double]
            psigauss: [60x51 double]
          sheargauss: [60x51 double]
            uplgauss: [60x51 double]
           psiNgauss: [60x51 double]
           psipgauss: [60x51 double]
                  Ip: [16x51 double]
                jtor: [16x51 double]
                jpar: [16x51 double]
                  q0: [1x51 double]
                 q95: [1x51 double]
                qmin: [1x51 double]
            jbsgauss: [60x51 double]
            jecgauss: [60x51 double]
            jnbgauss: [60x51 double]
            jlhgauss: [60x51 double]
            jicgauss: [60x51 double]
           jauxgauss: [60x51 double]
         signeogauss: [60x51 double]
           jpargauss: [60x51 double]
            johgauss: [60x51 double]
                 jbs: [16x51 double]
                jaux: [16x51 double]
                 joh: [16x51 double]
                 jec: [16x51 double]
                 jnb: [16x51 double]
                 jlh: [16x51 double]
                 jic: [16x51 double]
              signeo: [16x51 double]
            pohgauss: [60x51 double]
            pecgauss: [60x51 double]
           pnbegauss: [60x51 double]
           pnbigauss: [60x51 double]
            plhgauss: [60x51 double]
           picegauss: [60x51 double]
           picigauss: [60x51 double]
          pauxegauss: [60x51 double]
          pauxigauss: [60x51 double]
        palphaegauss: [60x51 double]
        palphaigauss: [60x51 double]
          pbremgauss: [60x51 double]
           pradgauss: [60x51 double]
            peigauss: [60x51 double]
           chiegauss: [60x51 double]
           chiigauss: [60x51 double]
                 poh: [16x51 double]
                 pec: [16x51 double]
                pnbe: [16x51 double]
                pnbi: [16x51 double]
                 plh: [16x51 double]
                pice: [16x51 double]
                pici: [16x51 double]
               pauxe: [16x51 double]
               pauxi: [16x51 double]
             palphae: [16x51 double]
             palphai: [16x51 double]
               pbrem: [16x51 double]
                prad: [16x51 double]
                 pei: [16x51 double]
                chie: [16x51 double]
                chii: [16x51 double]
                 dne: [16x51 double]
                 vne: [16x51 double]
                 sne: [16x51 double]
                pein: [16x51 double]
               peout: [16x51 double]
               penet: [16x51 double]
                piin: [16x51 double]
               piout: [16x51 double]
               pinet: [16x51 double]
             tegauss: [60x51 double]
             negauss: [60x51 double]
             tigauss: [60x51 double]
             nigauss: [60x51 double]
                  We: [1x51 double]
                  Wi: [1x51 double]
                 Wth: [1x51 double]
                Wpol: [1x51 double]
                 Ibs: [16x51 double]
                 Iec: [16x51 double]
                 Ilh: [16x51 double]
                 Iic: [16x51 double]
                 Inb: [16x51 double]
                Iaux: [16x51 double]
                 Poh: [1x51 double]
               Pauxe: [1x51 double]
               Pauxi: [1x51 double]
             Pauxtot: [1x51 double]
           Palphatot: [1x51 double]
             Palphai: [1x51 double]
             Palphae: [1x51 double]
               Pbrem: [1x51 double]
                Prad: [1x51 double]
               Pcycl: [1x51 double]
                 Pei: [1x51 double]
                 Pec: [1x51 double]
                Pnbe: [1x51 double]
                Pnbi: [1x51 double]
                Pice: [1x51 double]
                Pici: [1x51 double]
              Volume: [16x51 double]
             ne20avg: [1x51 double]
                neGR: [1x51 double]
              fne_gr: [1x51 double]
                 nel: [1x51 double]
              pprime: [60x51 double]
             ttprime: [60x51 double]
               rbphi: [16x51 double]
                beta: [1x51 double]
               betaN: [1x51 double]
             betapol: [1x51 double]
             rhoq1_1: [1x51 double]
             rhoq2_1: [1x51 double]
             rhoq3_1: [1x51 double]
             rhoq1_2: [1x51 double]
             rhoq2_2: [1x51 double]
             rhoq3_2: [1x51 double]
                  Li: [1x51 double]
                 li3: [1x51 double]
                 pav: [1x51 double]
                 Ini: [16x51 double]
                 Ioh: [16x51 double]
                 fbs: [1x51 double]
                faux: [1x51 double]
                 fni: [1x51 double]
                 foh: [1x51 double]
                 PLH: [1x51 double]
               Pecin: [1x51 double]
               Pnbin: [1x51 double]
               Picin: [1x51 double]
                Pein: [1x51 double]
               Peout: [1x51 double]
               Penet: [1x51 double]
                Piin: [1x51 double]
               Piout: [1x51 double]
               Pinet: [1x51 double]
             Pfusion: [1x51 double]
              Palpha: [1x51 double]
               Pheat: [1x51 double]
                   Q: [1x51 double]
              Ptotin: [1x51 double]
               dWtdt: [1x51 double]
               Ploss: [1x51 double]
                tauE: [1x51 double]
            tauE_H98: [1x51 double]
                 H98: [1x51 double]
              rhodep: []
                  ST: [1x51 logical]
             STtimes: [1x0 double]
               STper: 0
               df_dx: {1x51 cell}
            df_dxdot: {1x51 cell}
              df_dxk: {1x51 cell}
             df_dxk1: {1x51 cell}
               df_du: {1x51 cell}

## Customized display of results

    % surf plots
    subplot(231); surf(out.time,out.rho,out.psi); xlabel('t [s]'); ylabel('\rho'); title('\psi[Wb]');
    subplot(232); surf(out.time,out.rho,out.q); xlabel('t [s]'); ylabel('\rho'); title('q');
    subplot(233); surf(out.time,out.rho,out.Ip); xlabel('t [s]'); ylabel('\rho'); title('I_p');
    subplot(234); surf(out.time,out.rho,out.upl); xlabel('t [s]'); ylabel('\rho'); title('U_{pl} [V]');
    subplot(235); surf(out.time,out.rho,out.te); xlabel('t [s]'); ylabel('\rho'); title('T_e [eV]');
    subplot(236); surf(out.time,out.rho,out.tep); xlabel('t [s]'); ylabel('\rho'); title('dT_e/d\rho [eV]');

note that some quantities in this output are also given on the grid of gauss points: out.rhogauss, though the normal grid version is also available. Beware, the normal grid point values are just a linear interpolation/extrapolation of the gauss point values.

    clf;
    hp=plot(out.rhogauss,out.jbsgauss(:,end),'x',out.rho,out.jbs(:,end),'o');
    grid on; ylabel('[A/m^2]'); xlabel('\rho'); title('j_{BS}');
    legend(hp,{'gauss grid points','main grid points'},'location','northeast')

## Interactive plotting tool

There exists an interactive plotting tool to plot several RAPTOR results at the same time, allowing one to adjust the time of the plotted profiles with the mouse

    outs{1} = out; % store old out structure

    % Do a new RAPTOR run with different input parameters
    params.debug.iterplot = 0; % do not plot this time
    U = bsxfun(@times,U,linspace(1,2,numel(params.tgrid))); % scale U
    simres = RAPTOR_predictive(x0,g,v,U,model,params);
    outs{2} = RAPTOR_out(simres,model,params); % different structure

         it  telaps newt     res   t[ms]  dt[ms]  Ip[kA] Icd[kA] Ibs[kA] Ioh[kA]      qe    qmin      q0   Vl[V] Te0[keV] Ti0[keV] ne0[e19]   f_ss
          1   0.048    4 5.7e-10       0       1      80       0    1.41    78.5    15.3    5.55    5.55 3.4e+00    0.21    0.21    1.00 1.1e+01
          6    0.18    3 3.9e-15       5       1     137       0    3.29     134    8.94    4.05    4.05 6.8e+00    0.32    0.32    1.00 2.6e+01
         11     0.3    3 8.7e-15      10       1     203       0    5.55     198    6.03    3.44    3.44 6.9e+00    0.48    0.48    1.00 3.0e+01
         16    0.42    2 9.7e-09      15       1     278       0    7.87     270    4.41    2.69    3.05 6.8e+00    0.70    0.70    1.00 3.2e+01
         21    0.51    3 2.1e-11      20       1     350       0    10.2     340    3.51    2.14    2.79 5.5e+00    0.96    0.96    1.00 2.6e+01
         26    0.63    2 1.0e-09      25       1     375       0      12     363    3.27    1.82    2.61 4.4e+00    1.23    1.23    1.00 2.0e+01
         31    0.72    2 5.2e-11      30       1     400       0    13.1     387    3.07    1.61    2.48 4.1e+00    1.47    1.47    1.00 1.8e+01
         36    0.82    2 6.9e-12      35       1     425       0    13.8     411    2.89    1.48    2.37 3.8e+00    1.66    1.66    1.00 1.7e+01
         41    0.91    2 1.3e-12      40       1     450       0    14.2     436    2.73    1.37    2.27 3.7e+00    1.83    1.83    1.00 1.6e+01
         46       1    2 4.0e-13      45       1     475       0    14.4     460    2.58    1.27    2.18 3.5e+00    1.98    1.98    1.00 1.5e+01
         51     1.1    2 1.6e-13      50       1     500       0    14.5     485    2.45     1.2    2.09 3.4e+00    2.12    2.12    1.00 1.5e+01

Plot 'overview'

    hf = RAPTOR_plot_GUI(outs,'overview');
    % drag your mouse over the time traces to see the time evolution
    figure(hf); % make visible to export

    Warning: could not plot
     Xdata:out.rho
     Ydata:[1]*ones(size(out.rho))

    Warning: could not plot
     Xdata:out.rho
     Ydata:[2]*ones(size(out.rho))

    Warning: could not plot
     Xdata:out.rho
     Ydata:[3]*ones(size(out.rho))

    Warning: could not plot
     Xdata:out.rho
     Ydata:[1]*ones(size(out.rho))

    Warning: could not plot
     Xdata:out.rho
     Ydata:[2]*ones(size(out.rho))

    Warning: could not plot
     Xdata:out.rho
     Ydata:[3]*ones(size(out.rho))

end of tutorial