%% Script 2: Merging and comparing age models of the Chew Bahir record
% 
% 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
%
%    1 = rrmarch2021
%    2 = mubawa
%    3 = mhtXXX
%    4 = mhtXXX
%    5 = mhtXXX
%    6 = rrmarch2021+mhtXXX
%    7 = rrmarch2021+mhtXXX
%    8 = rrmarch2021+mhtXXX

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

    clear, clc, close all

% Define all base ages.

   allbaseages = [550;560;570];

% Define tuned age models to be used.

   agemodeltune    = allbaseages;

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

    warning('off','MATLAB:MKDIR:DirectoryExists');

    mkdir DirFiguresMerging
    mkdir DirResultsTiePointsMerging

% Importing age models and data.

    % Import tiepoints of all age models.
    dirname = 'DirDataTiepoints/agemodeltiepoints_';
    tiepoints1 = load([dirname,'rrmarch2021.txt']);
    tiepoints1(:,2:3) = tiepoints1(:,2:3)/1000;  % Convert to kyr
    tiepoints1(:,5:8) = tiepoints1(:,5:8)/1000;  % Convert to kyr
    tiepoints1 = tiepoints1(1:1000:end,:);       % Decimate
    tiepoints2 = load([dirname,'mubawa.txt']);
    dirname = 'DirResultsTiePointsTuning/agemodeltiepoints_';
    tiepoints3 = load([dirname,agemodeltunestr(1,:),'.txt']); 
    tiepoints4 = load([dirname,agemodeltunestr(2,:),'.txt']);
    tiepoints5 = load([dirname,agemodeltunestr(3,:),'.txt']);
              
    % 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];
    
    % Define Marker.
    agemarker = 'none';
   
% Comparing age-depth models.

    figure('Color',[1 1 1],...
         'Position',[100 600 600 600]);
    axes('Position',[0.2 0.15 0.675 0.7],...
         'Box','On',...
         'FontSize',14,...
         'LineWidth',1,...
         'XLim',[0 700],...
         'YLim',[0 300],...
         'YDir','Reverse',...
         'XAxisLocation','Top'); hold on
    line(tiepoints1(:,2),tiepoints1(:,1),...
        'LineWidth',1,...
        'Marker',agemarker,...
        'LineStyle','-',...
        'Color',[0.8500   0.3250   0.0980]); % RRMarch2021
    line(tiepoints2(:,2),tiepoints2(:,1),...
        'LineWidth',1,...
        'Marker',agemarker,...
        'LineStyle','-',...    
        'Color',[0   0.4470   0.7410]);      % MUBAWA
    line(tiepoints3(:,2),tiepoints3(:,1),...
        'LineWidth',1,......
        'Marker',agemarker,...
        'LineStyle','-',...
        'Color',[0.9290   0.6940   0.1250]); % Tuned age model 1
    line(tiepoints4(:,2),tiepoints4(:,1),...
        'LineWidth',1,......
        'Marker',agemarker,...
        'LineStyle','-.',...
        'Color',[0.9290   0.6940   0.1250]); % Tuned age model 2
    line(tiepoints5(:,2),tiepoints5(:,1),...
        'LineWidth',1,......
        'Marker',agemarker,...
        'LineStyle',':',...
        'Color',[0.9290   0.6940   0.1250]); % Tuned age model 3
    errorbar(deino(:,2),deino(:,1),deino(:,3),...
        'Horizontal',...
        'LineStyle','None',...
        'LineWidth',1,......
        'Color','k');
    xlabel('Age (kyr BP)');
    ylabel('Composite Depth (m)');
    legend('RRMarch2021','MUBAWA',...
         agemodeltunestr(1,:),...
         agemodeltunestr(2,:),...
         agemodeltunestr(3,:),...
        'Box','off');

    printname = strcat("DirFiguresMerging/",...
        "trauth_agemodels_figure_05_agedepth_mht",...
        num2str(agemodeltune(1)),"_",...
        num2str(agemodeltune(2)),"_",...
        num2str(agemodeltune(3)),".png");
    print(printname,'-dpng','-r300')

