%%
% Diel profile and percetiles calculation
% Written by Milena Ponczek - Nov/2019
%
% This routine reads the mat file generated by the data_validation routine
% and calculates and plots diel profiles for biogenic and aromatic VOCs
% (separately)
% 
%%

load Data_Validation/RioBranco-noBG-data_val_diel.mat

var_median = NaN(24, min(size(vocs))); %inicializa os vetores
var_mean = NaN(24, min(size(vocs)));
var_p25 = NaN(24, min(size(vocs)));
var_p75 = NaN(24, min(size(vocs)));

    for i_diel_hour=1:24

    idx_diel=find(hour(time_PTRMS)==i_diel_hour-1); %Acha indices para dados de cada hora
        
        for i_diel_voc = 1:min(size(vocs))

        var_median(i_diel_hour,i_diel_voc) = nanmedian(vocs(idx_diel, i_diel_voc)); %Calcula a mediana destes valores

        var_mean(i_diel_hour, i_diel_voc) = nanmean(vocs(idx_diel, i_diel_voc)); %Calcula a media destes valores
        
        var_p25(i_diel_hour, i_diel_voc) = prctile(vocs(idx_diel, i_diel_voc)',25);
        var_p75(i_diel_hour, i_diel_voc) = prctile(vocs(idx_diel, i_diel_voc)',75);
        

        end
    end

 
time_diel = [0:23]';
%%
time_diel_LT = NaN(24, 1);
time_zone = 5; %time_zone = difference from GMT

for m=1:24
    if time_diel(m)<time_zone
     time_diel_LT(m) = time_diel(m)-time_zone+24;
    else
     time_diel_LT(m) = time_diel(m)-time_zone;
    end
    
end

%%
% ------------------------------------------------------------------
%                          PLOT SECTION
% ------------------------------------------------------------------
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%  FIGURE 1  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fig1=figure;

[rs,t,d,z,b,r]=suncycle(-9.97,-67.81);

set(fig1,'InvertHardcopy','on');
set(gca,'FontSize', 12, 'LineWidth', 2); 

axis('off');
axes1 = axes('Parent',fig1,...
'XTickLabel',{'0','2','4','6','8','10','12','14','16','18','20','22','24'},...
'XTick',[0 2 4 6 8 10 12 14 16 18 20 22 24]);

% axes1 = axes('Parent',fig1)
hold on
rectangle('Parent',axes1,'Position',[rs(1),-1000,rs(2)-rs(1),200000],'FaceColor',[1 1 0.7],'EdgeColor',[1 1 0.7])
%for LT: substract rs(1) - time_zone
plot(time_diel,var_median(:,[4,6,7,8,9]), 'o', 'LineWidth', 1.5)
title(['Diel Profile Biogenic VOCs'])
legend({'Acetaldehyde','Acetone','Isoprene','MVK+MAC', 'MEK'}, 'Location','northwest')
legend('boxoff')
xlabel('Time (LT)')
ylabel('Mixing Ratio (ppb)')
ylim([0 5])
xlim(axes1,[-0.5 24])


%ylim(axes1,[0 nanmax(nanmax([vocs(idx_st:idx_end,Idx45) vocs(idx_st:idx_end,Idx43) vocs(idx_st:idx_end,Idx59)]))])
box on

set(gcf,'units','inch','position',[0,0,8,4]); %
set(gcf,'PaperPositionMode','auto')

nome_diel_1=[fig_dir, 'Diel_Profile_LT'];

print(fig1,'-dpng',[nome_diel_1,'.png']);
hold off

%%
%comparing boxplot with my own calculate results
figure
boxplot(vocs(:,4),hour(time_PTRMS))
hold on
plot(time_diel+1,var_median(:,4), 'o-')
plot(time_diel+1, var_p25(:,4))
plot(time_diel+1, var_p75(:,4))
plot(time_diel+1, var_mean(:,4),'*-')

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%  FIGURE 2  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%using error bars for the diel profile
fig1=figure;
errorbar(time_diel,var_median(:,4), var_median(:,4)- var_p25(:,4),var_p75(:,4)-var_median(:,4), '.--', 'markersize', 23, 'LineWidth', 1.75)
hold on
%rectangle('Parent',axes1,'Position',[rs(1),-1000,rs(2)-rs(1),200000],'FaceColor',[1 1 0.7],'EdgeColor',[1 1 0.7])
errorbar(time_diel,var_median(:,6), var_median(:,6)- var_p25(:,6),var_p75(:,6)-var_median(:,6), '.--', 'markersize', 23, 'LineWidth', 1.5)
errorbar(time_diel,var_median(:,7), var_median(:,7)- var_p25(:,7),var_p75(:,7)-var_median(:,7),'.--', 'markersize', 23, 'LineWidth', 1.0)
% grid on
legend({'Acetaldehyde','Acetone','Isoprene'},'Location','northwest', 'FontSize', 10)
legend('boxoff')
xlim([-0.5 24])
title(['Diel Profile Biogenic VOCs'])
xlabel('Time (UTC)')
ylabel('Mixing Ratio (ppb)')
set(gca,'Xtick', 0:23)

nome_diel_2=['Diel_Profile_error'];

print(fig1,'-dpng',[nome_diel_2,'.png']);
hold off

%%
figure
%boxplot(vocs(:,4),hour(time_PTRMS),'boxstyle', 'filled','colors' , 'r','outliersize', 0)
boxplot(vocs(:,4),hour(time_PTRMS),'outliersize', 0)


%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%  FIGURE 3  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%using multiples subplots for diferent vocs boxplots
fig1=figure

subplot(3,1,1)
boxplot(vocs(:,4),hour(time_PTRMS),'boxstyle','outline','colors', 'm','outliersize', 1)
ylabel('Acetaldehyde (ppb)')
ylim([0 10])
title(['Diel Profile Biogenic VOCs'])
hold on
subplot(3,1,2)
boxplot(vocs(:,6),hour(time_PTRMS),'boxstyle','outline','colors', 'c', 'outliersize', 1)
%xlabel('Time (UTC)')
ylabel('Acetone (ppb)')
ylim([0 10])
subplot(3,1,3)
boxplot(vocs(:,7),hour(time_PTRMS),'boxstyle','outline','colors', 'g', 'outliersize', 1)
xlabel('Time (UTC)')
ylabel('Isoprene (ppb)')
ylim([0 10])

nome_diel_3=[fig_dir, 'Diel_Profile_Boxplot'];
print(fig1,'-dpng',[nome_diel_3,'.png']);
hold off

%%
%------------------------------------------------------------------------
% ANTHROPOGENIC VOCS
% figures for Benzene, toluene and other anthropogenic VOCs
%------------------------------------------------------------------------
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%  FIGURE 4  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fig1=figure;

[rs,t,d,z,b,r]=suncycle(-9.97,-67.81);

set(fig1,'InvertHardcopy','on');
set(gca,'FontSize', 12, 'LineWidth', 2); 

axis('off');
axes1 = axes('Parent',fig1,...
'XTickLabel',{'0','2','4','6','8','10','12','14','16','18','20','22','24'},...
'XTick',[0 2 4 6 8 10 12 14 16 18 20 22 24]);

% axes1 = axes('Parent',fig1)
hold on
rectangle('Parent',axes1,'Position',[rs(1),-1000,rs(2)-rs(1),200000],'FaceColor',[1 1 0.7],'EdgeColor',[1 1 0.7])

plot(time_diel,var_median(:,[10:12,16]), 'o-', 'LineWidth', 1.5)
title(['Diel Profile Aromatics'])
legend({'Benzene','Toluene','Frag aromatics','Naftalene'}, 'Location','northwest')
legend('boxoff')
xlabel('Time (UTC)')
ylabel('Mixing Ratio (ppb)')
ylim([0 1.5])
xlim(axes1,[0 24])


%ylim(axes1,[0 nanmax(nanmax([vocs(idx_st:idx_end,Idx45) vocs(idx_st:idx_end,Idx43) vocs(idx_st:idx_end,Idx59)]))])
box on

set(gcf,'units','inch','position',[0,0,8,4]); %
set(gcf,'PaperPositionMode','auto')

nome_diel_4=[fig_dir, 'Diel_Profile_aromat'];

print(fig1,'-dpng',[nome_diel_4,'.png']);
hold off

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%  FIGURE 5  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%using error bars for the diel profile

fig1=figure;
errorbar(time_diel,var_median(:,10), var_p25(:,10),var_p75(:,10), '.--', 'markersize', 23, 'LineWidth', 1.75)
hold on
%rectangle('Parent',axes1,'Position',[rs(1),-1000,rs(2)-rs(1),200000],'FaceColor',[1 1 0.7],'EdgeColor',[1 1 0.7])
errorbar(time_diel,var_median(:,11), var_p25(:,11),var_p75(:,11), '.--', 'markersize', 23, 'LineWidth', 1.5)
errorbar(time_diel,var_median(:,12), var_p25(:,12),var_p75(:,12),'.--', 'markersize', 23, 'LineWidth', 1.0)

title(['Diel Profile Aromatics'])
legend({'Benzene','Toluene','Frag aromatics'},'Location','northwest', 'FontSize', 10)
legend('boxoff')
xlim([-0.5 24])
ylim([0 1.5])
xlabel('Time (UTC)')
ylabel('Mixing Ratio (ppb)')
set(gca,'Xtick', 0:23)

nome_diel_5=[fig_dir, 'Diel_Profile_aromat_error'];

print(fig1,'-dpng',[nome_diel_5,'.png']);
hold off

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%  FIGURE 6  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%using multiples subplots for diferent vocs boxplots
fig1 = figure;

subplot(3,1,1)
boxplot(vocs(:,10),hour(time_PTRMS),'boxstyle','outline','colors', 'm','outliersize', 1)
ylabel('Benzene (ppb)')
ylim([0 1.5])
title(['Diel Profile Aromatics'])
hold on
subplot(3,1,2)
boxplot(vocs(:,11),hour(time_PTRMS),'boxstyle','outline','colors', 'c', 'outliersize', 1)
%xlabel('Time (UTC)')
ylabel('Toluene (ppb)')
ylim([0 1.5])
subplot(3,1,3)
boxplot(vocs(:,12),hour(time_PTRMS),'boxstyle','outline','colors', 'g', 'outliersize', 1)
xlabel('Time (UTC)')
ylabel('Frag Aromatics (ppb)')
ylim([0 1.5])

nome_diel_6=[fig_dir, 'Diel_Profile_Boxplot_Aromat'];
print(fig1,'-dpng',[nome_diel_6,'.png']);
hold off

%%
filename_mat_file = strcat('fig/',station,'data_val_diel', '.mat');
save (filename_mat_file);