# RAPTOR tutorial: ad-hoc transport model tuning

In this tutorial we will vary the parameters of the electron transport model to see their effect on the plasma.

We use the heating and current drive setup from tutorial 3 again, but modify some parmeters: in particular we choose a smaller time grid to handle the rapid rise of the pedestal.

## Contents

- Changing parameters for the default ad-hoc transport model
- Simple H-mode model
- Bohm Gyrobohm
- QualiKyz - Artificial Neural Network
- Manual chie
- plot some comparisons
- Time-varying Manual Chie
- Analytical formula for chi_e with caracteristic scalelength

    close all hidden;
    clear; run('../RAPTOR_path.m');

    [config] = RAPTOR_config; % load default config
    config.grid.tgrid = [0:0.001:0.2];
    config.grid.rhogrid = linspace(0,1,21); % spatial grid
    config.debug.iterdisp = 20; % display progress every N iterations

    % set ECH/ECCD parameters
    config.echcd.params.active = true;
    config.echcd.params.rdep = [0.4]; %
    config.echcd.params.wdep = [0.35]; % deposition width
    config.echcd.params.cd_eff = [1]; % CD efficiency
    config.echcd.params.uindices = [2]; % index in input vector

    % set NBI parameters
    config.nbhcd.params.active = true;
    config.nbhcd.params.rdep = [0]; %
    config.nbhcd.params.wdep = [0.8]; % broad heating
    config.nbhcd.params.cd_eff = [0]; % no current drive
    config.nbhcd.params.uindices = [3]; % index in input vector

    % rerun the config file with the new parameter set
    [model,params,init,g,v,U] = build_RAPTOR_model(config);

    % plasma current
    U(1,:)= 200e3*ones(size(params.tgrid)); % input Ip trace: ramp from 80kA to 200kA
    U(1,params.tgrid<0.02) = linspace(80e3,200e3,sum(params.tgrid<0.02));

    % first EC actuator power:  2MW start at 40ms
    U(2,:) = zeros(size(params.tgrid)); U(2,params.tgrid>0.04) = 2e6;
    U(3,:) = zeros(size(params.tgrid)); U(3,params.tgrid>0.04) = 1e6;

    % initial condition
    init.Ip0 = U(1,1);
    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);

         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.08    4 5.1e-10       0       1      80       0    1.41    78.6    15.4    5.59    5.59 3.4e+00    0.21    0.21    1.00 1.2e+01
         21    0.77    3 9.0e-13      20       1     200       0    6.24     194    6.16    2.32    2.32 4.0e+00    0.64    0.64    1.00 1.6e+01
         41     1.5    2 9.2e-13      40       1     200       0    6.67     193    6.16    1.38    1.38 2.4e+00    0.95    0.95    1.00 8.4e+00
         61     2.2    2 1.8e-13      60       1     200     252     113    -165    6.16    1.46    1.46 5.7e-02   15.90   15.90    1.00 1.7e+00
         81     2.9    2 5.6e-12      80       1     200     291     134    -225    6.16    1.56    1.56 2.9e-02   18.74   18.74    1.00 1.6e+00
        101     3.5    2 8.0e-13     100       1     200     376     178    -353    6.16    1.52    1.67 1.0e-02   24.72   24.72    1.00 1.7e+00
        121     4.2    2 9.0e-14     120       1     200     446     217    -463    6.16    1.44     1.8 -1.2e-02   28.99   28.99    1.00 1.5e+00
        141     4.8    2 4.6e-14     140       1     200     467     234    -501    6.16    1.38    1.99 -2.6e-02   30.10   30.10    1.00 1.3e+00
        161     5.5    2 6.5e-14     160       1     200     474     246    -520    6.16    1.35    2.22 -3.6e-02   30.53   30.53    1.00 1.1e+00
        181     6.1    2 5.5e-14     180       1     200     479     260    -540    6.16    1.32    2.47 -4.5e-02   31.05   31.05    1.00 1.1e+00
        201     6.8    2 2.8e-14     200       1     200     483     274    -556    6.16    1.31    2.57 -5.3e-02   31.19   31.19    1.00 1.0e+00

## Changing parameters for the default ad-hoc transport model

RAPTOR includes an ad-hoc expression for \chi_e

where

and is an extra factor to reduce the transport at the edge to create an H-mode pedestal.

This model is written as

