%% RAPTOR tuturial on the use of the stationary state solver
% REGULAR RAPTOR RUN -----APPROACHING STATIONARY FINAL STATE --------------
clc
close all hidden;
run(fullfile(pwd, '..', 'RAPTOR_path.m')); % add RAPTOR path

tic
[config] = RAPTOR_config; % load default params
config.grid.rhogrid = linspace(0, 1, 21); % 21 knots (22 splines)
config.grid.tgrid = 0:1e-3:1; % simulate up to 1s (>>tau_crt, tau_crt~150ms for heated TCV plasmas)

config.echcd = RAPTORmodule('echcd_gaussian'); % add ECHCD
config.echcd.Configuration.n_units = 1;
config.echcd.Configuration.n_rdep_external = 0;

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

u(1, :) = 200e3*ones(size(params.tgrid)); 
u(2, :) = 1e6*ones(size(params.tgrid));
init.Ip0 = u(1, 1);

params.echcd.active = true; 
params.echcd.rdep = 0.4;
params.echcd.wdep = 0.2;  
params.echcd.cd_eff = 1;  

x0 = RAPTOR_initial_conditions(model, init, g, v);
simres = RAPTOR_predictive(x0, g, v, u, model, params, 'verbosity', 0);
out = RAPTOR_out(simres, model, params);
toc

itime(1) = find(out.time>=.15,1,'first');
itime(2) = find(out.time>=.30,1,'first');
itime(3) = find(out.time>=.45,1,'first');
itime(4) = numel(out.time);
figure;
subplot(121)
plot(out.rho,out.jpar(:,itime))
xlabel('\rho')
ylabel('j_{par} [A/m^2]')
legend({'t=0.15s','0.30s','0.45s','1s'})
title('time evolution radial j_{par} profile')
subplot(122)
plot(out.rho,out.q(:,itime))
xlabel('\rho')
ylabel('q')
title('time evolution radial q profile')

%% STATIONARY STATE SOLVER ------------------------------------------------
tic
[config] = RAPTOR_config;
config.grid.rhogrid = linspace(0, 1, 21);
config.grid.tgrid = 0; % no time evolution

config.echcd = RAPTORmodule('echcd_gaussian');
config.echcd.Configuration.n_units = 1;
config.echcd.Configuration.n_rdep_external = 0;

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

u(1) = 200e3;   % overwritten if solving for Ip
u(2) = 1e6;     % overwritten if solving for Pec

% --- initial conditions for Newton iterative solver ---
init.Ip0 = u(1); % overwritten if solving for Ip and inputs.Ip provided
init.te0 = 2e3;
init.jpow = 4;
% --> used to construct initial guess x0
% -------------------------------------------------------------------------
% new feature: an initial guess x0 can also be provided as an input to
% stationary_state_solver, as a varargin input, e.g. 
% [x, xdot, u, outs, res_vec] = stationary_state_solver(model, params, init, g, v, u, ...
%                                                       parameter, inputs, verbosity, x0);
% -------------------------------------------------------------------------                                          
params.echcd.active = true; 
params.echcd.rdep = 0.4;
params.echcd.wdep = 0.2;   
params.echcd.cd_eff = 1;

verbosity = 2; 
% 0: no report at all
% 1: global information reported
% 2: iteration step information reported
% 3: plots

scalar_unknown = 'Upl';
% Specify the additional scalar unknown (in addition to the stationary profiles)
% 1:	'Upl'	- initial guess inputs.Upl needed
% 2:	'Ip'	- initial guess inputs.Ip needed (if not provided: init.Ip0 provides initial guess)
% 3:	'Act' - initial guess inputs.Act needed
%             can solve for any element of the u vector, user-specified
%             (u(2) by default) as a varargin argument, e.g.
%             [x, xdot, u, outs, res_vec] = stationary_state_solver(model, params, init, g, v, u, ...
%                                           parameter, inputs, verbosity, 3)
%             to solve for u(3)
% for options 2, 3: inputs.Upl user provides the demanded loop voltage

% inputs structure is used for conditioning of the set of equations and
% provides the initial guess for the Newton solver
inputs.Upl = 1;                     % initial condition should be provided for unknown
% inputs.Ip = 200e3;
% inputs.Act = 1e6;  

[x, xdot, u, outs, res_vec] = stationary_state_solver(model, params, init, g, v, u, ...
                                                      scalar_unknown, inputs, verbosity);
toc
 
fprintf('plasma current: %.3eA\n', u(1));
fprintf('electron cyclotron power: %.3eW\n', u(2));
fprintf('plasma loop voltage: %.4fV\n', xdot(1));