% Comparing differences from (0-600 kyr)/(0-290 m) linear model vs. depth.

    deltatiepoints1 = tiepoints1(:,2) - 600/290 * tiepoints1(:,1);
    deltatiepoints2 = tiepoints2(:,2) - 600/290 * tiepoints2(:,1);
    deltatiepoints3 = tiepoints3(:,2) - 600/290 * tiepoints3(:,1);
    deltatiepoints4 = tiepoints4(:,2) - 600/290 * tiepoints4(:,1);
    deltatiepoints5 = tiepoints5(:,2) - 600/290 * tiepoints5(:,1);

    figure('Position',[150 550 1000 400],...
        'Color',[1 1 1]);
    axes('Position',[0.15 0.2 0.7 0.7],...
         'Box','On',...
         'FontSize',14,...
         'LineWidth',1,...
         'XLim',[0 300],...
         'YLim',[-100 100],...
         'YTick',-100:20:100); hold on
    line([0 700],[0 0],...
        'LineWidth',1,...
        'LineStyle',':',....
        'Color',[0 0 0]);                    % Linear
    line(tiepoints1(:,1),deltatiepoints1(:,1),...
        'LineWidth',1,...
        'LineStyle','-',....
        'Color',[0.8500   0.3250   0.0980]); % RRmarch2021
    line(tiepoints2(:,1),deltatiepoints2(:,1),...
        'LineWidth',1,...
        'LineStyle','-',....
        'Color',[0   0.4470   0.7410]);      % MUBAWA
    line(tiepoints3(:,1),deltatiepoints3(:,1),...
        'LineWidth',1,...
        'LineStyle','-',....
        'Color',[0.9290   0.6940   0.1250]); % MHT530
    line(tiepoints4(:,1),deltatiepoints4(:,1),...
        'LineWidth',1,...
        'LineStyle','-.',....
        'Color',[0.9290   0.6940   0.1250]); % MHT570
    line(tiepoints5(:,1),deltatiepoints5(:,1),...
        'LineWidth',1,...
        'LineStyle',':',...
        'Color',[0.9290   0.6940   0.1250]); % MHT610
    xlabel('Composite Depth (m)');
    ylabel('Delta-Age (kyr BP)');
    legend('Linear','RRMarch2021','MUBAWA',...
         agemodeltunestr(1,:),...
         agemodeltunestr(2,:),...
         agemodeltunestr(3,:),...
        'Box','off');

    printname = strcat("DirFiguresMerging/",...
        "trauth_agemodels_supplfigs_02_depthdeltaage_mht",...
        num2str(agemodeltune(1)),"_",...
        num2str(agemodeltune(2)),"_",...
        num2str(agemodeltune(3)),".png");
    print(printname,'-dpng','-r300')

% Comparing differences from (0-600 kyr)/(0-290 m) linear model vs. age.

    figure('Position',[200 500 1000 400],...
        'Color',[1 1 1]);
    axes('Position',[0.15 0.2 0.7 0.7],...
         'Box','On',...
         'FontSize',14,...
         'LineWidth',1,...
         'XLim',[0 700],...
         'YLim',[-100 100],...
         'YTick',-100:20:100); hold on
    line([0 700],[0 0],...
        'LineWidth',1,...
        'LineStyle',':',....
        'Color',[0 0 0]);                    % Linear
    line(tiepoints1(:,2),deltatiepoints1(:,1),...
        'LineWidth',1,...
        'LineStyle','-',....
        'Color',[0.8500   0.3250   0.0980]); % RRMarch2021
    line(tiepoints2(:,2),deltatiepoints2(:,1),...
        'LineWidth',1,...
        'LineStyle','-',....
        'Color',[0   0.4470   0.7410]);      % Mubawa
    line(tiepoints3(:,2),deltatiepoints3(:,1),...
        'LineWidth',1,...
        'LineStyle','-',....
        'Color',[0.9290   0.6940   0.1250]); % MHT530
    line(tiepoints4(:,2),deltatiepoints4(:,1),...
        'LineWidth',1,...
        'LineStyle','-.',....
        'Color',[0.9290   0.6940   0.1250]); % MHT570
    line(tiepoints5(:,2),deltatiepoints5(:,1),...
        'LineWidth',1,...
        'LineStyle',':',...
        'Color',[0.9290   0.6940   0.1250]); % MHT610
    xlabel('Age (kyr BP)');
    ylabel('Delta-Age (kyr BP)');
    legend('Linear','RRMarch2021','MUBAWA',...
         agemodeltunestr(1,:),...
         agemodeltunestr(2,:),...
         agemodeltunestr(3,:),...
        'Box','off');

    printname = strcat("DirFiguresMerging/",...
        "trauth_agemodels_supplfigs_03_agedeltaage_mht",...
        num2str(agemodeltune(1)),"_",...
        num2str(agemodeltune(2)),"_",...
        num2str(agemodeltune(3)),".png");
    print(printname,'-dpng','-r300')

