function [P,p]=bayes_lag_prob_discont1(ts,T_min,T_max,t1,y1,e1,t2,y2,e2,samplingrate,varargin) 
% This code was written by Giulia Sinnl and released 10/05/2022. Please, contact
% giulia.sinnl@nbi.ku.dk for any more info.
% If ussing this function, please cite the original creators of this algorithm [1], its first application to ice cores [2]
% and the paper accompanying this code [3]. 

% References:
% [1] Bronk Ramsey, C. B., van der Plicht, J., & Weninger, B. ‘Wiggle matching’radiocarbon dates. Radiocarbon, 43(2A), 381-389. doi: 10.1017/S0033822200038248 ,  2001
% [2] Adolphi, F., & Muscheler, R. Synchronizing the Greenland ice core and radiocarbon timescales over the Holocene–Bayesian wiggle-matching of cosmogenic radionuclide records. Climate of the Past, 12(1), 15-30. https://doi.org/10.5194/cp-12-15-2016 ,  2016
% [3] Sinnl G., Adolphi F., Christl M., Welten K.C.,  Woodruff T.,  Caffee M.,  Svensson A.,  Muscheler R.,   Rasmussen S.O. (2022) Synchronizing ice-core and U/Th time scales in the Last Glacial Maximum using Hulu Cave 14C and  new 10Be measurements from Greenland and Antarctica  Clim. Past. [submitted Aug 2022]

% INPUTS:
% ts=vector of offsets e.g. [-1000:1:2000]
% T_min and T_max are limits of time window being investigated e.g. T_min=23000 years b2k T_max=24300 years b2k
% t1,y1,e1 are ice-core time (e.g. GICC05 ages), modelled D14C, and D14C uncertainty
% t2,y2,e2 are Hulu time, D14C, and uncertainty (Hulu timescale already already converted to years b2k=BP+50)
% r=sampling rate (e.g. 60 years)
% varargin if plots are needed e.g. 'plots',true

if length(varargin)>1 
    if strcmp(varargin{1},'plots')
    flag_plot=varargin{2};
    else
        'wrong arg: should be >'
        'plots'
        return
    end
end

P=(zeros(length(ts),1))*nan; %initialize output probability
log_P=P;

format long g;
digits(4);

% mask hulu in the window 
mask2=t2>=T_min & t2<=T_max;
r=y2(mask2); % cut the signal
dr=e2(mask2);
t=t2(mask2);
r=detrend(r,'omitnan'); % detrend hulu data in window

% ice core data
R=y1(~isnan(y1)); 
dR=e1(~isnan(y1));
T_gicc=t1(~isnan(y1));

%some offsets are not allowed, if there is not enough ice core data
high_lim_ts=t(1)-t1(1);
low_lim_ts=t(end)-t1(end);
mask_ts_allowed=(ts>=low_lim_ts & ts<=high_lim_ts);
idx_ts_allowed=find(mask_ts_allowed);

P(~mask_ts_allowed)=nan;

if flag_plot
   figure;hold on;
   plot(t,r,'linewidth',1,'displayname','Hulu');
   plot(T_gicc,R,'linewidth',1,'displayname','Ice core');
   legend()
end

p=zeros(length(ts),length(t))*nan; % initialize matrix p: local probability to be multiplied at the end
for is=idx_ts_allowed % iterate only over the offsets actually allowed
    
    T_shift_gicc=T_gicc+ts(is); %shift ice core data        
    R_f=griddedInterpolant(T_shift_gicc,R,'linear','none');
    dR_f=griddedInterpolant(T_shift_gicc,dR,'linear','none');
    
    R_det=detrend(R_f(t),'omitnan'); %resample at the Hulu ages "t" and detrend ice core data
	dR_det=dR_f(t);
    
    if flag_plot &  ismember(ts(is),[-400,0,400])
                plot(t,R_det,'linewidth',1,'displayname',num2str(ts(is)))
    end
    
    for i=1:length(t)
        
        d=abs(r(i)-R_det(i));
        s=dr(i)^2+dR_det(i)^2;
        
        p(is,i)= exp(-0.5*d^2/s)/sqrt(s);
    
    end

end

for is=1:length(ts)
    log_P(is)=sum((log(p(is,:))),'omitnan'); %take log of prod(p) 
end

log_P(log_P==0)=nan;%just making sure

P=double(exp(vpa(log_P)));

P=P.^(1/samplingrate);

P(~mask_ts_allowed)=nan;

P=P/max(P); % highlight the maximum of probability to be 1

if std(P)<3e-12; %remove noise
   P(:)=0;
   P(find(ts==0))=1;
end
if flag_plot
    figure;
   plot(ts,P);
   xlabel('Offset [years]');
   ylabel('P/max(P)');
   title(['Offset probability for window ' num2str([T_min T_max],'%d - %d') ' years b2k']);
end
end