% Badness function for shooting to model solar coronal loops using loop.m
function badness = looperror(parameters, plot_title)
% If plot_title exists, then the loop will be plotted.



% Globals to be set by the main program
global desired_length      % Upper BC
global desired_temperature % Upper BC
global Tchromosphere       % Lower BC
global gradTchromosphere   % Lower BC
global kappa

% Declare global parameters (to share with loop.m)
global T0
global ne0
global EH

% Translate input parameters into cgs loop parameters
P0 = parameters(1) * 1e15;  % Pressure is nT rather than 2nkT, so 1e15 
   % corresponds to T = 1e6 K, and n = 1e9 cm^-3.
EH = parameters(2) * 1e-4;  % erg cm^-3 sec^-1

T0 = 1e6;
ne0 = P0/T0; % Again, pressure is nT rather than 2nkT


% Position array
Npoints = 1000;
s1 = 1e5; % Minimum value of coordinate axis (1e5 cm = 1 km)
sN = desired_length + s1; % Maximum value of coordinate axis
s = exp(linspace(log(s1), log(sN), Npoints));


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

% Evaluate the merit function (badness)
badness = ( (T(Npoints) - desired_temperature) / desired_temperature )^2 +...
          ( gradT(Npoints) * desired_length / desired_temperature )^2

if exist('plot_title','var')
   plot(s,T)
   xlabel('s (cm)','fontsize',18)
   ylabel('T (K)','fontsize',18)
   title(plot_title,'fontsize',20)
end