% SMPS-TSI-3082 (version3) data processing from txt files (LRizzo)
% Interpolated to a commom diameter sequence
% Funcoes utilizadas: fsmps_interp, fsmps_param, fmediaqualquercoisa2.m
%
% Output: alldata(level0) 
% 1 - matdate(utc) 
% 2:7 - Dpg(nm), Mode(nm), GeoStd, TotalConc(cm-3), SurfArea(um2/cm3), Vol(um3/cm3) 
% 8:112 - dN/dlogDp
% 113:119 - SampleTemp(C), SamplePress(kPa), SampleRH(%), SheathFlow(lpm), AerosolFlow(lpm), LowerSize(nm), UpperSize(nm)
%
% 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:
caminho='C:\A\Diadema2016\data_trailer_2016\smps3082\txt\';
caminho1='C:\A\Diadema2016\processados\smps\';
caminho2='C:\A\Diadema2016\processados\smps\plots\';
tit='Diadema-2016';
dens=1.6; % particle density, in g/cm3
ano=2016;
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	25.9	26.9	27.9	28.9	30	31.1	32.2	33.4	34.6	35.9	37.2	38.5	40	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	73.7	76.4	79.1	82	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	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	278.8	289	299.6	310.6	322	333.8	346	358.7	371.8	385.4	399.5	414.2	429.4];
filtmanut=true; % 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',42,'Delimiter',',');
    if(strcmp(config{1,1}(26),'FALSE'))
        disp('Warning: no multiple charge correction!!!')
    end
    aux=textscan(fid,'%s',31,'Delimiter',',');
    % verifica se chegou no início dos diâmetros
    if~(strcmp(aux{1,1}(31),'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,'dd/mm/yyyyHH: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))
%
alldata=cat(2,tempo,totconc,maisparam(:,4:5),sizedist,param);
alldata=sortrows(alldata,1); % sort data based on day of year
alldata(1,8:length(diam)+7)=diam;
save(strcat(caminho1,tit,'_smps_level0.mat'),'alldata','-mat')
csvwrite(strcat(caminho1,tit,'_smps_level0.csv'),[alldata,datevec(alldata(:,1))])
%
% 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')
    %csvwrite(strcat(caminho1,tit,'_smps_clean.csv'),[alldata,datevec(alldata(:,1))])
end


