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

% Set the global parameters
global T0
global ne0
global EH
T0 = 1e6
ne0 = 1e9
EH = 1e-3

% Position array
Npoints = 1000
L = 7.5807e8 % Maximum value of coordinate axis
s = [1:Npoints]*L/(Npoints);


% Integrate from loop top to footpoint
q0 = [T0, 0];
q = lsode(@loop, q0, s);
T = q(:,1);

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

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

figure(2)
s_backwards = s(Npoints:-1:1);
loglog(s_backwards, T, 'b-', s_backwards, T_backwards, 'g-')
legend('top to footpoint','footpoint to top')
legend('location','southeast')
axis([s_backwards(1), s_backwards(Npoints), 1e4,1e6])