%% Age-Elevation Relationship (AER)
% This script to plot(U-Th-Sm)/He thermochronology data trends in age, elevation,
% eU, and ESR.
%
% Input:
% A plain text file with 5 columns in order:
%   1. FT corrected age (Ma)
%   2. Error (Ma)
%   3. Elevation (m)
%   4. eU (ppm)
%   5. ESR (mm)
%
% See "..Example/NP.txt" for a template dataset (based on Macaulay et al., 2013
% and unpublished data by Lingxiao Gong).
%
% Output:
% A bubble chart with error bars, saved under "Example/" as both PDF and PNG.
%
% Lingxiao Gong (gong@uni-potsdam.de) 
% First release: April 2024
% Update: September 2025


clear all
close all

%% Paths
% Automatically set path relative to script location
% Script is inside AER_eU_ESR/Code
% Example data and output figures go to ../Example/
script_path  = fileparts(mfilename('fullpath'));
main_path    = fileparts(script_path);
example_path = fullfile(main_path, 'Example');

if ~exist(example_path, 'dir')
    mkdir(example_path);
end
% addpath(genpath('D:\MATLAB\MATLAB\toolbox\Crameri'));% colormap created by Fabio Crameri, need to download.

%% USER INPUT %% 
data_file = fullfile(example_path, 'NP.txt'); % Example input file, change to your file name
P = readmatrix(data_file);

A   = P(:,1); % Age [Ma]
Er  = P(:,2); % Error [Ma]
ELE = P(:,3); % Elevation [m]
eU  = P(:,4); % eU [ppm]
ESR = P(:,5); % ESR [mm]

%% Make Figure
figure(1)%(x,y,size,color)
bc = bubblechart(A,ELE,ESR,eU); % (x,y,size,color)
% Uncomment below if Fabio Crameri colormap is available
% crameri('hawaii');

b = bubblelegend('ESR [mm]','location','northeast');%add bubble legend

c = colorbar % add colorbar
c.Label.FontSize = 10;
%c =('Ticks',[10:10:60])% change the spacing of colorbar if you dislike the default one
c.Label.String = 'eU [ppm]'% add colorbar legend
xlabel('Age [Ma]')
ylabel('Elevation [m]')

% if you want to set the range of the axises
% ylim([2100 2700])
% xlim([0 30])
% ylim([2400 3700])
% xticks(0:5:30)
% yticks(2000:100:2700)
grid on;
bubblesize([5 25])% bubble size range

% Add horizontal error bars
hold on;
eb = errorbar(A, ELE, Er, 'horizontal', 'LineStyle','none');
set(eb, 'color', 'k', 'LineWidth', 1.2);

title('AER - North Profile')% NAME YOUR PLOTTING

%% Save plots
out_pdf = fullfile(example_path, 'NP_AER.pdf');
out_png = fullfile(example_path, 'NP_AER.png');

set(gcf, 'Units', 'inches');
screenposition = get(gcf, 'Position');
set(gcf, ...
    'PaperPosition',[0 0 screenposition(3:4)], ...
    'PaperSize',    [screenposition(3:4)]);

print(gcf, '-dpdf', '-painters', out_pdf);
saveas(gcf, out_png);