function varargout = RAPTOR_config(varargin)
% RAPTOR main configuration script
% Called to get RAPTOR configuration structure,
% or to create a RAPTOR model and parameter structure from a given
% configuration
%
% Typical usage:
% [config] = RAPTOR_config; % call with no inputs to get default config structure
% % then make modifications to parameters

% F. Felici CRPP 2010,2011, TU/e 2012,2017

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Init Main Variables
%%%%%%%%%%%%%%%%%%%%%%%
%default param struct

% start with some secondary structures

%%% debug options
debug = struct(...
    'verbosity',0,... % verbosity level
    'iterdisp',10,... % iteration display every N time steps
    'iterplot',0,... % plot display every N time steps
    'dispFailedNewtonStep',false,... % warn upon in case of failed Newton step
    'checkstateequation',false,... % Optional check of state equation jacobian
    'checkprofilejacobians',false,...  % Optional check of all gradients of sub-functions
    'plotprofilejacobians',true,... % Optional plot of gradients of sub-functions
    'stapProfilesToCheck','all',... % Optional list of states profiles to check (e.g. 'ne','q') instead of 'all'
    'trapProfilesToCheck','all',... % Optional list of profiles to check (e.g. 'jbsB','pec') instead of 'all'
    'checkfwsens',false,... % optional check of forward sensitivities vs. finite difference
    'pauseEachNewtonStep',false,... % optional plot of profiles at each Newton step individually
    'plotEachNewtonStep',false,... % optional plot of profiles at each Newton step individually
    'plotFailedNewtonStep',false,... % plot of prifiles at failed Newton steps
    'plotJacobians',false, ... % optional plot of Jacobian SV's and Netwon convergence daata
    'checkgeom',false, ... % optional check of geom mapping
    'keyboard_at',[],... % optional keyboard given time point
    'eq_tran_iter', 1, ... % store equilibrium and transp iteration convergence test
    'eq_tran_wtcmdl', 1 ... % every how many equilibrium and transport iteration write messaged to cmd. = 0 for no outputs
    );

%%% numerics parameters
numerics = struct(...
    'usemex',false,...  % option to call mex instead of RAPTOR_predictive_step
    'calc_sens',false, ...   % do not compute sensitivities by default
    'nmax' , uint32(50) ,... % Maximum Newton iterations
    'sporder',3,... %spline order
    'ngauss',4,... % gauss points
    'restol', 1e-8 ,... % residual tolerance
    'tau' , 1, ...   % Newton step size
    'Upl0', 1, ...  % initial guess of (flat) Upl profile. Important for first step convergence.
    'extrap',1, ... % extrapolation factor for first guess of next xk1
    'nextrap',5, ... % number of iterations to ramp use of extrap
    'theta',1,...  % implicit-explicit in time, set 1 for fully implicit
    'eq_tran_iter', struct() ... % 1 will be just compute/take from external file the equilibrium and run RAPTOR
    );   % implicit-explicit in time, set 1 for fully implicit
numerics.eq_tran_iter.N_max_iteration = 2; % =2 --> 2 equlibrium and 2 transport at each time iteration.
numerics.eq_tran_iter.convergence_threshold = 1e-4; % 1e-4 is faster, 1e-5 more precision
numerics.eq_tran_iter.relax_alpha = 1.; % Relaxation parameter between iterations. 1 for no relaxation

% code building options
buildopts = struct('mexify',false); % option to mexify RAPTOR_predictive_step

%%% geometry (equilibrium)
equi = struct(...
    'tokamak','TCV',...
    'source','chease',... % select source, can be 'chease',  'liuqe'. In future 'ESCOT', 'CRONOS', 'IDS'
    'filenames' , {{'TCV_41083_cheasedata'}}, ...  % Reference equilibrium .mat file (chease or astra)
    'times', [],... % default time corresponding to the equilibrium data on RAPTOR time grid
    'offset',0.,...% defaul offset time: starting point for RAPTOR simulation on an experimental time grid  [s]
    'eq_call', false,... % Call eq_solver during time iteration
    'epsilon',0.28,... % default inverse aspect ratio
    'kappa',1.5,... % typical kappa
    'delta',0.3,... % typical delta
    'B0'  , 1.44,... % Toroidal field [T]
    'R0'  , 0.88,... % Major radius [m]
    'Lext', 1e-6,... % external inductance (guess) [Henry]
    'gsplineorder', 3, ...% polynomial spline order for geometry interpolation
    'gspline_axis_tens', 0.5e-3, ... % Spline tension for geometrical quantities interpolation at the axis
    'gspline_edge_tens', 0.5e-3, ... % spline tension for geometrical quantities interpolation at the edge
    'gknots', []... 
    );


