% CPC TSI data processing - Diadema
% Flex data files (Labview DAQ)
% Read raw data files, apply corrections, calculate averages, make plots
% Functions used: frepeated; export_LFA_GoAmazon_jan2015
% Outputs: cpc_Level0 (mattime(utc);conc(cm-3))--> all raw data, without filters
%          cpc_Level1_5min (dd-mmm-yyyy HH:MM; sec_since_1970; mean; std; median; p25; p75; max; min; num)
%          cpc_Level1_30min
% Time stamps in the output files refer to the start of the period  
%
% Level 1: 1-min conc must be in the range 50-10^5 cm-3. No instrument error should be reported.
% Level 2: Apply STP corrections and coincidence corrections.
%          Os dados flex vêm direto da saída serial do CPC, utilizando o comando RALL.
%          Do manual do CPC, entendi que a concentração que sai via serial já está corrigida para erros de coincidência. 
% Level 3: Intercomparison with other CPC if available. Particle counts comparison must be within 20%.
%
% Parametros de entrada --------------------------------------------
caminho1='C:\A\Diadema2016\flex_data_files\CPC3772\csv\'; % arquivos flex (Labview)
caminho2='C:\A\Diadema2016\processados\CPC3772\';
tit='Diadema-2016';
pref=1013.25; % Reference pressure for STP correction (mbar)
tref=293.15; % Reference temp for STP correction (K)
corrigirtempo=false;
filtmanut=true; % logical para filtrar dados em períodos de manutenção ou com problemas
readprev=false; % logical for reading previous processed data 
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','cpc','A4:M13');
    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));    
end
%
% Program start -----------------------------------------------------
files=dir(caminho1);
alldata=zeros(1,3);
if (readprev)
    prevdata=load([caminho2,tit,'_cpc3772_level0.mat']);
    lerdaqui=floor(prevdata.alldata(length(prevdata.alldata),1));    
else
    lerdaqui=datenum(1900,1,1);
end
disp('-----Reading flexdata files')
for ifile=3:size(files,1)
    nome=getfield(files, {ifile,1}, 'name');
    aux1=strfind(nome,'3772-')+5;
    aux2=strfind(nome,'.csv')-1;
    filedia=datenum(nome(aux1:aux2),'mm_dd_yyyy');
    if(filedia>=lerdaqui)
        disp(nome)
        fid=fopen(strcat(caminho1,nome));
        C=cell(floor(getfield(files, {ifile,1}, 'bytes')/92),5);
        i=1;
        while (~feof(fid))
            C(i,:)=textscan(fid,'%*s%s%s%s %f%f %*s%*s%*s%*s%*s%*s%*s%*s%*s%*s',1,'Delimiter',',');
            i=i+1;
        end
        fclose(fid);
        %
        % removendo linhas com defeito ou cabeçalhos
        rem=0;
        for j=4:5
            aux=find(cellfun(@isempty,C(:,j)));
            rem=cat(1,rem,aux);
        end
        rem(1)=[];
        C(rem,:)=[];
        rem=find(cellfun(@isnan,C(:,4)));
        C(rem,:)=[];
        %
        data_(:,2:3)=cell2mat(C(:,4:5));
        for i=1:length(C)
            data_(i,1)=datenum(strcat(C{i,1},C{i,2},C{i,3}),'mm/dd/yyyyHH:MM:SSAM');
        end
        alldata=cat(1,alldata,data_);
        clear data_ C
    end    
end
alldata(1,:)=[];
if (readprev)
    aux=find(prevdata.alldata(:,1)<min(alldata(:,1)));
    alldata=cat(1,prevdata.alldata(aux,:),alldata);
