% SO2 monitor (Thermo 43i) data processing 
% Read raw data files, remove repetitions, calculate averages, save files
% according to LFA protocols
% Functions used: frepeated; export_LFA_GoAmazon_jan2015
% Outputs: SO2_Level0 (mat_time(utc); so2; hiso2; intt; rctt; pgast; pres; smplfl; pmtv; lmpv; lmpi; flag(0/1))--> all raw data, without filters
%          SO2_Level1_5min (dd-mmm-yyyy HH:MM; sec_since_1970; mean_so2; std_so2; median_so2; p25_so2; p75_so2; max_so2; min_so2; num_so2)
%          SO2_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: (Exclude records when) Flow < 0.38 LPM; rctt < 43oC; lmpi<88% (and the subsequent point)
%                      -10<ConcO3<100 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\SO2_1406960668sn\raw_dat\';% folder containing raw data
caminho1='C:\A\GOAMAZON\T1\SO2_1406960668sn\';% folder where outputs will be saved
fuso=0; % fuso horario para calcular hora local
tit='T1_2015-2016_SO2';
lv2=false; % logical para gerar arquivos level 2
graficos=false;
% inicio e termino de periodos de manutencao (horario em utc): removidos das analises
% manut=xlsread([caminho1,'manut_SO2_TIWA.xlsx'],'tabela','A2:E41'); % doy(utc)|type|ax|b
% manut(:,1)=manut(:,1)+datenum(2014,1,1)-1; % doy para matlab time
% manut(:,2)=manut(:,2)+datenum(2014,1,1)-1;
%
% ---------------- Fim dos parametros --------------------
%
data=zeros(1,12);
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)));
    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 Matlab date number
outdata(1,:)=[];
clear data
%
disp('----Checking for repeated data')
[newdata,erro]=frepeated(outdata,2);
if (strcmp(erro,'none')~=1) 
    disp(erro)
end
clear outdata
outdata=newdata;
clear newdata erro 
% 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
save([caminho1,tit,'_SO2_Level0.txt'],'outdata','-ASCII')
%
% ----------- Level0 plots --------------
if(graficos)
    figure
    subplot(2,3,1),plot(outdata(:,1)+(fuso/24),outdata(:,4),'b.'),hold on
    plot(outdata(:,1)+(fuso/24),outdata(:,5),'r.'),hold off,dynamicDateTicks
    ylabel('Temp (oC)','FontSize',12),grid on, axis tight
    subplot(2,3,4),plot(outdata(:,1)+(fuso/24),outdata(:,11),'c.'),dynamicDateTicks
    ylabel('Lamp Intensity (%)','FontSize',12),xlabel('UTC'),grid on, axis tight
    subplot(2,3,2),plot(outdata(:,1)+(fuso/24),outdata(:,8),'b.'),dynamicDateTicks
    title(tit,'FontSize',12)
    ylabel('Sample Flow','FontSize',12),grid on, axis tight
    subplot(2,3,5),plot(outdata(:,1)+(fuso/24),outdata(:,9),'m.'),dynamicDateTicks
    ylabel('PMT Voltage (V)','FontSize',12),xlabel('UTC'),grid on, axis tight
    subplot(2,3,3),plot(outdata(:,1)+(fuso/24),outdata(:,7),'g.'),dynamicDateTicks
    ylabel('Pressure','FontSize',12),grid on, axis tight
    subplot(2,3,6),plot(outdata(:,1)+(fuso/24),outdata(:,12),'r.'), axis tight,dynamicDateTicks
    ylabel('Flag','FontSize',12),xlabel('UTC'),grid on
    saveas(gcf,strcat(caminho1,tit,'-diag'),'jpg')
%   Time series and diurnal cycle 
    figure
    subplot(2,1,1),plot(outdata(:,1),outdata(:,2),'b.')
    title(tit,'FontSize',14)
    ylabel('SO2 [ppb]','FontSize',14),grid on, axis tight,dynamicDateTicks
    subplot(2,1,2),plot(mod(outdata(:,1)-4/24,1),outdata(:,2),'b.')
    ylabel('SO2 [ppb]','FontSize',14),grid on,xlabel('Local Time','FontSize',14), axis tight,dynamicDateTicks
    set(gca,'XLim',[0 1])
    set(gca,'XTick',0:0.25:1)
    saveas(gcf,strcat(caminho1,tit,'-conc'),'jpg')    
    % Zoom 
    aux=find(outdata(:,1)>=zoom(1) & outdata(:,1)<=zoom(2));
    figure
    subplot(2,1,1),plot(outdata(aux,1),outdata(aux,2),'b.')
    title(tit,'FontSize',14)
    ylabel('SO2 [ppb]','FontSize',14),grid on, axis tight,dynamicDateTicks
    subplot(2,1,2),plot(mod(outdata(aux,1)-4/24,1),outdata(aux,2),'b.')
    ylabel('SO2 [ppb]','FontSize',14),grid on,xlabel('Local Time','FontSize',14), axis tight,dynamicDateTicks
    set(gca,'XLim',[0 1])
    set(gca,'XTick',0:0.25:1)
    saveas(gcf,strcat(caminho1,tit,'-zoom'),'jpg')
