% Work out the correct parameters for a loop of length L, temperature Ttop,
% and (optionally) footpoint conductive flux FCfoot. The conductive flux
% is normally zero or negative.

function shootloop(L, Ttop, FCfoot)

% Global variables to be shared with looperror.m
global desired_length      % Upper BC
global desired_temperature % Upper BC
global Tchromosphere       % Lower BC
global gradTchromosphere   % Lower BC
global kappa

% Set BCs at the loop top
desired_temperature = Ttop
desired_length = L

% Set BCs at the footpoint.
kappa = 1e-6            % Plasma thermal conductivity (normally 1e-6)
Tchromosphere = 2e4  
if !exist('FCfoot', 'var')
   FCfoot = 0 % Default
end
gradTchromosphere = -FCfoot/(kappa*Ttop^(5/2))
   % positive FCfoot <=> negative gradT
   
% ODE solver options
lsode_options('integration method','adams');  % default solver is 'bdf'.
lsode_options('maximum order',-1);            % -1 means unconstrained.

% Generate first guess using Rosner et al. scaling laws (1978)
Pguess = desired_temperature^3 * desired_length / 1.6e12
EHguess = Pguess^(7/6) * desired_length^(-5/6) / 1.5e14
prams_guess = [Pguess/1e15, EHguess/1e-4]

% Plot results of first guess
figure(1)
badness_initial_guess = looperror(prams_guess, 'First Guess')


% Optmize the parameters for the desired length & temperature

bestfit_prams = amoeba('looperror', prams_guess, 100)

% Plot results of optimization
figure(2)
badness_bestfit = looperror(bestfit_prams, 'Best Fit')


% Translate optimized parameters into cgs units
P0_best = bestfit_prams(1) * 1e15  % Pressure is nT rather than 2nkT, so 1e15 
   % corresponds to T = 1e6 K, and n = 1e9 cm^-3.
EH_best = bestfit_prams(2) * 1e-4  % erg cm^-3 sec^-1