The relevant parameters can be set from the params structure

    params.chi_e

    ans =

                 cneo: 0.5000
                 cano: 7
                 witb: 3
                 aitb: 1
                 ditb: 0
        chi_csawtooth: 0
        chi_wsawtooth: 10
          te_constant: 0
              aTeterm: 0
             arhoterm: 1
             apsiterm: 1
          chiecentral: 10
         deltacentral: 0.1500
             chiiscal: 1
             implicit: 1
                check: 0

For example, to remove the confinement enhancement at low shear we can set

    params.chi_e.aitb = 0;

    % simulate
    [simres_noenhance] = RAPTOR_predictive(x0,g,v,U,model,params);
    out_noenhance = RAPTOR_out(simres_noenhance,model,params);
    RAPTOR_plot_GUI({out,out_noenhance},'overview')

         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.064    4 4.4e-10       0       1      80       0    1.41    78.6    15.4    5.57    5.57 3.3e+00    0.19    0.19    1.00 1.1e+01
         21    0.64    3 8.1e-13      20       1     200       0    5.66     194    6.16    2.42    2.42 4.1e+00    0.56    0.56    1.00 1.6e+01
         41     1.1    2 8.3e-13      40       1     200       0    6.27     194    6.16    1.42    1.42 2.4e+00    0.83    0.83    1.00 7.8e+00
         61     1.7    2 3.3e-14      60       1     200     215    93.8    -108    6.16    1.52    1.52 7.0e-02   11.77   11.77    1.00 1.7e+00
         81     2.3    2 4.4e-14      80       1     200     220    95.6    -115    6.16    1.62    1.62 3.7e-02   11.78   11.78    1.00 1.4e+00
        101     2.8    2 3.7e-14     100       1     200     223    97.2    -120    6.16    1.65    1.73 1.5e-02   11.74   11.74    1.00 1.2e+00
        121     3.3    2 4.9e-14     120       1     200     226    98.9    -125    6.16     1.6    1.87 -3.1e-04   11.69   11.69    1.00 1.0e+00
        141     3.8    2 4.4e-14     140       1     200     228     101    -129    6.16    1.57    2.06 -1.1e-02   11.62   11.62    1.00 8.9e-01
        161     4.3    1 9.1e-09     160       1     200     230     103    -132    6.16    1.54    2.35 -1.9e-02   11.55   11.55    1.00 8.0e-01
        181     4.7    1 5.7e-09     180       1     200     231     105    -136    6.16    1.52    2.81 -2.6e-02   11.47   11.47    1.00 7.2e-01
        201       5    1 3.8e-09     200       1     200     232     108    -140    6.16     1.5    3.64 -3.1e-02   11.38   11.38    1.00 6.5e-01

    ans =

      Figure (1: RAPTOR output) with properties:

          Number: 1
            Name: 'RAPTOR output'
           Color: [0.9400 0.9400 0.9400]
        Position: [0.0490 0.1983 0.9010 0.7058]
           Units: 'normalized'

      Use GET to show all properties

