%% Script 1: Traditional tuning 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

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

    clear, clc, close all

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

    mkdir DirFiguresTuning
    mkdir DirResultsTiePointsTuning
    mkdir DirResultsDataTuning

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

  allbaseages = 450 : 10 : 650;

% Figure labels
  abclabels = 'A':'Z';

% Loop for baseages

for loopi = 1 : length(allbaseages)

% Setting up tuning experiment.

    % Clear variables except for loopi for runn and tiepointscell.
    clearvars -except loopi allbaseages tiepointscell abclabels
    
    % Define essential parameters for age-depth modeling.
    proxy =  6;                   % Define proxy you wish to load (e.g. 6).
    printpng = 1;                 % 0=do not save, 1=save to file.
    baseage = allbaseages(loopi); % Initial age of the core (in kyrs).
    XLimLaskar = baseage + 300;   % XLim of Laskar data (in kyrs).
    XLimPlot = [0 700];           % XLim of final plots (in kyrs).
    YLimPlot = [-8 6];          % YLim of final plots (in kyrs).
    intv = 0.8;                   % Initial sampling interval (in kyrs).
    cutoff = [29 17];             % Passband of the filter (in kyrs).
    
% Importing data.

    % 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.
    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)),:) = [];
    
    % 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];

% Start tuning.

    % Defining the initial age of the base of the core.
    agerange = [0 0; max(data(:,1)) baseage];
    
    % We create an initial age axis running from 0 to baseage kyrs BP, with
    % 0.1 kyrs intervals.
    interv = 0.1;
    agemodel(:,1) = 0 : interv : baseage;
    agemodel(:,2) = ...
        interp1(agerange(:,1),agerange(:,2),agemodel(:,1),'pchip');
    dataii(:,1) = ...
        interp1(agemodel(:,1),agemodel(:,2),datai(:,1),'pchip');
    dataii(:,2) = datai(:,2);
    
    deinoi(:,1) = interp1(agemodel(:,1),agemodel(:,2),...
        deino(:,1),'pchip');
    deinoi(:,2) = deino(:,2);
    deinoi(:,3) = deino(:,3);
    
    % Interpolate data upon an evenly spaced time axes.
    clear data2i
    dataiii(:,1) = min(dataii(:,1)) : interv : max(dataii(:,1));
    dataiii(:,2) = ...
        interp1(dataii(:,1),dataii(:,2),dataiii(:,1),'pchip');
    
    % Tuning starts here! Load Laskar 2004 precession values.
    close all
    orbital = load('DirDataOther/insolaout6HoGWM.txt');
    orbital(:,1) = -orbital(:,1);
    orbital = flipud(orbital);
    
    % Decimate both proxy and precession values to a sampling interval intv
    % and display the resulting proxy and precession records. The limit for
    % Laskar data must be larger (e.g. +300 kyr) than the baseage.
    dataiiid(:,1) = 1:intv:baseage;
    dataiiid(:,2) = ...
        interp1(dataiii(:,1),dataiii(:,2),dataiiid(:,1),'pchip');
    
    orbitaliii(:,1) = 1:intv:XLimLaskar;       
    for i = 2:4
    orbitaliii(:,i) = ...
        interp1(orbital(:,1),orbital(:,i),orbitaliii(:,1),'pchip');
    end
    
    % Filter design. Laskar has a nice 23 and 19 kyr cycle.
    [b,a] = butter(6,1./cutoff/(0.5/intv),'bandpass');
    [h,w] = freqz(b,a,1024);
    f = (1/intv)*w/(2*pi);
    clear dataiiid_filt
    
    dataiiid_filt(:,1) = dataiiid(:,1);
    dataiiid_filt(:,2) = filtfilt(b,a,dataiiid(:,2));
    
    % Find maxima (= peaks) in both proxy and precession values. The
    % minimum distance of two subsequent peaks is 10 data points. Display
    % (1) the proxy record, (2) the filtered proxy and the (3) Laskar 2004
    % precession record.
    [pks1,locs1] = findpeaks(dataiiid_filt(:,2),'MinPeakDistance',10);
    [pks2,locs2] = findpeaks(orbitaliii(:,3),'MinPeakDistance',10);
    
    % Tuning peak by peak using by using the locations of the maxima of the
    % Laskar record and used them to interpolate the proxy record. This
    % then puts the maxima of the proxy record exactly at the same
    % locations as in the Laskar record.
    tiepoints(:,1) = dataiiid_filt(locs1,1);
    tiepoints(:,2) = orbitaliii(locs2(1:length(pks1)),1);
    
    % Interpolate filtered tuned data
    dataiiid_tune(:,1) = ...
        interp1(tiepoints(:,1),tiepoints(:,2),dataiiid_filt(:,1),'pchip');
    dataiiid_tune(:,2) = dataiiid_filt(:,2);
    
    % Interpolate tuned data.
    dataiii_tune(:,1) = ...
        interp1(tiepoints(:,1),tiepoints(:,2),dataii(:,1),'pchip');
    dataiii_tune(:,2) = dataii(:,2);
    
    % Interpolate tune ages.
    deinoi_tune(:,1) = ...
        interp1(tiepoints(:,1),tiepoints(:,2),deinoi(:,1),'pchip');
    deinoi_tune(:,2) = deinoi(:,2);
    deinoi_tune(:,3) = deinoi(:,3);
    
    % Interpolate everything to the same evenly-spaced time axis.
    dataiii_tuned(:,1) = 1:intv:max(dataiii_tune(:,1));
    dataiii_tuned(:,2) = interp1(dataiii_tune(:,1),dataiii_tune(:,2),...
        dataiii_tuned(:,1),'pchip');
    
    dataiiid_tuned(:,1) = 1:intv:max(dataiii_tune(:,1));
    dataiiid_tuned(:,2) = interp1(dataiiid_tune(:,1),dataiiid_tune(:,2),...
        dataiiid_tuned(:,1),'pchip');
    
    dataiii_tunedt(:,2) = detrend(dataiii_tuned(:,2));
    dataiiid_tunedt(:,2) = detrend(dataiiid_tuned(:,2));
    
    [Pxx,~] = periodogram(dataiii_tunedt(:,2),[],1024,1/intv);
    [Pxxd,f] = periodogram(dataiiid_tunedt(:,2),[],1024,1/intv);
    
    % Extracting age-depth model.
    tiepoints(:,1) = interp1(agemodel(:,2),agemodel(:,1),...
        tiepoints(:,1),'pchip');
    
    % Tuning peak by peak using by using the locations of the maxima of the
    % Laskar record and used them to interpolate the proxy record. This
    % then puts the maxima of the proxy record exactly at the same
    % locations as in the Laskar record.
    tiepoints(:,1) = dataiiid_filt(locs1,1);
    tiepoints(:,2) = orbitaliii(locs2(1:length(pks1)),1);
    
    % Interpolate filtered tuned data.
    dataiiid_tune(:,1) = interp1(tiepoints(:,1),tiepoints(:,2),...
        dataiiid_filt(:,1),'pchip');
    dataiiid_tune(:,2) = dataiiid_filt(:,2);
    
    % Interpolate tuned data.
    dataiii_tune(:,1) = interp1(tiepoints(:,1),tiepoints(:,2),...
        dataii(:,1),'pchip');
    dataiii_tune(:,2) = dataii(:,2);
    
    % Interpolate tuned ages.
    deinoi_tune(:,1) = interp1(tiepoints(:,1),tiepoints(:,2),...
        deinoi(:,1),'pchip');
    deinoi_tune(:,2) = deinoi(:,2);
    deinoi_tune(:,3) = deinoi(:,3);
    
    % Interpolate original data and calculate powerspectrum.
    interv_orig = 0.1;       % in meters
    datai_itp(:,1) = ...
        min(datai(:,1)) : interv_orig : max(datai(:,1));
    datai_itp(:,2) = ...
        interp1(datai(:,1),datai(:,2),datai_itp(:,1),'pchip');
    [Pxx_orig ,f_orig ] = periodogram(datai_itp(:,2),...
        [],32768,1/interv_orig);
    
    % Interpolate tuned data and calculate powerspectrum.
    interv_tune = 0.1;       % in kyr
    dataiii_tune_itp(:,1) = ...
        min(dataiii_tune(:,1)) : interv_tune : max(dataiii_tune(:,1));
    dataiii_tune_itp(:,2) = interp1(dataiii_tune(:,1),dataiii_tune(:,2),...
        dataiii_tune_itp(:,1),'pchip');
    [Pxx_tuned,f_tuned] = periodogram(dataiii_tune_itp(:,2),...
        [],32768,1/interv_tune);

% Display the results.

script_display_tuning

% Save data.

    % Save tiepoints to ASCII files.
    tiepoints_export(:,1) = interp1(agemodel(:,2),agemodel(:,1),...
        tiepoints(:,1),'pchip');
    tiepoints_export(:,2) = tiepoints(:,2);
    
    savename_tiepoints = [
        'DirResultsTiePointsTuning/agemodeltiepoints_mht',...
        num2str(baseage),'.txt'];
    save(savename_tiepoints,'tiepoints_export','-ascii')
    
    % Save tune data to ASCII files.
    savename_datatune = [
        'DirResultsDataTuning/data_chb14_2_mht',...
        num2str(baseage),'.txt'];
    save(savename_datatune,'dataiii_tune','-ascii')
    
    % Save graphics to file.
    if printpng == 1
        printname = ['DirFiguresTuning/trauth_agemodels_supplfigs_01_mht',...
           num2str(baseage),...
           '.png'];
        print(printname,'-dpng','-r300')
    end

    if printpng == 1 && loopi == 12
        printname = ...
           'DirFiguresTuning/trauth_agemodels_figure_04_tuning.png';
        print(printname,'-dpng','-r300')
    end

% End of for loop for different base ages and close all figures.

end

%close all
