figure;
subplot(121)
plot(out.rho,out.jpar(:,itime), ...
     outs.rho, outs.jpar, 'k*')
xlabel('\rho')
ylabel('j_{par} [A/m^2]')
legend({'t=0.15s','0.30s','0.45s','1s','stationary'})
title('time evolution radial j_{par} profile')
subplot(122)
plot(out.rho,out.q(:,itime), ...
     outs.rho, outs.q, 'k*')
xlabel('\rho')
ylabel('q')
title('time evolution radial q profile')

figure;
plot(out.rho,out.upl(:,itime), ...
     outs.rho, outs.upl, 'k*')
xlabel('\rho')
ylabel('U_{pl} [V]')
legend({'t=0.15s','0.30s','0.45s','1s','stationary'},'Location','SouthEast')
title('time evolution radial U_{pl} profile')
% flat loop voltage: no more drive for jpar diffusion
%%
u(2) = 0; % let's compare to stationary state resulting when no ECHCD is applied
[x2, xdot2, u2, outs2] = stationary_state_solver(model, params, init, g, v, u, ...
                         scalar_unknown, inputs, verbosity);

fprintf('plasma current: %.3eA\n', u2(1));
fprintf('electron cyclotron power: %.3eW\n', u2(2));
fprintf('plasma loop voltage: %.4fV\n', xdot2(1));

figure;
subplot(311)
plot(outs.rho,outs.q, ...
     outs2.rho,outs2.q)
legend({'1MW ECHCD','no ECHCD'})
xlabel('\rho')
ylabel('q')
title('radial q profile')
subplot(312)
plot(outs.rho,outs.upl, ...
     outs2.rho,outs2.upl)
xlabel('\rho')
ylabel('U_{pl}')
title('radial loop voltage profile')
subplot(313)
plot(outs.rho,outs.joh, ...
     outs2.rho,outs2.joh)
xlabel('\rho')
ylabel('j_{oh}')
title('radial ohmic current profile')
% EC current drive -> less ohmic current to be driven inductively and
% higher conductivty -> lower loop voltage
%%
% finally let's calculate how much ECHCD we need to drive all plasma
% current non-inductively, i.e. loop voltage Upl=0V
scalar_unknown = 'Act';
inputs.Upl = 0; % demanding zero loop voltage
inputs.Act = 1e6;  % initial guess
[x3, xdot3, u3, outs3] = stationary_state_solver(model, params, init, g, v, u, ...
                                           scalar_unknown, inputs, verbosity);

fprintf('plasma current: %.3eA\n', u3(1));
fprintf('electron cyclotron power: %.3eW\n', u3(2));
fprintf('plasma loop voltage: %.4fV\n', xdot3(1));   

figure;
plot(outs3.rho, outs3.upl)
title('radial loop voltage profile'); 
xlabel('\rho')
ylabel('U_{pl} [V]')
%% CALCULATE INITIAL CONDITIONS -------------------------------------------
% the idea here is to find a set of consistent initial profiles to initiate
% this illustrates how the non-linear solver can be used generically, i.e.
% not only to calculate stationary states with flat loop voltage
% a RAPTOR run
%       'init state'    - state with an xdot mimicing initial conditions
%                       - Upl sigmoid, inward diffusion 
%                       - radially flat dTe_dt profile of unknown magnitude 
%                         --> initial guess xdotprofiles.dTe_dt needed
[config] = RAPTOR_config;

config.grid.rhogrid = linspace(0, 1, 21);
config.grid.tgrid = 0;

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

u(1) = 80e3;  

% initial conditions for Newton solver
init.Ip0 = u(1);                        

verbosity = 2; 

Upl_core = 3;
Upl_edge = 5;
xdot_profiles.Upl = @(rho) Upl_core + (Upl_edge - Upl_core) ./ (1 + exp(-10*rho + 5));
xdot_profiles.dTe_dt = 100;

[~, xdot_init, u_init, out_init, ~] = initial_state_solver(model, params, init, g, v, u, ...
                                             xdot_profiles, verbosity);
fprintf('plasma current: %.3eA\n', u_init(1));
fprintf('electron cyclotron power: %.3eW\n', u_init(2));
fprintf('dTe_dt on axis: %.3eeV/s\n', xdot_init(end));

figure(1)
subplot(311); plot(out_init.rho, out_init.te); xlabel('\rho'); ylabel('T_e [eV]'); title('radial temperature profile'); 
subplot(312); plot(out_init.rho, out_init.q); xlabel('\rho'); ylabel('q'); title('radial q profile'); 
subplot(313); plot(out_init.rho, out_init.upl); xlabel('\rho'); ylabel('U_{pl} [V]'); ylim([0 7]); title('radial loop voltage profile'); 