% Ozone monitor (Thermo 49i) data processing 
%
% 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 -----------------------
caminho='C:\A\Diadema2016\data_trailer_2016\o3_lfa\';
caminho1='C:\A\Diadema2016\processados\O3\';
fuso=0; % fuso horario para calcular hora local
tit='Diadema-2016';
filtflag=false; % logical para filtrar dados com flag
filtmanut=true; % logical para filtrar dados em períodos de manutenção ou com problemas
lertxt=false; % logical para ler arquivos txt com concentrações
basicplots=true;
lev1=false;
lev2=false;
%
% inicio e termino de periodos de manutencao (horario em utc): removidos das analises
if(filtmanut)   
    C1=xlsread('C:\A\Diadema2016\logbook_trailer_Diadema2016.xlsx','todos','A4:M135');
    C2=xlsread('C:\A\Diadema2016\logbook_trailer_Diadema2016.xlsx','o3','A4:O34');
    manut(:,1)=cat(1,datenum(C1(:,1:6)),datenum(C2(:,1:6)));
    manut(:,2)=cat(1,datenum(C1(:,7:12)),datenum(C2(:,7:12)));
    manut(:,3)=cat(1,C1(:,13),C2(:,13));
    manut(:,4:5)=cat(1,zeros(length(C1),2),C2(:,14:15));
end
%
% inicio e termino de periodo de zoom para fazer grafico
zoom=[datenum(2016,11,22,0,0,0) datenum(2017,3,15,0,0,0)];
zoom1=[datenum(2017,03,05,0,0,0) datenum(2017,3,15,0,0,0)];%
% ---------------- Fim dos parametros --------------------
%
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);%-datenum(ano,1,1)+1; % matlab time
    aux2=find(strcmpi('--------',C{1,3})==0); % look for flaged data
    nflag=zeros(length(mattime),1);
    nflag(aux2)=1;
    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
%
if(lertxt)
   caminhotxt='C:\A\GOAMAZON\T1\O3_1109148077sn\raw_txt\processar\';
   datatxt=zeros(1,2);
   files=dir(caminhotxt);
   disp('-----Reading txt data files')
   for ifile=3:size(files,1)
       nome=getfield(files, {ifile,1}, 'name');
       disp(nome)
       C=importdata(strcat(caminhotxt,nome));
       aux1=strcat((C.textdata(:,1)),(C.textdata(:,2)));
       datatxt1(:,1)=datenum(aux1,'yyyy-mm-ddHH:MM:SS');
       datatxt1(:,2)=C.data;
       datatxt=cat(1,datatxt,datatxt1);
       clear aux aux1 C datatxt1 nome 
   end
   datatxt(1,:)=[];
   datatxt=sortrows(datatxt,1); % sort data based on day of year
   first =floor(datatxt(1,1));   % first averaging time
   last = ceil(datatxt(length(datatxt),1));  % last averaging time
   delta = 5/60/24;     % 5min
   tol=0;
   stat=fmedia_qualquer_coisa2(datatxt,first,last,delta,tol);
   aux=find(stat(:,1)<=outdata(length(outdata),1)); % procura pontos que já estão em outdata (arquivos do iport)
   stat(aux,:)=[];
   stat=cat(2,stat,zeros(length(stat),9));
   stat(:,3:18)=NaN;   
   stat(:,8)=stat(:,2);
   stat(:,2:7)=datevec(stat(:,1));   
   outdata=cat(1,outdata,stat);
   clear datatxt stat
end
%
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
if(filtmanut)
    outdataorig=outdata;
    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
    figure,plot(outdataorig(:,1),outdataorig(:,8),'r.'),hold on,grid on
    plot(outdata(:,1),outdata(:,8),'b.')
    dynamicDateTicks
    clear aux
end
%
% Remove flagged periods from dataset
if(filtflag)
    aux=find(outdata(:,18)==1);
    outdata(aux,:)=[];
    clear aux
end
%
% Corrigindo o time stamp do monitor para reportar o início do período
% O sampling time variou no tempo: 5min até 11/jan e 1min depois disso. 
% figure,plot(outdata(2:length(outdata),1),diff(outdata(:,1))*24*60,'b.'),grid on, dynamicDateTicks
deltat1=datenum(0,0,0,0,5,0);
deltat2=datenum(0,0,0,0,1,0);
aux=find(outdata(:,1)<datenum(2017,1,11,12,59,0));
outdata(aux,1)=outdata(aux,1)-deltat1;
aux=find(outdataorig(:,1)<datenum(2017,1,11,12,59,0));
outdataorig(aux,1)=outdataorig(aux,1)-deltat1;
aux=find(outdata(:,1)>=datenum(2017,1,11,12,59,0));
outdata(aux,1)=outdata(aux,1)-deltat2;
aux=find(outdataorig(:,1)>=datenum(2017,1,11,12,59,0));
outdataorig(aux,1)=outdataorig(aux,1)-deltat2;
%
% Saving output files (level0 and clean)
save([caminho1,tit,'_ozone_level0.mat'],'outdataorig','-mat')
if(filtflag || filtmanut)
    save([caminho1,tit,'_ozone_clean.mat'],'outdata','-mat')    