end
alldata=sortrows(alldata,1); % sort data based on day of year
%
if(corrigirtempo) % correção segundo o logbook
end
% 
% Obs: é comum a repetição de horários nos arquivos flexdata. Nada grave,
% dura alguns segundos por dia, com valores de concentração ligeiramente diferentes.
% Tem a ver com um descompasso entre o instante de coleta do LabView e a
% frequência de medidas do CPC. Assim, vou dispensar a verificação de 
% dados repetidos aqui, e simplesmente remover linhas com mesmo time stamp.
disp('----Checking for repeated data')
aux=find(diff(alldata(:,1))==0);
disp(['Número de linhas com timestamp repetido: ',num2str(length(aux))])
alldata(aux,:)=[];
clear aux
% [newdata,erro]=frepeated(alldata,2);
% if (strcmp(erro,'none')~=1) 
%     disp(erro)
% else
%     alldata=newdata;
% end
% clear aux erro 
%
save([caminho2,tit,'_cpc3772_level0.mat'],'alldata','-mat') 
% 
% Remove períodos de manutencao/contaminacao de acordo com o logbook
if(filtmanut)
   alldataorig=alldata;
   for i=1:size(manut,1)
        aux=find(alldata(:,1)>=manut(i,1) & alldata(:,1)<=manut(i,2));
        if(manut(i,3)==0)            
            alldata(aux,:)=[]; % remove maintanance periods from dataset
        end    
   end
   save([caminho2,tit,'_cpc3772_clean.mat'],'alldata','-mat') 
   csvwrite([caminho2,tit,'_cpc3772_clean.csv'],[alldata(:,1:2),datevec(alldata(:,1))]) 
   figure,plot(alldataorig(:,1),alldataorig(:,2),'r.'),hold on,grid on
   plot(alldata(:,1),alldata(:,2),'c.')
   ylabel('CPC3772 (cm^-^3)'),title(tit)
   dynamicDateTicks   
   saveas(gcf,strcat(caminho2,tit,'_cpc3772_level0'),'jpg')
else   
    figure
    plot(alldata(:,1),alldata(:,2),'c.'),grid on
    ylabel('CPC3772 (cm^-^3)'),title(tit)
    %set(gca,'XLim',[datenum(2016,11,22) datenum(2016,12,11)])
    dynamicDateTicks
    saveas(gcf,strcat(caminho2,tit,'_cpc3772_level0'),'jpg')
end
%
%
% ------------ Level 1 ------------------------------------------------
% Remove local contamination (clean data)
%
if(lev1)
    C=load('C:\A\Diadema2016\processados\CPC3772\Diadema-2016_cpc3772_clean.mat');
    alldata=C.alldata;
    C=load('C:\A\Diadema2016\processados\CPC3772\Diadema-2016_cpc3772_level0.mat');
    alldataorig=C.alldata;
    tit='Diadema-2016';
    %
    origsize=size(alldataorig,1);
    rem=find(alldata(:,2)<50 | alldata(:,2)>10^5);
    alldata(rem,:)=[];    
    figure,plot(alldataorig(:,1),alldataorig(:,2),'r.'),hold on,grid on
    plot(alldata(:,1),alldata(:,2),'c.'),ylabel('CPC3772 (cm^-^3)'),title(tit),dynamicDateTicks
    %
    % Save Level 1 data
    Instrument='CPC';
    Model='3772';
    serial_num=' ';
    Site=tit;
    Val_level='1';
    Proc_vers='1.0';
    labels{1}='conc';
    units{1}='(cm-3)';
    export_LFA_GoAmazon_jan2015(labels,units,Instrument,Model,serial_num,Site,Val_level,Proc_vers,origsize,alldata(:,1),alldata(:,2))
end

%
% ------------ Level 2 ------------------------------------------------
% Level 2: Apply STP corrections and coincidence corrections.
% Os dados flex vêm direto da saída serial do CPC, utilizando o comando RALL.
% Do manual do CPC, entendi que a concentração que sai via serial já está
% corrigida para erros de coincidência. 
if(lev2)        
    % Lê os arquivos Level 1 que ja salvou
    aux=importdata('C:\A\Diadema2016\processados\LFA_standard\Diadema-2016_CPC_3772_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';
    %
    % Lê dados de meteorologia e calcula fator stp
    C=importdata('C:\A\Diadema2016\meteorologia_IAG\Diadema-2016_Meteo_5min.csv');
    met(:,1)=datenum(C.textdata(3:length(C.textdata)),'dd-mmm-yyyy HH:MM:SS ');
    met(:,2:3)=C.data;
    met(:,1)=met(:,1)+3/24; % local time to utc time
    met(:,3)=met(:,3)+273.15; % oC to K
    met(met(:,1)<level1_5min(1,1),:)=[];
    met(met(:,1)>level1_5min(length(level1_5min),1),:)=[];
    if(length(met)~=length(level1_5min))
        disp('check length of meteorological data file')
        stop
    end    
    pref=1013.25; tref=293.15;
    stp=(pref*met(:,3))./(tref*met(:,2));
    level1_5min(:,2)=level1_5min(:,2).*stp;
    %
    % Save Level 2 data 
    Instrument='CPC';
    Model='3772';
    serial_num=' ';
    Site=tit;
    Val_level='2';
    Proc_vers='1.0';
    labels{1}='conc';
    units{1}='(cm-3)';
    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_CPC_3772_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_CPC_3772_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.')