%%% Atomic numbers
atom = struct(...
    'Zi'  , 1 ,...% Deuterium plasma, main ion charge
    'Ai'  , 2 ,... % Deuterium atomic weight
    'Z1'  , 6 ,...% first impurity
    'A1'  , 12, ...% first impurity A
    'Z2'  , 7 ,...% second impurity A
    'A2'  , 14, ...% second impurity A
    'Z3'  , 4 ,...% second impurity A
    'A3'  , 8 ...% second impurity A
    );

%%% Grid
grid = struct(...
    'rhogrid'       ,linspace(0,1,16),... % radial grid
    'tgrid'         ,[0:0.001:0.05]... % time grid
    );

% disturbance model default
dist = struct('ndist',0,'noisemodel','default','disttype','equidist','distgausswidth',0,'distgaussrpos',0);
% ndist: number of disturbance states for this equation (only for state observer)
% disttype: type of disturbances (equidistant or gaussian)
% distgausswidth: width of gaussian when disttype is set to gaussian
% distgaussrpos: center of gaussian when disttype is set to gaussian

% boundary condition

%%% Settings for the various equations
psi =  struct('method','state',  'BC', struct('type','Ip','rhoSource','Constant','rho',1,'valueSource','u'),   'dist',dist); % options for psi equation
te  =  struct('method','state',  'BC', define_BC('Dirichlet',50), 'dist',dist); % options for Te equation
ti  =  struct('method','tescal', 'BC', define_BC('Dirichlet',50), 'dist',dist);
ne  =  struct('method','direct', 'BC', define_BC('Dirichlet',1e18), 'dist',dist);
ni  =  struct('method','direct', 'BC', define_BC('Dirichlet',1e18), 'dist',dist);
ze  =  struct('method','direct', 'BC', define_BC('Dirichlet',1), 'dist',dist);
n1  =  struct('method','nescal', 'BC', define_BC('Dirichlet',0), 'dist',dist); 
n2  =  struct('method','nescal', 'BC', define_BC('Dirichlet',0), 'dist',dist); 
n3  =  struct('method','nescal', 'BC', define_BC('Dirichlet',0), 'dist',dist); 
vt  =  struct('method','direct', 'BC', define_BC('Dirichlet',0), 'dist',dist);

%% initialization
init = struct(...
    'te0'    , 200,... % Te at axis
    'tewidth', 0.6,... % Te gaussian witdh
    'ti0'    , 200,... % Ti at axis
    'tiwidth', 0.6,... % Ti gaussian width
    'ne0'    , 1e19,...% etc etc
    'newidth', 0.8,...
    'ni0'    , 1e19,...
    'niwidth', 0.8,...
    'n1scal'    ,0,...
    'n2scal'    ,0,...
    'n3scal'    ,0,...
    'ze0'    , 3.5,... % Zeff at axis
    'zee'    , 3.5,... % Zeff at edge
    'vt0'    , 0,...
    'vte'    , 0,...
    'jpow'   , 4, ...   % polynomial initial current density profile: j0 =  1-(rho).^jpow
    'Ip0'    , 80e3,...
    'psib0'  , 0.0... % initial psi at edge
    );


%%% environment options
envopts = struct(...
    'RAPTOR_path'   , RAPTOR_get_path , ...     % determine RAPTOR root path
    'equil_path'    , '', ... % Path to equilibrium data
    'data_path'     , fullfile( RAPTOR_get_path, 'personal'), ... % Path to external data, except equilibrium data
    'call_path'     , pwd);
  