%
% Basic plots
% 1 - matdate(utc) 
% 2:7 - Dpg(nm), Mode(nm), GeoStd, TotalConc(cm-3), SurfArea(um2/cm3), Vol(um3/cm3) 
% 8:112 - dN/dlogDp
% 113:119 - SampleTemp(C), SamplePress(kPa), SampleRH(%), SheathFlow(lpm), AerosolFlow(lpm), LowerSize(nm), UpperSize(nm)
lin=2:size(alldata,1);
figure
subplot(3,1,1),plot(alldata(lin,1),alldata(lin,115),'c*'),grid on
set(gca,'YLim',[20 60]),ylabel('SampleRH(%)')
title(tit)
dynamicDateTicks
subplot(3,1,2),plot(alldata(lin,1),alldata(lin,116),'m*'),grid on
set(gca,'YLim',[5 10]),ylabel('SheathFlow(lpm)')
dynamicDateTicks
subplot(3,1,3),plot(alldata(lin,1),alldata(lin,117),'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,8:112);
[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(2,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),xlabel('Local Time','FontSize',12),ylabel('Diameter (nm)','FontSize',12)
set(gca,'XLim',[datenum(2016,11,22) floor(alldata(size(alldata,1),1))])
set(gca,'XTick',datenum(2016,11,22):7:floor(alldata(size(alldata,1),1)))
set(gca,'XTickLabel', datestr(datenum(2016,11,22):7:floor(alldata(size(alldata,1),1)),'mm/dd/yy'),'FontSize',12)
subplot(2,1,2),plot(alldata(lin,1),alldata(lin,5),'ko'),grid on,hold on
addaxis(alldata(lin,1),alldata(lin,7)*dens,'b.')
set(gca,'XLim',[datenum(2016,11,22) floor(alldata(size(alldata,1),1))])
set(gca,'XTick',datenum(2016,11,22):7:floor(alldata(size(alldata,1),1)))
set(gca,'XTickLabel', datestr(datenum(2016,11,22):7:floor(alldata(size(alldata,1),1)),'mm/dd/yy'),'FontSize',12)
xlabel('Local Time','FontSize',12)
addaxislabel(1,'Total Conc (cm^-^3)');
addaxislabel(2,'Mass (ug/m^-^3)');
saveas(gcf,strcat(caminho2,tit,'_contour_tseries'),'jpg')
%
mediana=zeros(3,length(8:112));    
for ii=8:112
    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')
%
% ------------ 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:7 - Dpg(nm), Mode(nm), GeoStd, TotalConc(cm-3), SurfArea(um2/cm3), Vol(um3/cm3) 
% 8:112 - dN/dlogDp
% 113:119 - SampleTemp(C), SamplePress(kPa), SampleRH(%), SheathFlow(lpm), AerosolFlow(lpm), LowerSize(nm), UpperSize(nm)
%
if(lev1)       
    C=load('C:\A\Diadema2016\processados\smps\Diadema-2016_smps_clean.mat');
    alldata=C.alldata;
    C=load('C:\A\Diadema2016\processados\smps\Diadema-2016_smps_level0.mat');
    alldata_level0=C.alldata;
    diam=C.alldata(1,8:112);
    tit='Diadema-2016';
    %
    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 
    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,8:112));    
end
%
% ------------ Level 2 ------------------------------------------------
% Level2: STP correction (1013.25; 293.15).
if(lev2)        
    % Lê os arquivos Level 1 que ja salvou
    aux=importdata('C:\A\Diadema2016\processados\LFA_standard\Diadema-2016_SMPS_TSI_3082_Level1_5min.csv');
    diam=zeros(1,105);
    for j=24:128
        diam(j-23)=str2num(char(aux.textdata(2,j)));
    end
    level1_5min=[aux.data(:,1:2) aux.data(:,9) aux.data(:,16) aux.data(:,23:127)];        
    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));
    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.    