## Simple H-mode model

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

    Hmode_indicator = rampfun(params.tgrid,0.08,0,0.1,1) - rampfun(params.tgrid,0.18,0,0.19,1);
    plot(params.tgrid,Hmode_indicator);
    nt = numel(params.tgrid);
    x0 = RAPTOR_initial_conditions(model,init,g(:,1),v(:,1));
    [model,params,init,g,v,~] = build_RAPTOR_model(config);
    v = v*ones(size(params.tgrid));
    v(model.hmode.vind,:) = Hmode_indicator;
    params.hmode.active = true;

    [simres_Hmode] = RAPTOR_predictive(x0,g,v,U,model,params);
    out_Hmode = RAPTOR_out(simres_Hmode,model,params);
    RAPTOR_plot_GUI({out_noenhance,out_Hmode},'overview')

         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.052    4 5.1e-10       0       1      80       0    1.41    78.6    15.4    5.59    5.59 3.4e+00    0.21    0.21    1.00 1.2e+01
         21     0.6    3 9.0e-13      20       1     200       0    6.24     194    6.16    2.32    2.32 4.0e+00    0.64    0.64    1.00 1.6e+01
         41     1.1    2 9.2e-13      40       1     200       0    6.67     193    6.16    1.38    1.38 2.4e+00    0.95    0.95    1.00 8.4e+00
         61     1.7    2 1.8e-13      60       1     200     252     113    -165    6.16    1.46    1.46 5.7e-02   15.90   15.90    1.00 1.7e+00
         81     2.3    2 5.6e-12      80       1     200     291     134    -225    6.16    1.56    1.56 2.9e-02   18.74   18.74    1.00 1.6e+00
        101     2.9    2 5.1e-12     100       1     200     382     181    -363    6.16    1.53    1.67 3.8e-03   24.75   24.75    1.00 1.7e+00
        121     3.5    2 7.9e-14     120       1     200     453     220    -473    6.16    1.45     1.8 -1.6e-02   29.10   29.10    1.00 1.5e+00
        141       4    2 6.1e-14     140       1     200     475     238    -513    6.16    1.39    1.98 -2.7e-02   30.28   30.28    1.00 1.3e+00
        161     4.5    2 7.2e-14     160       1     200     482     250    -533    6.16    1.36    2.21 -3.6e-02   30.71   30.71    1.00 1.1e+00
        181     5.1    2 5.8e-14     180       1     200     487     265    -553    6.16    1.33    2.46 -4.3e-02   31.23   31.23    1.00 1.1e+00
        201     5.7    2 5.6e-14     200       1     200     485     276    -561    6.16    1.32     2.6 -4.0e-02   31.40   31.40    1.00 1.0e+00

    ans =

      Figure (1: RAPTOR output) with properties:

          Number: 1
            Name: 'RAPTOR output'
           Color: [0.9400 0.9400 0.9400]
        Position: [0.0490 0.1983 0.9010 0.7058]
           Units: 'normalized'

      Use GET to show all properties

## Bohm Gyrobohm

We can also use the Bohm-gyroBohm model first initialize new default parameters

    config.chi_e = chi_e('BgB');
    %
    % we now have a new set of parameters to change, but leave them fixed for now
    config.chi_e.params

    % make new model and parameters
    [model,params,init,g,v,U] = build_RAPTOR_model(config);

    % run simulation with new model
    params.debug.iterdisp = 10; % plot progress to see pedestal growing
    [simres_Hmode] = RAPTOR_predictive(x0,g,v,U,model,params);
    out_BgB = RAPTOR_out(simres_Hmode,model,params);
    RAPTOR_plot_GUI({out_noenhance,out_BgB},'overview')

    ans =

                  aeb: 8.0000e-05
                 aegb: 0.0700
                  aib: 1.6000e-04
                 aigb: 0.0175
                 cneo: 1.0000e-03
        chi_csawtooth: 0
        chi_wsawtooth: 5
               chimax: 20
                   a1: 2
                   z1: 1
             implicit: 1
                check: 0

         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.072    5 1.5e-09       0       1      80       0    1.41    78.6    15.4    5.67    5.67 2.3e+00    0.25    0.25    1.00 6.3e+00
         11    0.39    2 7.0e-12      10       1      80       0    6.76    73.2    15.4    3.99    3.99 1.4e+00    0.48    0.48    1.00 3.6e+00
         21    0.66    2 3.8e-14      20       1      80       0    7.27    72.7    15.4    3.02    3.02 1.3e+00    0.57    0.57    1.00 2.5e+00
         31    0.94    2 6.2e-15      30       1      80       0    7.51    72.5    15.4    2.48    2.48 1.2e+00    0.63    0.63    1.00 1.8e+00
         41     1.2    2 5.6e-15      40       1      80       0    7.66    72.3    15.4    2.17    2.17 1.1e+00    0.68    0.68    1.00 1.3e+00
         51     1.5    2 6.9e-15      50       1      80       0    7.76    72.2    15.4    1.98    1.98 1.0e+00    0.71    0.71    1.00 9.4e-01
         61     1.7    2 6.3e-15      60       1      80       0    7.82    72.2    15.4    1.85    1.85 1.0e+00    0.74    0.74    1.00 6.8e-01
         71     1.9    1 7.7e-09      70       1      80       0    7.87    72.1    15.4    1.77    1.77 1.0e+00    0.75    0.75    1.00 5.0e-01
         81     2.1    1 4.0e-09      80       1      80       0     7.9    72.1    15.4    1.71    1.71 9.8e-01    0.77    0.77    1.00 3.6e-01
         91     2.3    1 2.1e-09      90       1      80       0    7.92    72.1    15.4    1.67    1.67 9.7e-01    0.77    0.77    1.00 2.6e-01
        101     2.5    1 1.1e-09     100       1      80       0    7.94      72    15.4    1.64    1.64 9.6e-01    0.78    0.78    1.00 1.9e-01
        111     2.6    1 5.9e-10     110       1      80       0    7.95      72    15.4    1.62    1.62 9.6e-01    0.79    0.79    1.00 1.4e-01
        121     2.8    1 3.1e-10     120       1      80       0    7.95      72    15.4    1.61    1.61 9.5e-01    0.79    0.79    1.00 1.0e-01
        131       3    1 1.7e-10     130       1      80       0    7.96      72    15.4     1.6     1.6 9.5e-01    0.79    0.79    1.00 7.5e-02
        141     3.2    1 8.9e-11     140       1      80       0    7.96      72    15.4    1.59    1.59 9.5e-01    0.79    0.79    1.00 5.5e-02
        151     3.4    1 4.7e-11     150       1      80       0    7.97      72    15.4    1.59    1.59 9.4e-01    0.80    0.80    1.00 4.0e-02
        161     3.5    1 2.5e-11     160       1      80       0    7.97      72    15.4    1.58    1.58 9.4e-01    0.80    0.80    1.00 2.9e-02
        171     3.7    1 1.3e-11     170       1      80       0    7.97      72    15.4    1.58    1.58 9.4e-01    0.80    0.80    1.00 2.1e-02
        181     3.9    1 7.2e-12     180       1      80       0    7.97      72    15.4    1.58    1.58 9.4e-01    0.80    0.80    1.00 1.6e-02
        191     4.1    1 3.8e-12     190       1      80       0    7.97      72    15.4    1.58    1.58 9.4e-01    0.80    0.80    1.00 1.1e-02
         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
        201     4.2    1 2.0e-12     200       1      80       0    7.97      72    15.4    1.57    1.57 9.4e-01    0.80    0.80    1.00 8.3e-03
    Error using matlab.graphics.axis.Axes/get
    Invalid or deleted object.

    Error in RAPTOR_plot_GUI>gui_init/wbmf (line 291)
                haxposs = get(hax,'position');

    Error using drawnow
    Error while evaluating Figure WindowButtonMotionFcn

    ans =

      Figure (1: RAPTOR output) with properties:

          Number: 1
            Name: 'RAPTOR output'
           Color: [0.9400 0.9400 0.9400]
        Position: [0.0490 0.1983 0.9010 0.7058]
           Units: 'normalized'

      Use GET to show all properties

