% Try out the coronal loop model (loop.m)
function looptest2

% ODE solver options
lsode_options('integration method','bdf')  % default solver is 'bdf'.
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 = 5e16 % Pressure expressed as ne*T
T0 = 1e6
ne0 = P0/T0
EH = 2e-2
kappa = 1e-6; % Plasma thermal conductivity (normally 1e-6)

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


% Integrate from loop footpoint to top
Tchromosphere = 10^4.5  % Radiative loss model is invalid below 10^4.5 K.
gradTchromosphere = 0.00
FCchromosphere = kappa * Tchromosphere^(5/2) * gradTchromosphere
   % Conductive flux at the base temperature. Everything conducted
   % through the 1e5K isotherm should be radiated in H Ly alpha.
   % Vernazza & Reeves (1978) reported H Ly alpha values up to
   % 2.6e5 erg cm^-2 s^-1 sr^-1. If I assume that emission is optically
   % thin, then it would come out over 2pi steradians, making it
   % 1.6e6 erg cm^-2. This is probably too large for FC, as I think
   % limb brightening is not generally observed in H Ly alpha.

q0 = [Tchromosphere, gradTchromosphere];
q = lsode(@loop, q0, s);
T = q(:, 1);

% Integrate the other way.
q0_backwards = q(Npoints,:);
s_backwards = s(Npoints:-1:1);
q_backwards = lsode(@loop, q0_backwards, s_backwards);
T_backwards = q_backwards(:, 1);

% Compare the results.
figure(1)
plot(s, T, 'b-', s_backwards, T_backwards,'g-')
legend('footpoint to top', 'top to footpoint')
legend('location','southeast')

figure(2)
loglog(s, T, 'b-', s_backwards, T_backwards, 'g-')
legend('footpoint to top', 'top to footpoint')
legend('location','southeast')

figure(3)
plot(s, T, 'b-', s_backwards, T_backwards,'g-')
legend('footpoint to top', 'top to footpoint')
legend('location','southeast')
axis([s1, 10*s1, Tchromosphere, 3*Tchromosphere])