% Monthly "PC" as median GEP when PAR ~ 1000 and VPD ~ 1
% !! contains all gap-filled GEP. 
% Monthly surface conductance, "Gs_m", as median Gs when PAR ~ 1000 and VPD ~ 1
% !! contains all gap-filled gs, no rain filter. 


NDVIdatenum = datenum(DateTime_NDVI);
[xData, yData] = prepareCurveData( NDVIdatenum, NDVI );
% Set up fittype and options.
ft = fittype( 'smoothingspline' );
opts = fitoptions( 'Method', 'SmoothingSpline' );
opts.SmoothingParam = 1e-04;
% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );
s_NDVI = feval(fitresult,datenum(DateTime_CUP));

start_d = datenum(DateTime_LAI(1));
end_d = datenum(DateTime_LAI(297));
st_LAI = find(datenum(DateTime_CUP) == start_d);
en_LAI = find(datenum(DateTime_CUP) == end_d);
hh_num = datenum(DateTime_CUP(2)) - datenum(DateTime_CUP(1));
xq = start_d:hh_num:end_d;
% smooth LAI
LAIdatenum = datenum(DateTime_LAI);
[xData, yData] = prepareCurveData( LAIdatenum, LAI );
% Set up fittype and options.
ft = fittype( 'smoothingspline' );
excludedPoints = excludedata( xData, yData, 'Indices', [8 129 131 132] );
opts = fitoptions( 'Method', 'SmoothingSpline' );
opts.SmoothingParam = 4.46646226829198e-05;
opts.Exclude = excludedPoints;
% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );
s_LAI = feval(fitresult,xq);
PC = nan(12*3,3);
gs_m = nan(12*3,3);
date_m = datetime;
LAI_m = nan(36,1);
NDVI_m = nan(36,1);
i = 1;
for y = 1:3
    for m = 1:12
        LAI_m(i) = mean(s_LAI(Year_t(st_LAI:en_LAI) == 2013+y & Month_t(st_LAI:en_LAI) == m & Day_t(st_LAI:en_LAI) == 15));
        NDVI_m(i) = mean(s_NDVI(Year_t == 2013+y & Month_t == m & Day_t == 15));
        use = find(Year_t == 2013+y & Month_t == m & PAR > 800 & PAR < 1200 & VPD > 1 & VPD < 1.5 ...
            & qc == 0 & qc_Sc == 0 & AGC_c == 0 & isnan(NEE_c) == 0 & u > 0.2);
        use_gsmax = find(Year_t == 2013+y & Month_t == m & Precip_2daysago < 1 & Precip_1daysago < 0.5 ...
            & Precip_halfdayago < 0.2 & VPD > 0.9 & VPD < 1.3 & qc_h2o_flux == 0 & AGC_c == 0 & u > 0.2 & PAR > 800 & PAR < 1200);
        PC(i,1) = quantile(-GPP(use),0.25);
        PC(i,2) = quantile(-GPP(use),0.5);
        PC(i,3) = quantile(-GPP(use),0.75);
        gs_m(i,1) = quantile(gsmol(use_gsmax),0.25);
        gs_m(i,2) = quantile(gsmol(use_gsmax),0.5);
        gs_m(i,3) = quantile(gsmol(use_gsmax),0.75);
        date_m(i) = datetime(2013+y,m,15);
        i = i+1;
    end
end
% 
% figure; subplot(4,1,1); plot(DateTime_LAI,LAI,'LineStyle','none','Marker','o', ...
%     'MarkerFaceColor',[0.8 0.8 0.8],'MarkerEdgeColor','none','MarkerSize',2); 
% hold on; plot(xq,s_LAI,'LineWidth',2,'Color','k');
% subplot(4,1,2); plot(date_m,PC(:,2),'LineWidth',2,'Color','k');
% hold on; plot(date_m,PC(:,1),'LineWidth',1,'Color',[0.8 0.8 0.8]);
% plot(date_m,PC(:,3),'LineWidth',1,'Color',[0.8 0.8 0.8]);
% subplot(4,1,3); plot(date_m,gs_m(:,2),'LineWidth',2,'Color','k');
% hold on; plot(date_m,gs_m(:,1),'LineWidth',1,'Color',[0.8 0.8 0.8]);
% plot(date_m,gs_m(:,3),'LineWidth',1,'Color',[0.8 0.8 0.8]);
% subplot(4,1,1); ylim([0.6 1.1]);
% xlim([datenum(datetime(2014,01,01)) datenum(datetime(2017,01,01))]);
% set(gca,'FontSize',12); ylabel('LAI (m^2 m^-^2)');
% subplot(4,1,2); ylim([5 20]);
% xlim([datenum(datetime(2014,01,01)) datenum(datetime(2017,01,01))]);
% set(gca,'FontSize',12); ylabel('PC (\mumol m^-^2 s^-^1)');
% subplot(4,1,3); ylim([0.002 0.008]);
% xlim([datenum(datetime(2014,01,01)) datenum(datetime(2017,01,01))]);
% set(gca,'FontSize',12); ylabel('Gs,max (m s^-^1)');
% subplot(4,1,4); plot(DateTime_CUP,SWC_s,'LineWidth',2,'Color','k');
% xlim([datenum(datetime(2014,01,01)) datenum(datetime(2017,01,01))]);
% set(gca,'FontSize',12); ylabel('SWC _0_-_8_c_m');