%%%% Collect all in main params structure
config = struct( ...
    'grid',grid,...
    'equi',equi,... % equilibrium settings
    'atom',atom,... % atomic physics settings
    'numerics', numerics, ... % numerics
    'buildopts',buildopts, ... % build options
    'controllers', RAPTOR_controllers, ...
    'diagnostics',RAPTOR_diagnostics,...
    'psi',psi,...
    'te',te,...
    'ti',ti,...
    'ne',ne,...
    'ni',ni,...
    'ze',ze,...
    'n1',n1,...
    'n2',n2,...
    'n3',n3,...
    'vt',vt,...
    'echcd',echcd,...
    'ichcd',ichcd,...
    'nbhcd',nbhcd,...
    'lhhcd',lhhcd,...
    'runaways',runaways,...
    'cvp',RAPTOR_cvp,...
    'chi_e', chi_e, ...
    'vpdn_e', vpdn_e,...`
    'sn_e', sn_e,...
    'hmode', hmode, ...
    'ped','',...
    'saw', saw, ...
    'ntm', ntm, ...
    'neos', neos,...
    'palpha',palpha, ...
    'pbrem',pbrem, ...
    'prad',prad, ...
    'pei',pei, ...
    'envopts',envopts,...
    'debug', debug,...
    'realtime',false,...
    'init',init...
    );

%%%%%%%%%%%%%%%%%%%%%%%
%Default output data structure

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Input processing

if nargout>1
    error('This use is deprecated, please see the RAPTOR_tutorial_1 to learn how prepare RAPTOR runs')
end

if nargin == 0
    if config.debug.verbosity
        disp('calling default RAPTOR_config(''TCV'')');
    end
    config = RAPTOR_config('TCV'); % default
    varargout{1} = config;
    return %empty call, probably to get default structures
elseif nargin==1 && ischar(varargin{1})
    % call to get a customized configuration for a given tokamak
    config = device_config(config,varargin{1});
    varargout{1} = config;
    return
else
    error('This use is deprecated, please see the RAPTOR_tutorial_1 to learn how prepare RAPTOR runs')
end

end

function ppp = RAPTOR_get_path

pp = mfilename('fullpath'); % path of this RAPTOR_config file
ii=strfind(pp,['code' filesep 'initialize']);
if isempty(ii) 
    error('script run outside RAPTOR path');
else
    ppp = [pp(1:ii-1)];
end

end


function c = device_config(config,devicestr)

c = config;

switch devicestr
    case 'ITER'
        c.equi.tokamak = 'ITER';
        c.equi.filenames = {'ITER_hybrid_citrin_equil_cheasedata'}; % Reference chease equilibrium .mat file
        c.equi.B0 = 5.3; % Toroidal field [T]
        c.equi.R0 = 6.2; % Major radius [m]
        c.equi.epsilon = 0.32; % inverse aspect ratio
        c.equi.kappa = 1.87; % default scenario kappa
        c.equi.delta = 0.3; % default scenario delta
        c.equi.Lext =  1e-6; % external inductance (guess) [Henry]
        
        c.grid.rhogrid = linspace(0,1,31); % radial grid
        
        c.te.BC.defaultValue = 400;
        
        c.init.newidth = 2;
        c.init.niwidth = 2;
        c.init.te0 = 4e3;
        c.init.ti0 = 4e3;
        c.init.tewidth = 0.6;
        c.init.ze0 = 1.6;
        c.init.zee = 1.6;
        c.init.ne0 = 2e19;
        c.init.ni0 = 2e19;
        c.init.jpow = 1;
        c.init.Ip0 = 4e6;
        
        c.grid.tgrid = [20:0.05:21-0.05,21:0.5:100];
        c.chi_e = chi_e('BgB');
        
        c.numerics.restol = 1e-8;
        
        c.palpha.active = true;
        c.pbrem.active = true;
        c.pei.active = true;
        c.prad.params.active = false;
        
        c.saw.params.active=false;
        
        c.neos.jbsB_fact = 1.15;
        
        c.echcd = echcd('echcd_gaussian');
        c.echcd.params.active = true;
        c.echcd.params.rdep = [0.4 0.55 .2]; % two actuators, set rdep=0 for first, rdep=0.4 for second
        c.echcd.params.wdep = [.05 .05 .05];  % wdep =0.35 for each
        c.echcd.params.cd_eff = [10 11 10];  % current drive efficiency factor: pure ECH for first, co-ECCD for second
        c.echcd.params.uindices = uint32([2 3 4]); % index of power for each actuator in input vector
        
        c.nbhcd = nbhcd('nbhcd_gaussian');
        c.nbhcd.params.active = true;
        c.nbhcd.params.rdep = [0.26]; %
        c.nbhcd.params.wdep = [0.15]; % broad heating
        c.nbhcd.params.wdep_out = [0.75]; % broad heating
        c.nbhcd.params.cd_eff = [4]; % current drive
        c.nbhcd.params.uindices = uint32([5]);% index in input vector
        
    case 'TCV'
        c.nbhcd.params.uindices = uint32([4]);% index in input vector to avoid conflict w ECH
    case 'JET'
        c.equi.tokamak = 'JET';
        c.equi.filenames = {'JET_83224_t46_cheasedata'}; % Reference chease equilibrium .mat file
        c.equi.B0 = 2.4; % Toroidal field [T]
        c.equi.R0 = 2.875; % Major radius [m]
        c.equi.epsilon = 1.25/c.equi.R0; % inverse aspect ratio
        c.equi.kappa = 2.10/1.25; % default scenario kappa
        c.equi.delta = 0.3; % default scenario delta
        c.equi.Lext =  1e-6; % external inductance (guess) [Henry]
        
        c.grid.rhogrid = linspace(0,1,21); % radial grid
        
        c.te.BC.defaultValue = 120;
        
        c.init.te0 = 0.5e3;
        c.init.tewidth = 0.6;
        c.init.ti0 = 0.5e3;
        c.init.tiwidth = 0.6;
        c.init.ne0 = 1.6e19;
        c.init.newidth = 1;
        c.init.ni0 = 1.6e19;
        c.init.niwidth = 1;
        c.init.ze0 = 1.1;
        c.init.zee = 1.1;
        c.init.jpow = 2;
        c.init.Ip0 = 2e6;
        
        c.grid.tgrid = [41.5:0.005:43-0.005,43:0.5:52];
        c.chi_e = chi_e('BgB');
        c.numerics.restol = 1e-8;
        
        c.palpha.active = false;
        c.pbrem.active = true;
        c.pei.active = true;
        c.prad.params.active = false;
                
        c.nbhcd = nbhcd('nbhcd_gaussian');
        c.nbhcd.params.rdep = [0]; % deposition in centre
        c.nbhcd.params.wdep = [0.4]; % broad heating
        c.nbhcd.params.wdep_out = [0.25]; % broad heating
        c.nbhcd.params.cd_eff = [0]; % no current drive
        c.nbhcd.params.uindices = uint32([2]);% index in input vector
        
        c.echcd = echcd('none');
        
    case {'WEST'}
        c.equi.tokamak = 'WEST';
        c.equi.filenames = {'WEST_cheasedata.mat'};
        c.equi.B0 = 3.6; % Toroidal field [T]
        c.equi.R0 = 2.45; % Major radius [m]
        c.equi.epsilon = 0.45/c.equi.R0; % inverse aspect ratio
        c.equi.kappa = 1.3635; % default scenario kappa
        c.equi.delta = 0.3205; % default scenario delta
        c.equi.Lext =  1e-6; % external inductance (guess) [Henry]
        c.grid.rhogrid = linspace(0,1,11); % radial grid
        c.grid.tgrid = [0:0.01:10];
        c.te.BC.defaultValue = 120;
        c.init.te0 = 500;
        c.init.tewidth = 0.6;
        c.init.ti0 = 200;
        c.init.tiwidth = 0.6;
        c.init.ne0 = 2e19;
        c.init.newidth = 2;
        c.init.ni0 = 2e19;
        c.init.niwidth = 2;
        c.init.ze0 = 2.5;
        c.init.zee = 2.5;
        c.init.jpow = 3;
        c.init.Ip0 = 5e5;
        

        c.chi_e = chi_e('BgB');
        c.numerics.restol = 1e-8;
        
        c.palpha.active = false;
        c.pbrem.active = true;
        c.pei.active = true;
        c.prad.params.active = true;
        
        c.echcd = echcd('echcd_gaussian');
        c.echcd.params.active = true;
        c.echcd.params.rdep = [0.0 0.0 -1 -1];
        c.echcd.params.wdep = [0.3 0.3 0.3 0.3 ];
        c.echcd.params.cd_eff =[0  0   0   0   ];
        c.echcd.params.uindices = uint32(1 + (1:4)); % (8 heating + 4 rhodep)
        
        c.nbhcd = nbhcd('nbhcd_gaussian'); % params fit to transp results Ben Geiger for shot 30387
        c.nbhcd.params.active = true;
        c.nbhcd.params.rdep     = [0.05  0.1  0.05  0.05 0.17 0.43 0.37 0.08];
        c.nbhcd.params.wdep     = 2*[0.2 0.3  0.3   0.3  0.3  0.18  0.20 0.25]; % 2*half width
        c.nbhcd.params.wdep_out = 2*[0.4 0.35 0.35   0.35  0.35  0.29 0.30 0.35];
        c.nbhcd.params.cd_eff =   [0.5119    0.8123    0.8670    0.5539    0.8647    1.3500    1.5853    0.9539];
        c.nbhcd.params.eta_nbi0 = 1/5; % to get about 100kA per 5MW beam power
        c.nbhcd.params.frac_to_electrons  = [1.6926    1.6760    1.6396    1.6527    1.9848    2.4017    2.4039    1.9844]/5; % MW to electrons per MW input power
        c.nbhcd.params.uindices = uint32(5 + (1:8));
                
%         c.lhcd = lhcd('lhcd_gaussian');
%         c.lhcd.params.active = true;
%         c.lhcd.params.uindices = 13+[1:2]; % 2 heatings
%              
    case {'AUG'}
        c.equi.tokamak = 'AUG';
        c.equi.filenames = {'AUG_30594_t4_cheasedata.mat'};
        c.equi.B0 = 2.5; % Toroidal field [T]
        c.equi.R0 = 1.64; % Major radius [m]
        c.equi.epsilon = 0.5/c.equi.R0; % inverse aspect ratio
        c.equi.kappa = 1.485; % default scenario kappa
        c.equi.delta = 0.0804; % default scenario delta
        c.equi.Lext =  1e-6; % external inductance (guess) [Henry]
        
        c.grid.rhogrid = linspace(0,1,21); % radial grid
        
        c.te.BC.defaultValue = 120;
        
        c.init.te0 = 500;
        c.init.tewidth = 0.6;
        c.init.ti0 = 500;
        c.init.tiwidth = 0.6;
        c.init.ne0 = 1e19;
        c.init.newidth = 0.8;
        c.init.ni0 = 1e19;
        c.init.niwidth = 0.8;
        c.init.ze0 = 2;
        c.init.zee = 2;
        c.init.jpow = 4;
        c.init.Ip0 = 0.5e6;
        
        c.grid.tgrid = [0:0.01:10];
        c.chi_e = chi_e('BgB');
        c.numerics.restol = 1e-8;
        
        c.palpha.active = false;
        c.pbrem.active = true;
        c.pei.active = true;
        c.prad.params.active = false;
               
        c.echcd = echcd('echcd_gaussian');
        c.echcd.params.active = true;
        c.echcd.params.rdep = [0.0 0.0 0.0 0.0 -1 -1 -1 -1]; %use rhodep, not rdep
        c.echcd.params.wdep = [0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3];
        c.echcd.params.cd_eff =[0  0   0   0   0   0   0   0  ];
        c.echcd.params.uindices = uint32(1 + (1:12)); % (8 heating + 4 rhodep)
        
        c.nbhcd = nbhcd('nbhcd_gaussian'); % params fit to transp results Ben Geiger for shot 30387
        c.nbhcd.params.active = true;
        c.nbhcd.params.rdep     = [0.05  0.1  0.05  0.05 0.17 0.43 0.37 0.08];
        c.nbhcd.params.wdep     = 2*[0.2 0.3  0.3   0.3  0.3  0.18  0.20 0.25]; % 2*half width
        c.nbhcd.params.wdep_out = 2*[0.4 0.35 0.35   0.35  0.35  0.29 0.30 0.35];
        c.nbhcd.params.cd_eff =   [0.5119    0.8123    0.8670    0.5539    0.8647    1.3500    1.5853    0.9539];
        c.nbhcd.params.eta_nbi0 = 1/5; % to get about 100kA per 5MW beam power
        c.nbhcd.params.frac_to_electrons  = [1.6926    1.6760    1.6396    1.6527    1.9848    2.4017    2.4039    1.9844]/5; % MW to electrons per MW input power
        c.nbhcd.params.uindices = uint32(13 + (1:8));
        
       case {'RFX'}
        c.equi.tokamak = 'RFX';
        c.equi.filenames = {'RFX_11111_cheasedata.mat'};
        c.equi.B0 = 0.5; % Toroidal field [T]
        c.equi.R0 = 2; % Major radius [m]
        c.equi.epsilon = 0.5/c.equi.R0; % inverse aspect ratio 
        c.equi.kappa = 1.0; % default scenario kappa
        c.equi.delta = 0.0; % default scenario delta
        c.equi.Lext =  1e-6; % external inductance (guess) [Henry]
        
        c.grid.rhogrid = linspace(0,1,11); % radial grid
        
        c.te.BC.defaultValue = 50;
        
        c.init.te0 = 300;
        c.init.tewidth = 0.65;
        c.init.ti0 = 300;
        c.init.tiwidth = 0.65;
        c.init.ne0 = 2.7e19;
        c.init.newidth = 0.8;
        c.init.ni0 = 2.7e19;
        c.init.niwidth = 0.8;
        c.init.ze0 = 2;
        c.init.zee = 2;
        c.init.jpow = 4;
        c.init.Ip0 = 160e3;

        c.grid.tgrid = [0:0.001:0.5];
        c.chi_e = chi_e('FF');
        c.numerics.restol = 1e-8;

        c.palpha.active = false;
        c.pbrem.active = true;
        c.pei.active = false;
        c.prad.params.active = true;

        c.echcd = echcd('none');
        c.nbhcd = nbhcd('none'); % params fit to transp results Ben Geiger for shot 30387

case {'DEMO'}
        c.equi.tokamak = 'DEMO';
        c.equi.filenames = {'DEMO/chease_mat_SOF.mat'};
        % infos from DEMO1 reference 12 May 2015 ppt METIS
        c.equi.B0 = 5.667; % Toroidal field [T]
        c.equi.R0 = 9.073; % Major radius [m]
        c.equi.epsilon = 2.926/c.equi.R0; % inverse aspect ratio
        c.equi.kappa = 1.6994; % default scenario kappa
        c.equi.delta = 0.33; % default scenario delta
        c.equi.Lext =  1e-6; % external inductance (guess) [Henry]
        
        c.grid.rhogrid = linspace(0,1,11); % radial grid
        
        c.te.BC.defaultValue = 100;
        
        c.init.newidth = 2;
        c.init.niwidth = 2;
        c.init.te0 = 25e3;
        %         c.init.ti0 = 0.95*25e3;
        c.init.tewidth = 0.6;
        c.init.ze0 = 1.984;
        c.init.zee = 1.984;
        c.init.ne0 = 10e19;
        c.init.ni0 = 9.15e19;
        c.init.jpow = 1;
        c.init.Ip0 = 1e6;
        %         c.init.psib0 = 10;
        %         c.grid.tgrid = [20:0.05:21-0.05,21:0.5:100];
        c.chi_e = chi_e('BgB');
        
        c.numerics.restol = 1e-8;
        
        c.palpha.active = true;
        c.pbrem.active = true;
        c.pei.active = true;
        c.prad.params.active = true;
        
        c.saw.params.active=false;
              
        c.echcd = echcd('echcd_gaussian');
        c.echcd.params.active = true;
        c.echcd.params.rdep = [0.4]; 
        c.echcd.params.wdep = [.1];  
        c.echcd.params.cd_eff = [1];  
        c.echcd.params.uindices = uint32([2]); 
        
        c.nbhcd = nbhcd('nbhcd_gaussian');
        c.nbhcd.params.active = true;
        c.nbhcd.params.rdep = [0.26]; 
        c.nbhcd.params.wdep = [0.15]; 
        c.nbhcd.params.wdep_out = [0.75]; 
        c.nbhcd.params.cd_eff = [1]; 
        c.nbhcd.params.uindices = uint32([3]);
     
    otherwise
        error('unknown tokamak option ''%s''',devicestr)
end
end

function BC = define_BC(type,defaultValue)
% default
BC  =  struct('type','Dirichlet',...
    'rhoSource','Constant',...
    'rho',1,...
    'valueSource','v',...
    'defaultValue',0); % default BC setting for kinetic equations

if ~isempty(strfind(lower(type),'dirichlet'))
    BC.type = 'Dirichlet';
elseif ~isempty(strfind(lower(type),'neumann'))
    BC.type = 'Neumann';
    warning('not supported yet')
else
    BC.type = BC.type;
end

BC.defaultValue = defaultValue;

end
