function Apatite_model_non_ideal_V1
clear all; %Clear all variables before running the script
close all; %Close all figures before running the script

%This script runs the Apatite model
% This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 
% research and innovation programme (grant agreement No 864923) through a Consolidator grant awarded to PI M. Humphreys. 
% The model was developed as part of the STEMMS project (STorage and Eruption of Mushy Magma Systems).
% 
% How to cite this model:
% 
% Lormand, C., Humphreys, M.C.S., Coumans, J.P., Colby, D. J., Chelle-Michou, C. & Li, W. 
% Volatile budgets and evolution in porphyry-related magma systems, determined using apatite. Submitted to Lithos (2023)
% 
%  The model is built on that published by Humphreys et al. (2021), 
% Rapid pre-eruptive mush reorganisation and atmospheric volatile emissions from the 12.9 ka Laacher See eruption, determined using apatite.
% Earth and Planetary Science Letters 576, 117198. https://doi.org/j.epsl.2021.117198


%Determine if PC, MAC or UNIX machine for the delimeter
if ismac
    % Code to run on Mac plaform
    dlm = '/';
elseif isunix
    % Code to run on Linux plaform
elseif ispc
    dlm = '\';
else
    disp('Platform not supported')
end

%====================================================
%Gets the matlab filename
mfilename;
%Gets the folder (local computer) where the script is saved
mpath = strrep(which(mfilename),[mfilename '.m'],'');
mFolder = mpath;
%Adds the folder to the path so that local functions can be called
addpath(mFolder);
%===================================================


%+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!
%Choose which operation to run.
%+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!

%================================!!!
operation = 'Explore model and/or create multistart target';
%operation = 'Run multistart to a target';
%operation = 'Visualize multistart outputs';

%================================!!!


%+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!
%Choose which dataset to open and define the names of the data folder
%(where the input data is located) and the save folder (target and
%multistart outputs)
%+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!
DataFolder = [];
SaveFolder = 'TestOutputTF';
Dataset = 'Test';
Model = 'Model1';

%+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!
% Set Temperature
%+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!
% Set constant temperature
T = 1135; %Set a constant Temperature in oC 

%Read in TF curve
%TF_File = 'Paricutin_Aster_TF.csv';
%T = readtable(TF_File);

%+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!
%Operation switches
%+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!
%The switch which either runs the parameter exploration or the model
switch operation
    case 'Explore model and/or create multistart target'

        %Determine if reading a FT datafile or using a hardwired temperature
        %Determine the temperature vector based on the input

        %If the F-T file cannot be found, set a fixed temperature
        if istable(T) == 0
            T =T+273.15; %T in Kelvins
            F_start_Multi = 1; %Fraction of melt starting
            F_end_Multi = 0.1; %Fraction of melt ending
            TF_File = 'none';
            

        %Else read the file and get the start and final melt fraction
        elseif istable(T) ==1
            %T = readtable(char([mFolder , Model, '.csv']));
            [C,ia,ic] = unique(T.Var2, 'stable');
            T = T(ia,:);

            F_start_Multi = max(T.Var2); %Fraction of melt starting
            F_end_Multi = max([0.1, min(T.Var2)]); %Fraction of melt ending
        end

        %==================================================================
        %Set the target model properties
        %==================================================================
        num_frac = 100; %Number of fractionation steps for the model


        %Initial volatile concentrations
        Cl_init = 0.12;
        F_init = 0.1;
        H2O_init = 1.7;
        
        %Initial volatile crystal-melt distribution coefficients
        D_xl_m_Cl = 0.05;
        D_xl_m_F = 0.1;
        D_xl_m_H2O = 0.12;
        
        %Initial volatile fluid-melt distribution coefficients
        D_fl_m_Cl = 40;
        D_fl_m_F = 2;
        
        %Water saturation in wt. %
        H2O_Sat = 5.6;

        %Speciation for OH-H2O
        K_speciation = 0.363;

        %Primary Boiling e.g. for every 1% (0.01) xtal remove 0.25
        %(0.0025) water
        %For D_fl_m_PrimBoilH2O = 0.25, remove 0.25 wt% H2O
        D_fl_m_PrimBoilH2O=0.0; %Primary boiling off

        %Fractionation information
        X = [F_start_Multi, F_end_Multi, num_frac];


        %The parameters to input to the model, note the K_xx_yy are removed
        %from the input parameters because they are determined by the
        %non-ideal model in the apatite-crystallization section. The
        %calculated equilibrium coefficients are output with the rest of
        %the data.
        InputParam = [H2O_init Cl_init F_init D_xl_m_Cl D_xl_m_F...
            D_xl_m_H2O D_fl_m_Cl D_fl_m_F H2O_Sat];

        %Call the CreateTargets function which will run the model once for
        %the parameters used and will save the output
        CreateTargets(mpath,DataFolder,SaveFolder, Dataset,Model,...
            InputParam,X,D_fl_m_PrimBoilH2O,K_speciation,dlm, T,TF_File)

    case 'Run multistart to a target'
        %==================================================================
        %Set the multistart model properties
        %==================================================================
        num_frac = 100; %Number of fractionation steps for the model
        numStart=  10; %Number of start points

        %Whether to use parallel processing or not
        UseParallelProcessing = true;
        %UseParallelProcessing = false;

        %Numerical tolerance for RMSE
        Tolerance = 1e-2;

        RunMultistartModel(mpath, dlm, SaveFolder, Dataset, Model,...
            num_frac, Tolerance, numStart, UseParallelProcessing,T)

    case 'Visualize multistart outputs'
        %==================================================================
        %Set the visualization properties
        %==================================================================
        %Name of the file to visualize that is within the dataset output
        %folder
        %Fill follow the format of Multiouput_DatafileName_ModelName
        FileName = 'Multioutput_Test_Model1';

        %Run the plotting script
        MultistartFigures(mpath, dlm, SaveFolder, DataFolder, Dataset,...
            FileName, Model)

end
%==========================================================================
%==========================================================================
%vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
%Below are the functions that perform the main numerical operations:
%Explore model and create targets and run the multistart.
%vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
%==========================================================================
%==========================================================================


%This function explores a single model run for the set inputs
function CreateTargets(mpath,DataFolder,SaveFolder, Dataset,Model,...
    InputParam,X,D_fl_m_PrimBoilH2O,K_speciation,dlm, T,TF_File)
%Set the flag for a single run using parameters passed through
OperatorFlag = 'SingleRun';

[Output, XH2Ot_m, XCl_m, XF_m, XOH_m, F,NaCl_agg_wt] = ...
    Apatite_Crystallization_Multistart(InputParam,X, D_fl_m_PrimBoilH2O,...
    [],[],[],[],K_speciation,OperatorFlag, T);

if isnan(NaCl_agg_wt)
    NaCl_agg_wt = 0
else
    NaCl_agg_wt
end

XClap_out = Output(:,1);
XOHap_out = Output(:,2);
XFap_out = Output(:,3);
K_Cl_OH_out = Output(:,4);
K_F_OH_out = Output(:,5);
K_Cl_F_out = Output(:,6);

XCl_div_XOH = Output(:,7);
XF_div_XOH = Output(:,8);
XF_div_XCl = Output(:,9);