end
%
% ------ Filters for Level1 data ---------------------------
origsize=size(outdata,1);
orig=outdata;
aux=find(outdata(:,8)<0.38 | outdata(:,5)<43 | outdata(:,11)<88);
aux1=aux(diff(aux)>1)+1; % subsequent point
aux2=find(outdata(:,2)<=-10 | outdata(:,2)>=100);
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,04,00,00,00);...
    datenum(2015,12,05,21,55,00) datenum(2015,12,06,22,44,00);... %
    datenum(2015,12,06,19,00,00) datenum(2015,12,06,19,12,00);... %
    datenum(2015,12,08,12,53,00) datenum(2015,12,08,15,00,00);... % zero check (-0.25ppb) 
    datenum(2015,12,09,20,48,00) datenum(2015,12,09,21,00,00);... %
    datenum(2015,12,12,19,50,00) datenum(2015,12,12,20,52,00);... %
    datenum(2015,12,16,02,00,00) datenum(2015,12,16,02,40,00);... % inicializando
    datenum(2015,12,24,00,30,00) datenum(2015,12,24,01,06,00);... % zero check (pouco tempo para estabilizar) --> parece que o zero foi ajustado errado, pois tem um degrau na timeseries
    %datenum(2016,01,02,21,04,00) datenum(2016,01,03,00,52,00);... % fluxo aumentou (zero check?)
    datenum(2016,01,11,23,52,00) datenum(2016,01,12,00,57,00);... % zero check (-0.58ppb)
    datenum(2016,01,18,21,42,00) datenum(2016,01,18,22,40,00);... % zero check (-0.145ppb)  
    datenum(2016,01,25,15,05,00) datenum(2016,01,25,15,34,00);... % zero check (0.089ppb)
    datenum(2016,02,01,14,48,00) datenum(2016,02,01,15,25,00);... % zero check (0.041ppb)
    datenum(2016,02,08,14,36,00) datenum(2016,02,08,15,26,00);... % zero check (0.0ppb)
    datenum(2016,02,15,14,53,00) datenum(2016,02,15,15,39,00);... % zero check (0.036ppb)
    datenum(2016,02,18,12,00,00) datenum(2016,02,18,14,30,00);... % inicializando
    datenum(2016,02,22,12,34,00) datenum(2016,02,22,13,15,00);... % zero check (0.173ppb)
    datenum(2016,02,23,12,55,00) datenum(2016,02,23,15,05,00);... % inicializando
    datenum(2016,03,04,21,55,00) datenum(2016,03,04,22,20,00);... % zero check (-0.011ppb)
    datenum(2016,03,10,11,30,00) datenum(2016,03,10,19,00,00);... % inicializando
    datenum(2016,03,14,13,30,00) datenum(2016,03,14,13,50,00);... % zero check (0.024ppb)
    datenum(2016,03,22,13,06,00) datenum(2016,03,22,13,37,00);... % zero check (-0.031ppb)
    datenum(2016,03,29,14,00,00) datenum(2016,03,29,15,10,00);... % zero check (-0.017ppb)
    datenum(2016,04,06,00,00,00) datenum(2016,04,06,00,39,00);... % zero check (0.086ppb)
    datenum(2016,04,12,22,54,00) datenum(2016,04,12,23,15,00);... % zero check (pouco tempo para estabilizar)
    datenum(2016,04,19,12,30,00) datenum(2016,04,19,13,25,00)]; % zero check (0.0ppb)
for i=1:length(rem)
    outdata(outdata(:,1)>=rem(i,1) & outdata(:,1)<=rem(i,2),:)=[]; 
end
%
figure
plot(orig(:,1),orig(:,2),'-r.'),grid on, hold on,dynamicDateTicks
plot(outdata(:,1),outdata(:,2),'b.')
%
% ------------ Save Level 1 data -----------------------------------------
Instrument='Thermo';
Model='43iTLE';
serial_num='1406960668';
Site=tit;
Val_level='1';
Proc_vers='1.0';
labels{1}='SO2';
units{1}='(ppb)';
export_LFA_GoAmazon_jan2015(labels,units,Instrument,Model,serial_num,Site,Val_level,Proc_vers,origsize,outdata(:,1),outdata(:,2))
%
%
% Level2 -------------------------------------------------------------- 
% Specific to T1 data: during the zero checks, SO2 conc was ~0.02 ppb
% higher than the ambient measurements. Either the zero air was not "zero"
% for SO2, or/and the ambient concentrations were virtually zero most of the time.
% The timeseries show a diurnal cycle. For level2, so2 data was subtracted
% by a constant factor of 0.425ppb, in a way that 25% of data falls below
% 0.0ppb. 
if(lv2)
    % Lê os arquivos Level 1 
    aux=importdata('C:\A\GOAMAZON\T1\LFA_standard\T1_2015-2016_SO2_Thermo_43iTLE_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)));
    %
    figure,plot(level1_5min(:,1),level1_5min(:,2),'-r.'),grid on, hold on,dynamicDateTicks
    %
    % Specific correction for T1 data:
    aux=find(level1_5min(:,1)<=datenum(2015,12,08,15,00,00));
    level1_5min(aux,2)=level1_5min(aux,2)-(-0.25); 
    aux=find(level1_5min(:,1)>=datenum(2015,12,24,00,30,00) & level1_5min(:,1)<=datenum(2016,01,12,00,57,00));
    level1_5min(aux,2)=level1_5min(aux,2)-(0.58); % o zero parece ter sido ajustado errado neste período, pois há um degrau na série temporal      
    %    
    plot(level1_5min(:,1),level1_5min(:,2),'b.')
%
%     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);    
%     level1_5min(level1_5min(:,2)<-0.5,2)=NaN;
%     for i=1:size(manut,1) % 5min level 1 data
%         aux=find(level1_5min(:,1)>=manut(i,1) & level1_5min(:,1)<=manut(i,2));
%          if(manut(i,3)==0)
%                 level1_5min(aux,2:9)=NaN; % remove maintanance periods from dataset
%          else
%              % 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         
%     end
%     clear aux
    %
    Instrument='Thermo';
    Model='43iTLE';
    serial_num='1406960668';
    Site=tit;
    Val_level='2';
    Proc_vers='1.0';
    labels{1}='SO2';
    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
%