% Merging RRMarch2021 and MHTXXX age models. The RRMarch2021 and some tuned
% age models MHTXXX intersect at ~70 m and ~240 m.

    tiepoints6 = tiepoints1(1:15,1:2);
    tiepoints6(16:15+length(tiepoints3(7:end,1:2)),:) = ...
        tiepoints3(7:end,1:2);
    tiepoints6(30:29+length(tiepoints1(49:61,1:2)),:) = ...
        tiepoints1(49:61,1:2);

    tiepoints7 = tiepoints1(1:15,1:2);
    tiepoints7(16:15+length(tiepoints4(7:end,1:2)),:) = ...
        tiepoints4(7:end,1:2);
    tiepoints7(30:29+length(tiepoints1(49:61,1:2)),:) = ...
        tiepoints1(49:61,1:2);

    tiepoints8 = tiepoints1(1:15,1:2);
    tiepoints8(16:15+length(tiepoints5(7:end,1:2)),:) = ...
        tiepoints5(7:end,1:2);
    tiepoints8(30:29+length(tiepoints1(49:61,1:2)),:) = ...
        tiepoints1(49:61,1:2);

    agemodeltunemergestr = ['RRMarch2021+',agemodeltunestr(1,:);
                            'RRMarch2021+',agemodeltunestr(2,:);
                            'RRMarch2021+',agemodeltunestr(3,:)];

% Comparing merged models with MUBAWA and RRMarch2021 with 2s conf bounds
% stored in columns 7 and 8 of tiepoints1.

    tpatch = vertcat(tiepoints1(:,1),flipud(tiepoints1(:,1)));
    xpatch = vertcat(tiepoints1(:,7),flipud(tiepoints1(:,8)));

    figure('Color',[1 1 1],...
         'Position',[250 450 600 600]);
    axes('Position',[0.2 0.15 0.675 0.7],...
         'Box','On',...
         'FontSize',14,...
         'LineWidth',1,...
         'XLim',[0 700],...
         'YLim',[0 300],...
         'YDir','Reverse',...
         'XAxisLocation','Top'); hold on
    line(tiepoints1(:,2),tiepoints1(:,1),...
        'LineWidth',1,...
        'Marker',agemarker,...
        'LineStyle','-',...
        'Color',[0.8500   0.3250   0.0980]); % RRMarch2021
    line(tiepoints1(:,7),tiepoints1(:,1),...
        'LineWidth',1,...
        'Marker',agemarker,...
        'LineStyle',':',...
        'Color',[0.8500   0.3250   0.0980]); % RRMarch2021 lower 2s
    line(tiepoints1(:,8),tiepoints1(:,1),...
        'LineWidth',1,...
        'Marker',agemarker,...
        'LineStyle',':',...
        'Color',[0.8500   0.3250   0.0980]); % RRMarch2021 lower 2s
    patch(xpatch,tpatch,...
        [0.8500   0.3250   0.0980],...
        'FaceAlpha',0.05,...
        'LineStyle','None')                  % Fill area between +/- 2s
    line(tiepoints2(:,2),tiepoints2(:,1),...
        'LineWidth',1,...
        'Marker',agemarker,...
        'LineStyle','-',...    
        'Color',[0   0.4470   0.7410]);      % MUBAWA
    line(tiepoints7(:,2),tiepoints7(:,1),...
        'LineWidth',1,......
        'Marker',agemarker,...
        'LineStyle','-',...
        'Color',[0.9290   0.6940   0.1250]); % Merged RRMarch2021+MHT560
    errorbar(deino(:,2),deino(:,1),deino(:,3),...
        'Horizontal',...
        'LineStyle','None',...
        'LineWidth',1,......
        'Color','k');
    xlabel('Age (kyr BP)');
    ylabel('Composite Depth (m)');
    legend('RRMarch2021',...
        '','','',...
        'MUBAWA',...
         agemodeltunemergestr(2,:),...
        'Box','off');

    printname = strcat("DirFiguresMerging/",...
        "trauth_agemodels_figure_06_agedepthmergedconf_mht",...
        num2str(agemodeltune(2)),".png");
    print(printname,'-dpng','-r300')

% Exporting merged age models.

    filename = ['DirResultsTiePointsMerging/',...
        'agemodeltiepoints_rrmarch2021_mht',num2str(agemodeltune(1)),'.txt'];
    save(filename,'tiepoints6','-ascii')
    filename = ['DirResultsTiePointsMerging/',...
        'agemodeltiepoints_rrmarch2021_mht',num2str(agemodeltune(2)),'.txt'];
    save(filename,'tiepoints7','-ascii')
    filename = ['DirResultsTiePointsMerging/',...
        'agemodeltiepoints_rrmarch2021_mht',num2str(agemodeltune(3)),'.txt'];
    save(filename,'tiepoints8','-ascii')

    %close all

