%% Script 3: Compare merged K record from CHB14-2 with ODP 967 data
% 
% 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.

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

    mkdir DirFiguresCorrelation

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

  clear, clc, close all

% Define all base ages.

   allbaseages = 560;

% 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 all age models.
    agemodeltunemergestr = strcat(...
        "RRMarch2021_MHT",num2str(agemodeltune(1)));

% Chew Bahir CHB14-2 - age model 1 RRMarch2021.

    % Load long core data.
    data_cbh14_1 = ...
       load('DirResultsLongCores/data_chb14_2_age_orig_rrmarch2021.txt');
    data_cbh14_1(:,1) = data_cbh14_1(:,1)/1000;
    data_cbh14_1(:,2) = data_cbh14_1(:,2) - mean(data_cbh14_1(:,2),'omitnan');
    data_cbh14_1(:,2) = data_cbh14_1(:,2) / std(data_cbh14_1(:,2),'omitnan');

% Chew Bahir CHB14-2 - age model 2 MUBAWA.

    % Load long core data.
    data_cbh14_2 = ...
       load('DirResultsLongCores/data_chb14_2_age_orig_mubawa.txt');
    data_cbh14_2(:,1) = data_cbh14_2(:,1);
    data_cbh14_2(:,2) = data_cbh14_2(:,2) - mean(data_cbh14_2(:,2),'omitnan');
    data_cbh14_2(:,2) = data_cbh14_2(:,2) / std(data_cbh14_2(:,2),'omitnan');

% Chew Bahir CHB14-2 - age model 7 RRMarch2021+MHTXXX.
    
    % Load long core data.
    data_cbh14_3 = ...
        load(strcat("DirResultsLongcores/data_chb14_2_age_orig_",...
        agemodeltunemergestr(1,:),".txt"));

    % Normalize data.
    data_cbh14_3(:,1) = data_cbh14_3(:,1);
    data_cbh14_3(:,2) = data_cbh14_3(:,2) - mean(data_cbh14_3(:,2),'omitnan');
    data_cbh14_3(:,2) = data_cbh14_3(:,2) / std(data_cbh14_3(:,2),'omitnan');
    
    % Load short core data, now with age model based on IntCal20 14C ages. 
    data_cb01 = load('DirDataOther/data_cb01_cb03_InCal20.txt');
    data_cb01(:,1) = data_cb01(:,1)/1000;
    data_cb01(:,2) = data_cb01(:,2) - mean(data_cb01(:,2),'omitnan');
    data_cb01(:,2) = data_cb01(:,2) / std(data_cb01(:,2),'omitnan');
    
    % Merge data.
    data_cbh14_3_merged = vertcat(data_cb01(data_cbh14_3(:,1)<20,:),...
        data_cbh14_3(data_cbh14_3(:,1)>20,:));
    data_cbh14_3_merged(data_cbh14_3_merged(:,1)<20,2) = ...
        data_cb01(data_cb01(:,1)<20,2) + 0.85;   % Offset

% ODP Site 967 (Grant et al. 2017)
    %     1 = age (kyrs BP)
    %     2 = PC1 sapropel
    %     3 = PC2 sapropel
    %     4 = dust proxy
    %     5 = wetness index

    fid = fopen('DirDataOther/data_odp_967_grant.txt');
    formstrg = [repmat('%s ',1,4),repmat('%f ',1,23)];
    A1 = textscan(fid,formstrg,'Delimiter','\t',...
        'Headerlines',9);
    fclose(fid);
    for i = 1:23
        datagrant(:,i) = A1{i+4};
    end
    dataodp967(:,1) = datagrant(:,22);       % Age
    dataodp967(:,2) = datagrant(:,18);       % PC1 sapropel
    dataodp967(:,3) = datagrant(:,19);       % PC2 detritus
    dataodp967(:,4) = datagrant(:,21);       % Dust proxy
    dataodp967(:,5) = datagrant(:,23);       % Wetness index

% Remove gaps.

    data_cbh14_1(isnan(data_cbh14_1(:,1)),:) = [];
    data_cbh14_2(isnan(data_cbh14_2(:,1)),:) = [];
    data_cbh14_3(isnan(data_cbh14_3(:,1)),:) = [];
    data_cbh14_3_merged(isnan(data_cbh14_3_merged(:,1)),:) = [];
    dataodp967(isnan(dataodp967(:,1)),:) = [];