Fluid_inc_mass= Output(:,10);
            XNaClEq_fluid= Output(:,11);
            XH2O_Fluid= Output(:,12);
            NaClagg = Output(:,13);

%=========================
%Save the model output as a target file to run
%=========================
WriteData(F,XF_div_XOH,XCl_div_XOH,XF_div_XCl,XFap_out, XClap_out, XOHap_out,...
    XH2Ot_m, XCl_m, XF_m, XOH_m, InputParam, D_fl_m_PrimBoilH2O,...
    K_Cl_OH_out, K_F_OH_out, K_Cl_F_out, K_speciation,...
    mpath, SaveFolder, Dataset, Model, dlm,T,TF_File,Fluid_inc_mass,XNaClEq_fluid,XH2O_Fluid,NaClagg)

%Read the data from the .csv file saved
[XCl_div_XOH_Data,XF_div_XOH_Data,XF_div_XCl_Data]=ReadGeoData(mpath,...
    DataFolder,Dataset,dlm);



%=============================================
%Plot results from the single run
%figure 1 plots the T-dependent KDs versus the melt fraction
%figure 2 plots the data and the target trend in the XCl/XOH versus XF/XOH and XCl/XOH versus XF/XCl 
%=============================================
figure(1)
hold on
grid on
box on
plot(F,log10(K_Cl_OH_out),'r-','linewidth',2);
plot(F,log10(K_F_OH_out),'g-','linewidth',2);
plot(F,log10(K_Cl_F_out),'b-','linewidth',2);

%Create a legend
h(1) = plot(NaN,NaN,'linestyle','-','color','r','linewidth',2);
h(2) = plot(NaN,NaN,'linestyle','-','color','g','linewidth',2);
h(3) = plot(NaN,NaN,'linestyle','-','color','b','linewidth',2);
% lbl = {'Data';'Model'};
lbl = {'K_{Cl-OH}'; 'K_{F-OH}'; 'K_{Cl-F}'};
%Set positions of the Legend
legend(h,lbl,'Location','northeast','FontName','Times New Roman', ...
    'FontSize',10);
xlabel('Melt Fraction');
ylabel('log_{10}Kd');


figure(2)
subplot(1,2,1)
hold on
grid on
box on
plot(XF_div_XOH_Data,XCl_div_XOH_Data,'ko');
plot(XF_div_XOH,XCl_div_XOH,'b-','linewidth',2);

xlabel('XF/XOH');
ylabel('XCl/XOH');

subplot(1,2,2)
hold on
grid on
box on
plot(XF_div_XCl_Data,XCl_div_XOH_Data,'ko');
plot(XF_div_XCl,XCl_div_XOH,'b-','linewidth',2);

%Create a legend

h(1) = plot(NaN,NaN,'linestyle','-','color','b','linewidth',2);
h(2) = plot(NaN,NaN,'ko');
lbl = {'Data';'Model'};
% lbl = {'Model'};
%Set positions of the Legend
legend(h,lbl,'Location','northeast','FontName','Times New Roman', ...
    'FontSize',10);


xlabel('XF/XCl');
ylabel('XCl/XOH');


%This function does a parameter sweep comparing the model results against a
%target curve. The idea is to investigate what parameters may produce
%similar results

function RunMultistartModel(mpath, dlm, SaveFolder, Dataset, Model1,...
    num_frac, Tolerance, numStart, UseParallelProcessing,T)

     %Determine if reading a FT datafile or using a hardwired temperature
    %Determine the temperature vector based on the input

    %If the F-T file cannot be found, set a fixed temperature
        if istable(T) == 0
            T =T+273.15; %T in Kelvins
            F_start_Multi = 1; %Fraction of melt starting
            F_end_Multi = 0.1; %Fraction of melt ending

        %Else read the file and get the start and final melt fraction
        elseif istable(T) ==1
            %T = readtable(char([mFolder , Model, '.csv']));
            [C,ia,ic] = unique(T.Var2, 'stable');
            T = T(ia,:);

            F_start_Multi = max(T.Var2); %Fraction of melt starting
            F_end_Multi = max([0.1, min(T.Var2)]); %Fraction of melt ending
        end


OperatorFlag = 'Multistart';
%==========================================================================
%Get the data from the target file
%==========================================================================
%Where are the output files going to be saved?
SaveDirectory = [mpath SaveFolder dlm];

%Write the model outputs to a .csv file
Model_savePath = [SaveDirectory Dataset '_'  Model1 '_Tgt' '.csv'];

%Read the target curve points
%Read in the data as a table (Mixed number and strings)
Target_data = readtable(Model_savePath);

%Store the volatile ratios
F_OBS = Target_data.Frac_Melt;
XCl_div_XOH_TGT = Target_data.XCl_div_XOH;
XF_div_XOH_TGT = Target_data.XF_div_XOH;
XF_div_XCl_TGT = Target_data.XF_div_XCl;

%==========================================================================
%Get the constant parameters from the log file
%==========================================================================

log_savePath = [SaveDirectory Dataset '_' Model1 '_log' '.csv'];
Parameter_data = readtable(log_savePath);
%Define the values that will remain constant during the multistart
D_fl_m_PrimBoilH2O = Parameter_data.D_fl_m_PrimBoilH2O;

K_speciation = Parameter_data.K_speciation;
F_start = Parameter_data.Frac_0;
F_end = Parameter_data.Frac_f;

%==========================================================================
%Define the variable parameters for the crystallization model
%==========================================================================
%Initial melt compositions

H2O_init_low = 0.1;
H2O_init_high = 2;

Cl_init_low= 0.01;
Cl_init_high= 1;

F_init_low= 0.01;
F_init_high= 1;

%Initial volatile=melt distribution coefficients
D_xl_m_Cl_low = 0.01;
D_xl_m_Cl_high = 0.5;
D_xl_m_F_low = 0.01;
D_xl_m_F_high = 0.5;
D_xl_m_H2O_low = 0.01;
D_xl_m_H2O_high = 0.75;

%Initial volatile=fluid distribution coefficients
D_fl_m_Cl_low = 35;
D_fl_m_Cl_high =45;
D_fl_m_F_low = 1;
D_fl_m_F_high = 4;

H2O_Sat_low = 3;
H2O_Sat_high = 6;

%==========================================================================
%Create bound vectors for the variable parameters
%==========================================================================

%Place variables that vary as inputs as a coefficient vector that is then
%unpacked, then constants afterwards.
lb = [H2O_init_low, Cl_init_low, F_init_low, D_xl_m_Cl_low,...
    D_xl_m_F_low, D_xl_m_H2O_low, D_fl_m_Cl_low, D_fl_m_F_low,H2O_Sat_low];

ub = [H2O_init_high, Cl_init_high, F_init_high, D_xl_m_Cl_high,...
    D_xl_m_F_high, D_xl_m_H2O_high, D_fl_m_Cl_high, D_fl_m_F_high,H2O_Sat_high];


%Initial random starting guess
% InputCoeff = lb + rand(size(lb)).*(ub - lb);

InputCoeff = [Parameter_data.H2O_init, Parameter_data.Cl_init,...
    Parameter_data.F_init, Parameter_data.D_xl_m_Cl, ...
    Parameter_data.D_xl_m_F, Parameter_data.D_xl_m_H2O, ...
    Parameter_data.D_fl_m_Cl, Parameter_data.D_fl_m_F,Parameter_data.H2O_Sat];