## QualiKyz - Artificial Neural Network

Use ANN model first initialize new default parameters

    config.chi_e = chi_e('qlkANN');
    %
    % we now have a new set of parameters to change, but leave them fixed for now
    config.chi_e.params

    % make new model and parameters
    [model,params,init,g,v,U] = build_RAPTOR_model(config);
    % run simulation with new model
    params.debug.iterdisp = 10; % plot progress to see pedestal growing
    [simres_qlkANN] = RAPTOR_predictive(x0,g,v,U,model,params);
    out_qlkANN = RAPTOR_out(simres_qlkANN,model,params);
    RAPTOR_plot_GUI({out_noenhance,out_qlkANN},'overview')

    ans =

           constrains: 1
           constrainq: 1
             shearmin: 0.1000
             shearmax: 3
                 qmin: 1
                 qmax: 5
               IW_eef: [20x4 double]
              L1W_eef: [20x20 double]
              L2W_eef: [1x20 double]
               b1_eef: [20x1 double]
               b2_eef: [20x1 double]
               b3_eef: -1.8186
          max_qlk_eef: 4.9989
          min_qlk_eef: 5.9621e-04
          prepros_eef: -1
               IW_ief: [20x4 double]
              L1W_ief: [20x20 double]
              L2W_ief: [1x20 double]
               b1_ief: [20x1 double]
               b2_ief: [20x1 double]
               b3_ief: 1.4270
          max_qlk_ief: 19.9961
          min_qlk_ief: 0.0021
          prepros_ief: -1
               IW_dfe: [20x4 double]
              L1W_dfe: [20x20 double]
              L2W_dfe: [1x20 double]
               b1_dfe: [20x1 double]
               b2_dfe: [20x1 double]
               b3_dfe: -0.1042
          max_qlk_dfe: 49.8982
          min_qlk_dfe: 0.0026
          prepros_dfe: -1
               IW_vce: [20x4 double]
              L1W_vce: [20x20 double]
              L2W_vce: [1x20 double]
               b1_vce: [20x1 double]
               b2_vce: [20x1 double]
               b3_vce: 0.3464
          max_qlk_vce: 29.9993
          min_qlk_vce: -8.7546
          prepros_vce: -1
               IW_vte: [20x4 double]
              L1W_vte: [20x20 double]
              L2W_vte: [1x20 double]
               b1_vte: [20x1 double]
               b2_vte: [20x1 double]
               b3_vte: 0.4158
          max_qlk_vte: 0.1151
          min_qlk_vte: -29.9764
          prepros_vte: -1
        constrainrlti: 1
        constraintite: 1
          shift_chi_e: 0.7500
              chiemin: 0.4000
              chiimin: 0.4000
                 Dmin: 0
        allowpospinch: 1
        thresholdscan: 0
             implicit: 1
                check: 0
              rhocore: 0.3000
               rhoped: 0.8000
             chieedge: 5
             chiecore: 2
             chiiedge: 5
             chiicore: 2
              rltimin: 2
              rltimax: 12
              titemin: 0.3000
              titemax: 3
        out_sm_kernel: 1

         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.09    5 6.0e-09       0       1      80       0    1.41    78.6    15.4    5.79    5.79 2.0e+00    0.24    0.24    1.00 4.5e+00
    Warning: Exceeded maximum newton iterations
    Warning: exceeded maximum newton iterations
          3    0.65   50 5.9e-03       2       1      80       0    5.11    74.9    15.4    5.43    5.43 1.3e+00    0.32    0.32    1.00 1.4e+00
    Warning: NANs in X, RAPTOR simulation did not complete successfully
    Warning: could not plot
     Xdata:out.time
     Ydata:out.tauE

    ans =

      Figure (1: RAPTOR output) with properties:

          Number: 1
            Name: 'RAPTOR output'
           Color: [0.9400 0.9400 0.9400]
        Position: [0.0490 0.1983 0.9010 0.7058]
           Units: 'normalized'

      Use GET to show all properties

