% Ozone monitor (Thermo 49i) data processing 
% Read raw data files, remove repetitions, calculate averages, make plots
% Functions used: frepeated; export_LFA_GoAmazon_jan2015
% Outputs: O3_Level0 (mat_time(utc); 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
% Criteria for Level2: remove calibration and maintanace periods and eventually correct the baseline
%                      
% Obs: trace gas Thermo monitors report mixing ratios normalized to 1013.25 hPa and 0 oC
% Para converter de unix time para matlabtime: datestr(level1_5min(1,1)/(24*60*60)+datenum(1970,1,1))
% ---------------- Parametros -----------------------
caminho='C:\A\GOAMAZON\T1\O3_1109148077sn\raw_dat\';
caminho1='C:\A\GOAMAZON\T1\O3_1109148077sn\';
tit='T1_2015-2016_O3';
level2=false; % logical to produce level2 data
ano=2015;
if(level2)
    % inicio e termino de periodos de manutencao (horario em utc): removidos das analises no level 2
%     manut=xlsread('C:\A\GOAMAZON\ZF2\O3_49i\manut_O3_ZF2.xlsx','tabela','A2:E25'); % doy(utc)|type|ax|b
%     %manut=xlsread([caminho1,'manut_O3_TIWA.xlsx'],'tabela','A2:E26'); % doy(utc)|type|ax|b
%     manut(:,1)=manut(:,1)+datenum(ano,1,1)-1; % doy para matlab time
%     manut(:,2)=manut(:,2)+datenum(ano,1,1)-1;
end
% ---------------- Fim dos parametros --------------------
%
data=zeros(1,12);
%
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)));
    aux2=datevec(aux1,'mm-dd-yyHH:MM'); % date vector
    tmat=datenum(aux2); % date number in Matlab
    aux3=find(strcmpi('--------',C{1,3})==0); % look for flaged data
    nflag=zeros(length(tmat),1);
    nflag(aux3)=1;
    data1=cat(2,tmat,aux,nflag);
    data=cat(1,data,data1);
    flag=cat(1,flag,C{1,3});
    clear aux aux1 aux2 aux3 nflag C data1 nome icol
    fclose(fid);
end
outdata=sortrows(data,1); % sort data based on day of year
outdata(1,:)=[];
clear data
outdata=cat(2,outdata(:,1),datevec(outdata(:,1)),outdata(:,2:12));
%
disp('----Checking for repeated data')
[newdata,erro,linha]=frepeated(outdata,8);
if (strcmp(erro,'none')~=1) 
    disp(erro)
    disp(linha)
else
    clear outdata
    outdata=newdata;
    clear newdata erro
