% Map out the phase space of a loop with a given pressure and heating rate.
% CCK 2009-Apr-02  Based on looptest2.m
% function loopPhaseSpace

% ODE solver options
lsode_options('integration method','adams')% default solver is 'bdf'='stiff'.
lsode_options('maximum order',-1)          % -1 means unconstrained.
lsode_options                              % display all the settings.

% Set the global parameters
global T0
global ne0
global EH
global kappa
P0 = 1e16 % Pressure expressed as ne*T
T0 = 1e6
ne0 = P0/T0
EH = 1e-3
kappa = 1e-6; % Plasma thermal conductivity (normally 1e-6)

% Position array
Npoints = 100
s1 = 1e5 % Minimum value of coordinate axis (1e5 cm = 1 km)
sN = 1e9 % 1.2e9 % Maximum value of coordinate axis (1e8 cm = 1 Mm)
s = exp(linspace(log(s1), log(sN), Npoints));


% Integrate from loop top to footpoint
% FootpointTemperatures = ((1:6) - 1/2)*1e5 ;
FootpointTemperatures = ((2:7) / 5)*1e5 ;
colors = ['k-';'r-';'g-';'c-';'b-';'m-'];
gradTfoot = 0.00
hold off
for i=1:6
   Tfoot = FootpointTemperatures(i)
   q0 = [Tfoot, gradTfoot];
   q = lsode(@loop, q0, s);
   T = q(:, 1);
   gradT = q(:, 2);
   plot(T/1e6, gradT*1e8/1e6, colors(i,:), 'linewidth', 2)
   hold on
   plot(T/1e6, -gradT*1e8/1e6, colors(i,:), 'linewidth', 2)
end
text(0.2, 3, "<---")
text(0.27, 3.3, "Trajectories diverge with decreasing T.\nAdvice: Don't integrate that direction!")
xlabel('T (MK)')
ylabel('dT/ds (MK/Mm)')
title('Loop Models: E_H = 10^{-3}, n_eT = 10^{16}')
axis([0,1.1,-1,15])
print('loopPhaseSpace.pdf','-dpdf')