% Example illustrating how each of the predictors comprising the
% leaf index and bloom index are generated using the functions
% in ../si_funcs/.

%% Clear workspace and load appropriate data (select6.mat, in this case):
clear
close all
%% Select a station and year
stn_nums=[114442]; % Station #
yrq=47; % index of year (47 = 1896 AD)

% Create vectors of Tmin/Tmax
for i =1:length(stn_nums);
stn_id=num2str(stn_nums(i));
eval(['tmin1year(:,:,i)=convert_temp(USC00' stn_id '.TMIN.data(yrq,:),' char(39) 'C' char(39) ',' char(39) 'F' char(39) ');']);
eval(['tmax1year(:,:,i)=convert_temp(USC00' stn_id '.TMAX.data(yrq,:),' char(39) 'C' char(39) ',' char(39) 'F' char(39) ');']);
eval(['lat(i)=USC00' stn_id '.lat;']);

if i ==1
eval(['stn_time=USC00' stn_id '.time;']);
end

end

% Calculate SI
[LFMTX_1year,BLMTX_1year]=calc_si(tmin1year,tmax1year,lat);

%% Now show each of the individual predictor variables
DAYLEN=eval(['calc_daylen(240,USC00' num2str(stn_nums)  '.lat)']);
BASET=31;
T=nan(24,1);
AGDH=zeros(8,1);
ENUM=1;
SYNTOT=0;
SYNTOT_bl=0;
SYNFLAG=0;
SYNFLAG_bl=0;

for i=2:240;
MN=tmin1year(i);
MX=tmax1year(i);

[GDH,D_temp]=growdh(MX,MN,DAYLEN(i),BASET,T);
GDH_ts(i)=GDH;
if i>8
LAG=(GDH_ts(i-7:i));
[DDE2, DD57, SYNFLAG]=synval(LAG,GDH,ENUM);
DDE2_ts(i)=DDE2;
DD56_ts(i)=DD57;
AGDH(i)=AGDH(i-1)+GDH;

[DDE2, DD57, SYNFLAG_bl]=synval(LAG,GDH,3);
end

SYNTOT(i)=SYNTOT(i-1)+SYNFLAG;
SYNTOT_bl(i)=SYNTOT_bl(i-1)+SYNFLAG_bl;

end

%%
figure(1)
clf

Tmax=tmax1year;
Tmin=tmin1year;
daymax=240;
GDH=GDH_ts;
SYNOP=SYNTOT;
DDE2=DDE2_ts;
DD57=DD56_ts;
mnLFraw=LFMTX_1year(1);
mnBLraw=BLMTX_1year(1);

ax(1)=subplot(511);
p1=plot([Tmax(1:daymax)' Tmin(1:daymax)']);
set(gca,'ylim',[-22 90])
ylabel 'Temp (^oF)'
lg=legend('T_{max}','T_{min}','location','northwest');
legend boxon
set(lg,'pos',[    0.15    0.9    0.1351    0.0562])

% subplot(512)
%     p2=plot(1:daymax,1:daymax);
%     ylabel 'DOY'
%     xlabel 'Time (DOY)'

ax(2)=subplot(512);
p2=plot(1:daymax,GDH);
ylabel 'GDH'

ax(3)=subplot(513);
p3=stairs(1:daymax,SYNOP);
ylabel 'SYNOP'

ax(4)=subplot(514);
p4=plot(1:daymax,DDE2,1:daymax,ones(1,daymax)*637,'--');
ylabel 'DDE2'

ax(5)=subplot(515);
p5=plot(1:daymax,DD57);
ylabel 'DD57'
xlabel 'Time (DOY)'

for i =1:5
set(eval(['p' num2str(i)]),'color','k','linewidth',2);
set(gcf,'currentax',ax(i))
hold on
ylims=get(gca,'ylim');

text(mnLFraw,ylims(2)+ylims(2)*.1,'Leaf Index');

axpos=get(gca,'pos');

set(gca,axset{:},'box','on','pos',[axpos(1)-.005 axpos(2:end)],'xlim',[0 115]);
plot([mnLFraw mnLFraw],ylims,'linestyle','--','color',darkgreen_rgb,'linewidth',2)
if i ==4
text(daymax+1,637,{'Event';'Threshold'},axset{1:4})
end
end

set(p1(1),'color',darkred_rgb);
set(p1(2),'color','b');

h=7.2; w=6.5; res=300; filename='../figs/predictor_fig_lf';printertype='-painters';
set(gcf,'paperposition',[1 1 w h])
PrevFig(1)
print('-depsc2',['-r' num2str(res)],'-cmyk','-loose',[filename '.eps'],printertype)
saveas(gcf,[filename '.fig'],'fig');

%%
figure(2)
clf
Tmax=tmax1year;
Tmin=tmin1year;
GDH=GDH_ts;
SYNOP=SYNTOT_bl;
DDE2=DDE2_ts;
DD57=DD56_ts;
mnLFraw=LFMTX_1year(1);

ax(1)=subplot(311);
p1=plot([Tmax(1:daymax)' Tmin(1:daymax)']);
set(gca,'ylim',[-10 110])
ylabel 'Temp (^oF)'
lg=legend('T_{max}','T_{min}','location','northwest');
legend boxon
set(lg,'pos',[    0.15    0.9    0.1351    0.0562])

% subplot(512)
%     p2=plot(1:daymax,1:daymax);
%     ylabel 'DOY'
%     xlabel 'Time (DOY)'

ax(2)=subplot(312);
p2=plot(1:daymax,GDH);
ylabel 'GDH'

ax(3)=subplot(313);
p3=stairs(1:daymax,AGDH./100);
ylabel 'AGDH \cdot 10^{-2}'
set(gca,'ylim',[0 600]);

for i =1:3
set(eval(['p' num2str(i)]),'color','k','linewidth',2);
set(gcf,'currentax',ax(i))
hold on
ylims=get(gca,'ylim');

text(mnLFraw,ylims(2)+ylims(2)*.1,'Leaf Index');
text(mnBLraw,ylims(2)+ylims(2)*.1,'Bloom Index');

axpos=get(gca,'pos');

set(gca,axset{:},'box','on','pos',[axpos(1)-.005 axpos(2:end)],'xlim',[60 140]);
plot([mnBLraw mnBLraw],ylims,'linestyle','--','color','m','linewidth',2)
plot([mnLFraw mnLFraw],ylims,'linestyle','--','color',darkgreen_rgb,'linewidth',2)
if i ==4
text(daymax+1,2001,{'Event';'Threshold'},axset{1:4})
end
end

set(p1(1),'color',darkred_rgb);
set(p1(2),'color','b');

% Output figure settings:
h=7.2*(3/5); w=6.5; res=300; filename='../figs/predictor_fig_bl';printertype='-painters';
set(gcf,'paperposition',[1 1 w h])
PrevFig(2)
print('-depsc2',['-r' num2str(res)],'-cmyk','-loose',[filename '.eps'],printertype)
saveas(gcf,[filename '.fig'],'fig');