end
%
% -------------------------------------------------------------------------
%   Basic Plots
if(basicplots)
    figure
    subplot(2,3,1),plot(outdata(:,1)+(fuso/24),outdata(:,10),'b.')
    dynamicDateTicks
    ylabel('Intensity Cell A (Hz)','FontSize',12),grid on, axis tight
    subplot(2,3,4),plot(outdata(:,1)+(fuso/24),outdata(:,11),'c.')
    dynamicDateTicks
    ylabel('Intensity Cell B (Hz)','FontSize',12),grid on, axis tight
    subplot(2,3,2),plot(outdata(:,1)+(fuso/24),outdata(:,15),'b.')
    dynamicDateTicks
    title(tit,'FontSize',12)
    ylabel('Flow Cell A','FontSize',12),grid on, axis tight
    subplot(2,3,5),plot(outdata(:,1)+(fuso/24),outdata(:,16),'m.')
    dynamicDateTicks
    ylabel('Flow Cell B','FontSize',12),grid on, axis tight
    subplot(2,3,3),plot(outdata(:,1)+(fuso/24),outdata(:,17),'g.')
    dynamicDateTicks
    ylabel('Pressure','FontSize',12),grid on, axis tight
    subplot(2,3,6),plot(outdata(:,1)+(fuso/24),outdata(:,18),'r.'), axis tight
    dynamicDateTicks
    ylabel('Flag','FontSize',12),grid on
    saveas(gcf,strcat(caminho1,tit,'_O3_diag'),'jpg')
    %
    figure
    subplot(2,1,1),plot(outdata(:,1)+(fuso/24),outdata(:,8),'b.')
    title(tit,'FontSize',14)
    ylabel('O_3 (ppb)','FontSize',14),grid on, axis tight
    xlabel('UTC Time')
    dynamicDateTicks
    subplot(2,1,2),plot(mod(outdata(:,1)+(fuso/24),1),outdata(:,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,tit,'_O3_conc'),'jpg')
    %
    aux=find(outdata(:,1)>=zoom(1) & outdata(:,1)<=zoom(2));
    figure
    subplot(2,1,1),plot(outdata(aux,1)+(fuso/24),outdata(aux,8),'b.')
    title(tit,'FontSize',14)
    ylabel('O_3 (ppb)','FontSize',14),grid on, axis tight
    %xlabel('UTC Time')
    %dynamicDateTicks
    set(gca,'XTick',zoom(1):2:zoom(2))
    set(gca,'XTickLabel',datestr(zoom(1):2:zoom(2),'dd/mm'))
    set(gca,'XLim',[zoom(1) zoom(2)])
    subplot(2,1,2),plot(mod(outdata(aux,1)+(fuso/24),1),outdata(aux,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'})
    saveas(gcf,strcat(caminho1,tit,'O3_zoom1_conc'),'jpg')
    %
    aux=find(outdata(:,1)>=zoom1(1) & outdata(:,1)<=zoom1(2));
    figure
    subplot(2,1,1),plot(outdata(aux,1)+(fuso/24),outdata(aux,8),'b.')
    title(tit,'FontSize',14)
    ylabel('O_3 (ppb)','FontSize',14),grid on, axis tight
    xlabel('UTC Time')
    set(gca,'XTick',zoom1(1):2:zoom1(2))
    set(gca,'XTickLabel',datestr(zoom1(1):2:zoom1(2),'dd/mm'))
    set(gca,'XLim',[zoom1(1) zoom1(2)])
    subplot(2,1,2),plot(mod(outdata(aux,1)+(fuso/24),1),outdata(aux,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'})
    saveas(gcf,strcat(caminho1,tit,'O3_zoom2_conc'),'jpg')
    %
    % Ciclo diurno médio
    first = 0;   % first averaging time
    last = 23.5/24;  % last averaging time
    delta = 30/60/24;     % 30min
    mdata=cat(2,mod(outdata(:,1),1),outdata(:,8));
    tol=0;
    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('Average Diurnal Cycle','FontSize',14)
    hgsave(gcf,strcat(caminho1,tit,'_O3_dcycle_level0'))
    saveas(gcf,strcat(caminho1,tit,'_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)
    C=load('C:\A\Diadema2016\processados\O3\Diadema-2016_ozone_clean.mat');
    outdata=C.outdata;
    C=load('C:\A\Diadema2016\processados\O3\Diadema-2016_ozone_level0.mat');
    outdataorig=C.outdataorig;
    tit='Diadema-2016';
    %
    origsize=size(outdataorig,1);
    aux=find(outdata(:,15)<0.55 | outdata(:,16)<0.55 | outdata(:,10)<75000 | outdata(:,11)<75000);
    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('C:\A\Diadema2016\processados\LFA_standard\Diadema-2016_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='Diadema-2016';
    %
    % 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('C:\A\Diadema2016\processados\LFA_standard\Diadema-2016_Thermo_49i_Level1_5min.csv');
% level1_5min=aux.data(:,1:2);
% level1_5min(:,1)=level1_5min(:,1)/(24*60*60)+datenum(1970,1,1); % unix time to matlab time
% aux=importdata('C:\A\Diadema2016\processados\LFA_standard\Diadema-2016_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(level1_5min(:,1),level1_5min(:,2),'b.'),grid on, hold on,dynamicDateTicks
% plot(level2_5min(:,1),level2_5min(:,2),'r.')