% linear fit monthly PC / NDVI / gs,max figures

mdl = fitlm(LAI_m,PC(:,2));
figure; subplot(2,3,1);
h = plot(mdl);
Dots_h = h(1); % h(1) is dots handle
Dots_h.Color = [0.6 0.6 0.6];
Dots_h.Marker = 'o';
Line_h = h(2); % h(2) is line handle
Line_h.Color = [0 0 0];
Line_h.LineWidth = 2;
CB_h1 = h(3); % h(3) is confidence bound handle (bottom)
CB_h1.Color = [0.7 0.7 0.7];
CB_h1.LineStyle = '-';
CB_h2 = h(4); % h(3) is confidence bound handle (above)
CB_h2.Color = [0.7 0.7 0.7];
CB_h2.LineStyle = '-';
set(get(gca,'XLabel'),'Interpreter','tex'); % The LinearModel object turns this interpreter off, so things like my_variable won't be displayed with a subscript v.
set(get(gca,'YLabel'),'Interpreter','tex');
xlabel('LAI (m^2 m^-^2)'); ylabel('PC (\mumol m^-^2 s^-^1)');
box off; title(''); 
legend('off');
% h = legend; set(h,'FontSize',10);
% set(h,'Location','northwest');
h = text(0.08,0.95,'PC = 10.9*LAI + 1.8','Units', 'Normalized', 'VerticalAlignment', 'Top'); 
h.FontSize = 10;
h = text(0.08,0.85,'r^2 = 0.24','Units', 'Normalized', 'VerticalAlignment', 'Top');
h.FontSize = 10;
h = text(0.08,0.75,'p = 0.02','Units', 'Normalized', 'VerticalAlignment', 'Top');
h.FontSize = 10;
h = text(0.90,0.98,'a','Units', 'Normalized', 'VerticalAlignment', 'Top');
h.FontSize = 10;
axis([0.6 1.1 4 18]);

mdl = fitlm(NDVI_m,PC(:,2));
subplot(2,3,2);
h = plot(mdl);
Dots_h = h(1); % h(1) is dots handle
Dots_h.Color = [0.6 0.6 0.6];
Dots_h.Marker = 'o';
Line_h = h(2); % h(2) is line handle
Line_h.Color = [0 0 0];
Line_h.LineWidth = 2;
CB_h1 = h(3); % h(3) is confidence bound handle (bottom)
CB_h1.Color = [0.7 0.7 0.7];
CB_h1.LineStyle = '-';
CB_h2 = h(4); % h(3) is confidence bound handle (above)
CB_h2.Color = [0.7 0.7 0.7];
CB_h2.LineStyle = '-';
set(get(gca,'XLabel'),'Interpreter','tex'); % The LinearModel object turns this interpreter off, so things like my_variable won't be displayed with a subscript v.
set(get(gca,'YLabel'),'Interpreter','tex');
xlabel('NDVI'); ylabel(''); set(gca,'yticklabel',[]);
box off; title(''); 
legend('off');
% h = legend; set(h,'FontSize',10);
% set(h,'Location','northwest');
h = text(0.08,0.95,'PC = 2*NDVI + 10','Units', 'Normalized', 'VerticalAlignment', 'Top'); 
h.FontSize = 10;
h = text(0.08,0.85,'r^2 = 0.003','Units', 'Normalized', 'VerticalAlignment', 'Top');
h.FontSize = 10;
h = text(0.08,0.75,'p = 0.8','Units', 'Normalized', 'VerticalAlignment', 'Top');
h.FontSize = 10;
h = text(0.90,0.98,'b','Units', 'Normalized', 'VerticalAlignment', 'Top');
h.FontSize = 10;
%axis([0.6 1.1 4 18]);

