%% Exercise 5.1 - Simulating the plasma 0D energy balance
clear
close all

%%
% Define constants and wrap in struct
p.n = 1E20; % density [#/m^3]
p.R_0 = 8;
p.a = 2;
p.kappa = 2;
p.B_0 = 7;
p.A = 2;
p.Z_eff = 1.5;
p.f_DT = 0;
p.I_p = 15; % plasma current [MA]

% Define time span, input signal and initial condition
t_span = [0 200]; % time span [s]
P_aux = [10E6*ones(1,800) 25E6*ones(1,800) 50E6*ones(1,801)]; % apply 20 MW input power
time = linspace(t_span(1),t_span(end),size(P_aux,2)); % time grid for input power
T_0 = 1E3; % initial temperature [eV]

% Simulate thermal energy balance and compute power density sources and pressure
[~,T] = ode45(@(t,T) thermal_energy_balance(t,T,time,P_aux,p),time,T_0);
S = sources(T,P_aux',p);

%%
% Plot temperature evolution
fig = figure;
subplot(2,2,1)
plot(time,T)
grid on; xlabel('Time [s]'); ylabel('Temperature [eV]');
subplot(2,2,3)
plot(time,S.pressure)
grid on; xlabel('Time [s]'); ylabel('Pressure');
subplot(2,2,[2 4])
hold on
plot(time,S.S_aux,'b')
plot(time,S.S_ohm,'k')
plot(time,S.S_alpha,'m')
plot(time,S.S_rad,'g')
plot(time,S.S_cond,'r')
legend('S_{aux}','S_{ohm}','S_{alpha}','S_{rad}','S_{cond}','Location','Northwest')
grid on; xlabel('Time [s]'); ylabel('Power density [W/m^3]');
set(fig,'PaperPositionMode','auto');
print -depsc exercise_5_1_simulation