end    
% Corrigindo o time stamp do monitor para reportar o início do período
aux=diff(outdata(1:2,1));
outdata(2:length(outdata),1)=outdata(2:length(outdata),1)-diff(outdata(:,1));
outdata(1,1)=outdata(1,1)-aux;
clear aux
outdata(outdata(:,1)<datenum(ano,1,1),:)=[];
save([caminho1,tit,'_Level0.txt'],'outdata','-ASCII')
%
orig=outdata;
%
figure
subplot(2,3,1),plot(outdata(:,1),outdata(:,10),'b.')
d1=outdata(1,1); d2=outdata(length(outdata),1);
set(gca,'XTick',d1:floor((d2-d1)/2):d2)
set(gca,'XTickLabel',datestr(d1:floor((d2-d1)/2):d2,'dd-mmm'),'FontSize',8)
ylabel('Intensity Cell A (Hz)','FontSize',12),grid on, axis tight
subplot(2,3,4),plot(outdata(:,1),outdata(:,11),'c.')
set(gca,'XTick',d1:floor((d2-d1)/2):d2)
set(gca,'XTickLabel',datestr(d1:floor((d2-d1)/2):d2,'dd-mmm'),'FontSize',8)
ylabel('Intensity Cell B (Hz)','FontSize',12),grid on, axis tight
subplot(2,3,2),plot(outdata(:,1),outdata(:,15),'b.')
set(gca,'XTick',d1:floor((d2-d1)/2):d2)
set(gca,'XTickLabel',datestr(d1:floor((d2-d1)/2):d2,'dd-mmm'),'FontSize',8)
title(tit,'FontSize',12)
ylabel('Flow Cell A','FontSize',12),grid on, axis tight
subplot(2,3,5),plot(outdata(:,1),outdata(:,16),'m.')
set(gca,'XTick',d1:floor((d2-d1)/2):d2)
set(gca,'XTickLabel',datestr(d1:floor((d2-d1)/2):d2,'dd-mmm'),'FontSize',8)
ylabel('Flow Cell B','FontSize',12),grid on, axis tight
subplot(2,3,3),plot(outdata(:,1),outdata(:,17),'g.')
set(gca,'XTick',d1:floor((d2-d1)/2):d2)
set(gca,'XTickLabel',datestr(d1:floor((d2-d1)/2):d2,'dd-mmm'),'FontSize',8)
ylabel('Pressure','FontSize',12),grid on, axis tight
subplot(2,3,6),plot(outdata(:,1),outdata(:,18),'r.'), axis tight
set(gca,'XTick',d1:floor((d2-d1)/2):d2)
set(gca,'XTickLabel',datestr(d1:floor((d2-d1)/2):d2,'dd-mmm'),'FontSize',8)
ylabel('Flag','FontSize',12),grid on
% %saveas(gcf,strcat(caminho1,'O3_diag'),'jpg')
%
figure
plot(orig(:,1),orig(:,8),'r.'),grid on, hold on,dynamicDateTicks
%
% ------ Filters for Level1 data ---------------------------
origsize=size(outdata,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
%
% ------- Additional periods to remove (specific for T1) ---------------
rem=[datenum(2015,12,01,18,00,00) datenum(2015,12,13,00,00,00);...
    datenum(2015,12,23,22,20,00) datenum(2015,12,23,23,30,00);... % zero check 
    datenum(2016,01,11,23,53,00) datenum(2016,01,12,00,05,00);... % zero check
    datenum(2016,01,25,15,00,00) datenum(2016,01,25,15,30,00);... % zero check (0ppb)
    datenum(2016,02,01,14,49,00) datenum(2016,02,01,15,24,00);... % zero check (0ppb)
    datenum(2016,02,08,14,32,00) datenum(2016,02,08,15,29,00);... % zero check (0ppb)
    datenum(2016,02,15,14,48,00) datenum(2016,02,15,15,26,00);... % zero check (0ppb)
    datenum(2016,02,18,12,00,00) datenum(2016,02,18,15,00,00);... % inicializando
    datenum(2016,02,22,12,35,00) datenum(2016,02,22,13,07,00);... % zero check (0ppb)
    datenum(2016,02,23,12,00,00) datenum(2016,02,23,18,00,00);... % 
    datenum(2016,03,10,18,00,00) datenum(2016,03,10,18,55,00);... % inicializando
    datenum(2016,03,14,13,27,00) datenum(2016,03,14,13,40,00);... % se foi zero check, pouco tempo para estabilizar
    datenum(2016,03,22,13,00,00) datenum(2016,03,22,13,20,00);... % se foi zero check, pouco tempo para estabilizar 
    datenum(2016,03,29,14,28,00) datenum(2016,03,29,14,40,00);... % zero check (0ppb)
    datenum(2016,04,06,00,00,00) datenum(2016,04,06,00,15,00);... % se foi zero check, pouco tempo para estabilizar 
    datenum(2016,04,12,22,50,00) datenum(2016,04,12,23,03,00);... % se foi zero check, pouco tempo para estabilizar 
    datenum(2016,04,19,12,45,00) datenum(2016,04,19,12,52,00)]; % se foi zero check, pouco tempo para estabilizar 
for i=1:length(rem)
    outdata(outdata(:,1)>=rem(i,1) & outdata(:,1)<=rem(i,2),:)=[]; 
end
%
figure
plot(orig(:,1),orig(:,8),'r.'),grid on, hold on,dynamicDateTicks
plot(outdata(:,1),outdata(:,8),'b.')
plot([outdata(1,1) outdata(length(outdata),1)],[0 0],'-y')
%
% ------------ Save Level 1 data -----------------------------------------
Instrument='Thermo';
Model='49i';
serial_num='1406960636';
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))
%
if(level2)
    % Lê os arquivos Level 1 
    tit='T1_2015-2016_O3';
    aux=importdata('C:\A\GOAMAZON\T1\LFA_standard\T1_2015-2016_O3_Thermo_49i_Level1_5min.csv');
    level1_5min=aux.data;
    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)));
    orig=level1_5min;
    %
    % Detrending baseline accoring to zero checks
%     manut(:,1)=manut(:,1)-mod(manut(:,1),5/60/24); % arredondando horários de manutenção para 5 minutos inteiros
%     manut(:,2)=manut(:,2)-mod(manut(:,2),5/60/24);
%     for i=1:size(manut,1)
%         aux=find(level1_5min(:,1)>=manut(i,1) & level1_5min(:,1)<=manut(i,2));
%         if(manut(i,3)==0 && ~isempty(aux))
%             level1_5min(aux,2:size(level1_5min,2))=NaN; % remove maintanance periods from dataset
%         elseif(~isempty(aux))
%             % detrend baseline (ozone_detrend.m)
%             if(manut(i,4)==0)
%                 level1_5min(aux,2)=level1_5min(aux,2)-manut(i,5);
%             else
%                 level1_5min(aux,2)=level1_5min(aux,2)+(polyval(manut(i,4:5),level1_5min(aux(1),1)))+(-polyval(manut(i,4:5),level1_5min(aux,1)));           
%             end
%         end
%         clear aux aux1
%     end      
    figure
    plot(orig(:,1),orig(:,2),'r.'),grid on, hold on
    plot(level1_5min(:,1),level1_5min(:,2),'b.'),dynamicDateTicks
    %    
    % Save Level 2 data
    Instrument='Thermo';
    Model='49i';
    serial_num='1406960636';
    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
%
%




