% SMPS-TSI-3082 (version3) data processing from txt files (LRizzo)
% Interpolated to a commom diameter sequence
% Funcoes utilizadas: fsmps_interp, fsmps_param, fmedia2.m
%
% Obs: às vezes o software AIM gera uma linha de amostra com caracteres especiais e falha nos dados. 
% Por exemplo: "3.55,1.00e+03,-1.#IO,-1.#IO,-1.#IO,0.000,"
% O script do matlab falha ao ler essas linhas, pois elas têm um número de colunas diferente do resto do arquivo. 
% Geralmente excluo essas linhas na mão no arquivo txt.
%
% Output: alldata(level0) 
% 1: matdate(utc) 
% 2:8 - Dpg(nm), Mode(nm), GeoStd, TotalConc(cm-3), N50(cm-3), SurfArea(um2/cm3), Vol(um3/cm3) 
% 9:113 - dN/dlogDp
% 113:118 - SampleTemp(C), SamplePress(kPa), SampleRH(%), SheathFlow(lpm), AerosolFlow(lpm)
%
% Level1: 1-min aerosol concentration must be >50 cm-3 and <100 000 cm-3. 
%         Sheath flow must be > 3 lpm and < 12 lpm. 
%         Sample flow must be > 0.3 lpm and < 2 lpm. 
%         Status flag must be “Normal Scan”. 
%         Multiple charge correction must have been applied.
% Level2: STP correction (1013.25; 293.15).
% Level3: Integrated aerosol number must be within 30% CPC measurement. 
%         Considering an aerosol density of 1.4 g cm-3, integrated mass must be within 50%-120% of TEOM PM2.5.    
%  
% Parametros de entrada:
user='LRizzo';%'LRizzo'; % 'luvar'
caminho=['C:\Users\',user,'\Dropbox\A\FSP2023\raw\smps3082\txt\'];
caminho1=['C:\Users\',user,'\Dropbox\A\FSP2023\processados\'];
caminho2=['C:\Users\',user,'\Dropbox\A\FSP2023\processados\plots\'];
fuso=0; % fuso horario para calcular hora local
tit='FSP2023';
dens=1.0; % particle density, in g/cm3
ano=2023;
lastday=[31 28 31 30 31 30 31 31 30 31 30 31];
timezone=0;
% diameter sequence to which you want the data to be interpolated (this code does not extrapolate!!)
diam=[10.2, 10.6, 10.9, 11.3, 11.8, 12.2, 12.6, 13.1, 13.6, 14.1, 14.6, 15.1, 15.7, 16.3, 16.8, 17.5, 18.1, 18.8, 19.5, 20.2, 20.9, 21.7, 22.5, 23.3, 24.1, 25.0, 25.9, 26.9, 27.9, 28.9, 30.0, 31.1, 32.2, 33.4, 34.6, 35.9, 37.2, 38.5, 40.0, 41.4, 42.9, 44.5, 46.1, 47.8, 49.6, 51.4, 53.3, 55.2, 57.3, 59.4, 61.5, 63.8, 66.1, 68.5, 71.0, 73.7, 76.4, 79.1, 82.0, 85.1, 88.2, 91.4, 94.7, 98.2,101.8,105.5,109.4,113.4,117.6,121.9,126.3,131.0,135.8,140.7,145.9,151.2,156.8,162.5,168.5,174.7,181.1,187.7,194.6,201.7,209.1,216.7,224.7,232.9,241.4,250.3,259.5,269.0,278.8,289.0,299.6,310.6,322.0,333.8,346.0,358.7,371.8,385.4,399.5,414.2,429.4];
filtmanut=false; % logical para filtrar dados em períodos de manutenção ou com problemas
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');
%     manut(:,1)=datenum(C1(:,1:6));
%     manut(:,2)=datenum(C1(:,7:12));
%     manut(:,3)=C1(:,13);
end
%
files=dir(caminho);
if(size(files,1)==0)
    disp('Error: raw data files path might be wrong')
    stop
end
tempo=zeros(1,1);param=zeros(1,7);totconc=zeros(1,4);sizedist=zeros(1,length(diam));
disp('-----Reading raw data files')
for ifile=3:size(files,1)
    nome=getfield(files, {ifile,1}, 'name');
    disp(nome)
    fid=fopen(strcat(caminho,nome));
    config=textscan(fid,'%s',63,'Delimiter',',');
    if(strcmp(config{1,1}(55),'FALSE'))
        disp('Warning: no multiple charge correction!!!')
    end
    aux=textscan(fid,'%s',9,'Delimiter',',');
    % verifica se chegou no início dos diâmetros
    if~(strcmp(aux{1,1}(9),'Diameter Midpoint (nm)'))
        disp('Did not find diameter start in raw data file.')
        stop
    end
    aux=textscan(fid,'%f',1,'Delimiter',',');
    idiam=1;
    diam_=nan(1,100);
    while(cell2mat(aux)>0)
        diam_(idiam)=cell2mat(aux);
        idiam=idiam+1;
        aux=textscan(fid,'%f',1,'Delimiter',',');
    end
    aux=textscan(fid,'%s',29,'Delimiter',',');
    % verifica se chegou no final do cabeçalho
    if~(strcmp(aux{1,1}(29),'Comment'))
        disp('Did not find end of header.')
        stop
    end
    %
    f1='%*f%s%s %f%f%f%*f%*f%*s';
    f3='%*f%*f%*f%*f %f%f %*f%*f%*f %f%f %*f%*f%*f%*f%*f%*f %f%f%f%f %*s%*s%*s%*s%*s%*s%*s%*s';
    f2='';
    for i=1:length(diam_)
        f2=strcat(f2,'%f');
    end
    format=strcat(f1,f2,f3);
    C=textscan(fid,format,'Delimiter',',');
    fclose(fid);
    %
    % Retrieve normalized size distribution
    sizedist_=zeros(size(C{1,1},1),length(diam_));        
    for j=6:length(diam_)+5  
        sizedist_(:,j-5)=C{1,j};
    end
    % Retrieve control parameters [1-SampleTemp(C), 2-SamplePress(kPa), 3-SampleRH(%), ...
    ... 4-SheathFlow(lpm), 5-AerosolFlow(lpm), 6-LowerSize(nm), 7-UpperSize(nm)]
    param_=zeros(size(C{1,1},1),7);
    for j=1:3
        param_(:,j)=C{1,j+2};
    end
    jj=4;
    for j=length(diam_)+6:size(C,2)-4
        param_(:,jj)=C{1,j};
        jj=jj+1;
    end
    % Retrieve [Geo.Mean(nm) Mode(nm) Geo.Std.Dev. TotalConc.(cm-3)]
    totconc_=zeros(size(C{1,1},1),4);
    for j=size(C,2)-3:size(C,2)
        totconc_(:,j-size(C,2)+4)=C{1,j};
    end
    %disp('Convertendo data e hora para matriz de tempo')
    aux1=strcat((C{1,1}(:,1)),'-',(C{1,2}(:,1)));
    tempo_=datenum(aux1,'yyyy/mm/dd-HH:MM:SS'); % matlab date 
    %
    interp=false;
    if(length(diam_)==length(diam))          
        if(sum(diam_==diam)~=length(diam))
            interp=true;
        end
    else
        interp=true;
    end    
    if(interp) % Interpolate to match a diameter sequence defined by the user, if necessary 
        disp('Interpolating to a commom diameter sequence')
        [newsizedist]=fsmps_interp(diam,diam_,sizedist_,tempo_);
        sizedist_=newsizedist(:,2:size(newsizedist,2));        
    end
    %
    tempo=cat(1,tempo,tempo_);
    param=cat(1,param,param_);
    totconc=cat(1,totconc,totconc_);
    sizedist=cat(1,sizedist,sizedist_);
    %
    clear C aux aux1 aux2 param_ totconc_ diam_ tempo_ sizedist_
end
%
% Calculate size dist parameters
[maisparam]=fsmps_param(diam,9,sizedist,false); % N(cm-3), Dpg(nm), sg, SurfaceArea(um2/cm3), Volume(um3/cm3), SMD(nm), VMD(nm), Nlt50(cm-3), Ngt100(cm-3), Nlt70(cm-3)
%

% Saving Level 0 data --------------------------------------------------
% 1: matdate(utc) 
% 2:8 - Dpg(nm), Mode(nm), GeoStd, TotalConc(cm-3), N50(cm-3), SurfArea(um2/cm3), Vol(um3/cm3) 
% 9:113 - dN/dlogDp
% 113:118 - SampleTemp(C), SamplePress(kPa), SampleRH(%), SheathFlow(lpm), AerosolFlow(lpm)
nomes = ["datetime(utc)", "Dpg(nm)", "Mode(nm)", "GeoStd", "TotConc(cm-3)", "N50(cm-3)", "SurfArea(um2/cm3)", "Vol(um3/cm3)", ...
    string(diam), "SampleTemp(C)", "SamplePress(kPa)", "SampleRH(%)", "SheathFlow(lpm)", "AerosolFlow(lpm)"];
alldata=cat(2,tempo,totconc,maisparam(:,8),maisparam(:,4:5),sizedist,param(:,1:5));
alldata=sortrows(alldata,1); % sort data based on day of year
alldata(1,9:length(diam)+8)=diam;
level0 = array2table(alldata);
level0 = renamevars(level0,1:width(level0),nomes);
level0(1,:)=[];
level0.("datetime(utc)") = datetime(datevec(level0.("datetime(utc)")));
save(strcat(caminho1,tit,'_smps_level0.mat'),'alldata','-mat')
writetable(level0,strcat(caminho1,tit,'_smps_level0.csv'))
%
% Remove maintanance periods from dataset
if(filtmanut)
    for i=1:size(manut,1)
        aux=find(alldata(:,1)>=manut(i,1) & alldata(:,1)<=manut(i,2));
        alldata(aux,:)=[];
    end
    clear aux
    save(strcat(caminho1,tit,'_smps_clean.mat'),'alldata','-mat')
end

%
% Basic plots --------------------------------------------------
% 1: matdate(utc) 
% 2:8 - Dpg(nm), Mode(nm), GeoStd, TotalConc(cm-3), N50(cm-3), SurfArea(um2/cm3), Vol(um3/cm3) 
% 9:113 - dN/dlogDp
% 113:118 - SampleTemp(C), SamplePress(kPa), SampleRH(%), SheathFlow(lpm), AerosolFlow(lpm)
lin=2:size(alldata,1);
figure
subplot(3,1,1),plot(alldata(lin,1),alldata(lin,116),'c*'),grid on
set(gca,'YLim',[20 60]),ylabel('SampleRH(%)')
title(tit)
dynamicDateTicks
subplot(3,1,2),plot(alldata(lin,1),alldata(lin,117),'m*'),grid on
set(gca,'YLim',[5 10]),ylabel('SheathFlow(lpm)')
dynamicDateTicks
subplot(3,1,3),plot(alldata(lin,1),alldata(lin,118),'g*'),grid on
set(gca,'YLim',[0.5 1.7]),ylabel('AerosolFlow(lpm)')
dynamicDateTicks
saveas(gcf,[caminho2,tit,'_smps_diag'],'jpg')
%
disp('Calculations for contourplots')
aux=find(diff(alldata(:,1))>=0.5);
alldata2(1,:)=alldata(1,:);
for i=2:length(aux) % insere linhas em branco em períodos com mais de 12h sem dados
    alldata2=cat(1,alldata2,alldata(aux(i-1)+1:aux(i),:),nan(1,size(alldata,2)));
end
alldata2=cat(1,alldata2,alldata(aux(length(aux))+1:size(alldata,1),:));
aux=find(isnan(alldata2(:,1)));
alldata2(aux,1)=alldata2(aux-1,1)+0.4;
timesec=alldata2(lin,1);
Z=alldata2(lin,9:113);
aux=find(diff(timesec)<=0);
timesec(aux)=[];
Z(aux,:)=[];
[r,c]=find(Z<0);  % replace negative values per 0
for k=1:length(r)
    Z(r(k),c(k))=0;
end
clear r c
Znorm=Z;
for ii=1:size(Z,1)
    Znorm(ii,:)=Z(ii,:)./max(Z(ii,:));
end
figure,subplot(3,1,1)
contourf(timesec,diam,Znorm','LineStyle','none'),grid on
set(gca,'YScale','log')
set(gca,'YTick',[20 50 100 200 400])
% imagesc(timesec,log10(diam),Znorm')
% set(gca,'YDir','normal')
% set(gca,'YTick',[1.0 log10(30) log10(50) 2.0 log10(200) log10(400)])
% set(gca,'YTickLabel',[10 30 50 100 200 400],'FontSize',12)
legend('normalized')%,colorbar 
title(tit,'FontSize',15),ylabel('Diameter (nm)','FontSize',12)
%dynamicDateTicks
set(gca,'XLim',[datenum(2023,09,11) floor(alldata(size(alldata,1),1))])
set(gca,'XTick',datenum(2023,09,11):7:floor(alldata(size(alldata,1),1)))
set(gca,'XTickLabel', datestr(datenum(2023,09,11):7:floor(alldata(size(alldata,1),1)),'mm/dd/yy'),'FontSize',12)
subplot(3,1,2),plot(alldata(lin,1),alldata(lin,5),'ko'),grid on
ylabel('Total Conc (cm^-^3)');
set(gca,'YLim',[0 10^5]);
set(gca,'XLim',[datenum(2023,09,11) floor(alldata(size(alldata,1),1))])
set(gca,'XTick',datenum(2023,09,11):7:floor(alldata(size(alldata,1),1)))
set(gca,'XTickLabel', datestr(datenum(2023,09,11):7:floor(alldata(size(alldata,1),1)),'mm/dd/yy'),'FontSize',12)
subplot(3,1,3),plot(alldata(lin,1),alldata(lin,8),'b.'), grid on
hold on, plot(alldata(lin,1),alldata(lin,2),'m.')
ylabel('Vol (um^3/cm^-^3); Dpg (nm)');
legend('Vol','Dpg')
set(gca,'XLim',[datenum(2023,09,11) floor(alldata(size(alldata,1),1))])
set(gca,'XTick',datenum(2023,09,11):7:floor(alldata(size(alldata,1),1)))
set(gca,'XTickLabel', datestr(datenum(2023,09,11):7:floor(alldata(size(alldata,1),1)),'mm/dd/yy'),'FontSize',12)
xlabel('Local Time','FontSize',12)
saveas(gcf,strcat(caminho2,tit,'_contour_tseries'),'jpg')
clear alldata2
%
figure
subplot(2,1,1),plot(alldata(lin,1),alldata(lin,5),'ko'),grid on
ylabel('Total Conc (cm^-^3)');
set(gca,'YLim',[0 10^5]);
set(gca,'XLim',[datenum(2023,09,11) floor(alldata(size(alldata,1),1))])
set(gca,'XTick',datenum(2023,09,11):7:floor(alldata(size(alldata,1),1)))
set(gca,'XTickLabel', datestr(datenum(2023,09,11):7:floor(alldata(size(alldata,1),1)),'mm/dd/yy'),'FontSize',12)
subplot(2,1,2),plot(alldata(lin,1),alldata(lin,2),'mo'), grid on
ylabel('Dpg (nm)');
set(gca,'XLim',[datenum(2023,09,11) floor(alldata(size(alldata,1),1))])
set(gca,'XTick',datenum(2023,09,11):7:floor(alldata(size(alldata,1),1)))
set(gca,'XTickLabel', datestr(datenum(2023,09,11):7:floor(alldata(size(alldata,1),1)),'mm/dd/yy'),'FontSize',12)
xlabel('Local Time','FontSize',12)
saveas(gcf,strcat(caminho2,tit,'_N_Dpg_tseries'),'jpg')
%
mediana=zeros(3,length(9:113));    
for ii=9:113
    mediana(1,ii-7)=prctile(alldata(:,ii),25);
    mediana(2,ii-7)=prctile(alldata(:,ii),50);
    mediana(3,ii-7)=prctile(alldata(:,ii),75);
end
figure
errorbar(diam,mediana(2,1:105),mediana(1,1:105),mediana(3,1:105),'-b.')
title(strcat(tit,'Median Size Distribution'),'FontSize',12,'fontweight','b'),xlabel('Dp (nm)','FontSize',12,'fontweight','b'),ylabel('dN/dlogDp (cm^-^3)','FontSize',12,'fontweight','b')
set(gca,'XScale','log'), grid on
set(gca,'XTick',[20 50 100 200 400],'FontSize',12,'fontweight','b')
set(gca,'XLim',[10 500])
saveas(gcf,strcat(caminho2,'\median_sizedist'),'jpg')
%
% Ciclo diurno médio N e Dpg
first = 0;
last = 1;
delta = datenum(0,0,0,0,30,0);
tol = 0;
dados = cat(2,mod(alldata(:,1),1),alldata(:,2),alldata(:,5:6));
dados(1,:)=[];
[~,statorder] = fmedia2(dados,first,last,delta,tol);
figure
subplot(3,1,1) 
errorbar(statorder(:,1)*24,statorder(:,9),statorder(:,9)-statorder(:,12),statorder(:,15)-statorder(:,9),'-b.')
set(gca,'XLim',[0 24],'XTick',0:3:24), grid on
ylabel('Total Conc (cm-3)')
title([tit, ' Median Diurnal Cycle'])
subplot(3,1,2) 
errorbar(statorder(:,1)*24,statorder(:,10),statorder(:,10)-statorder(:,13),statorder(:,16)-statorder(:,10),'-c.')
set(gca,'XLim',[0 24],'XTick',0:3:24), grid on
ylabel('N50 (cm-3)')
subplot(3,1,3) 
errorbar(statorder(:,1)*24,statorder(:,8),statorder(:,8)-statorder(:,11),statorder(:,14)-statorder(:,8),'-m.')
set(gca,'XLim',[0 24],'XTick',0:3:24), grid on
ylabel('Dpg (nm)'), xlabel('UTC Time')
saveas(gcf,strcat(caminho2,'\median_diurnalcycle'),'jpg')

%
% ------------ Level 1 ------------------------------------------------
% Level1: 1-min aerosol concentration must be >50 cm-3 and <100 000 cm-3. 
%         Sheath flow must be > 3 lpm and < 12 lpm. 
%         Sample flow must be > 0.3 lpm and < 2 lpm. 
%         Status flag must be “Normal Scan”. 
%         Multiple charge correction must have been applied.
% alldata:
% 1: matdate(utc) 
% 2:8 - Dpg(nm), Mode(nm), GeoStd, TotalConc(cm-3), N50(cm-3), SurfArea(um2/cm3), Vol(um3/cm3) 
% 9:113 - dN/dlogDp
% 113:118 - SampleTemp(C), SamplePress(kPa), SampleRH(%), SheathFlow(lpm), AerosolFlow(lpm)
%
if(lev1)   % Level1 --------------------------------------------------    
    alldata_level0=alldata;
    diam=alldata(1,9:113);
    %
    origsize=size(alldata_level0,1);
    alldata(alldata(:,5)<50,2:size(alldata,2))=NaN; % remove lines with totconc<50cm-3
    alldata(alldata(:,5)>10^5,2:size(alldata,2))=NaN; % remove lines with totconc>10^5cm-3
    alldata(alldata(:,116)<3,2:size(alldata,2))=NaN; % remove lines with sheath<3 lpm
    alldata(alldata(:,116)>12,2:size(alldata,2))=NaN; % remove lines with sheath>12 lpm
    alldata(alldata(:,117)<0.3,2:size(alldata,2))=NaN; % remove lines with sampleflow<0.3 lpm
    alldata(alldata(:,117)>2,2:size(alldata,2))=NaN; % remove lines with sampleflow>2 lpm
    %
    aux=2:length(alldata);
    figure,plot(alldata_level0(aux,1),alldata_level0(aux,5),'r.'),grid on, hold on
    plot(alldata(aux,1),alldata(aux,5),'b.')
    dynamicDateTicks
    %   
    % Save Level 1 data 
    % Averages 5 and 30 min 
    first = datenum(2023,9,11);
    last = datenum(2024,01,26);
    delta = datenum(0,0,0,0,5,0);
    tol = 0;
    [~,statorder] = fmedia2(alldata,first,last,delta,tol);
    dados5 = statorder(:,1:size(alldata,2));
    delta = datenum(0,0,0,0,30,0);
    [~,statorder] = fmedia2(alldata,first,last,delta,tol);
    dados30 = statorder(:,1:size(alldata,2));

    nomes = ["datetime(utc)", "Dpg(nm)", "Mode(nm)", "GeoStd", "TotConc(cm-3)", "N50(cm-3)", "SurfArea(um2/cm3)", "Vol(um3/cm3)", ...
    string(diam), "SampleTemp(C)", "SamplePress(kPa)", "SampleRH(%)", "SheathFlow(lpm)", "AerosolFlow(lpm)"];
    level1_5min = array2table(dados5);
    level1_30min = array2table(dados30);
    level1_5min = renamevars(level1_5min,1:width(level1_5min),nomes);
    level1_30min = renamevars(level1_30min,1:width(level1_30min),nomes);
    level1_5min.("datetime(utc)") = datetime(datevec(level1_5min.("datetime(utc)")));
    level1_30min.("datetime(utc)") = datetime(datevec(level1_30min.("datetime(utc)")));
    writetable(level1_5min,strcat(caminho1,tit,'_smps_level1_5min.csv'))
    writetable(level1_30min,strcat(caminho1,tit,'_smps_level1_30min.csv'))

%     Instrument='SMPS_TSI';
%     Model='3082';
%     serial_num='3082001629001';
%     Site=tit;
%     Val_level='1';
%     Proc_vers='1.0';
%     labels{1}='Total concentration';
%     labels{2}='Total volume';
%     labels{3}='Mean diameter';
%     units{1}='(cm-3)';
%     units{2}='(um3.cm-3)';
%     units{3}='(nm)';
%     export_LFA_GoAmazon_jan2015(labels,units,Instrument,Model,serial_num,Site,Val_level,Proc_vers,origsize,alldata(aux,1),[alldata(aux,5) alldata(aux,7) alldata(aux,2)],diam,alldata(aux,9:113));    
end
%
% ------------ Level 2 ------------------------------------------------
% Level2: Apply STP correction (1013.25; 293.15).
if(lev2)        
% % %     %
% % %     % 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));
% % %     for j=2:109
% % %         if(j~=4) % skip Dpg                 
% % %             level1_5min(:,j)=level1_5min(:,j).*stp;
% % %         end
% % %     end
% % %     %   
% % %     % Save Level 2 data 
% % %     Instrument='SMPS_TSI';
% % %     Model='3082';
% % %     serial_num='3082001629001';
% % %     Site=tit;
% % %     Val_level='2';
% % %     Proc_vers='1.0';
% % %     labels{1}='Total concentration';
% % %     labels{2}='Total volume';
% % %     labels{3}='Mean diameter';
% % %     units{1}='(cm-3)';
% % %     units{2}='(um3.cm-3)';
% % %     units{3}='(nm)';
% % %     export_LFA_GoAmazon_jan2015(labels,units,Instrument,Model,serial_num,Site,Val_level,Proc_vers,origsize,level1_5min(:,1),level1_5min(:,2:4),diam,level1_5min(:,5:109));    
end
%
%
% ------------ Level 3 ------------------------------------------------
% Level3: Integrated aerosol number must be within 30% CPC measurement. 
%         Considering an aerosol density of 1.4 g cm-3, integrated mass must be within 50%-120% of TEOM PM2.5.    
%
% if(lev3)        
%     
% end
% -------------------------------------------------------------------------
