% SMPS-TSI-3082 (version3) data processing from txt files 
% Written by Luciana Rizzo
% 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\GOAMAZON\T1\SMPS_71331078sn_CPC_71025292sn\txt\';
caminho1='C:\A\GOAMAZON\T1\SMPS_71331078sn_CPC_71025292sn\';
caminho2='C:\A\GOAMAZON\T1\SMPS_71331078sn_CPC_71025292sn\plots\';
tit='T1_2015-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=[9.47,9.82, 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,445.1];
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,6);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));
    if ifile<6
        config=textscan(fid,'%s',59,'Delimiter',','); % ajustar num de linhas de acordo com o modelo do SMPS
        if(strcmp(config{1,1}(51),'FALSE')) % ajustar posição de acordo com o modelo do SMPS
            disp('Warning: no multiple charge correction!!!')
            stop
        end
    else
        config=textscan(fid,'%s',57,'Delimiter',','); 
        if(strcmp(config{1,1}(49),'FALSE')) % ajustar posição de acordo com o modelo do SMPS
            disp('Warning: no multiple charge correction!!!')
            stop
        end
    end    
    aux=textscan(fid,'%s',8,'Delimiter',',');
    % verifica se chegou no início dos diâmetros
    if~(strcmp(aux{1,1}(8),'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,90);
    while(cell2mat(aux)>0)
        diam_(idiam)=cell2mat(aux);
        idiam=idiam+1;
        aux=textscan(fid,'%f',1,'Delimiter',',');
    end
    aux=textscan(fid,'%s',26,'Delimiter',','); % ajustar posição de acordo com o modelo do SMPS
    % verifica se chegou no final do cabeçalho
    if~(strcmp(aux{1,1}(26),'Comment'))
        disp('Did not find end of header.')
        stop
    end
    %
    % Adaptar formato de acordo com o modelo do SMPS
    f1='%*f%s%s%f%f%*f%*f%*s';
    f3='%*f	%*f	%*f	%*s	%*f	%f	%f	%*f	%*f	%*f	%f	%f	%*f	%*f	%*f	%*f	%*f	%*f	%f	%f	%f	%f	%*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',',');
    for j=1:size(C,2)
        if(length(C{1,j})<length(C{1,1}))
            disp('Check raw data file for # strings')
            stop
        end
    end
    fclose(fid);
    %
    % Retrieve normalized size distribution
    sizedist_=zeros(size(C{1,1},1),length(diam_));    
    % Adaptar índice da coluna de acordo com o modelo do SMPS
    for j=5:length(diam_)+4  
        sizedist_(:,j-4)=C{1,j};
    end
    % Retrieve control parameters [1-SampleTemp(C), 2-SamplePress(kPa), ...
    ... 3-SheathFlow(lpm), 4-AerosolFlow(lpm), 5-LowerSize(nm), 6-UpperSize(nm)]
    param_=zeros(size(C{1,1},1),6);
    for j=1:2
        param_(:,j)=C{1,j+2};
    end
    jj=3;
    for j=length(diam_)+5: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,'mm/dd/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), Nlt50(cm-3), Ngt100(cm-3), Nlt70)
%
tempo=tempo+3/24; % changing to utc !!! txt files were exported from aim files in a computer with Sao Paulo time zone (UTC-3h)
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
dynamicDateTicks
subplot(2,1,1),plot(alldata(lin,1),alldata(lin,118),'m*'),grid on
set(gca,'YLim',[5 10]),ylabel('SheathFlow(lpm)')
dynamicDateTicks
title(tit)
subplot(2,1,2),plot(alldata(lin,1),alldata(lin,119),'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:115);
[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('UTC','FontSize',12),ylabel('Diameter (nm)','FontSize',12)
set(gca,'XLim',[datenum(2015,12,02) datenum(2016,4,26)])
set(gca,'XTick',datenum(2015,12,02):10:datenum(2016,4,26))
set(gca,'XTickLabel',datestr(datenum(2015,12,02):10:datenum(2016,4,26),'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(2015,12,02) datenum(2016,4,26)])
set(gca,'XTick',datenum(2015,12,02):10:datenum(2016,4,26))
set(gca,'XTickLabel',datestr(datenum(2015,12,02):10:datenum(2016,4,26),'mm/dd/yy'),'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:115));    
for ii=8:115
    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:108),mediana(1,1:108),mediana(3,1:108),'-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:115 - dN/dlogDp
% 116:121 - SampleTemp(C), SamplePress(kPa), SheathFlow(lpm), AerosolFlow(lpm), LowerSize(nm), UpperSize(nm)
%
if(lev1)       
    origsize=size(alldata,1);
    orig=alldata;
    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(:,118)<3,2:size(alldata,2))=NaN; % remove lines with sheath<3 lpm
    alldata(alldata(:,118)>12,2:size(alldata,2))=NaN; % remove lines with sheath>12 lpm
    alldata(alldata(:,119)<0.3,2:size(alldata,2))=NaN; % remove lines with sampleflow<0.3 lpm
    alldata(alldata(:,119)>2,2:size(alldata,2))=NaN; % remove lines with sampleflow>2 lpm
    %
    aux=2:length(alldata);
    figure,plot(orig(aux,1),orig(aux,5),'-r.'),grid on, hold on
    plot(alldata(aux,1),alldata(aux,5),'b.')
    dynamicDateTicks
    %   
    % Save Level 1 data 
    Instrument='SMPS_TSI';
    Model='3081';
    serial_num='71331078, with CPC3772,WCPC3787,CPC3010';
    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:115));    
end
%
% ------------ Level 2 ------------------------------------------------
% Level2: STP correction (1013.25; 293.15).
if(lev2)        
    % Lê os arquivos Level 1 
    aux=importdata('C:\A\GOAMAZON\T1\LFA_standard\T1_2015-2016_SMPS_TSI_3081_Level1_5min.csv');
    diam=zeros(1,108);
    for j=24:131
        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:130)];        
    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;
    tit='T1_2015-2016';
    %
    %    
    % Lê fator stp
    C=importdata('C:\A\GOAMAZON\T1\LFA_standard\T1_2015-2016_AWS_Gill_Level2_5min.csv');
    stp=[C.data(:,1) C.data(:,10)];  
    stp(:,1)=stp(:,1)/(24*60*60)+datenum(1970,1,1); 
    if(length(stp)~=length(level1_5min))
        disp('check length of meteorological data file')
        stop
    end     
    for j=2:size(level1_5min,2)
        if(j~=4) % skip Dpg                 
            level1_5min(:,j)=level1_5min(:,j).*stp(:,2);
        end
    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='SMPS_TSI';
    Model='3081';
    serial_num='71331078, with CPC3772,WCPC3787,CPC3010';
    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:size(level1_5min,2)));    
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)        
    %dens=1.6; %g/cm3
    % Lê os arquivos Level 2 
    aux=importdata('C:\A\GOAMAZON\T1\LFA_standard\T1_2015-2016_SMPS_TSI_3081_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)));
    orig=param;
    tit='T1_2015-2016';   
    caminho2='C:\A\GOAMAZON\T1\SMPS_71331078sn_CPC_71025292sn\plots\';
    %    
    aux=importdata('C:\A\GOAMAZON\T1\LFA_standard\T1_2015-2016_CPC_3776_Level2_5min.csv');
    cpc=aux.data;
    cpc(:,1)=cpc(:,1)/(24*60*60)+datenum(1970,1,1); % unix time to matlab time
    aux=importdata('C:\A\GOAMAZON\T1\LFA_standard\T1_2015-2016_BC_Aethalometer_AE33_Level2_5min.csv');
    bc=[aux.data(:,1) aux.data(:,7)]; % bc880
    bc(:,1)=bc(:,1)/(24*60*60)+datenum(1970,1,1); % unix time to matlab time
    aux=importdata('C:\A\GOAMAZON\T1\LFA_standard\T1_2015-2016_Nephelometer_Aurora_Ecotech_2000_Level3_5min.csv');   
    scat=aux.data(:,1:2); % scat450
    scat(:,1)=scat(:,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
%     msdist=param(:,3)*dens; % ug/m3
%     macsmmaap=sum(acsm(:,2:6),2)+bc(:,2); % ug/m3    
%
    figure
    subplot(4,1,1),plot(param(:,1),param(:,2),'b.'),grid on,hold on,dynamicDateTicks
    plot(cpc(:,1),cpc(:,2),'r.')
    ylabel('cm^-^3','FontSize',12),title(tit,'FontSize',12)
    legend('SMPS','CPC')%,'ACSM+MAAP/SMPS')
    subplot(4,1,2),plot(param(:,1),param(:,4),'m.'),grid on,dynamicDateTicks
    ylabel('Dpg(nm)','FontSize',12)
    subplot(4,1,3),plot(bc(:,1),bc(:,2),'k.'),grid on,dynamicDateTicks
    ylabel('BC880(ug/m^3)','FontSize',12)
    subplot(4,1,4),plot(scat(:,1),scat(:,2),'c.'),grid on,dynamicDateTicks
    ylabel('scat450(Mm^-^1','FontSize',12)
    % 
    % Correções e exclusões no level3 ----------------------------
    % O período 19mar-07abr foi excluído porque as concentrações estavam relativamente
    % baixas (média~1500cm-3), apesar dos diâmetros parecerem ok. Não há
    % dados do cpc. Scat e BC não mostram diminuições significativas nesse
    % periodo. Como o sizing parece correto, talvez seja possível corrigir
    % o smps fazendo um estudo de closure com espalhamento e absorção
    % utilizando o código Mie.
    paramorig=param;
    rem=[datenum(2015,12,27,00,00,00) datenum(2015,12,30,03,00,00);...
        datenum(2015,12,19,19,00,00) datenum(2015,12,19,21,00,00);...
        datenum(2016,02,17,00,25,00) datenum(2016,02,17,15,00,00);...
        datenum(2016,03,19,00,00,00) datenum(2016,04,07,09,00,00)];
    for i=1:size(rem,1)
        aux=find(param(:,1)>=rem(i,1) & param(:,1)<=rem(i,2));
        param(aux,2:4)=NaN;
        sdist(aux,:)=NaN;
    end
    % Ajuste de relógio, com base em smps, cpc, neph e aeth
    aux=find(param(:,1)<=datenum(2016,02,21,07,45,00));
    param(aux,1)=param(aux,1)-datenum(0,0,0,1,0,0); % subtraindo 1h do smps    
    %
    % Colocando cpc e smos na mesma base de tempo após ajuste de relógio
    paramplot(:,1)=cpc(:,1); paramplot(:,2:4)=nan;
    for i=1:length(cpc)
        aux=find(param(:,1)==cpc(i,1));
        if(~isempty(aux))
            paramplot(i,2:4)=param(aux,2:4);
        end
    end
    figure,plot(cpc(:,1),paramplot(:,2)./cpc(:,2),'-b.'),grid on,dynamicDateTicks
    ylabel('SMPS/CPC')
    saveas(gcf,[caminho2,tit,'_ratio_SMPSCPC_tseries.jpg'],'jpg')
    %
    figure
    scatter(cpc(:,2),paramplot(:,2),4,cpc(:,1)),grid on,hold on
    [a(1),b(1),rsquare(1)]=flinreg(cpc(:,2),paramplot(:,2));
    x=[min(cpc(:,2)) max(cpc(:,2))];
    plot(x,a(1)*x+b(1),'-k','LineWidth',3)
    aux=find(cpc(:,1)<datenum(2015,12,8));
    [a(2),b(2),rsquare(2)]=flinreg(cpc(aux,2),paramplot(aux,2));
    x=[min(cpc(:,2)) max(cpc(:,2))];
    plot(x,a(2)*x+b(2),'-k')
    aux=find(cpc(:,1)>datenum(2015,12,8) & cpc(:,1)<datenum(2015,12,14));
    [a(3),b(3),rsquare(3)]=flinreg(cpc(aux,2),paramplot(aux,2));
    x=[min(cpc(:,2)) max(cpc(:,2))];
    plot(x,a(3)*x+b(3),'-b')
    aux=find(cpc(:,1)>datenum(2016,1,12) & cpc(:,1)<datenum(2016,1,29));
    [a(4),b(4),rsquare(4)]=flinreg(cpc(aux,2),paramplot(aux,2));
    x=[min(cpc(:,2)) max(cpc(:,2))];
    plot(x,a(4)*x+b(4),'-c')
    aux=find(cpc(:,1)>datenum(2016,1,29) & cpc(:,1)<datenum(2016,2,1));
    [a(5),b(5),rsquare(5)]=flinreg(cpc(aux,2),paramplot(aux,2));
    x=[min(cpc(:,2)) max(cpc(:,2))];
    plot(x,a(5)*x+b(5),'-g')
    aux=find(cpc(:,1)>datenum(2016,2,2));
    [a(6),b(6),rsquare(6)]=flinreg(cpc(aux,2),paramplot(aux,2));
    x=[min(cpc(:,2)) max(cpc(:,2))];
    plot(x,a(6)*x+b(6),'-r')
    legend('data',...
        ['Alldata ',num2str(a(1),'%1.2f'),'*x+',num2str(b(1),'%3.0f'),' R^2=',num2str(rsquare(1),'%1.2f')],...
        ['Dec04-08 ',num2str(a(2),'%1.2f'),'*x+',num2str(b(2),'%3.0f'),' R^2=',num2str(rsquare(2),'%1.2f')],...
        ['Dec08-14 ',num2str(a(3),'%1.2f'),'*x+',num2str(b(3),'%3.0f'),' R^2=',num2str(rsquare(3),'%1.2f')],...
        ['Jan12-29 ',num2str(a(4),'%1.2f'),'*x+',num2str(b(4),'%3.0f'),' R^2=',num2str(rsquare(4),'%1.2f')],...
        ['Jan29-Feb01 ',num2str(a(5),'%1.2f'),'*x+',num2str(b(5),'%3.0f'),' R^2=',num2str(rsquare(5),'%1.2f')],...
        ['Feb02-Apr25 ',num2str(a(6),'%1.2f'),'*x+',num2str(b(6),'%3.0f'),' R^2=',num2str(rsquare(6),'%1.2f')])
    ylabel('ConcSMPS (cm^-^3)','FontSize',12),xlabel('CPC (cm^-^3)','FontSize',12)
    saveas(gcf,[caminho2,tit,'_SMPSxCPC.jpg'],'jpg')    
    %
    % Save Level 3 data 
    Instrument='SMPS_TSI';
    Model='3081';
    serial_num='71331078, with CPC3772,WCPC3787,CPC3010. Data between 19mar-07apr2016 was excluded due to low concentrations, but might be corrected through a closure between smps, neph and aeth, using a Mie code';
    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
% -------------------------------------------------------------------------
