%% Script 4: Analyzing the Chew Bahir tuning ensemble
% 
% Trauth, M.H., Asrat, A., Fischer, M., Foerster, V., Kaboth-Bahr, S.,
% Lamb, H.F., Marwan, N., Roberts, H.M., Schaebitz, F., 2024, Combining
% orbital tuning and direct dating approaches to age-depth model
% development for Chew Bahir, Ethiopia. Quaternary Science Advances
% 
% Martin H Trauth 7 June 2024

% Clear workspace, clear command window and close figure windows.

    clear, clc, close all
    
% Create directories.

    mkdir DirFiguresEnsemble

% For loop to run the experiment several times for different base ages
% which are defined here:

    allbaseages = 450 : 10 : 650;

% Loop.

   for loopi = 1 : length(allbaseages)/3

% Define tuned age models to be used.

   agemodeltune    = allbaseages(1,1+3*(loopi-1):3*(loopi))

   [Lia,Locb] = ismember(agemodeltune,allbaseages);

% Create character string for all age models.
    agemodeltunestr = ['MHT',num2str(agemodeltune(1));
                       'MHT',num2str(agemodeltune(2));
                       'MHT',num2str(agemodeltune(3))];

% Importing age models and data.

    % Import all tuned age models.
    dirname = 'DirResultsTiePointsTuning';
     
    tiepointstune = NaN(50,2,26);

    for loopi = 1 : length(allbaseages)
        baseage = allbaseages(loopi);
        filename = [dirname,...
            '/agemodeltiepoints_mht',num2str(baseage),'.txt'];
        tiepointsraw = load(filename);
        tiepointstune(1:length(tiepointsraw),:,loopi) = tiepointsraw;
    end

    % Argon and tephrochronological ages.
    % 2B-035E-1,  17- 24 cm = 68.855 mcd mid-point
    % 2B-101Q-1,  11-120 cm = 234.480 mcd mid-point
    % 2A-097Q-1,  23- 52 cm = 234.067 mcd mid-point
    % 2A-111Q-1,  50- 63 cm = 276.725 mcd mid-point
    % 2A-117Q-2,  25- 30 cm = 292.118 mcd mid-point
    deino = [ 68.855  172    12 
              75.000  155     7
             234.480  428.2   7.7 
             234.067  438    37
             276.725  561    14
             292.118  628    11];

% Display results.
    
    figure('Color',[1 1 1],...
         'Position',[100 600 600 600]);
    axes('Position',[0.2 0.15 0.675 0.7],...
         'Box','On',...
         'FontSize',12,...
         'LineWidth',1,...
         'XLim',[0 700],...
         'YLim',[0 300],...
         'YDir','Reverse',...
         'XAxisLocation','Top'); hold on
    for loopi = 1 : length(allbaseages)
    line(tiepointstune(:,2,loopi),tiepointstune(:,1,loopi),...
        'LineWidth',1,...
        'LineStyle','-',....
        'Color',[0.9290   0.6940   0.1250]);
    end
    axes('Position',[0.2 0.15 0.675 0.7],...
         'Box','On',...
         'FontSize',12,...
         'LineWidth',1,...
         'XLim',[0 700],...
         'YLim',[0 300],...
         'YDir','Reverse',...
         'XAxisLocation','Top',...
         'Color','none'); hold on
    line(tiepointstune(:,2,Locb(1)),tiepointstune(:,1,Locb(1)),...
        'LineWidth',1,...
        'LineStyle','-',....
        'Color',[0 0 0]);
    line(tiepointstune(:,2,Locb(2)),tiepointstune(:,1,Locb(2)),...
        'LineWidth',1.5,...
        'LineStyle','-.',....
        'Color',[0 0 0]);
    line(tiepointstune(:,2,Locb(3)),tiepointstune(:,1,Locb(3)),...
        'LineWidth',1.5,...
        'LineStyle',':',....
        'Color',[0 0 0]);
    line(tiepointstune(:,2,length(allbaseages)),...
        tiepointstune(:,1,length(allbaseages)),...
        'LineWidth',1.5,...
        'LineStyle','-',....
        'Color',[0.9290   0.6940   0.1250]);
    errorbar(deino(:,2),deino(:,1),deino(:,3),...
        'Horizontal',...
        'LineStyle','None',...
        'LineWidth',1,......
        'Color','k');
    xlabel('Age (kyr BP)');
    ylabel('Composite Depth (m)');
    legend(agemodeltunestr(1,:),...
        agemodeltunestr(2,:),...
        agemodeltunestr(3,:),...
        'ALL',...
        'Box','off');

    printname = ...
        'DirFiguresEnsemble/trauth_agemodels_supplfigs_0X_ensemble.png';
    print(printname,'-dpng','-r300')

    end

    close all