%==========================================================================
%Targets to minimize
%==========================================================================
X = [F_start, F_end, num_frac];

%define the function to be minimized
FunAp = @(VariableCoeff)Apatite_Crystallization_Multistart(VariableCoeff,...
    X, D_fl_m_PrimBoilH2O,XF_div_XOH_TGT,XCl_div_XOH_TGT,...
    XF_div_XCl_TGT,F_OBS,K_speciation,OperatorFlag, T);

%Use the 'fmincon' local solver
opts = optimoptions('fmincon','MaxIterations',2000,'FunctionTolerance', ...
    Tolerance,'StepTolerance',Tolerance,'display','off');

%Setup the multistart problem
ms = MultiStart('UseParallel',UseParallelProcessing,'Display','iter', ...
    'StartPointsToRun','bounds','XTolerance',0);

problem = createOptimProblem('fmincon','x0',InputCoeff,'objective',...
    FunAp, 'lb',lb,'ub',ub,'options',opts);

%run the minimizer
%40 seconds for 100 iterations with no parallel
tic
[BestFitCoeff,fval,eflag,output,manymins] = run(ms,problem,numStart);
toc

%Define the file to save the multistart output
filename = char([SaveDirectory dlm 'Multioutput_' Dataset '_'  Model1 '.mat']);
[OutCoeff, Xinit,Fval] = Allbestfit(manymins);
numBestFit = size(OutCoeff,1);

XCl_div_XOH_out = cell(numBestFit,1);
XF_div_XOH_out = cell(numBestFit,1);
XF_div_XCl_out = cell(numBestFit,1);
F_out = cell(numBestFit,1);

%Run the model for the best fit conditions

FunAp = @(VariableCoeff)Apatite_Crystallization_Multistart(VariableCoeff,...
    X, D_fl_m_PrimBoilH2O,XF_div_XOH_TGT,XCl_div_XOH_TGT,...
    XF_div_XCl_TGT,F_OBS,K_speciation,'DoesntMatter', T);

for W = 1:1:numBestFit

%     Output = [XClap_out, XOHap_out,XFap_out,...
%     K_Cl_OH_out, K_F_OH_out, K_Cl_F_out, ...
%     XCl_div_XOH_out, XF_div_XOH_out, XF_div_XCl_out]

    [Output_out, XCl_div_XOH, XF_div_XOH, XF_div_XCl, F] = FunAp(OutCoeff(W,:));
%     [~, XCl_div_XOH, XF_div_XOH, XF_div_XCl, F] = FunAp(OutCoeff(W,:));

    
    XCl_div_XOH = Output_out(:,7);
    XF_div_XOH = Output_out(:,8);
    XF_div_XCl = Output_out(:,9);

    disp(XCl_div_XOH);
   

    XCl_div_XOH_out{W} = XCl_div_XOH;
    XF_div_XOH_out{W} = XF_div_XOH;
    XF_div_XCl_out{W} = XF_div_XCl;
    F_out{W} = F;
end


%Save the best fit runs (that are within the tolerances). The initial
%guess, which is the target parameters, is removed from this output.
OutputVariable = {XCl_div_XOH_out(2:end); XF_div_XOH_out(2:end); ...
    XF_div_XCl_out(2:end); F_out(2:end); ...
    OutCoeff(2:end,:); Fval(2:end)};
%Save all minimizations then load the coefficients
save(filename,'OutputVariable');

A = struct2cell(load(filename));
A = A{1,1};
Tbl = [A{5}, A{6}];
Variables = {'H2O','Cl','F','Dxtl-mCl','Dxtl-mF','Dxtl-mH2O','Dfl-mCl','Dfl-mF','H2OSat','RMSE'};
Tbl = array2table(Tbl);
Tbl.Properties.VariableNames = Variables;
filename2 = erase(filename,'.mat');
filename2 = [filename2 '.csv'];
writetable(Tbl,filename2);

% Uncomment to display the best fit values in the command window
% disp('H2O_init');
% disp(BestFitCoeff(1));
% disp('Cl_init');
% disp(BestFitCoeff(2));
% disp('F_init');
% disp(BestFitCoeff(3));
% disp('D_xl_m_Cl');
% disp(BestFitCoeff(4));
% disp('D_xl_m_F');
% disp(BestFitCoeff(5));
% disp('D_xl_m_H2O');
% disp(BestFitCoeff(6));
% disp('D_fl_m_Cl');
% disp(BestFitCoeff(7));
% disp('D_fl_m_F');
% disp(BestFitCoeff(8));
% 
% 
% disp(numBestFit);

%==========================================================================
%==========================================================================
%vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
%Below are the functions that are used to perform the model calculations,
%output results...etc
%vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
%==========================================================================
%==========================================================================


%==========================================================================
%Apatite Model Function
%==========================================================================

function [Output, XH2Ot_m, XCl_m, XF_m, XOH_m, F,NaCl_agg_wt] = ...
    Apatite_Crystallization_Multistart(VariableCoeff, X,...
    D_fl_m_PrimBoilH2O,XF_div_XOH_OBS,XCl_div_XOH_OBS,...
    XF_div_XCl_OBS,F_OBS,K_speciation, OperatorFlag, T)

NaCl_agg_wt=0;

%Extract the individual variables
H2O_init = VariableCoeff(1);
Cl_init = VariableCoeff(2);
F_init = VariableCoeff(3);
D_xl_m_Cl = VariableCoeff(4);
D_xl_m_F = VariableCoeff(5);
D_xl_m_H2O = VariableCoeff(6);
D_fl_m_Cl = VariableCoeff(7);
D_fl_m_F = VariableCoeff(8);
H2O_Sat = VariableCoeff(9);

F_start = X(1);
F_end = X(2);
num_frac = X(3);

%Declare the vectors for melt fraction
F = linspace(F_start,F_end,num_frac);
dF_vec = diff(F);
dF = abs(dF_vec(1));

Melt_H2Ot_frac = zeros(length(F),1);
Melt_Cl_frac= zeros(length(F),1);
Melt_F_frac = zeros(length(F),1);

Melt_H2O_Remaining = zeros(length(F),1);
Melt_Cl_Remaining= zeros(length(F),1);
Melt_F_Remaining= zeros(length(F),1);

Fluid_H2O_TotMass = zeros(length(F),1);
Fluid_H2O_IncMass = zeros(length(F),1);

%Variables for salinity calculations
Fluid_inc_mass = zeros(length(F),1);
Mass_NaCl_fluid_ist = zeros(length(F),1);
Cl_fluid_ist_wt = zeros(length(F),1);
NaCleq_fluid_inst = zeros(length(F),1);
NaCleq_fluid_wt = zeros(length(F),1);

%Set the concentration of volatile species (Cl,F) in the fluid = 0
Fluid_Cl_frac = zeros(length(F),1);
Fluid_F_frac = zeros(length(F),1);

%!!!!!!!NEW ADDITIONS BELOW!!!!

%define the output vectors for the volatile mole fractions in apatite
XClap_out = zeros(length(F),1);
XOHap_out = zeros(length(F),1);
XFap_out = zeros(length(F),1);

%define the output vectors for the KD's for each volatile
K_Cl_OH_out = zeros(length(F),1);
K_F_OH_out = zeros(length(F),1);
K_Cl_F_out = zeros(length(F),1);


