%%
% Data validation filters (for further processing)
% Written by Milena Ponczek - Nov/2019
%
%
% This routine reads the mat file generated by the quick_plot routine and
% perform some data validation to remove points that should not be consider
% for level 3 results
% 
%
%%
clear all

% USER INPUT:  Name of mat file we want to load
load mat-files/RioBranco-noBG.mat;

orig_size=max(size(time_cps));

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Primary Ion criteria
% filter data when PrimIon intensity is <= 7E6 cps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%diagnosis
idx_Prim_criteria = find(Prim <= 7E6);
disp([num2str(100.*(max(size(idx_Prim_criteria)./orig_size)),3),'% of the raw data is out of Prim intensity criteria'])

%% Removing data...
% remove data that doesn't match the criteria - PrimIon intensity is <= 7E6 cps

% Prim_criteria = Prim;
% Prim_criteria(idx_Prim_criteria)=[];
% Prim_criteria_time = time_PTRMS;
% Prim_criteria_time(idx_Prim_criteria)=[];

time_PTRMS(idx_Prim_criteria)=[];
time_cps(idx_Prim_criteria)=[];
cps(idx_Prim_criteria,:) = [];
vocs(idx_Prim_criteria,:)=[];
Prim(idx_Prim_criteria)=[];
m37(idx_Prim_criteria)=[];
m32(idx_Prim_criteria)=[];
m30(idx_Prim_criteria)=[];
ratioCluster = (m37./Prim) *100;

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Water cluster criteria
% filter data when the ratio mz37/Prim >= 10%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%diagnosis
ratioCluster = (m37./Prim) *100;
idx_mz37 = find(ratioCluster >= 10);
disp([num2str(100.*(max(size(idx_mz37)./orig_size)),3),'% of the raw data is out of water cluster criteria'])


%% Removing data...
% data when the ratio mz37/Prim >= 10%

time_PTRMS(idx_mz37)=[];
time_cps(idx_mz37)=[];
cps(idx_mz37,:) = [];
vocs(idx_mz37,:)=[];
Prim(idx_mz37)=[];
m37(idx_mz37)=[];
m32(idx_mz37)=[];
m30(idx_mz37)=[];
ratioCluster = (m37./Prim) *100;
ratioO2 = ((m32)./Prim) *100;
ratioNO = (m30./Prim) *100;
ratioImpurities = ratioO2 + ratioNO;

meanCluster = mean(ratioCluster);
meanImpurities = mean(ratioImpurities);

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Impurity criteria
% filter data when the ratio mz32/Prim >= 10% *
% *Obs: this criteria can be changed to 5% depending on the case!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%diagnosis
idx_mz32 = find(ratioO2 >= 4);
disp([num2str(100.*(max(size(idx_mz32)./orig_size)),3),'% of the raw data is out of water cluster criteria'])

%% Removing data...
% data when the ratio mz32/Prim >= 10%

time_PTRMS(idx_mz32)=[];
time_cps(idx_mz32)=[];
cps(idx_mz32,:) = [];
vocs(idx_mz32,:)=[];
Prim(idx_mz32)=[];
m37(idx_mz32)=[];
m32(idx_mz32)=[];
m30(idx_mz32)=[];

%%
%removing specific points or days that  have weird results

time_cps(time_cps>datenum(2018,09,20,0,0,0)&time_cps<datenum(2018,09,26,00,0,0))=[];
time_PTRMS(time_cps>datenum(2018,09,20,0,0,0)&time_cps<datenum(2018,09,26,00,0,0))=[];
cps((time_cps>datenum(2018,09,20,0,0,0)&time_cps<datenum(2018,09,26,00,0,0)),:) = [];
vocs((time_cps>datenum(2018,09,20,0,0,0)&time_cps<datenum(2018,09,26,00,0,0)),:)=[];
Prim(time_cps>datenum(2018,09,20,0,0,0)&time_cps<datenum(2018,09,26,00,0,0))=[];
m37(time_cps>datenum(2018,09,20,0,0,0)&time_cps<datenum(2018,09,26,00,0,0))=[];
m32(time_cps>datenum(2018,09,20,0,0,0)&time_cps<datenum(2018,09,26,00,0,0))=[];
m30(time_cps>datenum(2018,09,20,0,0,0)&time_cps<datenum(2018,09,26,00,0,0))=[];

