%% Script 3: Analyzing the Chew Bahir age models
% 
% 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 = rrmarch2021+mhtXXX

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

  clear, clc, close all

% Define all base ages.

   allbaseages = [550;560;570];

% Loop.

   for loopi = 1 : length(allbaseages)

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

   clearvars -except allbaseages loopi
   clc, close all
    
% Define tuned age model to be used.

   agemodeltune = allbaseages(loopi);

% Create character string for age models.

   agemodeltunemergestr = ['RRMarch2021+MHT',num2str(agemodeltune(1))];

% Create directories.

    mkdir DirFiguresAnalyzing

% Importing age models and data.

    % Import tiepoints of all age models.
    dirname = 'DirDataTiepoints/agemodeltiepoints_';
    tiepoints1 = load([dirname,'rrmarch2021.txt']);
    tiepoints1(:,2) = tiepoints1(:,2)/1000;
    tiepoints2 = load([dirname,'mubawa.txt']);
    dirname = 'DirResultsTiePointsMerging/agemodeltiepoints_';
    tiepoints3 = load([dirname,'rrmarch2021_mht',...
        num2str(agemodeltune(1)),'.txt']);
    
    % Load PC2 of color reflectance values. Since high PC2 scores reflect
    % wet conditions during precession minima we are using change the sign
    % of the PC2 scores.
    filename = 'DirDataOther/pc2_colref.txt';
    data = load(filename);
    data(:,2) = -data(:,2);
    
    % Remove identical values in the first column.
    clear data2b
    data = sortrows(data,1);
    [datai,IA,IC] = unique(data(:,1));
    datai(:,2) = data(IA,2);
    
    % Remove NaNs in the 2nd column.
    datai(isnan(datai(:,2))==1,:) = [];
    
    % Load Laskar 2004 orbital parameters.
    close all
    orbital = load('DirDataOther/insolaout6HoGWM.txt');
    orbital(:,1) = -orbital(:,1);
    orbital = flipud(orbital);
                
    % 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];

% Interpolating color data and argon/argon ages to age models.

    % Interpolate PC2 of color reflectance values to age models.
    dataage1(:,1) = ...
        interp1(tiepoints1(:,1),tiepoints1(:,2),datai(:,1),'pchip');
    dataage1(:,2) = datai(:,2);
    dataage2(:,1) = ...
        interp1(tiepoints2(:,1),tiepoints2(:,2),datai(:,1),'pchip');
    dataage2(:,2) = datai(:,2);
    dataage3(:,1) = ...
        interp1(tiepoints3(:,1),tiepoints3(:,2),datai(:,1),'pchip');
    dataage3(:,2) = datai(:,2);
    
    % Interpolate the argon/argon ages to age models.
    deinoi_tune1(:,1) = interp1(tiepoints1(:,1),tiepoints1(:,2),...
        deino(:,1),'pchip');
    deinoi_tune1(:,2) = deino(:,2);
    deinoi_tune1(:,3) = deino(:,3);
    deinoi_tune2(:,1) = interp1(tiepoints2(:,1),tiepoints2(:,2),...
        deino(:,1),'pchip');
    deinoi_tune2(:,2) = deino(:,2);
    deinoi_tune2(:,3) = deino(:,3);
    deinoi_tune3(:,1) = interp1(tiepoints3(:,1),tiepoints3(:,2),...
        deino(:,1),'pchip');
    deinoi_tune3(:,2) = deino(:,2);
    deinoi_tune3(:,3) = deino(:,3);

% Spectral analysis.

    % We can calculate the normalized Lomb-Scargle power spectrum and false
    % alarm probabilities using the MATLAB function plomb. Annotate the
    % levels that correspond to false-alarm probabilities confid.
    dataage1 = sortrows(dataage1);
    dataage2 = sortrows(dataage2);
    dataage3 = sortrows(dataage3);
    [px1,f1] = plomb(dataage1(:,2),dataage1(:,1),'normalized');
    [px2,f2] = plomb(dataage2(:,2),dataage2(:,1),'normalized');
    [px3,f3] = plomb(dataage3(:,2),dataage3(:,1),'normalized');

% Display results.

    script_display_analyzing

    end

    close all