## Manual chie

We can also choose a fixed chi_e profile manually, by setting `'modeltype'` to `'manual'` and assigning the profile to `'cneo'`.

    config.chi_e = chi_e('manual');

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

    params.chi_e.chi_e = 10; % constant value
    params.debug.iterplot = 0; % no plots

    [simres_manual1] = RAPTOR_predictive(x0,g,v,U,model,params);
    out_manual1 = RAPTOR_out(simres_manual1,model,params);

    params.chi_e.chi_e = (5+10*model.rgrid.rhogauss.^2); % spatially varying profile
    params.debug.iterplot = 0; % no plots

    [simres_manual2] = RAPTOR_predictive(x0,g,v,U,model,params);
    out_manual2 = RAPTOR_out(simres_manual2,model,params);

         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.053    5 1.0e-14       0       1      80       0    1.41    78.6    15.4    5.78    5.78 2.1e+00    0.22    0.22    1.00 5.3e+00
         21    0.55    2 5.8e-15      20       1      80       0    5.22    74.8    15.4    4.13    4.13 1.4e+00    0.33    0.33    1.00 7.0e-01
         41    0.97    1 2.9e-09      40       1      80       0    5.24    74.7    15.4    3.84    3.84 1.4e+00    0.33    0.33    1.00 1.6e-01
         61     1.3    1 1.5e-10      60       1      80       0    5.24    74.7    15.4    3.78    3.78 1.4e+00    0.33    0.33    1.00 3.7e-02
         81     1.6    1 8.0e-12      80       1      80       0    5.24    74.7    15.4    3.77    3.77 1.4e+00    0.33    0.33    1.00 8.5e-03
        101     1.9    1 4.3e-13     100       1      80       0    5.24    74.7    15.4    3.76    3.76 1.4e+00    0.33    0.33    1.00 2.0e-03
        121     2.2    1 3.3e-14     120       1      80       0    5.24    74.7    15.4    3.76    3.76 1.4e+00    0.34    0.34    1.00 4.5e-04
        141     2.5    1 4.8e-14     140       1      80       0    5.24    74.7    15.4    3.76    3.76 1.4e+00    0.34    0.34    1.00 1.0e-04
        161     2.8    1 3.8e-14     160       1      80       0    5.24    74.7    15.4    3.76    3.76 1.4e+00    0.34    0.34    1.00 2.4e-05
        181     3.1    0 6.1e-09     180       1      80       0    5.24    74.7    15.4    3.76    3.76 1.4e+00    0.34    0.34    1.00 5.9e-06
        201     3.3    0 2.6e-09     200       1      80       0    5.24    74.7    15.4    3.76    3.76 1.4e+00    0.34    0.34    1.00 1.2e-06
         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.053    5 5.9e-15       0       1      80       0    1.41    78.6    15.4    5.74    5.74 2.3e+00    0.23    0.23    1.00 6.4e+00
         21    0.51    2 3.0e-14      20       1      80       0    5.05    74.9    15.4    3.53    3.53 1.5e+00    0.37    0.37    1.00 1.3e+00
         41    0.93    2 1.4e-14      40       1      80       0     5.1    74.9    15.4    3.08    3.08 1.5e+00    0.38    0.38    1.00 3.6e-01
         61     1.2    1 9.0e-10      60       1      80       0    5.11    74.9    15.4    2.97    2.97 1.4e+00    0.39    0.39    1.00 1.0e-01
         81     1.5    1 7.3e-11      80       1      80       0    5.11    74.9    15.4    2.94    2.94 1.4e+00    0.39    0.39    1.00 2.9e-02
        101     1.8    1 5.9e-12     100       1      80       0    5.11    74.9    15.4    2.93    2.93 1.4e+00    0.39    0.39    1.00 8.2e-03
        121     2.1    1 4.8e-13     120       1      80       0    5.11    74.9    15.4    2.93    2.93 1.4e+00    0.39    0.39    1.00 2.3e-03
        141     2.4    1 5.5e-14     140       1      80       0    5.11    74.9    15.4    2.93    2.93 1.4e+00    0.39    0.39    1.00 6.7e-04
        161     2.7    1 4.8e-14     160       1      80       0    5.11    74.9    15.4    2.93    2.93 1.4e+00    0.39    0.39    1.00 1.9e-04
        181       3    1 4.4e-14     180       1      80       0    5.11    74.9    15.4    2.93    2.93 1.4e+00    0.39    0.39    1.00 5.4e-05
        201     3.3    1 3.4e-14     200       1      80       0    5.11    74.9    15.4    2.93    2.93 1.4e+00    0.39    0.39    1.00 1.5e-05