% Interpolation after removing duplicate data points.

    t = 0 : 0.1 : 700;

    % 1

    data_cbh14_1 = sortrows(data_cbh14_1);

    data_cbh14_12 = data_cbh14_1;
    for i = 1:size(data_cbh14_12,1)-1
        if data_cbh14_12(i,1) == data_cbh14_12(i+1,1)
            data_cbh14_12(i,1) = NaN;
        end
    end
    data_cbh14_12(isnan(data_cbh14_12(:,1))==1,:) = [];
    data_cbh14_1 = data_cbh14_12;
    clear data_cbh14_12

    data_cbh14_1_intp(:,1) = t;
    data_cbh14_1_intp(:,2) = interp1(data_cbh14_1(:,1),...
        data_cbh14_1(:,2),t,'linear','extrap');

    % 2

    data_cbh14_2 = sortrows(data_cbh14_2);

    data_cbh14_22 = data_cbh14_2;
    for i = 1:size(data_cbh14_22,1)-1
        if data_cbh14_22(i,1) == data_cbh14_22(i+1,1)
            data_cbh14_22(i,1) = NaN;
        end
    end
    data_cbh14_22(isnan(data_cbh14_22(:,1))==1,:) = [];
    data_cbh14_2 = data_cbh14_22;
    clear data_cbh14_22

    data_cbh14_2_intp(:,1) = t;
    data_cbh14_2_intp(:,2) = interp1(data_cbh14_2(:,1),...
        data_cbh14_2(:,2),t,'linear','extrap');

    % 3

    data_cbh14_3 = sortrows(data_cbh14_3);

    data_cbh14_32 = data_cbh14_3;
    for i = 1:size(data_cbh14_32,1)-1
        if data_cbh14_32(i,1) == data_cbh14_32(i+1,1)
            data_cbh14_32(i,1) = NaN;
        end
    end
    data_cbh14_32(isnan(data_cbh14_32(:,1))==1,:) = [];
    data_cbh14_3 = data_cbh14_32;
    clear data_cbh14_32

    data_cbh14_3_intp(:,1) = t;
    data_cbh14_3_intp(:,2) = interp1(data_cbh14_3(:,1),...
        data_cbh14_3(:,2),t,'linear','extrap');

    % 3 merged

    data_cbh14_3_merged = sortrows(data_cbh14_3_merged);

    data_cbh14_3_merged2 = data_cbh14_3_merged;
    for i = 1:size(data_cbh14_3_merged2,1)-1
        if data_cbh14_3_merged2(i,1) == data_cbh14_3_merged2(i+1,1)
            data_cbh14_3_merged2(i,1) = NaN;
        end
    end
    data_cbh14_3_merged2(isnan(data_cbh14_3_merged2(:,1))==1,:) = [];
    data_cbh14_3_merged = data_cbh14_3_merged2;
    clear data_cbh14_3_merged2

    data_cbh14_3_merged_intp(:,1) = t;
    data_cbh14_3_merged_intp(:,2) = interp1(data_cbh14_3_merged(:,1),...
        data_cbh14_3_merged(:,2),t,'linear','extrap');

    % ODP 967

    dataodp967_int(:,1) = t;
    for i = 2:5
    dataodp967_int(:,i) = interp1(dataodp967(:,1),...
        dataodp967(:,i),t,'linear','extrap');
    end

% Piecewise cross correlation within window of size w. The window size of n
% corresponds to n times 100 years which is the resolution of the t vector.

    w = 500; 

    p_data_cbh14_1(:,1) = t;
    for i = w/2+1:length(t)-w/2
        clear p
        p = corr(data_cbh14_1_intp(i-w/2:i+w/2,2),...
            dataodp967_int(i-w/2:i+w/2,2),'Type','Spearman');
        p_data_cbh14_1(i,2) = p;
    end
    p_data_cbh14_1(1:w/2,2) = p_data_cbh14_1(w/2+1,2) * ones(w/2,1); 
    p_data_cbh14_1(length(t)-w/2+1:length(t),2) = ...
        p_data_cbh14_1(length(t)-w/2,2) * ones(w/2,1);
    p_data_cbh14_1(p_data_cbh14_1(:,1)>max(data_cbh14_1(:,1)),:) = [];

    p_data_cbh14_2(:,1) = t;
    for i = w/2+1:length(t)-w/2
        clear p
        p = corr(data_cbh14_2_intp(i-w/2:i+w/2,2),...
            dataodp967_int(i-w/2:i+w/2,2),'Type','Spearman');
        p_data_cbh14_2(i,2) = p;
    end
    p_data_cbh14_2(1:w/2,2) = p_data_cbh14_2(w/2+1,2) * ones(w/2,1); 
    p_data_cbh14_2(length(t)-w/2+1:length(t),2) = ...
        p_data_cbh14_2(length(t)-w/2,2) * ones(w/2,1);
    p_data_cbh14_2(p_data_cbh14_2(:,1)>max(data_cbh14_2(:,1)),:) = [];

    p_data_cbh14_3(:,1) = t;
    for i = w/2+1:length(t)-w/2
        clear p
        p = corr(data_cbh14_3_intp(i-w/2:i+w/2,2),...
            dataodp967_int(i-w/2:i+w/2,2),'Type','Spearman');
        p_data_cbh14_3(i,2) = p;
    end
    p_data_cbh14_3(1:w/2,2) = p_data_cbh14_3(w/2+1,2) * ones(w/2,1); 
    p_data_cbh14_3(length(t)-w/2+1:length(t),2) = ...
        p_data_cbh14_3(length(t)-w/2,2) * ones(w/2,1);
    p_data_cbh14_3(p_data_cbh14_3(:,1)>max(data_cbh14_3(:,1)),:) = [];

    p_data_cbh14_3_merged(:,1) = t;
    for i = w/2+1:length(t)-w/2
        clear p
        p = corr(data_cbh14_3_merged_intp(i-w/2:i+w/2,2),...
            dataodp967_int(i-w/2:i+w/2,2),'Type','Spearman');
        p_data_cbh14_3_merged(i,2) = p;
    end
    p_data_cbh14_3_merged(1:w/2,2) = p_data_cbh14_3_merged(w/2+1,2) * ones(w/2,1); 
    p_data_cbh14_3_merged(length(t)-w/2+1:length(t),2) = ...
        p_data_cbh14_3_merged(length(t)-w/2,2) * ones(w/2,1);
    p_data_cbh14_3_merged(p_data_cbh14_3_merged(:,1)> ...
        max(data_cbh14_3_merged(:,1)),:) = [];
    
% Display the data

    script_display_correlation_1_rrmarch2021
    script_display_correlation_2_mubawa
    script_display_correlation_3_rrmarch2021_mhtXXX
    script_display_correlation_3_rrmarch2021_mhtXXX_cb0103
    script_display_correlation_3_compare

    close all

    end