% time_cps(time_cps>datenum(2018,10,20,03,11,0)&time_cps<datenum(2018,10,21,16,33,0))=[];
% time_PTRMS(time_cps>datenum(2018,10,20,03,11,0)&time_cps<datenum(2018,10,21,16,33,0))=[];
% cps((time_cps>datenum(2018,10,20,03,11,0)&time_cps<datenum(2018,10,21,16,33,0)),:) = [];
% vocs((time_cps>datenum(2018,10,20,03,11,0)&time_cps<datenum(2018,10,21,16,33,0)),:)=[];
% Prim(time_cps>datenum(2018,10,20,03,11,0)&time_cps<datenum(2018,10,21,16,33,0))=[];
% m37(time_cps>datenum(2018,10,20,03,11,0)&time_cps<datenum(2018,10,21,16,33,0))=[];
% m32(time_cps>datenum(2018,10,20,03,11,0)&time_cps<datenum(2018,10,21,16,33,0))=[];
% m30(time_cps>datenum(2018,10,20,03,11,0)&time_cps<datenum(2018,10,21,16,33,0))=[];

%%
%diagnostic very hight voc values
% idx_outlier =  [];

[max_val,idx_time_max_val]= max(vocs);
%plot(day(time_PTRMS(idx_time_max_val)), max_val, 'o')

%não ta dando certo com esse 'for' assim...
% for i=1:max(size(header))
%      idx_outlier = find(vocs(:,i) > 10);
% end

%%  
% PLOT SECTION
%--------------------------------------------------
%plot mz21 corrected (Prim) and mz37
%--------------------------------------------------
fig1=figure;
plot(time_cps,cps(:,1),'m') %mz21 raw
hold on
plot(time_cps,(cps(:,1)*500),'r') %mz21 x 500 -> primary ion
plot(time_cps,cps(:,4),'b')       %mz37 -> water cluster
title(['Prim and Water Cluster'])
legend('mz21','Primary',' H_2O Cluster')
xlabel('Date')
ylabel('Intensity (cps)')

%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
dynamicDateTicks([], [], 'dd/mm');

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

%nome=['../RioBranco/fig/PTRMS_' mydir '_Time_series'];
nome2=[fig_dir, 'Time_series_Primary'];

print(fig1,'-dpng',[nome2,'.png']);

%%
%--------------------------------------------------
%plot water cluster and impurities analysis
%--------------------------------------------------
ratioCluster = (m37./Prim) *100;
ratioO2 = ((m32)./Prim) *100;
ratioNO = (m30./Prim) *100;
ratioImpurities = ratioO2 + ratioNO;
meanCluster = mean(ratioCluster);
meanImpurities = mean(ratioImpurities);

fig1=figure;
plot(time_PTRMS,ratioCluster,'c')
hold on
plot(time_PTRMS,ratioO2,'b')
plot(time_PTRMS,ratioNO,'r')
plot(time_PTRMS,ratioImpurities,'k')

title(['H_2O Cluster and Impurities'])
legend('H_2O Cluster','O_2^+','NO^+', 'Total Impurity')

xlabel('Date')
ylabel('Ratio')

%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
dynamicDateTicks([], [], 'dd/mm');

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

nome3 = [fig_dir, 'Time_series_water_cluster_ratio'];

print(fig1,'-dpng',[nome3 '.png']);


%%
%--------------------------------------------------------------------
% Plot full time series
%--------------------------------------------------------------------

disp('Plot: time series')

if isunix 
    fig1 = figure('visible','off');
else
    fig1=figure;
end;

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

plot(time_PTRMS,vocs(:,Idx45),'m')
hold on
%plot(time_PTRMS,vocs(:,Idx42),'b')
plot(time_PTRMS,vocs(:,Idx59),'c')
plot(time_PTRMS,vocs(:,Idx69),'g')

plot(time_PTRMS,vocs(:,Idx71),'r')
plot(time_PTRMS,vocs(:,Idx79),'k')
plot(time_PTRMS,vocs(:,Idx93),'y')
title([title_PTRMS,station])
legend('Acetaldehyde','Acetone','Isoprene','MVK+MACR','Benzene', 'Toluene')
xlabel('Date')
ylabel(label_PTRMS)

%ylim(axes1,[0 nanmax(nanmax([vocs(idx_st:idx_end,Idx45) vocs(idx_st:idx_end,Idx69) vocs(idx_st:idx_end,Idx59)]))])
ylim([0 50])
box on
dynamicDateTicks([], [], 'dd/mm');

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

%nome=['../RioBranco/fig/PTRMS_' mydir '_Time_series'];
nome=[fig_dir, 'Time_series'];

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

%%
filename_mat_file = strcat('mat-files/',station,'-data_val','.mat');
save (filename_mat_file);

%%