## plot some comparisons

    figure; clf;
    subplot(211)
    plot(out_manual1.rho,out_manual1.te(:,end)/1e3,'r',...
         out_manual2.rho,out_manual2.te(:,end)/1e3,'k'...
         ); ylabel('T_e [keV]'), xlabel('\rho')

    subplot(212)
    plot(out_manual1.rho,out_manual1.chie(:,end),'r',...
         out_manual2.rho,out_manual2.chie(:,end),'k'...
         ); ylabel('\chi_e'), xlabel('\rho')

## Time-varying Manual Chie

It is also possible to assign to `'cneo'` a time-varying cie_e profile. This profile should be provided as a matrix of dimension:$\rho_{gauss}$ by t. Here is an example of a profile which is based on a parabolic function with time dependent gradient. A default value for tgrid is used of dimension 2 by 1 so that the previous examples hold because these have dimension # of rhogauss by 1.

    % Create a matrix of the right size with a parabol which gradient increases
    % in time:

    [X,Y] = meshgrid(25*params.tgrid,model.rgrid.rhogauss);
    params.chi_e.chi_e = 5 + 2.^(X).*Y.^2;

    [simres_manual_time] = RAPTOR_predictive(x0,g,v,U,model,params);
    out_manual_time = RAPTOR_out(simres_manual_time,model,params);

    % Plot the result
    subplot(211)
    plot(out_manual_time.rho,out_manual_time.te(:,end)/1e3,'r');
        ylabel('T_e [keV]'), xlabel('\rho')

    subplot(212) % Chi at the beginning and the end
    plot(out_manual_time.rho,out_manual_time.chie(:,end),'r',...
         out_manual_time.rho,out_manual_time.chie(:,1),'b'...
        ); ylabel('\chi_e'), xlabel('\rho')

    return

         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.069    5 3.2e-14       0       1      80       0    1.41    78.6    15.4    5.78    5.78 2.0e+00    0.23    0.23    1.00 5.2e+00
         21    0.56    2 6.3e-15      20       1      80       0    7.26    72.7    15.4     4.1     4.1 1.1e+00    0.42    0.42    1.00 1.0e+00
         41       1    2 6.4e-15      40       1      80       0    7.06    72.9    15.4     3.6     3.6 1.1e+00    0.42    0.42    1.00 4.2e-01
         61     1.5    2 1.6e-14      60       1      80       0    6.78    73.2    15.4    3.41    3.41 1.1e+00    0.42    0.42    1.00 2.3e-01
         81     1.9    2 1.1e-14      80       1      80       0    6.42    73.6    15.4     3.3     3.3 1.2e+00    0.42    0.42    1.00 1.9e-01
        101     2.4    2 3.6e-14     100       1      80       0    5.99      74    15.4    3.21    3.21 1.3e+00    0.41    0.41    1.00 2.0e-01
        121     2.8    2 4.0e-14     120       1      80       0    5.51    74.5    15.4     3.1     3.1 1.4e+00    0.40    0.40    1.00 2.3e-01
        141     3.3    2 1.5e-14     140       1      80       0    4.98      75    15.4    2.98    2.98 1.5e+00    0.38    0.38    1.00 2.5e-01
        161     3.8    2 2.4e-14     160       1      80       0    4.43    75.6    15.4    2.85    2.85 1.6e+00    0.37    0.37    1.00 2.6e-01
        181     4.2    2 2.9e-14     180       1      80       0    3.88    76.1    15.4    2.72    2.72 1.8e+00    0.35    0.35    1.00 2.5e-01
        201     4.7    2 4.3e-14     200       1      80       0    3.35    76.6    15.4    2.59    2.59 2.0e+00    0.34    0.34    1.00 2.3e-01

