% Ozone monitor (Thermo 49i) data processing 
% Experimento FSP2023 (Faculdade de Saúde Pública da USP)
%
% Function used: fmedia_qualquer_coisa2; frepeated; export_LFA_GoAmazon_jan2015
%
% Outputs: O3_Level0 (mat_time(utc);y;m;d;h;m;s; o3; hio3; cellai; cellbi; bncht; lmpt; o3lt; flowa; flowb; pres; flag(0/1))--> all raw data, without filters
%          O3_Level1_5min (dd-mmm-yyyy HH:MM; sec_since_1970; mean_o3; std_o3; median_o3; p25_o3; p75_o3; max_o3; min_o3; num_o3)
%          O3_Level1_30min
%
% Time stamps in the output files refer to the start of the period (note that Thermo gas monitors time stamps refer to the end of the period) 
% Criteria for Level1: (Keep records when) Flow A and B >= 0.55 LPM; 
%                       Intensity A and B >= 75000 Hz (and the subsequent point)
%                      -50<ConcO3<1000 ppb
%                      Remove calibration periods according to the fiedl logbook.  
% Criteria for Level2: correct the baseline when needed
%                      
% Obs: trace gas Thermo monitors report mixing ratios normalized to 1013.25 hPa and 0 oC
%
% Flags (High/Low Alarms): Intensity B | Intensity A | Flow B | Flow A | Pressure | Ozone Lamp Temp | Bench Lamp Temp | Bech Temp
% Melhorar: identificar flag ----L--- (low intensity) e -----L-- (low flow)
%
% ---------------- Parametros -----------------------
user='luvar';%'LRizzo'; % 'luvar'
caminho=['C:\Users\',user,'\Dropbox\A\FSP2023-2024\raw\o3\'];
caminho1=['C:\Users\',user,'\Dropbox\A\FSP2023-2024\processados\'];
caminho2=['C:\Users\',user,'\Dropbox\A\FSP2023-2024\processados\plots\'];
fuso=0; % fuso horario para calcular hora local
tit='FSP2023-2024';
filtflag=true; % logical para filtrar dados com flag
filtmanut=false; % logical para filtrar dados em períodos de manutenção ou com problemas
basicplots=true;
lev1=false;
lev2=false;
inicio=datenum(2023,9,11); % início do experimento
final=datenum(2024,9,20); % final do experimento
%
% inicio e termino de periodos de manutencao (horario em utc): removidos das analises
if(filtmanut)   
    manut(1,1:2)=[datenum(2023,09,20,18,10,0) datenum(2023,09,20,18,37,0)];   
    manut(2,1:2)=[datenum(2023,09,29,20,05,0) datenum(2023,09,29,20,35,0)];
    manut(3,1:2)=[datenum(2023,10,20,18,20,0) datenum(2023,10,20,18,38,0)];
    manut(4,1:2)=[datenum(2023,10,31,11,05,0) datenum(2023,10,31,11,35,0)];
    manut(5,1:2)=[datenum(2023,11,16,11,30,0) datenum(2023,11,16,11,58,0)];
    manut(6,1:2)=[datenum(2023,11,28,18,35,0) datenum(2023,11,28,19,15,0)];
    manut(7,1:2)=[datenum(2023,12,21,19,35,0) datenum(2023,12,21,20,08,0)];
    manut(8,1:2)=[datenum(2024,01,10,20,15,0) datenum(2024,01,10,20,42,0)];
    manut(9,1:2)=[datenum(2024,01,31,19,15,0) datenum(2024,01,31,20,05,0)];
    manut(10,1:2)=[datenum(2024,08,14,10,55,0) datenum(2024,08,14,11,32,0)];
    manut(11,1:2)=[datenum(2024,09,13,12,55,0) datenum(2024,09,13,13,33,0)];
    manut(:,3)=0;
    manut(:,4:5)=0;
end
%
% inicio e termino de periodo de zoom para fazer grafico
%zoom=[datenum(2023,11,01,0,0,0) datenum(2023,11,20,0,0,0)];
zoom1=[datenum(2023,09,11,0,0,0) datenum(2024,02,20,0,0,0)];
zoom2=[datenum(2024,05,15,0,0,0) datenum(2024,09,30,0,0,0)];
% ---------------- Fim dos parametros --------------------
%
% ---------------- Correção no horário  --------------------
% Obs: mudaram o instrumento para horário local em 18:28  07-16-24 utc
% Retornei para utc em 09-24-24
% Separei os arquivos em utc e lt, e converti o resultado final para utc.
% --------------------------------------------------
%
data=zeros(1,18);
flag='--------';
files=dir(caminho);
disp('-----Reading raw data files')
for ifile=3:size(files,1)
    nome=getfield(files, {ifile,1}, 'name');
    disp(nome)
    fid=fopen(strcat(caminho,nome));
    textscan(fid,'%s',1);  
    textscan(fid,'%s %s %s %s %s %s %s %s %s',1);  
    textscan(fid,'%s',1);   
    textscan(fid,'%s %s %s %s %s %s %s',1);  
    textscan(fid,'%s',1);  
    label=textscan(fid,'%s %s %s %s %s %s %s %s %s %s %s %s %s',1);   
    C=textscan(fid,'%s %s %s %f %f %f %f %f %f %f %f %f %f');  
    for icol=1:size(C,2)-3
        aux(:,icol)=C{1,icol+3};
    end
    aux1=strcat((C{1,2}(:,1)),(C{1,1}(:,1)));
    vtempo=datevec(aux1,'mm-dd-yyHH:MM'); % date vector
    mattime=datenum(vtempo); % matlab time
    aux2=find(strcmpi('--------',C{1,3})==0); % look for flaged data
    nflag=zeros(length(mattime),1);
    nflag(aux2)=1;
    if (~isempty(regexp(nome,'_lt'))) % se os dados estiverem em local time, converte para utc
        mattime = mattime+datenum(0,0,0,3,0,0);
        vtempo = datevec(mattime);
    end
    data1=cat(2,mattime,vtempo,aux,nflag);
    data=cat(1,data,data1);
    flag=cat(1,flag,C{1,3});
    clear aux aux1 aux2 vtempo doy nflag C data1 nome icol mattime
    fclose(fid);
end
outdata=sortrows(data,1); % sort data based on day of year
outdata(1,:)=[];
clear data
%aux={data(:,1) flag}; % falta implementar flags --> dificil manter outdata
%e flags na mesma base temporal; talvez seja melhor transformar os diferentes flags em
%um codigo numerico
%
% Excluindo dados anteriores ao início do experimento
outdata(outdata(:,1)<inicio,:)=[];
%
disp('----Checking for repeated data')
[newdata,erro]=frepeated(outdata,8);
if (strcmp(erro,'none')~=1) 
    disp(erro)
else
    outdata=newdata;
end
clear aux erro 
%
% Remove maintanance periods from dataset
outdataorig=outdata;
if(filtmanut)
    for i=1:size(manut,1)
        aux=find(outdata(:,1)>=manut(i,1) & outdata(:,1)<=manut(i,2));
        if(manut(i,3)==0)            
            outdata(aux,:)=[]; % remove maintanance periods from dataset
        else
            % detrend baseline (ozone_detrend.m)
            if(manut(i,4)==0)
                outdata(aux,8)=outdata(aux,8)+(-polyval(manut(i,4:5),outdata(aux,1)));
            else
                outdata(aux,8)=outdata(aux,8)+(polyval(manut(i,4:5),outdata(aux(1),1))-polyval(manut(i,4:5),outdata(aux,1)));
            end
        end        
    end
 end
%
% Remove flagged periods from dataset
% Optei por não remover os pontos com flag em jan-fev-2024.
% Os flags são por causa de temperatura alta (ar condicionado quebrado).
% Mas os dados parecem ok.
if(filtflag)
    outdata(outdata(:,18)==1 & outdata(:,2)<2024,:)=[];
end
%
figure,plot(outdataorig(:,1),outdataorig(:,8),'r.'),hold on,grid on
plot(outdata(:,1),outdata(:,8),'b.')
plot(outdataorig(:,1),outdataorig(:,18)*100,'k*')
dynamicDateTicks
clear aux
%
% Corrigindo o time stamp do monitor para reportar o início do período
% O sampling time foi constante, a cada 1 min 
% figure,plot(outdata(2:length(outdata),1),diff(outdata(:,1))*24*60,'-b.'),grid on, dynamicDateTicks
deltat=datenum(0,0,0,0,1,0);
outdata(:,1)=outdata(:,1)-deltat;
%
% Saving output files (level0 and clean)
save([caminho1,'level0\',tit,'_ozone_level0.mat'],'outdataorig','-mat')
if(filtflag || filtmanut)
    save([caminho1,tit,'_ozone_clean.mat'],'outdata','-mat')    
end
%
% -------------------------------------------------------------------------
%   Basic Plots
if(basicplots)

    aux1=find(outdata(:,1)>=zoom1(1) & outdata(:,1)<=zoom1(2));
    aux2=find(outdata(:,1)>=zoom2(1) & outdata(:,1)<=zoom2(2));

    figure
    subplot(2,3,1),plot(outdata(aux1,1),outdata(aux1,10),'b.'), dynamicDateTicks
    ylabel('Intensity Cell A (Hz)','FontSize',12),grid on, axis tight
    subplot(2,3,4),plot(outdata(aux1,1),outdata(aux1,11),'c.'), dynamicDateTicks
    ylabel('Intensity Cell B (Hz)','FontSize',12),grid on, axis tight
    subplot(2,3,2),plot(outdata(aux1,1),outdata(aux1,15),'b.'), dynamicDateTicks
    title(tit,'FontSize',12)
    ylabel('Flow Cell A','FontSize',12),grid on, axis tight
    subplot(2,3,5),plot(outdata(aux1,1),outdata(aux1,16),'m.'), dynamicDateTicks
    ylabel('Flow Cell B','FontSize',12),grid on, axis tight
    subplot(2,3,3),plot(outdata(aux1,1),outdata(aux1,17),'g.'), dynamicDateTicks
    ylabel('Pressure','FontSize',12),grid on, axis tight
    subplot(2,3,6),plot(outdata(aux1,1),outdata(aux1,18),'r.'), axis tight
    dynamicDateTicks, ylabel('Flag','FontSize',12),grid on
    saveas(gcf,strcat(caminho1,'level0\',tit,'set2023-fev2024_O3_diag'),'jpg')
    %
    figure
    subplot(2,3,1),plot(outdata(aux2,1),outdata(aux2,10),'b.'), dynamicDateTicks
    ylabel('Intensity Cell A (Hz)','FontSize',12),grid on, axis tight
    subplot(2,3,4),plot(outdata(aux2,1),outdata(aux2,11),'c.'), dynamicDateTicks
    ylabel('Intensity Cell B (Hz)','FontSize',12),grid on, axis tight
    subplot(2,3,2),plot(outdata(aux2,1),outdata(aux2,15),'b.'), dynamicDateTicks
    title(tit,'FontSize',12)
    ylabel('Flow Cell A','FontSize',12),grid on, axis tight
    subplot(2,3,5),plot(outdata(aux2,1),outdata(aux2,16),'m.'), dynamicDateTicks
    ylabel('Flow Cell B','FontSize',12),grid on, axis tight
    subplot(2,3,3),plot(outdata(aux2,1),outdata(aux2,17),'g.'), dynamicDateTicks
    ylabel('Pressure','FontSize',12),grid on, axis tight
    subplot(2,3,6),plot(outdata(aux2,1),outdata(aux2,18),'r.'), axis tight
    dynamicDateTicks, ylabel('Flag','FontSize',12),grid on
    saveas(gcf,strcat(caminho1,'level0\',tit,'jul-set2024_O3_diag'),'jpg')
    %
    figure
    plot(outdata(:,1),outdata(:,12),'r.')
    dynamicDateTicks
    ylabel('Bench Temp (oC)','FontSize',12),grid on, axis tight
    saveas(gcf,strcat(caminho1,'level0\',tit,'_O3_temp'),'jpg')
    %
    figure
    subplot(2,1,1),plot(outdata(aux1,1)+(fuso/24),outdata(aux1,8),'b.')
    title(tit,'FontSize',14), ylabel('O_3 (ppb)','FontSize',14), xlabel('UTC Time')
    grid on, axis tight
    dynamicDateTicks
    subplot(2,1,2),plot(mod(outdata(aux1,1)+(fuso/24),1),outdata(aux1,8),'b.')
    ylabel('O_3 (ppb)','FontSize',14),grid on,xlabel('UTC Time','FontSize',14), axis tight
    set(gca,'XTick',[0.125 0.25 0.375 0.5 0.625 0.75 0.875 1.0])
    set(gca,'XTickLabel',{'3:00';'6:00';'9:00';'12:00';'15:00';'18:00';'21:00';'0:00'})
    %hgsave(gcf,strcat(caminho1,'O3_conc'))
    saveas(gcf,strcat(caminho1,'\final\plots\',tit,'set2023-fev2024_O3_conc'),'jpg')
    %
    figure
    subplot(2,1,1),plot(outdata(aux2,1)+(fuso/24),outdata(aux2,8),'b.')
    title(tit,'FontSize',14), ylabel('O_3 (ppb)','FontSize',14), xlabel('UTC Time')
    grid on, axis tight
    dynamicDateTicks
    subplot(2,1,2),plot(mod(outdata(aux2,1)+(fuso/24),1),outdata(aux2,8),'b.')
    ylabel('O_3 (ppb)','FontSize',14),grid on,xlabel('UTC Time','FontSize',14), axis tight
    set(gca,'XTick',[0.125 0.25 0.375 0.5 0.625 0.75 0.875 1.0])
    set(gca,'XTickLabel',{'3:00';'6:00';'9:00';'12:00';'15:00';'18:00';'21:00';'0:00'})
    %hgsave(gcf,strcat(caminho1,'O3_conc'))
    saveas(gcf,strcat(caminho1,'\final\plots\',tit,'jul-set2024_O3_conc'),'jpg')
    %
    %
    % Ciclo diurno médio
    first = 0;   % first averaging time
    last = 23.5/24;  % last averaging time
    delta = datenum(0,0,0,0,30,0);     % 30min
    tol=0;
    
    mdata=cat(2,mod(outdata(aux1,1),1),outdata(aux1,8));
    stat=fmedia_qualquer_coisa2(mdata,first,last,delta,tol);
    figure
    errorbar(stat(:,1),stat(:,2),stat(:,3),'b.')
    dynamicDateTicks, set(gca,'XLim',[0 1])
    ylabel('O_3 (ppb)','FontSize',14),grid on,xlabel('UTC Time','FontSize',14)
    title('Mean Diurnal Cycle set2023-fev2024','FontSize',14)
    saveas(gcf,strcat(caminho1,'\final\plots\',tit,'set2023-fev2024_O3_dcycle_level0'),'jpg')    

    mdata=cat(2,mod(outdata(aux2,1),1),outdata(aux2,8));
    stat=fmedia_qualquer_coisa2(mdata,first,last,delta,tol);
    figure
    errorbar(stat(:,1),stat(:,2),stat(:,3),'b.')
    dynamicDateTicks, set(gca,'XLim',[0 1])
    ylabel('O_3 (ppb)','FontSize',14),grid on,xlabel('UTC Time','FontSize',14)
    title('Mean Diurnal Cycle jul-set2024','FontSize',14)
    saveas(gcf,strcat(caminho1,'\final\plots\',tit,'jul-set2024_O3_dcycle_level0'),'jpg')    
end
%
% -------------------------------------------------------------------------
% Level 1
% Keep records when: 
% Flow A and B >= 0.55LPM; Intens A and B >= 75000Hz (and subsequent point); -50<ConcO3<1000 ppb
% Remove calibration periods according to the field logbook.  
%
if(lev1)    
    %
    origsize=size(outdataorig,1);
    aux=find(outdata(:,15)<0.5 | outdata(:,16)<0.5 | outdata(:,10)<60000 | outdata(:,11)<60000);
    aux1=aux(diff(aux)>1)+1; % subsequent point
    aux2=find(outdata(:,8)<=-50 | outdata(:,8)>=1000);
    rem=cat(1,aux,aux1,aux2);
    outdata(rem,:)=[];
    clear aux aux1 aux2 rem
    figure,plot(outdataorig(:,1),outdataorig(:,8),'r.'),grid on, hold on
    plot(outdata(:,1),outdata(:,8),'b.'),dynamicDateTicks
    %
    % Save Level 1 data
    Instrument='Thermo';
    Model='49i';
    serial_num='0830332761';
    Site=tit;
    Val_level='1';
    Proc_vers='1.0';
    labels{1}='O3';
    units{1}='(ppb)';
    export_LFA_GoAmazon_jan2015(labels,units,Instrument,Model,serial_num,Site,Val_level,Proc_vers,origsize,outdata(:,1),outdata(:,8))
end
%
%
% -------------------------------------------------------------------------
% Level 2
% Criteria for Level2: correct the baseline when needed
if(lv2)
    % Lê os arquivos Level 1 que já salvou
    aux=importdata([caminho1,'FSP2023_Thermo_49i_Level1_5min.csv']);
    level1_5min=aux.data(:,1:2);
    clear aux 
    level1_5min(:,1)=level1_5min(:,1)/(24*60*60)+datenum(1970,1,1); % unix time to matlab time
    origsize=sum(~isnan(level1_5min(:,2)));
    tit='FSP2023';
    %
    % Save Level 2 data
    Instrument='Thermo';
    Model='49i';
    serial_num='0830332761';
    Site=tit;
    Val_level='2';
    Proc_vers='1.0';
    labels{1}='O3';
    units{1}='(ppb)';
    export_LFA_GoAmazon_jan2015(labels,units,Instrument,Model,serial_num,Site,Val_level,Proc_vers,origsize,level1_5min(:,1),level1_5min(:,2))
end
%
% Verificação final
% aux=importdata([caminho1,'FSP2023_Thermo_49i_Level2_5min.csv']);
% level2_5min=aux.data(:,1:2);
% level2_5min(:,1)=level2_5min(:,1)/(24*60*60)+datenum(1970,1,1); % unix time to matlab time
% figure, plot(level2_5min(:,1),level2_5min(:,2),'r.'),grid on, hold on,dynamicDateTicks
% 