%Set counter for primary boiling modification. The counter is set to -1
%because primary boiling begins on the second iteration.
count = -1;


for i = 1:1:length(F)

    %Calculate the fraction of water in the melt
    Melt_H2Ot_frac(i,1) =(H2O_init*F(i)^(D_xl_m_H2O-1))*0.01;

    %If the water concentration in the melt is undersaturated
    if Melt_H2Ot_frac(i,1) <= (H2O_Sat*0.01)

        %Compute melt trace element concentrations via Raleigh
        %fractionation
        Melt_Cl_frac(i,1) =(Cl_init*F(i)^(D_xl_m_Cl-1))*0.01;
        Melt_F_frac(i,1) =(F_init*F(i)^(D_xl_m_F-1))*0.01;

        %Because there is no fluid exsolved, the remaining Cl and F will be
        %the calculated  Cl and F wt fractions
        Melt_H2O_Remaining(i,1) = Melt_H2Ot_frac(i,1) - ...
            Fluid_H2O_TotMass(i,1);
        Melt_Cl_Remaining(i,1) = Melt_Cl_frac(i,1) - Fluid_Cl_frac(i,1);
        Melt_F_Remaining(i,1) =Melt_F_frac(i,1) - Fluid_F_frac(i,1);
    end


    %If the water concentration in the melt is oversaturated
    if Melt_H2Ot_frac(i,1) > (H2O_Sat*0.01)
        %Add a count for primary boiling
        count = count+1;


        %Calculate the concentration of Cl and in the melt
        Melt_Cl_frac(i,1) =(Melt_Cl_Remaining(i-1,1)*(F(i) / ...
            F(i-1))^(D_xl_m_Cl-1));

        Melt_F_frac(i,1) =(Melt_F_Remaining(i-1,1)*(F(i) / ...
            F(i-1))^(D_xl_m_F-1));


        %Calculate the mass of water exsolved from the melt
        Fluid_H2O_TotMass(i,1) = (H2O_init*F(i)^(D_xl_m_H2O-1))-H2O_Sat;
        %Calculate the incremental mass of water exsolved from each step
        %and convert into wt. fraction (/100)
        Fluid_H2O_IncMass(i,1) = (Fluid_H2O_TotMass(i,1) - ...
            Fluid_H2O_TotMass(i-1,1))/100;

        %Calculate the concentrations of Cl and F
        Fluid_Cl_mass = D_fl_m_Cl*Melt_Cl_frac(i,1);
        Fluid_F_mass = D_fl_m_F*Melt_F_frac(i,1);


        %Calculate the mass of F and Cl entering the fluid by taking the
        %concentration of F and Cl in the fluid and multiplying this by the
        %wt fraction of H2O liquid that has exsolved from the melt.
        Fluid_Cl_frac(i,1) = Fluid_Cl_mass* Fluid_H2O_IncMass(i,1);
        Fluid_F_frac(i,1) = Fluid_F_mass* Fluid_H2O_IncMass(i,1);

         %%% Salinity Calculations %%%
        % Incremental mass of the fluid
        Fluid_inc_mass(i,1) = Fluid_Cl_frac(i,1) + Fluid_H2O_IncMass(i,1);
        Mass_NaCl_fluid_ist(i,1) = Fluid_Cl_frac(i,1)*58.44/35.453;
        Cl_fluid_ist_wt(i,1) = Fluid_Cl_frac(i,1)/( Fluid_H2O_IncMass(i,1)+Fluid_Cl_frac(i,1));
        NaCleq_fluid_inst(i,1) = Cl_fluid_ist_wt(i,1)*(58.44/35.453);
        NaCleq_fluid_wt(i,1) = NaCleq_fluid_inst(i,1)*100;

        %Mole fraction
        XNaCleq_fluid(i,1) = (Mass_NaCl_fluid_ist(i,1)/58.44)/((Mass_NaCl_fluid_ist(i,1)/58.44)+((Fluid_H2O_IncMass(i,1)/18.015)));
        XH2O_fluid(i,1) = (Fluid_H2O_IncMass(i,1)/18.015)/((Mass_NaCl_fluid_ist(i,1)/58.44)+(Fluid_H2O_IncMass(i,1)/18.015));
        NaClagg(i,1) = NaCleq_fluid_inst(i,1)*Fluid_inc_mass(i,1);


        Melt_H2O_Remaining(i,1) = Melt_H2Ot_frac(i,1) - ...
            (Fluid_H2O_TotMass(i,1)/100)-(dF*D_fl_m_PrimBoilH2O);

        H2O_Sat = H2O_Sat - (100*dF*D_fl_m_PrimBoilH2O);

        Melt_Cl_Remaining(i,1) = Melt_Cl_frac(i,1) - Fluid_Cl_frac(i,1);
        Melt_F_Remaining(i,1) =Melt_F_frac(i,1) - Fluid_F_frac(i,1);
    end

end

% Salinity calculations to work out the aggregate amount of NaCleq in the
% fluid
TotalFluid_mass = sum(Fluid_inc_mass);
Sum_NaCl = sum(NaClagg);
NaCleq_Agg = Sum_NaCl/TotalFluid_mass;
NaCl_agg_wt = NaCleq_Agg*100;

%Get the mole fraction of the volatiles in the melt
XH2Ot_m = (Melt_H2O_Remaining/18.015)./ ...
    ((Melt_H2O_Remaining/18.015)+((1-Melt_H2O_Remaining)/32));
XCl_m =  (Melt_Cl_Remaining/35.5)./ ...
    ((Melt_H2O_Remaining/18.015)+((1-Melt_H2O_Remaining)/32));
XF_m =  (Melt_F_Remaining/19)./ ...
    ((Melt_H2O_Remaining/18.015)+((1-Melt_H2O_Remaining)/32));
XOH_m = equilibrium(XH2Ot_m,K_speciation);


%%%NEW!!!!%%====================
%Define the constants used in the non-ideal model
Wg_ClOH = 5;
Wg_FOH = 7;
Wg_FCl = 16;

R = 8.314; %gas constant in J/K/mol

% Difference between two Wgs (to be used for Kd calculation
%if non-ideal)
Wg_diff1 = Wg_FOH-Wg_FCl;
Wg_diff2 = Wg_ClOH-Wg_FCl;


%Loop through each melt composition (function of F and computed in the
%previous loop) and determine the KD's and the associated melt composition