## Analytical formula for chi_e with caracteristic scalelength

L-mode simulation

    config.chi_e = chi_e('MS');
    % A new set of parameters for chi_e module
    config.chi_e.params

    % Time grid
    config.grid.tgrid = [0:0.001:0.15];
    % Define reference He factor(==tauEe/tauE_H98) on RAPTOR time grid
    config.chi_e.params.hefactor = 0.2*ones(size(config.grid.tgrid));
    % Make new model and parameters
    [model,params,init,g,v,U] = build_RAPTOR_model(config);

    % Run simulation with new model
    simres_MS = RAPTOR_predictive(x0,g,v,U,model,params);
    out_MS = RAPTOR_out(simres_MS,model,params);

    % Profile of electron temperature
    plot(out_MS.rho,out_MS.te(:,end),'m','LineWidth',2); hold on
    % Define inded for rho_inv on gauss grid
    [~,ind_rinv] = min((out_MS.rhogauss - params.chi_e.rhoinv).^2);
    % rho_ped index
    [~,ind_rped] = min((out_MS.rhogauss - params.chi_e.rhoped(1)).^2);
    % Te on gauss grid
    tegauss = eval_te(simres_MS.X(:,end),simres_MS.G(:,1),simres_MS.V(:,1),model,'gauss');
    % dlnTe/drho
    dlnTe_drhogauss = bsxfun(@rdivide,diff(log(tegauss(:,end))),diff(out_MS.rhogauss(:,end)));
    dlnTe_drhogauss = [dlnTe_drhogauss' dlnTe_drhogauss(end)]';
    % Get lambdaTe value from simulation
    slambdaTe = -dlnTe_drhogauss(ind_rinv:ind_rped,end);
    % Get muTe value from simulation
    smuTe = -tegauss(ind_rped:end,end).*dlnTe_drhogauss(ind_rped:end,end);
    muTe = mean(smuTe);

    % Add core part
    coeffl = tegauss(ind_rped)/exp(-params.chi_e.lambdaTe(1)*out_MS.rhogauss(ind_rped));
    te_core = coeffl*exp(-params.chi_e.lambdaTe(1)*out_MS.rhogauss(ind_rinv:ind_rped));
    plot(out_MS.rhogauss(ind_rinv:ind_rped),te_core,'b--','LineWidth',2);

    % Add pedestal part
    coeffm = tegauss(end) + muTe*out_MS.rhogauss(end);
    te_edge = coeffm - muTe*out_MS.rhogauss(ind_rped:end);
    plot(out_MS.rhogauss(ind_rped:end),te_edge,'r--','LineWidth',2);

    xlabel('rho');legend('Te\_prof','labmdaTe','muTe','Location','SouthWest');
    title('Te profile');