mdl = fitlm(gs_m(:,3),PC(:,2));
subplot(2,3,3);
h = plot(mdl);
Dots_h = h(1); % h(1) is dots handle
Dots_h.Color = [0.6 0.6 0.6];
Dots_h.Marker = 'o';
Line_h = h(2); % h(2) is line handle
Line_h.Color = [0 0 0];
Line_h.LineWidth = 2;
CB_h1 = h(3); % h(3) is confidence bound handle (bottom)
CB_h1.Color = [0.7 0.7 0.7];
CB_h1.LineStyle = '-';
CB_h2 = h(4); % h(3) is confidence bound handle (above)
CB_h2.Color = [0.7 0.7 0.7];
CB_h2.LineStyle = '-';
set(get(gca,'XLabel'),'Interpreter','tex'); % The LinearModel object turns this interpreter off, so things like my_variable won't be displayed with a subscript v.
set(get(gca,'YLabel'),'Interpreter','tex');
xlabel('G_s_,_m_a_x (mmol m^-^2 s^-^1)'); ylabel(''); set(gca,'yticklabel',[]);
box off; title(''); 
legend('off');
% h = legend; set(h,'FontSize',10);
% set(h,'Location','northwest');
h = text(0.08,0.95,'PC = 8.8*G_s_,_m_a_x + 9.3','Units', 'Normalized', 'VerticalAlignment', 'Top'); 
h.FontSize = 10;
h = text(0.08,0.85,'r^2 = 0.08','Units', 'Normalized', 'VerticalAlignment', 'Top');
h.FontSize = 10;
h = text(0.08,0.75,'p = 0.13','Units', 'Normalized', 'VerticalAlignment', 'Top');
h.FontSize = 10;
h = text(0.90,0.98,'c','Units', 'Normalized', 'VerticalAlignment', 'Top');
h.FontSize = 10;
%axis([0.6 1.1 4 18]);

mdl = fitlm(LAI_m,gs_m(:,3));
subplot(2,3,4);
h = plot(mdl);
Dots_h = h(1); % h(1) is dots handle
Dots_h.Color = [0.6 0.6 0.6];
Dots_h.Marker = 'o';
Line_h = h(2); % h(2) is line handle
Line_h.Color = [0 0 0];
Line_h.LineWidth = 2;
CB_h1 = h(3); % h(3) is confidence bound handle (bottom)
CB_h1.Color = [0.7 0.7 0.7];
CB_h1.LineStyle = '-';
CB_h2 = h(4); % h(3) is confidence bound handle (above)
CB_h2.Color = [0.7 0.7 0.7];
CB_h2.LineStyle = '-';
set(get(gca,'XLabel'),'Interpreter','tex'); % The LinearModel object turns this interpreter off, so things like my_variable won't be displayed with a subscript v.
set(get(gca,'YLabel'),'Interpreter','tex');
xlabel('LAI (m^2 m^-^2)'); ylabel('G_s_,_m_a_x (mmol m^-^2 s^-^1)'); %set(gca,'yticklabel',[]);
box off; title(''); 
legend('off');
% h = legend; set(h,'FontSize',10);
% set(h,'Location','northwest');
h = text(0.08,0.95,'G_s_,_m_a_x = 0.46*LAI - 0.13','Units', 'Normalized', 'VerticalAlignment', 'Top'); 
h.FontSize = 10;
h = text(0.08,0.85,'r^2 = 0.28','Units', 'Normalized', 'VerticalAlignment', 'Top');
h.FontSize = 10;
h = text(0.08,0.75,'p = 0.004','Units', 'Normalized', 'VerticalAlignment', 'Top');
h.FontSize = 10;
h = text(0.90,0.98,'d','Units', 'Normalized', 'VerticalAlignment', 'Top');
h.FontSize = 10;
%axis([0.6 1.1 4 18]);

mdl = fitlm(NDVI_m,gs_m(:,3));
subplot(2,3,5);
h = plot(mdl);
Dots_h = h(1); % h(1) is dots handle
Dots_h.Color = [0.6 0.6 0.6];
Dots_h.Marker = 'o';
Line_h = h(2); % h(2) is line handle
Line_h.Color = [0 0 0];
Line_h.LineWidth = 2;
CB_h1 = h(3); % h(3) is confidence bound handle (bottom)
CB_h1.Color = [0.7 0.7 0.7];
CB_h1.LineStyle = '-';
CB_h2 = h(4); % h(3) is confidence bound handle (above)
CB_h2.Color = [0.7 0.7 0.7];
CB_h2.LineStyle = '-';
set(get(gca,'XLabel'),'Interpreter','tex'); % The LinearModel object turns this interpreter off, so things like my_variable won't be displayed with a subscript v.
set(get(gca,'YLabel'),'Interpreter','tex');
xlabel('NDVI'); ylabel(''); set(gca,'yticklabel',[]);
box off; title(''); 
legend('off');
% h = legend; set(h,'FontSize',10);
% set(h,'Location','northwest');
h = text(0.08,0.95,'G_s_,_m_a_x = 0.74*NDVI - 0.24','Units', 'Normalized', 'VerticalAlignment', 'Top'); 
h.FontSize = 10;
h = text(0.08,0.85,'r^2 = 0.24','Units', 'Normalized', 'VerticalAlignment', 'Top');
h.FontSize = 10;
h = text(0.08,0.75,'p = 0.002','Units', 'Normalized', 'VerticalAlignment', 'Top');
h.FontSize = 10;
h = text(0.90,0.98,'e','Units', 'Normalized', 'VerticalAlignment', 'Top');
h.FontSize = 10;
%axis([0.6 1.1 4 18]);