%% A matlab code to deconvolve the weathering rate from the delta18O of detrital zircon
clear;
tic

%% import data from excell file, run it for the first time
%[num,txt,raw]=xlsread('Datafile_2020.xlsx','Table S1_Data Compilation');
%zirconX(:,1)=num(:,7).*1000;
%zirconX(:,2)=num(:,5);
%zirconX(:,3)=num(:,1);
%idX=find(isnan(zirconX(:,2))==0);
%zirconNJU=zirconX(idX,:);
%save zirconNJU zirconNJU;

%% loading data
load zirconNJU;
zircon=zirconNJU;

%% Setting for model parameters
AR=200;                 %age resolution of the calculation (Ma)
Delta0=5.3;             %mantle delta18O
nMont=4;                %number of Mont Carlo sampling 4*datalenth
ninterp=2;              %number of interplot

%% global paramenters
nsample=fix(max(zircon(:,1))/AR)+1;     %number of sample in each simulation
ZirconL=length(zircon);
nSE=sqrt(nsample)/sqrt(ZirconL);

%% generating the Monte Carlo index, run it for the first time
%[bbb idm]=sort(random('Normal',1,1,nMont*ZirconL,nsample));
%idm=mod(idm,ZirconL)+1;
load idm;

%% Changing rate of delta18O of all zicorn
i=find((zircon(:,1)-zircon(:,3))<=AR);          %zircon with small delta age
AA(i,1)=(zircon(i,2)-Delta0)/AR;
i=find((zircon(:,1)-zircon(:,3))>AR);
AA(i,1)=(zircon(i,2)-Delta0)./(zircon(i,1)-zircon(i,3));

%% calculating the average alteraton rate
nOUT=1;
for j=1:ninterp
    clear BB;
    clear xrand;
    clear xtemp;
    %shifting age
    zircon(:,1)=zircon(:,1)-AR/ninterp*(j-1);
    zircon(:,3)=zircon(:,3)-AR/ninterp*(j-1);
    
    %correct for <0 ages
    i=find(zircon(:,3)<0);
    zircon(i,3)=0;
    
    i=find(zircon(:,1)<0);
    zircon(i,1)=0;
    
    %The starting and ending cell of matrix
    zircon(:,4)=fix(zircon(:,3)/AR)+1;
    zircon(:,5)=fix(zircon(:,1)/AR)+1;
    
 
    
    %age difference
    zircon(:,6)=zircon(:,1)-zircon(:,3);
    
    
    %U-Pb age larger than Hf age
    i=find(zircon(:,4)>zircon(:,5));
    for ii=1:length(i)
        BB(i(ii),zircon(i(ii),4))=1;
    end
    
    %age range within the same cell
    i=find(zircon(:,4)==zircon(:,5));
    for ii=1:length(i)
        BB(i(ii),zircon(i(ii),4))=1;
    end
    
    %age crosses two cells
    i=find((zircon(:,5)-zircon(:,4))==1);
    for ii=1:length(i)
        BB(i(ii),zircon(i(ii),4))=(zircon(i(ii),4)*AR-zircon(i(ii),3))/zircon(i(ii),6);
        BB(i(ii),zircon(i(ii),5))=(zircon(i(ii),1)-zircon(i(ii),4)*AR)/zircon(i(ii),6);
    end
    
    %age crosses more than two cells
    i=find((zircon(:,5)-zircon(:,4))>=2);
    for ii=1:length(i)
        BB(i(ii),zircon(i(ii),4))=(zircon(i(ii),4)*AR-zircon(i(ii),3))/zircon(i(ii),6);
        BB(i(ii),zircon(i(ii),5))=(zircon(i(ii),1)-(zircon(i(ii),5)-1)*AR)/zircon(i(ii),6);
        BB(i(ii),[zircon(i(ii),4)+1:zircon(i(ii),5)-1])=AR/zircon(i(ii),6);
    end
    
    %Mont Carlo method
    for ii=1:length(idm)
        BBB=BB(idm(ii,:),:);
        xtemp=lsqnonneg(BBB,AA(idm(ii,:)));     %positive solution for each random sampling
        xrand(ii,:)=xtemp';
        xrand(ii,find(sum(BBB)==0))=-1;
    end
    
    %Calculating average and SE of the simulations
    Age=((1:length(BB(1,:)))*AR-AR/2+(j-1)*AR/ninterp)/1000;
    for ii=1:length(BB(1,:))
        XOUT(nOUT,1)=Age(ii);
        idf=find(xrand(:,ii)>=0);
        XOUT(nOUT,2)=mean(xrand(idf,ii));
        XOUT(nOUT,3)=std(xrand(idf,ii))*nSE;
        nOUT=nOUT+1;
        %plot(Age(ii),XOUT(1,ii),'r.','markersize',50);
        %plot(Age(ii)*[1 1],XOUT(1,ii)+XOUT(2,ii)*[1 -1]*nSE,'r-','linewidth',3);
    end
    toc
    j/ninterp
end


%% data plot

subplot(2,2,1);hold on;box;xlim([0 4.5]);ylim([0 8e-3]);
[temp id]=sort(XOUT(:,1));
plot(XOUT(id,1),XOUT(id,2),'r.-','linewidth',3,'markersize',25);hold on;
plot(XOUT(id,1),XOUT(id,2)+XOUT(id,3),'r-');
plot(XOUT(id,1),XOUT(id,2)-XOUT(id,3),'r-');

%% save data
save XOUT XOUT;
%%