%% RAPTOR tutorial 6: Specifying boundary conditions

clear; close all hidden;
run ../RAPTOR_path.m

[config] = RAPTOR_config; % load default config

%% BC structure
% boundary conditions are set via the BC structure of each species
disp(config.ne.BC); 

%% Solve ne equation, and prescribe fixed boundary condition at pedestal top 
config.ne.method = 'state';
config.ne.BC.type = 'Dirichlet'; % VolumeAverage is also possible
config.ne.BC.defaultValue =  1e19;
config.ne.BC.rho = 0.9; % pedestal top at this rho value;

config.ni.method = 'nescal';

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

x0 = RAPTOR_initial_conditions(model,init,g(:,1),v(:,1)); 
simres = RAPTOR_predictive(x0,g,v,U,model,params);
out = RAPTOR_out(simres,model,params);

%% Time-varying boundary condition value and location
config.ne.method = 'state';
config.ne.BC.type = 'Dirichlet';
config.ne.BC.defaultValue =  1e18;
config.ne.BC.rhoSource = 'v';
config.ne.BC.valueSource = 'v';
config.ni.method = 'nescal';

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

% input time-varying boundary condition value and location via v vector
V = repmat(v,1,numel(params.tgrid));

rho_BC = 1 - 0.05*(1-cos(2*2*pi*params.tgrid/params.tgrid(end)));
val_BC = 1e18*ones(size(params.tgrid)); val_BC(params.tgrid>(params.tgrid(end)/2))=2e18;

V(model.ne.BC.vind_rho,:) = rho_BC;
V(model.ne.BC.vind_value,:) = val_BC; 

% initial condition
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);
%% plot contours including special ones for the boundary conditions
clf;
plot(out.rho,out.ne,'-',rho_BC,val_BC,'xk'); hold on;