%% Tutorial showing some options for how to specify kinetic profiles
% Kinetic profiles that are not solved for (i.e. are evolved as part of the
% state vector |x|, can be prescribed by other means. Either by scaling
% another profile, or by direct specification. 
% This tutorial provides some examples
%%
run(fullfile(pwd,'..','RAPTOR_path.m')); % add RAPTOR path
close all hidden; clear;
[config] = RAPTOR_config; % load default params
config.debug.iterdisp = 0; % nice and quiet

%% Basic default case
config.ti.method
config.ni.method

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

U(1,:) = 80e3*ones(size(params.tgrid)); % input Ip trace: constant 80kA
init.Ip0 = U(1); % initial condition for current
x0 = RAPTOR_initial_conditions(model,init,g(:,1),v(:,1)); % Define the initial condition
simres = RAPTOR_predictive(x0,g,v,U,model,params);
out = RAPTOR_out(simres,model,params);

%% ti defined as independent state
config.ti.method = 'state';
config.ti.BC.type = 'Dirichlet';
config.ti.BC.defaultValue = config.te.BC.defaultValue;
[model,params,init,g,v,U] = build_RAPTOR_model(config); % generate model structure for these params.

U(1,:) = 80e3*ones(size(params.tgrid)); % input Ip trace: constant 80kA
init.Ip0 = U(1); % initial condition for current
x0 = RAPTOR_initial_conditions(model,init,g(:,1),v(:,1)); % Define the initial condition
simres = RAPTOR_predictive(x0,g,v,U,model,params);
out = RAPTOR_out(simres,model,params);

%% ti scaled with respect to te
config.ti.method = 'tescal';
[model,params,init,g,v,U] = build_RAPTOR_model(config);
disp(init)

init.ti0scal = 0.9;
init.tiescal = 1;
% Construct kinetic profiles according to the new settings.
v = build_kinetic_profiles(model,init);

U(1,:) = 80e3*ones(size(params.tgrid)); % input Ip trace: constant 80kA
init.Ip0 = U(1); % initial condition for current
x0 = RAPTOR_initial_conditions(model,init,g(:,1),v(:,1)); % Define the initial condition
v(model.ti.vind)

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

plot(out.rho,[out.te(:,end),out.ti(:,end)]);
disp(out.ti(:,end)'./out.te(:,end)')

%% ni scaled w.r.t. ne
config.ni.method = 'nescal';
[model,params,init,g,v,U] = build_RAPTOR_model(config);
disp(init)
init.ni0scal = 0.9;
init.niescal = 1;
% Construct kinetic profiles according to the new settings.
v = build_kinetic_profiles(model,init);

U(1,:) = 80e3*ones(size(params.tgrid)); % input Ip trace: constant 80kA
init.Ip0 = U(1); % initial condition for current
x0 = RAPTOR_initial_conditions(model,init,g(:,1),v(:,1)); % Define the initial condition
simres = RAPTOR_predictive(x0,g,v,U,model,params);
out = RAPTOR_out(simres,model,params);

plot(out.rho,[out.ne(:,end),out.ni(:,end)]);
disp(out.ni(:,end)'./out.ne(:,end)')