%%NEW!!!!%%====================
%in this section we create a temperature vector
if istable(T) == 1
    %T_vec = interp1(F_data,T_data, F_query
    T_vec = interp1(T.Var2,T.Var1,F) + 273.15;
else
    %Create a vector of length F equal with all values equal to input
    T_vec = ones(length(F),1)*T;
end

%%STILL NEW!!!!%%====================
for j = 1:1:length(F)

    %Get the current molar fraction of volatiles in the melt
    XF_m_curr = XF_m(j);
    XCl_m_curr = XCl_m(j);
    XOH_m_curr = XOH_m(j);
    
    T_curr = T_vec(j);
    
    % Gibbs free energy of exchange reactions deltaG [kJ/mol]
    deltaG_ClOH = 72.9-0.034*T_curr;
    deltaG_FOH = 94.6-0.04*T_curr;

    %Generate anonymous functions to calculate the KD's
    K_OH_Cl = @(XClap_init, XOHap_init, XFap_init) exp((1000*(-deltaG_ClOH-((XClap_init-XOHap_init).*Wg_ClOH+XFap_init.*Wg_diff1)))./(R.*T_curr));
    K_OH_F = @(XClap_init, XOHap_init, XFap_init) exp(1000*(-deltaG_FOH-((XFap_init-XOHap_init).*Wg_FOH+XClap_init.*Wg_diff2))./(R.*T_curr));
    K_Cl_F = @(XClap_init, XOHap_init, XFap_init) K_OH_F(XClap_init, XOHap_init, XFap_init)./K_OH_Cl(XClap_init, XOHap_init, XFap_init);

    K_Cl_OH = @(XClap_init, XOHap_init, XFap_init) 1/K_OH_Cl(XClap_init, XOHap_init, XFap_init);
    K_F_OH = @(XClap_init, XOHap_init, XFap_init) 1/K_OH_F(XClap_init, XOHap_init, XFap_init);
    
    %linear inequality constraints for minimization
    A = [];
    b = [];

    %linear equality constraints:
    % 1*XClap_init + 1*XOHap_init + 1*XFap_init = 1

    %How this works: Aeq is [1,1,1] because the multiplier infront of the
    %of the molar fraction is 1
    %beq is 1 because the sums of all the mole fractions is = 1
    Aeq = [1,1,1];
    beq = 1;

    %Lower and upper bounds are between 0 and 1 for any given mole fraction
    lbA = [0,0,0];
    ubA = [1,1,1];

    %Create an minimization target function that uses our Li model function
    %defined lower (GetApMoleFractions). The reason we do this is because
    %the minimizer only accepts one variable (x). So we make an anonymous
    %function with variable x that calls GetApMoleFractions with the other
    %variables. This is how we can pass additional variables to a
    %minimization routine! a useful trick
    fun = @(x)GetApMoleFractions(x, XF_m_curr, XCl_m_curr, XOH_m_curr,...
        K_Cl_OH, K_F_OH, K_Cl_F);

    %Initial guess for [XClap, XOHap, XFap] (random, but if you have a good
    %guess can use that to potentially reduce iterations and speed up the
    %code)

    x0 = [0.7, 0.1, 0.2];


    %Run fmin constrained (fmincon) using:
    %upper and lower bound conditions (0-1)
    %Equality measurements (1*XClap_init + 1*XOHap_init + 1*XFap_init = 1)

    %Can set fmincon options here
    options = optimoptions('fmincon','Algorithm','active-set','Display', 'off', 'OptimalityTolerance', 1e-9, 'MaxIterations', 10000,...
        'StepTolerance', 1e-9, 'UseParallel', false) ;
    x_out = fmincon(fun,x0,A,b,Aeq,beq,lbA,ubA,[],options);

    XClap_curr = x_out(1);
    XOHap_curr = x_out(2);
    XFap_curr = x_out(3);

    %Subsititue the calculated molar fraction of volatiles in apatite back
    %into the KD in order to get the KD associated with this composition
    K_Cl_OH_curr = K_Cl_OH(XClap_curr, XOHap_curr, XFap_curr);
    K_F_OH_curr = K_F_OH(XClap_curr, XOHap_curr, XFap_curr);
    K_Cl_F_curr = K_Cl_F(XClap_curr, XOHap_curr, XFap_curr);

    %Store the calculated values to associated vectors for output
    XClap_out(j,1) = XClap_curr;
    XOHap_out(j,1) = XOHap_curr;
    XFap_out(j,1) = XFap_curr;

    %Store the associated KD's to vectors for output
    K_Cl_OH_out(j,1) = K_Cl_OH_curr;
    K_F_OH_out(j,1) = K_F_OH_curr;
    K_Cl_F_out(j,1) = K_Cl_F_curr;

    
end


%Compute the component ratios
XCl_div_XOH_out = XClap_out./XOHap_out;
XF_div_XOH_out = XFap_out./XOHap_out;
XF_div_XCl_out = XFap_out./XClap_out;

Output = [XClap_out, XOHap_out,XFap_out,...
    K_Cl_OH_out, K_F_OH_out, K_Cl_F_out, ...
    XCl_div_XOH_out, XF_div_XOH_out, XF_div_XCl_out,Fluid_inc_mass,XNaCleq_fluid,XH2O_fluid,NaClagg];

switch OperatorFlag
    case 'Multistart'

        %Looking at residuals between the Observations and the model
        XCl_div_XOH_interp = interp1(F,XCl_div_XOH_out,F_OBS);
        XF_div_XOH_interp = interp1(F,XF_div_XOH_out,F_OBS);
        XF_div_XCl_interp = interp1(F,XF_div_XCl_out,F_OBS);


        Res_XCl_div_XOH_interp = abs(XCl_div_XOH_OBS - XCl_div_XOH_interp);
        Res_XF_div_XOH_interp= abs(XF_div_XOH_OBS - XF_div_XOH_interp);
        %Res_XF_div_XCL_interp = abs(XF_div_XCL_OBS - XF_div_XCL_interp);

        %Res_all = [Res_XCl_div_XOH_interp ; Res_XF_div_XOH_interp ;...
        %Res_XF_div_XCL_interp];

        Res_all = [Res_XCl_div_XOH_interp; Res_XF_div_XOH_interp];
        Output = sqrt(mean(Res_all.^2));

    case 'SingleRun'
%         Output = [];
end


%==========================================================================
%Function which computes the equilibrium OH speciation
%==========================================================================

function XOH_m = equilibrium(XH2Ot,K_speciation)
XOH_m=(0.5-(0.25-(K_speciation-4).*(XH2Ot-XH2Ot.^2)/K_speciation).^(1/2))./...
    ((K_speciation-4)./(2*K_speciation));

function out = GetApMoleFractions(Coeff, XF_m_out, XCl_m_out,...
    XOH_m_out, Kd_ClOH_fun, Kd_FOH_fun, Kd_ClF_fun)

XClap_init = Coeff(1);
XOHap_init = Coeff(2);
XFap_init = Coeff(3);

%These functions are to be a zero, we are aiming to find the values of
%XFap_init, XClap_init, XOHap_init where x, y and z are equal to zero.
%The unknowns in the bottom equation are all in terms of those 3 variables


x = XFap_init - ( 1./(1+((XCl_m_out.*Kd_ClF_fun(XClap_init, XOHap_init, XFap_init))./XF_m_out) +...
    (XOH_m_out./(XF_m_out.*Kd_FOH_fun(XClap_init, XOHap_init, XFap_init)))) );

y = XClap_init - ( 1./(1+(XF_m_out./(XCl_m_out.*Kd_ClF_fun(XClap_init, XOHap_init, XFap_init))) +...
    (XOH_m_out./(XCl_m_out.*Kd_ClOH_fun(XClap_init, XOHap_init, XFap_init)))) );

z = XOHap_init - ( 1./(1+((XF_m_out.*Kd_FOH_fun(XClap_init, XOHap_init, XFap_init))./XOH_m_out) +...
    ((XCl_m_out.*Kd_ClOH_fun(XClap_init, XOHap_init, XFap_init))./XOH_m_out)) );

%The output is the sum of each of the above equations, use the absolute
%value incase x,y or z are negative. Fmincon wants to make this value as
%small as possible, or as close to 0.
out = abs(x) + abs(y) + abs(z);

%==========================================================================
%Function to extract all best fit coefficients
%==========================================================================

function [OutCoeff, Xinit, Fval] = Allbestfit(AllDat)
numSol = length(AllDat);

OutCoeff = NaN(numSol,9);
Xinit = NaN(numSol,9);
Fval = NaN(numSol,1);

for H = 1:1:numSol
    OutCoeff(H,:) = AllDat(H).X;
    Xinit(H,:) = cell2mat(AllDat(H).X0)';
    Fval(H,1) = AllDat(H).Fval;
end
%==========================================================================
%Function to write out the datafiles
%==========================================================================

function WriteData(F,XF_div_XOH,XCl_div_XOH,XF_div_XCl,XFap, XClap, XOHap,...
    XH2Ot_m, XCl_m, XF_m, XOH_m, InputParam, D_fl_m_PrimBoilH2O, ...
    K_Cl_OH_out, K_F_OH_out, K_Cl_F_out, K_speciation,...
    mpath, SaveFolder, Dataset, Model1, dlm,T,TF_File,Fluid_inc_mass,XNaClEq_fluid,XH2O_Fluid,NaClagg)

%Where are the output files going to be saved
SaveDirectory = [mpath SaveFolder dlm];

%Check if the directory exists, if not make it
if ~exist(SaveDirectory, 'dir')
    mkdir(SaveDirectory)
end

%Write the model outputs to a .csv file
Model_savePath = [SaveDirectory Dataset '_'  Model1 '_Tgt' '.csv'];
Frac_Melt = F';

T_output = table(Frac_Melt,K_F_OH_out,K_Cl_OH_out,K_Cl_F_out,...
    XF_div_XOH, XCl_div_XOH, XF_div_XCl,...
    XH2Ot_m, XCl_m, XF_m, XOH_m,...
    XClap, XFap, XOHap,Fluid_inc_mass,XNaClEq_fluid,XH2O_Fluid,NaClagg);

writetable(T_output, Model_savePath);

%Get the initial model parameters
%===
H2O_init = InputParam(1);Cl_init = InputParam(2); F_init = InputParam(3);
D_xl_m_Cl = InputParam(4); D_xl_m_F = InputParam(5); D_xl_m_H2O = InputParam(6);
D_fl_m_Cl = InputParam(7);D_fl_m_F = InputParam(8);H2O_Sat = InputParam(9);
%====
%Create a table with the initial parameters
Frac_0 = F(1);
Frac_f = F(end);

Directory = string(mpath);
OutputFile = string(SaveFolder);
TargetCurve = string(Model_savePath);
Data = string(Dataset);
Model = string(Model1);
if istable(T)
    T_val = string([mpath dlm TF_File]);
else
    T_val = T-273.15;
end
T_log = table(Frac_0, Frac_f, Cl_init, F_init, H2O_init,  D_xl_m_Cl, ...
    D_xl_m_F, D_xl_m_H2O, D_fl_m_Cl, D_fl_m_F, ...
    D_fl_m_PrimBoilH2O, H2O_Sat, K_speciation,Directory,OutputFile,TargetCurve,T_val,Model,Data);

%Save the initial parameters used for the model to a .csv file
log_savePath = [SaveDirectory Dataset '_' Model1 '_log' '.csv'];
writetable(T_log, log_savePath);

%==========================================================================
%Plotting
%==========================================================================

%==========================================================================
%Function to get the dataset of interest
%==========================================================================
function [XCl_div_XOH_Data,XF_div_XOH_Data,XF_div_XCL_Data] = ...
    ReadGeoData(mpath,DataFolder,Dataset,dlm)

%Define the path of the file to read in
CSVdataFilePath = [mpath DataFolder dlm Dataset '.csv'];

%Read in the data as a table (Mixed number and strings)
data = readtable(CSVdataFilePath);

%Store the volatile ratios
XCl_div_XOH_Data = data.Cl_OH;
XF_div_XOH_Data = data.F_OH;
XF_div_XCL_Data = data.F_Cl;


function MultistartFigures(mpath, dlm, SaveFolder, DataFolder, Dataset,...
    FileName,Model)

%Read the data from the .csv file saved
[XCl_div_XOH_Data,XF_div_XOH_Data,XF_div_XCl_Data]=ReadGeoData(mpath,...
    DataFolder,Dataset,dlm);


%Read in the multistart data
SaveDirectory = [mpath SaveFolder dlm];
filename = char([SaveDirectory dlm FileName '.mat']);


%Write the model outputs to a .csv file
Model_savePath = [SaveDirectory Dataset '_'  Model '_Tgt' '.csv'];
Target_data = readtable(Model_savePath);

%Store the volatile ratios
F_OBS = Target_data.Frac_Melt;
XCl_div_XOH_TGT = Target_data.XCl_div_XOH;
XF_div_XOH_TGT = Target_data.XF_div_XOH;
XF_div_XCl_TGT = Target_data.XF_div_XCl;


%Load the data into a cell matrix
ModelOutputs = struct2cell(load(filename));
ModelOutputs = ModelOutputs{1,1};

%Extract the individual parameters
XCl_div_XOH_out = ModelOutputs{1,1};
XF_div_XOH_out = ModelOutputs{2,1};
XF_div_XCl_out = ModelOutputs{3,1};
F_out = ModelOutputs{4,1};
OutCoeff = ModelOutputs{5,1};

%Define the markersize
sz = 2;

%How many fits are there?
numBestFit = size(XCl_div_XOH_out,1);

%WHat is the RMSE Value
Fval = ModelOutputs{6,1};

%Sort the RMSE
[~,srtIdx] = sort(log10(Fval),'descend');

%Apply the sorting index to the data so it is arranged based on RMSE
XF_div_XOH_out = XF_div_XOH_out(srtIdx,1);
XCl_div_XOH_out = XCl_div_XOH_out(srtIdx,1);
XF_div_XCl_out = XF_div_XCl_out(srtIdx,1);
Fval = Fval(srtIdx);


%========================================================================
%Figure 1: Bivariate plots of geochemical data
%========================================================================
figure(1)
%Create the tight subplot
% [ha, pos] = tight_subplot(Nh, Nw, gap, marg_h, marg_w)
ha = tight_subplot(1,2,[0.075 0.1],[.125 .025],[.1 .025]);
hold on
Colour = jet;
ClVal = linspace(min(log10(Fval)),max(log10(Fval)),length(Colour))';
%------------Panel 1------------

ha_Targ = ha(1);
axes(ha_Targ);
for W = 1:1:numBestFit
    Clr = [interp1(ClVal,Colour(:,1),log10(Fval(W)))...
        interp1(ClVal,Colour(:,2),log10(Fval(W)))...
        interp1(ClVal,Colour(:,3),log10(Fval(W)))];

    hold on
    plot(XF_div_XOH_out{W,1},XCl_div_XOH_out{W,1},'-','color',...
        Clr,'linewidth',1);
end



plot(XF_div_XOH_TGT,XCl_div_XOH_TGT,'k-','linewidth',2);
plot(XF_div_XOH_Data,XCl_div_XOH_Data,'ko','markersize',5);

box on
grid on
a = get(ha_Targ,'XTick');
set(ha_Targ,'XTickLabel','')
set(ha_Targ,'XTick',a,'FontName','Arial','fontsize',12)
set(ha_Targ, 'linewidth',1);
set(ha_Targ,'TickLength',[0.025, 0.025])
set(ha_Targ,'XMinorTick','on','YMinorTick','on')

xlabel('X_F/X_O_H','FontName','Arial','FontSize',12);
ylabel('X_C_l/X_O_H','FontName','Arial','FontSize',12);

Y_limits = [0 1.5];
ylim(Y_limits);
X_limits = [0 20];
xlim(X_limits)

set(ha_Targ, 'xtick', (min(X_limits):2.5:max(X_limits)));

xticklabels('auto');
yticklabels('auto');

%------------Panel 2------------
ha_Targ = ha(2);
axes(ha_Targ);
for W = 1:1:numBestFit
    Clr = [interp1(ClVal,Colour(:,1),log10(Fval(W)))...
        interp1(ClVal,Colour(:,2),log10(Fval(W)))...
        interp1(ClVal,Colour(:,3),log10(Fval(W)))];
    hold on
    plot(XF_div_XCl_out{W,1},XCl_div_XOH_out{W,1},'-','color', ...
        Clr,'linewidth',1);
end

hold on
plot(XF_div_XCl_TGT,XCl_div_XOH_TGT,'k-','linewidth',2);
plot(XF_div_XCl_Data,XCl_div_XOH_Data,'ko','markersize',5);

box on
grid on
a = get(ha_Targ,'XTick');
set(ha_Targ,'XTickLabel','')
set(ha_Targ,'XTick',a,'FontName','Arial','fontsize',12)
set(ha_Targ, 'linewidth',1);
set(ha_Targ,'TickLength',[0.025, 0.025])
set(ha_Targ,'XMinorTick','on','YMinorTick','on')

xlabel('X_F/X_C_l','FontName','Arial','FontSize',12);
ylabel('X_C_l/X_O_H','FontName','Arial','FontSize',12);

Y_limits = [0 1.5];
ylim(Y_limits);
X_limits = [0 50];
xlim(X_limits)
set(ha_Targ, 'xtick', (min(X_limits):10:max(X_limits)));

xticklabels('auto');
yticklabels('auto');


%Format the colorbar
hcb2 = colorbar;
colormap(ha_Targ,Colour)
caxis([ClVal(1),ClVal(end)]);
colorTitleHandle = get(hcb2,'Title');
titleString = 'log RMSE';
set(colorTitleHandle ,'String',titleString,'FontName', ...
    'Arial','FontSize',12);


%Set positions of the colorbar
Axespos = get(ha_Targ,'Position');
Cbarpos = get(hcb2,'Position');
set(hcb2,'position', ...
    [Axespos(1)+0.275 Axespos(2)+0.325 Cbarpos(3) 0.5*Cbarpos(4)])

%Set the size of the figure
set(gcf, 'Units', 'Inches', 'Position', [6, 3, 6.5, (6.5/(6/4))],...
    'PaperUnits', 'Inches', 'PaperSize', [8.5, 11]);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%==========================================================================
%Figure 2: Histogram plots of model runs
%==========================================================================
figure(2)
%Create the tight subplot
% [ha, pos] = tight_subplot(Nh, Nw, gap, marg_h, marg_w)
ha = tight_subplot(2,5,[0.1 0.05],[.125 .025],[.05 .025]);


LabelFontSize = 10;
TickFontSize = 8;
LabelVector = {'H_2O wt. %', 'Cl wt. %', 'F wt. %', 'D_x_t_l_-_m Cl', ...
    'D_x_t_l_-_m F', 'D_x_t_l_-_m H_2O', 'D_f_l_-_m Cl', ...
    'D_f_l_-_m F', 'H_2O Sat', 'log RMSE'};

for i = 1:1:10
    hold on
    box on
    ha_Targ = ha(i);
    axes(ha_Targ);
    if i ~= 10
        histogram(OutCoeff(:,i),'BinMethod','fd');
    else
        histogram(log10(Fval),'BinMethod','fd')
    end
    a = get(ha_Targ,'XTick');
    set(ha_Targ,'XTickLabel','')
    set(ha_Targ,'XTick',a,'FontName','Arial','fontsize',TickFontSize);
    set(ha_Targ, 'linewidth',1);
    set(ha_Targ,'TickLength',[0.025, 0.025])
    xlabel(LabelVector{i},'FontName','Arial','FontSize',LabelFontSize);
    
    xticklabels('auto');
    yticklabels('auto');
end


set(gcf, 'Units', 'Inches', 'Position', [6, 3, 6.5, (6.5/(6/4))],...
    'PaperUnits', 'Inches', 'PaperSize', [8.5, 11]);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%==========================================================================
%Figure 3: RMSE scatter diagrams
%==========================================================================
figure(3)
%Create the tight subplot
% [ha, pos] = tight_subplot(Nh, Nw, gap, marg_h, marg_w)
%Labels on the Y axis
%Labels on the X axis
ha = tight_subplot(3,3,[0.1 0.05],[.125 .025],[.05 .025]);

ha_Targ = ha(1);
axes(ha_Targ);
%plot(log10(Fval),OutCoeff(srtIdx,1),'bo','markersize',sz);
plot(OutCoeff(srtIdx,1),log10(Fval),'bo','markersize',sz);
a = get(ha_Targ,'XTick');
set(ha_Targ,'XTickLabel','')

set(ha_Targ,'XTick',a,'FontName','Arial','fontsize',TickFontSize);
set(ha_Targ, 'linewidth',1);
set(ha_Targ,'TickLength',[0.025, 0.025])

ylim([floor(min(log10(Fval))) 0]);
xlabel('H_2O wt. %','FontName','Arial','FontSize',LabelFontSize);
xticklabels('auto');
yticklabels('auto');


ha_Targ = ha(2);
axes(ha_Targ);
% plot(log10(Fval),OutCoeff(srtIdx,2),'bo','markersize',sz);
plot(OutCoeff(srtIdx,2),log10(Fval),'bo','markersize',sz);
a = get(ha_Targ,'XTick');
set(ha_Targ,'XTickLabel','')
set(ha_Targ,'XTick',a,'FontName','Arial','fontsize',TickFontSize);
set(ha_Targ, 'linewidth',1);
set(ha_Targ,'TickLength',[0.025, 0.025])

ylim([floor(min(log10(Fval))) 0]);
xlabel('Cl wt. %','FontName','Arial','FontSize',LabelFontSize);
xticklabels('auto');
yticklabels('auto');

ha_Targ = ha(3);
axes(ha_Targ);
%plot(log10(Fval),OutCoeff(srtIdx,3),'bo','markersize',sz);
plot(OutCoeff(srtIdx,3),log10(Fval),'bo','markersize',sz);
a = get(ha_Targ,'XTick');
set(ha_Targ,'XTickLabel','')
set(ha_Targ,'XTick',a,'FontName','Arial','fontsize',TickFontSize);
set(ha_Targ, 'linewidth',1);
set(ha_Targ,'TickLength',[0.025, 0.025])
ylim([floor(min(log10(Fval))) 0]);
xlabel('F wt. %','FontName','Arial','FontSize',LabelFontSize);
xticklabels('auto');
yticklabels('auto');

ha_Targ = ha(4);
axes(ha_Targ);
%plot(log10(Fval),OutCoeff(srtIdx,4),'bo','markersize',sz);
plot(OutCoeff(srtIdx,4),log10(Fval),'bo','markersize',sz);
a = get(ha_Targ,'XTick');
set(ha_Targ,'XTickLabel','')
set(ha_Targ,'XTick',a,'FontName','Arial','fontsize',TickFontSize);
set(ha_Targ, 'linewidth',1);
set(ha_Targ,'TickLength',[0.025, 0.025])
ylim([floor(min(log10(Fval))) 0]);
xlabel('D_x_t_l_-_m Cl','FontName','Arial','FontSize',LabelFontSize);
xticklabels('auto');
yticklabels('auto');

ha_Targ = ha(5);
axes(ha_Targ);
%plot(log10(Fval),OutCoeff(srtIdx,5),'bo','markersize',sz);
plot(OutCoeff(srtIdx,5),log10(Fval),'bo','markersize',sz);
a = get(ha_Targ,'XTick');
set(ha_Targ,'XTickLabel','')
set(ha_Targ,'XTick',a,'FontName','Arial','fontsize',TickFontSize);
set(ha_Targ, 'linewidth',1);
set(ha_Targ,'TickLength',[0.025, 0.025])
ylim([floor(min(log10(Fval))) 0]);
xlabel('D_x_t_l_-_m F','FontName','Arial','FontSize',LabelFontSize);
xticklabels('auto');
yticklabels('auto');

ha_Targ = ha(6);
axes(ha_Targ);
%plot(log10(Fval),OutCoeff(srtIdx,6),'bo','markersize',sz);
plot(OutCoeff(srtIdx,6),log10(Fval),'bo','markersize',sz);
a = get(ha_Targ,'XTick');
set(ha_Targ,'XTickLabel','')
set(ha_Targ,'XTick',a,'FontName','Arial','fontsize',TickFontSize);
set(ha_Targ, 'linewidth',1);
set(ha_Targ,'TickLength',[0.025, 0.025])
ylim([floor(min(log10(Fval))) 0]);
xlabel('D_x_t_l_-_m H_2O','FontName','Arial','FontSize',LabelFontSize);
xticklabels('auto');
yticklabels('auto');

ha_Targ = ha(7);
axes(ha_Targ);
%plot(log10(Fval),OutCoeff(srtIdx,7),'bo','markersize',sz);
plot(OutCoeff(srtIdx,7),log10(Fval),'bo','markersize',sz);
a = get(ha_Targ,'XTick');
set(ha_Targ,'XTickLabel','')
set(ha_Targ,'XTick',a,'FontName','Arial','fontsize',TickFontSize);
set(ha_Targ, 'linewidth',1);
set(ha_Targ,'TickLength',[0.025, 0.025])
ylim([floor(min(log10(Fval))) 0]);
xlabel('D_f_l_-_m Cl','FontName','Arial','FontSize',LabelFontSize);
xticklabels('auto');
yticklabels('auto');

ha_Targ = ha(8);
axes(ha_Targ);
%plot(log10(Fval),OutCoeff(srtIdx,8),'bo','markersize',sz);
plot(OutCoeff(srtIdx,8),log10(Fval),'bo','markersize',sz);
a = get(ha_Targ,'XTick');
set(ha_Targ,'XTickLabel','')
set(ha_Targ,'XTick',a,'FontName','Arial','fontsize',TickFontSize);
set(ha_Targ, 'linewidth',1);
set(ha_Targ,'TickLength',[0.025, 0.025])
ylim([floor(min(log10(Fval))) 0]);
xlabel('D_f_l_-_m F','FontName','Arial','FontSize',LabelFontSize);
xticklabels('auto');
yticklabels('auto');

ha_Targ = ha(9);
axes(ha_Targ);
plot(OutCoeff(srtIdx,9),log10(Fval),'bo','markersize',sz);
a = get(ha_Targ,'XTick');
set(ha_Targ,'XTickLabel','')
set(ha_Targ,'XTick',a,'FontName','Arial','fontsize',TickFontSize);
set(ha_Targ, 'linewidth',1);
set(ha_Targ,'TickLength',[0.025, 0.025])
ylim([floor(min(log10(Fval))) 0]);
xlabel('H_2O Sat','FontName','Arial','FontSize',LabelFontSize);
xticklabels('auto');
yticklabels('auto');

%set(get(ha_Targ,'Children'),'Visible','off');

% subplot(3,4,12);
% BinSize = 2*((iqr(Fval))/(numBestFit^(1/3)));
% bins = round((max(Fval)-min(Fval))/BinSize);
% histogram(Fval,bins)
% xlabel('Fval (RMSE)');

set(gcf, 'Units', 'Inches', 'Position', [6, 3, 6.5, (6.5/(6/4))],...
    'PaperUnits', 'Inches', 'PaperSize', [8.5, 11]);


%Function to create tight subplots
function [ha, pos] = tight_subplot(Nh, Nw, gap, marg_h, marg_w)

% tight_subplot creates "subplot" axes with adjustable gaps and margins
%
% [ha, pos] = tight_subplot(Nh, Nw, gap, marg_h, marg_w)
%
%   in:  Nh      number of axes in height (vertical direction)
%        Nw      number of axes in width (horizontaldirection)
%        gap     gaps between the axes in normalized units (0...1)
%                   or [gap_h gap_w] for different gaps in height and width
%        marg_h  margins in height in normalized units (0...1)
%                   or [lower upper] for different lower and upper margins
%        marg_w  margins in width in normalized units (0...1)
%                   or [left right] for different left and right margins
%
%  out:  ha     array of handles of the axes objects
%                   starting from upper left corner, going row-wise as in
%                   subplot
%        pos    positions of the axes objects
%
%  Example: ha = tight_subplot(3,2,[.01 .03],[.1 .01],[.01 .01])
%           for ii = 1:6; axes(ha(ii)); plot(randn(10,ii)); end
%           set(ha(1:4),'XTickLabel',''); set(ha,'YTickLabel','')

% Pekka Kumpulainen 21.5.2012   @tut.fi
% Tampere University of Technology / Automation Science and Engineering


if nargin<3; gap = .02; end
if nargin<4 || isempty(marg_h); marg_h = .05; end
if nargin<5; marg_w = .05; end

if numel(gap)==1;
    gap = [gap gap];
end
if numel(marg_w)==1;
    marg_w = [marg_w marg_w];
end
if numel(marg_h)==1;
    marg_h = [marg_h marg_h];
end

axh = (1-sum(marg_h)-(Nh-1)*gap(1))/Nh;
axw = (1-sum(marg_w)-(Nw-1)*gap(2))/Nw;

py = 1-marg_h(2)-axh;

% ha = zeros(Nh*Nw,1);
ii = 0;
for ih = 1:Nh
    px = marg_w(1);

    for ix = 1:Nw
        ii = ii+1;
        ha(ii) = axes('Units','normalized', ...
            'Position',[px py axw axh], ...
            'XTickLabel','', ...
            'YTickLabel','');
        px = px+axw+gap(2);
    end
    py = py-axh-gap(1);
end
if nargout > 1
    pos = get(ha,'Position');
end
ha = ha(:);