%         Obs:
%         Removido período entre 08fev 9:00 e 15fev 19:00 (high aerosol flow, low pressure, inlet entupido)  
%         Testes em laboratório com esferas de PSL após o experimento indicam que o 
%         sizing está correto, mas a contagem do SMPS está ~50% subestimada.
%         As distribuições dN/logDp foram divididas pelo fator constante de 0,58
%         (slope SMPSxCPC). Porém, mesmo após essa correção, as concentrações
%         do SMPS devem ser interpretadas apenas qualitativamente.
%
if(lev3)        
    dens=1.6; %g/cm3
    % Lê os arquivos Level 2 que ja salvou
    aux=importdata('C:\A\Diadema2016\processados\LFA_standard\Diadema-2016_SMPS_TSI_3082_Level2_30min.csv');
    diam=aux.data(1,:);
    aux.data(length(aux.data)+1:length(aux.textdata)-1,:)=NaN;
    sdist=aux.data(2:length(aux.data),:);  
    param=str2double(aux.textdata(3:length(aux.textdata),2:5)); % TotConc(cm-3); TotVol(cm3); Dpg(nm)
    param(:,1)=param(:,1)/(24*60*60)+datenum(1970,1,1); % unix time to matlab time
    origsize=sum(~isnan(param(:,2)));
    tit='Diadema-2016';
    %
    aux=load('C:\A\Diadema2016\processados\smps\Diadema-2016_smps_clean.mat');
    level0=[aux.alldata(:,1) aux.alldata(:,113:117)]; %SampleTemp(C), SamplePress(kPa), SampleRH(%), SheathFlow(lpm), AerosolFlow(lpm)
    level0(1,:)=[];
    %
    aux=importdata('C:\A\Diadema2016\processados\LFA_standard\Diadema-2016_CPC_3772_Level2_30min.csv');
    cpc=aux.data;
    cpc(:,1)=cpc(:,1)/(24*60*60)+datenum(1970,1,1); % unix time to matlab time
    %
    aux=importdata('C:\A\Diadema2016\processados\LFA_standard\Diadema-2016_ACSM_dados preliminares_Level1_30min.csv');
    acsm=aux.data; % {'Org','SO4','NO3','NH4','Chl','OOA','HOA','BBOA'}
    acsm(:,1)=acsm(:,1)/(24*60*60)+datenum(1970,1,1); % unix time to matlab time
    %
    aux=importdata('C:\A\Diadema2016\processados\LFA_standard\Diadema-2016_MAAP_5012_Level2_30min.csv');
    maap=aux.data(:,1:2); % BC
    maap(:,1)=maap(:,1)/(24*60*60)+datenum(1970,1,1); % unix time to matlab time
    %
    aux=importdata('C:\A\Diadema2016\processados\LFA_standard\Diadema-2016_OPS_TSI_3330_Level2_30min.csv');
    aux.data(length(aux.data)+1:length(aux.textdata)-1,:)=NaN;
    dops=aux.data(1,:)*1000; % nm
    ops(:,2:17)=aux.data(2:length(aux.data),:);
    for i=3:length(aux.textdata)
        ops(i-2,1)=str2double(char(aux.textdata{i,2}))/(24*60*60)+datenum(1970,1,1); % unix time to matlab time;
    end       
    clear aux
    %    
    msdist=param(:,3)*dens; % ug/m3
    macsmmaap=sum(acsm(:,2:6),2)+maap(:,2); % ug/m3    
    %
    rem=[datenum(2017,2,8,9,0,0) datenum(2017,2,15,19,0,0)]; 
    aux=find(param(:,1)>=rem(1,1) & param(:,1)<=rem(1,2));
    %
    figure
    subplot(2,1,1),plot(param(:,1),cpc(:,2)./param(:,2),'b.'),grid on,hold on,dynamicDateTicks
    ylabel('Ratio','FontSize',12),title(tit,'FontSize',12)
    plot(param(:,1),macsmmaap./msdist,'r.')
    legend('CPC/SMPS','ACSM+MAAP/SMPS')
    set(gca,'XLim',[datenum(2016,12,1) datenum(2017,3,16)])
    subplot(2,1,2),plot(level0(:,1),level0(:,2),'r.'),grid on,dynamicDateTicks
    ylabel('SampleTemp(C)','FontSize',12)
    set(gca,'XLim',[datenum(2016,12,1) datenum(2017,3,16)])
    figure
    subplot(2,1,1),plot(param(:,1),cpc(:,2)./param(:,2),'b.'),grid on,hold on,dynamicDateTicks
    ylabel('Ratio','FontSize',12),title(tit,'FontSize',12)
    plot(param(:,1),macsmmaap./msdist,'r.')
    legend('CPC/SMPS','ACSM+MAAP/SMPS')
    set(gca,'XLim',[datenum(2016,12,1) datenum(2017,3,16)])
    subplot(2,1,2),plot(level0(:,1),level0(:,3),'g.'),grid on,dynamicDateTicks
    ylabel('SamplePress(kPa)','FontSize',12)
    set(gca,'XLim',[datenum(2016,12,1) datenum(2017,3,16)])
    figure
    subplot(2,1,1),plot(param(:,1),cpc(:,2)./param(:,2),'b.'),grid on,hold on,dynamicDateTicks
    ylabel('Ratio','FontSize',12),title(tit,'FontSize',12)
    plot(param(:,1),macsmmaap./msdist,'r.')
    legend('CPC/SMPS','ACSM+MAAP/SMPS')
    set(gca,'XLim',[datenum(2016,12,1) datenum(2017,3,16)])
    subplot(2,1,2),plot(level0(:,1),level0(:,4),'c.'),grid on,dynamicDateTicks
    ylabel('SampleRH(%)','FontSize',12)
    set(gca,'XLim',[datenum(2016,12,1) datenum(2017,3,16)])
    figure
    subplot(2,1,1),plot(param(:,1),cpc(:,2)./param(:,2),'b.'),grid on,hold on,dynamicDateTicks
    ylabel('Ratio','FontSize',12),title(tit,'FontSize',12)
    plot(param(:,1),macsmmaap./msdist,'r.')
    legend('CPC/SMPS','ACSM+MAAP/SMPS')
    set(gca,'XLim',[datenum(2016,12,1) datenum(2017,3,16)])
    subplot(2,1,2),plot(level0(:,1),level0(:,5),'c.'),grid on,dynamicDateTicks
    ylabel('SheathFlow(lpm)','FontSize',12)
    set(gca,'XLim',[datenum(2016,12,1) datenum(2017,3,16)])
    set(gca,'YLim',[5 7])
    figure
    subplot(2,1,1),plot(param(:,1),cpc(:,2)./param(:,2),'b.'),grid on,hold on,dynamicDateTicks
    ylabel('Ratio','FontSize',12),title(tit,'FontSize',12)
    plot(param(:,1),macsmmaap./msdist,'r.')
    legend('CPC/SMPS','ACSM+MAAP/SMPS')
    set(gca,'XLim',[datenum(2016,12,1) datenum(2017,3,16)])
    subplot(2,1,2),plot(level0(:,1),level0(:,6),'b.'),grid on,dynamicDateTicks
    ylabel('AerosolFlow(lpm)','FontSize',12)
    set(gca,'XLim',[datenum(2016,12,1) datenum(2017,3,16)])
    figure
    subplot(2,1,1),plot(param(:,1),cpc(:,2)./param(:,2),'b.'),grid on,hold on,dynamicDateTicks
    ylabel('Ratio','FontSize',12),title(tit,'FontSize',12)
    plot(param(:,1),macsmmaap./msdist,'r.')
    legend('CPC/SMPS','ACSM+MAAP/SMPS')
    set(gca,'XLim',[datenum(2016,12,1) datenum(2017,3,16)])
    subplot(2,1,2),plot(param(:,1),param(:,4),'c.'),grid on,dynamicDateTicks
    hold on,  plot(param(aux,1),param(aux,4),'rx')
    ylabel('Dpg(nm)','FontSize',12)
    set(gca,'XLim',[datenum(2016,12,1) datenum(2017,3,16)])
    figure
    subplot(2,1,1),plot(param(:,1),cpc(:,2)./param(:,2),'b.'),grid on,hold on,dynamicDateTicks
    ylabel('Ratio','FontSize',12),title(tit,'FontSize',12)
    plot(param(:,1),macsmmaap./msdist,'r.')
    legend('CPC/SMPS','ACSM+MAAP/SMPS')
    set(gca,'XLim',[datenum(2016,12,1) datenum(2017,3,16)])
    subplot(2,1,2),plot(param(:,1),ops(:,3)./sdist(:,101),'k.'),grid on,dynamicDateTicks
    hold on,  plot(param(aux,1),ops(aux,3)./sdist(aux,101),'rx')
    ylabel('OPS/SMPS @370nm ','FontSize',12)
    set(gca,'XLim',[datenum(2016,12,1) datenum(2017,3,16)])
    %
    paramorig=param;
    param(param(:,1)>=rem(1,1) & param(:,1)<=rem(1,2),2:4)=NaN;
    msdist(param(:,1)>=rem(1,1) & param(:,1)<=rem(1,2))=NaN;
    param(:,2:3)=param(:,2:3)/0.58; % correção utilizando o slope com o cpc
    msdist=param(:,3)*dens; % ug/m3
    figure
    subplot(1,2,1),plot(cpc(:,2),param(:,2),'b.'),grid on,hold on
    [a,b,rsquare]=flinreg(cpc(:,2),param(:,2));
    x=[min(cpc(:,2)) max(cpc(:,2))];
    plot(x,a*x+b,'-r')
    legend('data (dens 1.6 g/cm3)',[num2str(a,'%1.2f'),'*x+',num2str(b,'%3.0f'),' R^2=',num2str(rsquare,'%1.2f')])
    ylabel('ConcSMPS (cm^-^3)','FontSize',12),xlabel('CPC (cm^-^3)','FontSize',12)
    title(tit,'FontSize',12)
    subplot(1,2,2),plot(macsmmaap,msdist,'c.'),grid on,hold on
    [a,b,rsquare]=flinreg(macsmmaap,msdist);
    x=[min(macsmmaap) max(macsmmaap)];
    plot(x,a*x+b,'-r')
    legend('data',[num2str(a,'%1.2f'),'*x+',num2str(b,'%3.0f'),' R^2=',num2str(rsquare,'%1.2f')])
    ylabel('ConcSMPS (cm^-^3)','FontSize',12),xlabel('ACSM+BC (ug/m^3)','FontSize',12)
    %
    % Lê os arquivos Level 2 5min para corrigir e salvar o Level3
    aux=importdata('C:\A\Diadema2016\processados\LFA_standard\Diadema-2016_SMPS_TSI_3082_Level2_5min.csv');
    diam=aux.data(1,:);
    aux.data(length(aux.data)+1:length(aux.textdata)-1,:)=NaN;
    sdist=aux.data(2:length(aux.data),:);  
    param=str2double(aux.textdata(3:length(aux.textdata),2:5)); % TotConc(cm-3); TotVol(cm3); Dpg(nm)
    param(:,1)=param(:,1)/(24*60*60)+datenum(1970,1,1); % unix time to matlab time
    origsize=sum(~isnan(param(:,2)));
    tit='Diadema-2016';
    rem=[datenum(2017,2,8,9,0,0) datenum(2017,2,15,19,0,0)];
    param(param(:,1)>=rem(1,1) & param(:,1)<=rem(1,2),2:4)=NaN;
    sdist(param(:,1)>=rem(1,1) & param(:,1)<=rem(1,2),:)=NaN;
    param(:,2:3)=param(:,2:3)/0.58; % correção utilizando o slope com o cpc
    sdist=sdist/0.58;
    figure
    subplot(2,1,1),plot(param(:,1),param(:,2),'-b.'),grid on,dynamicDateTicks
    ylabel('SMPS Conc','FontSize',12),title([tit,'Após correção'],'FontSize',12)
    set(gca,'XLim',[datenum(2016,12,1) datenum(2017,3,16)])
    subplot(2,1,2),plot(param(:,1),param(:,4),'-c.'),grid on,dynamicDateTicks
    ylabel('Dpg(nm)','FontSize',12)
    set(gca,'XLim',[datenum(2016,12,1) datenum(2017,3,16)])
    figure
    cor='brgkcmbrgkcmbrgkcm';
    delta=10;
    for i=31969:delta:length(sdist)
        if(sum(isnan(sdist(i:i+9,5)))<10)
            for ii=1:delta
                semilogx(diam,sdist(i+ii-1,:),['-',cor(ii),'.']),grid on,hold on
                set(gca,'XLim',[10 500])
                title(datestr(param(i,1)))
            end
            pause%(0.5)
        end
        hold off
    end
    %
    % Save Level 3 data 
    Instrument='SMPS_TSI';
    Model='3082';
    serial_num='3082001629001 - Sizing quantitative, Concentration qualitative';
    Site=tit;
    Val_level='3';
    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,param(:,1),param(:,2:4),diam,sdist);    
    %
