% Simplified equations for energy balance in a solar coronal loop. The
% model assumes a quasi-static, constant-width loop in zero gravity.
% The volumetric heating is constant. The radiative loss is simplified 
% (compare Martens, Kankelborg & Berger 2000). Units are cgs.
function dqds = loop(q,s)

% The following parameters need to be declared and initialized externally:
global ne0  % An arbitrary electron density value somewhere in the loop.
global T0   % Temperature of the loop corresponding to density ne0.
global EH    % Volumetric heating rate.
global kappa % Conduction constant (usually 1e-6)
T = q(1);
gradT = q(2);


% Radiative loss function
Lambda = 3e-21 * T^(-1/2);

% Differential equations (hermal conductivity of kappa*T^(5/2) divided out):
dqds(1) = gradT;
dqds(2) = (ne0 * T0)^2 * Lambda / (kappa * T^(9/2)) ... % Radiative loss term
          - (5/2)* gradT^2 / T ... % Thermal conduction term (this plus the LHS)
          - EH/(kappa * T^(5/2));   % Heating term
