% Fig 4
clear
close all

%% load
load('CL_TCCIP_c0_end_prec_transform_raw.mat');
Cevap_e = [nan(7,1);Cevap(1:end-7)];
Ctrans_e = [nan(7,1);Ctrans(1:end-7)];
CH2O_e = [nan(7,1);CH2O(1:end-7)];
clear Swnet Lwnet Swdw Swup U10 WIND QDRIP QINTR RSSUN RSSHA
clear Cevap Ctrans Gevap CH2O LE QBOT Q2M RH2M Prec SH TBOT T2M PBOT photo TV TG

load('CL_TCCIP_ctl_transform_raw.mat');
Cevap_c = [nan(7,1);Cevap(1:end-7)];
Ctrans_c = [nan(7,1);Ctrans(1:end-7)];
CH2O_c = [nan(7,1);CH2O(1:end-7)];
clear Swnet Lwnet Swdw Swup U10 WIND QDRIP QINTR RSSUN RSSHA
clear Cevap Ctrans Gevap CH2O LE QBOT Q2M RH2M Prec SH TBOT T2M PBOT photo TV TG

%% annual diurnal cycle (calculate & plot)
day_data = 24;
dCH2O = CH2O_e-CH2O_c; dCevap = Cevap_e-Cevap_c; dCtrans = Ctrans_e-Ctrans_c;
ano_CH2O = nan(1,24); ano_Cevap = nan(1,24); ano_Ctrans = nan(1,24);
for i = 1:day_data
    ano_CH2O(i) = mean(dCH2O(i:day_data:end),'omitnan');
    ano_Cevap(i) = mean(dCevap(i:day_data:end),'omitnan');
    ano_Ctrans(i) = mean(dCtrans(i:day_data:end),'omitnan');
end
load('Fig_4_6_7_8_9_10_11_quantile.mat')
time = 0:1:23;

%% plot
% plot
hfig = figure('units','inches','position',[0 0 16 5.5]);
% subplot 1
sf1 = subplot(1,2,1);
plot(time,ano_CH2O,'linewidth',3,'MarkerSize',5,'Color',[0 0.45 0.74])
hold on
grid on
upper_C = dCH2O_Q3(4,:);
lower_C = dCH2O_Q1(4,:);
p_C = patch([time fliplr(time)],[upper_C,fliplr(lower_C)],[0 0.45 0.74]);
p_C.EdgeColor = 'none';
p_C.FaceAlpha = 0.2; % transparency of patch
plot([0 23],[0 0],'k--','linewidth',2)
xlim([0 23])
ylim([-0.4 0.4])
set(sf1,'xTick',(0:2:23),'fontsize',20);
set(sf1,'yTick',[-0.4:0.2:0.4],'fontsize',20);
xlabel('hour','fontsize',20);
ylabel('(mm)','fontsize',20);
title({'(a) Diurnal Cycle of Canopy Water';'Anomaly'},'fontsize',20,'Interpreter','none');
% subplot 2
sf2 = subplot(1,2,2);
plot(0:1:23,ano_Cevap,'linewidth',3,'MarkerSize',5,'Color','b')
hold on
grid on
plot(0:1:23,ano_Ctrans,'linewidth',3,'MarkerSize',5,'Color',[0 0.5 0])
upper_Ce = dCevap_Q3(4,:);
lower_Ce = dCevap_Q1(4,:);
p_Ce = patch([time fliplr(time)],[upper_Ce,fliplr(lower_Ce)],[0 0.45 0.74]);
p_Ce.EdgeColor = 'none';
p_Ce.FaceAlpha = 0.2; % transparency of patch
upper_Ct = dCtrans_Q3(4,:);
lower_Ct = dCtrans_Q1(4,:);
p_Ct = patch([time fliplr(time)],[upper_Ct,fliplr(lower_Ct)],[0 0.5 0]);
p_Ct.EdgeColor = 'none';
p_Ct.FaceAlpha = 0.3; % transparency of patch
plot([0 23],[0 0],'k--','linewidth',2)
xlim([0 23])
ylim([-40 20])
set(sf2,'xTick',(0:2:23),'fontsize',20);
set(sf2,'yTick',[-40:20:20],'fontsize',20);
xlabel('hour','fontsize',20);
ylabel('(W/m^2)','fontsize',20);
legend({'Canopy E','Transpir.'},'location','NorthEast','fontsize',20)
title({'(b) Diurnal Cycle of Canopy Evapotranspiration';'Anomaly'},'fontsize',20,'Interpreter','none');
set(gcf,'PaperUnits','inches');
set(gcf,'PaperPosition',[0 0 16 5.5])
set(gcf,'PaperSize',[16 6])
img = gcf; 
img.PaperUnits = 'inches';  
img.PaperPosition = [0 0 16 6]; 
img.Units = 'inches'; 
img.PaperSize=[16 6]; 
print(img, '-dpdf', '-r600', './Fig_4.pdf')