end
% -------------------------------------------------------------------------
% %   Basic statistics 
%     conc=cat(2,alldata(:,1),alldata(:,8:11));
%     sizedist=alldata(:,8:size(alldata,2));
%     sizedist=cat(2,alldata(:,1),sizedist);
%     sizedist(1,:)=[];
% %
% %   Averages considering LAST 30min  
%     disp('Calculating 30min averages')
%     first = conc(2,1)-mod(conc(2,1),1) + (0.5/24);   % first averaging time
%     last = conc(length(conc),1)-mod(conc(length(conc),1),1) + 1 + (0/24);  % last averaging time
%     delta = 30/60/24;     % 30min
%     [aux1,aux2]=fmedia_qualquer_coisa(conc,first,last,delta);
%     avg30min=aux1;  % (time | avg values | std values | #)
%     aux3=find(avg30min(:,6)==0);
%     avg30min(aux3,2:6)=NaN;
%     clear aux1 aux2 aux3
%     save([caminho1,tit,'totconc_avg30min.txt'],'avg30min','-ASCII')
%     [aux1,aux2]=fmedia_qualquer_coisa(sizedist,first,last,delta);
%     sdavg30min=aux1;  % (time | avg values | std values | #)
%     aux3=find(avg30min(:,6)==0);
%     sdavg30min(aux3,2:6)=NaN;
%     clear aux1 aux2 aux3
%     sdavg30min=cat(1,zeros(1,size(sdavg30min,2)),sdavg30min);
%     sdavg30min(1,2:length(diam)+1)=diam;
%     save([caminho1,tit,'sizedist_avg30min.txt'],'sdavg30min','-ASCII')
%     %
%     mediana=zeros(6,size(sizedist,2)-1);    
%     for ii=2:size(sizedist,2)
%         mediana(1,ii-1)=prctile(sizedist(:,ii),25);
%         mediana(2,ii-1)=prctile(sizedist(:,ii),50);
%         mediana(3,ii-1)=prctile(sizedist(:,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')
%     %
%     baixa=find(alldata(:,11)<1000);
%     alta=find(alldata(:,11)>5000);
%     for ii=2:size(sizedist,2)
%         mediana(1,ii-1)=prctile(sizedist(baixa,ii),25);
%         mediana(2,ii-1)=prctile(sizedist(baixa,ii),50);
%         mediana(3,ii-1)=prctile(sizedist(baixa,ii),75);
%         mediana(4,ii-1)=prctile(sizedist(alta,ii),25);
%         mediana(5,ii-1)=prctile(sizedist(alta,ii),50);
%         mediana(6,ii-1)=prctile(sizedist(alta,ii),75);
%     end
%     figure
%     subplot(1,2,1),plot(diam,mediana(2,:),'-b.')
%     title(strcat(tit,'Median Size Dist. (N<1000cm^-^3)'),'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])
%     subplot(1,2,2),plot(diam,mediana(5,:),'-r.')
%     title(strcat(tit,'Median Size Distribution (N>5000cm^-^3)'),'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,'\alta_baixa_sizedist'),'jpg')
%     
%     %
%     figure
%     subplot(2,1,1),plot(conc(2:size(conc,1),1),conc(2:size(conc,1),5),'co')
%     title([tit],'FontSize',12),ylabel('SMPS total conc [cm^-^3]','FontSize',12),grid on,axis tight
%     d1=conc(2,1); d2=conc(length(conc),1);
%     set(gca,'XTick',floor(d1):floor((d2-d1)/2):floor(d2))
%     set(gca,'XTickLabel',datestr(floor(d1):floor((d2-d1)/2):floor(d2),'dd-mmm'))
%     subplot(2,1,2),semilogy(conc(2:size(conc,1),1),conc(2:size(conc,1),3),'ro')
%     ylabel('SMPS peak [nm]','FontSize',12),grid on,axis tight    
%     set(gca,'XTick',floor(d1):floor((d2-d1)/2):floor(d2))
%     set(gca,'XTickLabel',datestr(floor(d1):floor((d2-d1)/2):floor(d2),'dd-mmm'))
%     xlabel('UTC Time','FontSize',12)
%     saveas(gcf,strcat(caminho1,tit,'totconc_timeseries'),'jpg')
%     %
%     figure
%     subplot(2,1,1),plot(mod(conc(2:size(conc,1),1),1),conc(2:size(conc,1),5),'co')
%     title([tit],'FontSize',12),ylabel('SMPS total conc [cm^-^3]','FontSize',12),grid on,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'})
%     subplot(2,1,2),semilogy(mod(conc(2:size(conc,1),1),1),conc(2:size(conc,1),3),'ro')
%     ylabel('SMPS peak [nm]','FontSize',12),grid on,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'})
%     xlabel('UTC Time','FontSize',12)
%     saveas(gcf,strcat(caminho1,tit,'totconc_dcycle'),'jpg')
%     %
% % ---------------------------------------------------------------
% % Contourplots - daily
% year=alldata(2,2);month=alldata(2,3);iday=alldata(2,4); % Start date
% day=iday-1;
% clear tempo 
% tempo=alldata(2:size(alldata,1),1); 
% sizedist(:,1)=[];
% nday=tempo-mod(tempo,1);
% figure
% for i=nday(1):nday(length(nday))
%     aux=find(nday==i);
%     if (isempty(aux))
%         day=day+1;
%         continue
%     end
%     timesec=tempo(aux);
%     Z=sizedist(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
%     %
%     day=day+1;
%     if(day>lastday(month))
%         day=day-lastday(month);
%         month=month+1;
%     end
%     %
%     % figure
%     contourf(mod(timesec,1),diam,Z),shading flat
%     colorbar
%     title(strcat(tit,num2str(year),num2str(month,'%02d'),num2str(day,'%02d')),'FontSize',12),xlabel('UTC Time','FontSize',12),ylabel('Dp (um)','FontSize',12)
%     set(gca,'YScale','log'),grid on
%     set(gca,'YTick',[20 50 100 200 400])
%     set(gca,'XLim',[0 1])
%     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'},'FontSize',12)    
%     saveas(gcf,strcat(caminho2,'\',num2str(year),num2str(month,'%02d'),num2str(day,'%02d'),'_contour'),'jpg')
%     %
%     % figure
% %     contourf(mod(timesec,1),diam,log10(Z)),shading flat
% %     colorbar, legend('log10')
% %     title(strcat(num2str(year),num2str(month,'%02d'),num2str(day,'%02d'),tit),'FontSize',12),xlabel('UTC Time','FontSize',12),ylabel('Dp(um)','FontSize',12)
% %     set(gca,'YScale','log')
% %     set(gca,'YTick',[20 50 100 200 400])
% %     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'},'FontSize',12)
% %     saveas(gcf,strcat(caminho2,'\',num2str(year),num2str(month,'%02d'),num2str(day,'%02d'),'_log'),'jpg')
%     %
%     Znorm=Z;
%     for j=1:size(Z,2)
%         Znorm(:,j)=Z(:,j)./max(Z(:,j));
%     end
%     %  figure
%     contourf(mod(timesec,1),diam,Znorm),shading flat
%     colorbar, legend('normalized')
%     title(strcat(tit,num2str(year),num2str(month,'%02d'),num2str(day,'%02d')),'FontSize',12),xlabel('UTC Time','FontSize',12),ylabel('Dp (um)','FontSize',12)
%     set(gca,'YScale','log'), grid on
%     set(gca,'YTick',[20 50 100 200 400])
%     set(gca,'XLim',[0 1])
%     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'},'FontSize',12)
%     saveas(gcf,strcat(caminho2,'\',num2str(year),num2str(month,'%02d'),num2str(day,'%02d'),'_norm'),'jpg')
% end
