function OccurrenceExample
% Compute the  occurrence of compound dry and hot events based on precipitation and temperature
% Please check the paper below for detail: 
% XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
% Hao Z. et al. (2019). A monitoring and prediction system for compound dry and
% hot events. Environ. Res. Lett. , in press
% XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
 clc;clear all;close all
PD=rand(468,1); % Monthly precipitation (39 year)
TD=rand(468,1); % Monthly temperature(39 year)
% Time scale
ts=3;
% Compute accumulation/average
APD=zeros(size(PD)).*nan;
ATD=zeros(size(TD)).*nan;
for i=1:(length(PD)-ts+1)
    APD(ts+i-1)=mean(PD(i:(i+ts-1)));
    ATD(ts+i-1)=mean(TD(i:(i+ts-1)));
end

% Occurrence computation of each month based on P50/T50
% Occurrence computation of each month based on P30/T70
prcp=30; prct=70;

OCCU=zeros(size(PD)).*nan;

for k=1:12
    clear dp dt
    if k<ts
        dp=APD(k+12:12:end);
        dt=ATD(k+12:12:end);
    else
        dp=APD(k:12:end);
        dt=ATD(k:12:end);
    end
    
    if k<ts
        OCCU(k+12:12:end)=GetDH2(dp,dt,prcp,prct);
    else
        OCCU(k:12:end)=GetDH2(dp,dt,prcp,prct);
    end
end

n=length(OCCU);
figure(1);set(gcf,'position',[100,-100,900,500])


[I1,I]=find(OCCU>0);
stem(I1,OCCU(I1),'filled','Color','r','LineStyle','-.','MarkerFaceColor','red','MarkerEdgeColor','red');
% stem(1:n,OCCU,'filled','Color','r','LineStyle','-.','MarkerFaceColor','red','MarkerEdgeColor','red');
hold on; plot(0:(n+1),ones(n+2,1).*0.5,'k-.')
xlabel('Year','FontSize',12);ylabel('Occurrence','FontSize',12);
ylim([-0.05 1.2])
xk=[1980:3:2018];for i=1:length(xk) xl{i}=num2str(xk(i)); end
for i=1:length(xk) xt(i)=(xk(i)-1980).*12+1; end
set(gca,'xtick',xt,'xticklabel',xl)
title(['Occurrence of compound dry-hot events ','(P',num2str(prcp),'/','T',num2str(prct),')'],'FontSize',12)

 function y=GetDH2(mp,mt,prcp,prct)

n=length(mp);
y=zeros(n,1);

y01=zeros(n,1);
y02=zeros(n,1);

p0=prctile(mp,prcp);
t0=prctile(mt,prct);

y01(mp<=p0)=1;
y02(mt>t0)=1;
 
y=y01.*y02;

