%Data processing of AirPhoton Nephelometer. Adaptation by Marco A. M.
%Franco
clc; clear all; close all;

%cd 'D:\Documentos\Google Drive\Doutorado\Backup Rio Branco\Nefelometro'
cd 'D:\Documentos\Google Drive\Doutorado\Backup Rio Branco\Nefelometro'
%cd 'D:\Documentos\Google Drive\Doutorado\Backup Rio Branco\Nefelometro\teste'

file = dir('*.csv');

for i = 1:length(file)
    myfilename = file(i).name;
    data{i} = fnef_airphoton(myfilename);
end

data_nef = [];
for i = 1:length(data)
    
    date = data{i}(:,2);
    date_conv = datenum(date(:,1),'mm/dd/yy HH:MM:SS');
    
    %atmospheric quantities - temperature, pressure and RH%
    tprh = cell2mat(data{i}(:,15:17));
    
    %backscattering (R G B) 632nm, 532nm and 450nm
    bsck_r = cell2mat(data{i}(:,21));
    bsck_g = cell2mat(data{i}(:,22));
    bsck_b = cell2mat(data{i}(:,23));
    
    %totalscattering (R G B) 632nm, 532nm and 450nm
    totalscat_r = cell2mat(data{i}(:,24));
    totalscat_g = cell2mat(data{i}(:,25));
    totalscat_b = cell2mat(data{i}(:,26));
    
    %conditions for Level 1: scat for 450 nm between 2 Mm-1 and 900 Mm-1;
    %scat for 532 between 2 Mm-1 and 900 Mm-1; scat for 632 between 2 Mm-1 and 900 Mm-1
    %Zerocheck excluded;  Relative humidity must be smaller than 60%;
    
    totalscat(find(totalscat_r(:,1)<=2),1) = NaN;
    totalscat(find(totalscat_g(:,1)<=2),1) = NaN;
    totalscat(find(totalscat_b(:,1)<=2),1) = NaN;
    totalscat(find(totalscat_r(:,1)>=900)) = NaN;
    totalscat(find(totalscat_g(:,1)>=900)) = NaN;
    totalscat(find(totalscat_b(:,1)>=900)) = NaN;
    
    bsck(find(bsck_r(:,1)<=1),1) = NaN;
    bsck(find(bsck_g(:,1)<=1),1) = NaN;
    bsck(find(bsck_b(:,1)<=1),1) = NaN;
    bsck(find(bsck_r(:,1)>=900)) = NaN;
    bsck(find(bsck_g(:,1)>=900)) = NaN;
    bsck(find(bsck_b(:,1)>=900)) = NaN;
    
    data_nef = vertcat(data_nef,[date_conv,totalscat_b,totalscat_g,totalscat_r, bsck_b,bsck_g,bsck_r, tprh]); %Usefull data
    data_nef(find(tprh(:,3)>60),:) = []; %Excludes pressure higher than 60% of relative humidity
    
end

figure,plot(data_nef(:,1),data_nef(:,2),'.b',data_nef(:,1),data_nef(:,3),'.g',...
    data_nef(:,1),data_nef(:,4),'.r');
xlabel(''); ylabel('Scattering coef. (Mm-1)');
legend('450 nm','532 nm','632 nm')
dynamicDateTicks
%
% figure, plot(data_nef(:,1),data_nef(:,9),'.g'); %apenas 532 nm
% xlabel(' UTC '); ylabel('Scattering coef. (Mm-1)');
% legend('532 nm')
% dynamicDateTicks

%close all
cd 'D:\Google Drive\Doutorado\Campanha Rio Branco - AC\Dados\Processamento_LFA\Nef_AirPhoton'

[data_nef,N0,Nf,Delta] = correcao_brancos(data_nef);

figure, plot(data_nef(:,1),data_nef(:,2),'.b');
xlabel(''); ylabel('Scattering coef. (Mm-1)');
dynamicDateTicks

%%%%
lev1 = true; lev2 = true;
%%%%
%%  Level 1
if(lev1)
    maxrh = 60;
    origsize=size(data_nef,1);
    orig=data_nef;
    data_nef(data_nef(:,10)>maxrh,2:size(data_nef,2))=NaN; % remove lines with rh>60%
    data_nef(data_nef(:,10)<5,2:size(data_nef,2))=NaN; % remove lines with rh<5% (erro ocasional)
    data_nef(data_nef(:,8)>80 | data_nef(:,8)<10,2:size(data_nef,2))=NaN; % remove lines with temp>80oC or temp<10oC (erro ocasional)
    
    dif(:,1)=data_nef(1:length(data_nef)-1,1);
    dif(:,2)=diff(data_nef(:,2)); % Mm-1
    dif(:,3)=diff(data_nef(:,1))*24*60*60; % sec
    %      picos=find(abs(dif(:,2))./dif(:,3)>=1); % procura variações bruscas no scat450, geralmente associadas a troca de sílica
    %      for i=1:length(picos)
    %          intervalo=[dif(picos(i))-0.75/60/24 dif(picos(i))+0.75/60/24]; % +/- 4min em torno do pico (era 5,coloquei 2)
    %          data_nef(data_nef(:,1)>=intervalo(1) & data_nef(:,1)<=intervalo(2),2:size(data_nef,2))=NaN;
    %      end
    rem=[];
    for i=datenum(2018,8,21):datenum(2018,11,11)
        aux=find(data_nef(:,1)>=i+datenum(0,0,0,0,1,0) & data_nef(:,1)<=i+datenum(0,0,0,0,4,0));
        rem=cat(1,rem,aux);
    end
    
    data_nef(rem,2:size(data_nef,2))=NaN;
    %
    data_nef(isnan(data_nef(:,2)),:)=[];
    
    data_nef(data_nef(:,1)<datenum(2018,8,21,17,00,00)) = NaN;
    data_nef(data_nef(:,1)>=datenum(2018,8,21,18,41,00) & data_nef(:,1)<=datenum(2018,8,21,19,02,00),2:size(data_nef,2))=NaN;
    
    %Blanck corrections
    
    %[data_nef] = correcao_brancos(data_nef);
    
    figure,plot(data_nef(:,1),data_nef(:,2),'.b');
    xlabel(''); ylabel('Scattering coef. (Mm-1)');
    legend('450 nm','532 nm','632 nm')
    dynamicDateTicks
    
    Instrument='Nephelometer';
    Model='AirPhoton';
    serial_num='PM1 inlet with diffusion dryer';
    Site='T0rb';
    Val_level='1';
    Proc_vers='1.0';
    labels={'scat450';'scat532';'scat632';'bscat450';'bscat532';'bscat632';'RH'};
    units={'(Mm-1)';'(Mm-1)';'(Mm-1)';'(Mm-1)';'(Mm-1)';'(Mm-1)';'(prct)'};
    
    export_LFA_Leveis(labels,units,Instrument,Model,serial_num,Site,Val_level,Proc_vers,origsize,data_nef(:,1),cat(2,data_nef(:,2:7),data_nef(:,10)))
    
end

%% Level 2
clc; clear all;close all;
lambda=[450 532 632]; lev2 = true;
ano=2018;
% Level 2: Correct for truncation error according to Müller et al,2011b.
% Apply STP corrections. Ansgstrom must be within -1 and 10.

if(lev2)
    %     C=load('C:\A\ATTO\met\STP_InstTower_ATTO_filled_2014-2016_5min.mat');
    %     C = load([caminho,'T0aNeph2018C4_level0.mat']);
    %     stp=C.stp5(C.stp5(:,1)>=datenum(ano,1,1) & C.stp5(:,1)<=datenum(ano+1,1,1),2);
    %     %C=xlsread('C:\A\ATTO\Neph_Aurora\ATTO_silica_change_2014.xls','Sheet2');
    %     %silica(:,1)=datenum(C(:,1:6)); % start silica change (utc)
    %     %silica(:,2)=datenum(C(:,7:12)); % finish silica change
    %     %clear C
    %cd 'D:\Google Drive\Doutorado\Campanha Rio Branco - AC\Dados\dados_meteorologicos'
    cd 'D:\Documentos\Google Drive\Doutorado\Campanha Rio Branco - AC\Dados\dados_meteorologicos\'
    %meteo = import_meteo('RioBranco_AMES_TPR159_Level0_5min.csv');
    %meteo = importBDMEP('RioBrancoBDMEP_BDMEP Meteo_RioBranco_Level1_5min.csv');
    meteo = importMeteo('RioBranco-2018_AMES_TPR159_Level0_5min.csv');
    %temperatura = str2double(meteo.temperatura(2:end));
    %pressao = str2double(meteo.pressao(2:end));
    temperatura = meteo.Temperature_C_;
    pressao = meteo.Pressure_mbar_;
    
    stp = ((temperatura+293.15)./293.15).*(1013.25./pressao);
end

%caminho1 = 'D:\Google Drive\Doutorado\Campanha Rio Branco - AC\Dados\Processamento_LFA\Nef_AirPhoton\0_Ascii-files\Level1\';
caminho1 = 'D:\Documentos\Google Drive\Doutorado\Campanha Rio Branco - AC\Dados\Processamento_LFA\Nef_AirPhoton\0_Ascii-files\Level1\'


if(lev2)
    % Lê os arquivos Level 1 que acabou de salvar
    aux=importdata([caminho1,'T0rb_Nephelometer_AirPhoton_Level1_5min.csv']);
    level1_5min=cat(2,aux.data(:,1:2),aux.data(:,9),aux.data(:,16),aux.data(:,23),aux.data(:,30),aux.data(:,37),aux.data(:,44));
    %level1_5min=cat(2,aux.data(:,1:2),aux.data(:,3),aux.data(:,10),aux.data(:,17),aux.data(:,24),aux.data(:,31),aux.data(:,38));
    
    clear aux
    
    level1_5min(:,1)=level1_5min(:,1)/(24*60*60)+datenum(1970,1,1); % unix time to matlab time
    orig=level1_5min;
    origsize=sum(~isnan(level1_5min(:,2)));
    %
    % Calculate Angstrom B/G B/R G/R
    angstrom=NaN(length(level1_5min),5);
    level1_5min(level1_5min(:,2)<0,2:size(level1_5min,2))=NaN;
    level1_5min(level1_5min(:,3)<0,2:size(level1_5min,2))=NaN;
    level1_5min(level1_5min(:,4)<0,2:size(level1_5min,2))=NaN;
    
    for i=1:length(angstrom)
        if(mod(i,100)==0)
            
        end
        if(~isnan(level1_5min(i,2)))
            %angstrom(i,1)=-log10(level1_5min(i,2)./level1_5min(i,3))/log10(lambda(1)/lambda(2));
            angstrom(i,1)=-log(level1_5min(i,2)./level1_5min(i,3))/log(lambda(1)/lambda(2));
            %angstrom(i,2)=-log10(level1_5min(i,2)./level1_5min(i,4))/log10(lambda(1)/lambda(3));
            angstrom(i,2)=-log(level1_5min(i,2)./level1_5min(i,4))/log(lambda(1)/lambda(3));
            %angstrom(i,3)=-log10(level1_5min(i,3)./level1_5min(i,4))/log10(lambda(2)/lambda(3));
            angstrom(i,3)=-log(level1_5min(i,3)./level1_5min(i,4))/log(lambda(2)/lambda(3));
        end
    end
    angstrom(angstrom(:,2)<=-1 | angstrom(:,2)>=10,:)=NaN;
    
    figure, plot(level1_5min(:,1),angstrom)
    
    %
    % STP and truncation error correction (Muller et al., 2011)
    %aux=find(level1_5min(:,1)<datenum(2014,5,4)); % PM10 inlet
    %trunc(aux,1)=1.455-0.189*angstrom(aux,1); % Muller et al. 2011 truncation correction at 450 nm (no cut)
    %trunc(aux,2)=1.434-0.176*angstrom(aux,2); % same for 525 nm
    %trunc(aux,3)=1.403-0.156*angstrom(aux,3); % same for 635 nm
    %trunc(aux,4)=0.963; % backscatter correction 450 nm (no cut)
    %trunc(aux,5)=0.971; % same for 525 nm
    %trunc(aux,6)=0.968; % same for 635 nm
    %aux=find(level1_5min(:,1)>=datenum(2014,5,4)); % PM1 inlet
    %trunc(aux,1)=1.213-0.060*angstrom(aux,1); % Muller et al. 2011 truncation correction at 450 nm (sub-um)
    %trunc(aux,2)=1.207-0.061*angstrom(aux,2); % same for 525 nm
    %trunc(aux,3)=1.176-0.053*angstrom(aux,3); % same for 635 nm
    %trunc(aux,4)=0.932; % backscatter correction 450 nm (sub-um)
    %trunc(aux,5)=0.935; % same for 525 nm
    %trunc(aux,6)=0.935; % same for 635 nm
    %trunc(trunc(:,1:3)<1)=1; % truncation correction factor should be >1.0
    
    level1_5min(:,2)=level1_5min(:,2).*stp;%.*trunc(:,1); % corrected scat
    level1_5min(:,3)=level1_5min(:,3).*stp;%.*trunc(:,2);
    level1_5min(:,4)=level1_5min(:,4).*stp;%.*trunc(:,3);
    level1_5min(:,5)=level1_5min(:,5).*stp;%.*trunc(:,4);
    level1_5min(:,6)=level1_5min(:,6).*stp;%.*trunc(:,5);
    level1_5min(:,7)=level1_5min(:,7).*stp;%.*trunc(:,6);
    
    lambda=[450 532 632];
    disp('Calculating Angstrom and Scat637')
    for i=1:length(level1_5min)
        if(mod(i,100)==0)
            %disp(i)
        end
        if(level1_5min(i,2)*level1_5min(i,3)*level1_5min(i,4)>0)
            p=polyfit(log10(lambda),log10(level1_5min(i,2:4)),1);
            angstrom(i,4)=-p(1);   % Angstrom
            angstrom(i,5)=10^(polyval(p,log10(637))); % scat637
            clear p
        end
    end
    %
    
    %cd 'D:\Google Drive\Doutorado\Campanha Rio Branco - AC\Dados\Processamento_LFA\Nef_AirPhoton'
    cd 'D:\Documentos\Google Drive\Doutorado\Campanha Rio Branco - AC\Dados\Processamento_LFA\Nef_AirPhoton';
    
    
    figure, plot(orig(:,1),orig(:,2),'k.'),grid on, hold on
    plot(level1_5min(:,1),level1_5min(:,2),'r.')
    dynamicDateTicks
    figure,plot(level1_5min(:,1),angstrom(:,4),'g.'),grid on,dynamicDateTicks
    %
    % Save Level 2 data
    Instrument='Nephelometer';
    Model='AirPhoton';
    serial_num='PM1 inlet with diffusion dryer';
    Site='T0rb';
    Val_level='2';
    Proc_vers='1.0';
    labels={'scat450';'scat532';'scat632';'bscat450';'bscat532';'bscat632';'RH';'Angstrom';'scat637'};
    units={'(Mm-1)';'(Mm-1)';'(Mm-1)';'(Mm-1)';'(Mm-1)';'(Mm-1)';'(prct)';'(adim)';'(Mm-1)'};
    
    export_LFA_Leveis(labels,units,Instrument,Model,serial_num,Site,Val_level,Proc_vers,origsize,level1_5min(:,1),[level1_5min(:,2:8) angstrom(:,4:5)])
    
end

% Exportation - Raw data
% cd 'D:\Documentos\Google Drive\Doutorado\Campanha Rio Branco - AC\Dados\Processamento_LFA\Nef_AirPhoton'
%
%  data1 = num2cell(data_nef(:,2:end));
%  data2 = cellstr(datestr(data_nef(:,1)));
%  dado_total = [data2 data1];
%  isvarname cabeca
%  cabeca = {'Time_UTC','backscattering_450','backscattering_532',...
%      'backscattering_632','Scattering_450','Scattering_532',...
%      'Scattering_632'};
%  dado = cell2table(dado_total);
%  dado.Properties.VariableNames = cabeca;
%  %writetable(dado,'neph_raw.csv','Delimiter',';')
%  writetable(dado,'T0rb_Level0.csv','Delimiter',';')

%% Verification
close all
%cd 'D:\Google Drive\Doutorado\Campanha Rio Branco - AC\Dados\Processamento_LFA\Nef_AirPhoton\0_Ascii-files\Level2'

cd 'D:\Documentos\Google Drive\Doutorado\Campanha Rio Branco - AC\Dados\Processamento_LFA\Nef_AirPhoton\0_Ascii-files\Level2\'

file = dir('*.csv');
data = [];
for i = 1:length(file)
    myfilename = file(i).name;
    data{i} = importNef(myfilename);
end

data_nef = [];
for i = 1:length(data)
    
    date = data{i}(:,1);
    date_conv = datenum(date,'dd-mmm-yyyy HH:MM');
    
    %atmospheric quantities - temperature, pressure and RH%
    %tprh = cell2mat(data{i}(:,15:17));
    
    %backscattering (R G B) 450nm, 532nm and 632nm
    bsck = str2double(data{i}(:,6:8));
    
    %totalscattering (R G B) 450nm, 532nm and 632nm
    totalscat = str2double(data{i}(:,3:5));
    
    %Angstrom
    
    angs = str2double(data{i}(:,10));
    
    %scat 637nm
    scat = str2double(data{i}(:,11));
    
    data_nef = vertcat(data_nef,[date_conv, totalscat,bsck,angs,scat]); %Usefull data
    
end

figure,plot(data_nef(:,1),data_nef(:,2),'.b',data_nef(:,1),data_nef(:,3),'.g',...
    data_nef(:,1),data_nef(:,4),'.r');
xlabel(''); ylabel('Scattering coef. (Mm-1)');
legend('450 nm','532 nm','632 nm')
dynamicDateTicks


figure,plot(data_nef(:,1),data_nef(:,5),'.b',data_nef(:,1),data_nef(:,6),'.g',...
    data_nef(:,1),data_nef(:,7),'.r');
xlabel(''); ylabel('Backscattering coef. (Mm-1)');
legend('bscat450 nm','bscat532 nm','bscat632 nm')
dynamicDateTicks

figure,plot(data_nef(:,1),data_nef(:,8),'.k');
xlabel(''); ylabel('Angstrom exponent (mean)');
title('Scattering Angstron exponent (mean), PM 2.5, at Rio Branco - AC, Brazil')
dynamicDateTicks

%% Single Scattering albedo

cd 'D:\Google Drive\Doutorado\Campanha Rio Branco - AC\Dados\Processamento_LFA\MAAP\comparacao_level1'
%cd 'D:\Google Drive\Doutorado\Campanha Rio Branco - AC\Dados\Processamento_LFA\MAAP\comparacao_level1'

data = [];
file = dir('*.csv');
for i = 1:length(file)
    myfilename = file(i).name;
    data{i} = importBC(myfilename);
end

%%% adaptacoes feitas para comparacao entre MAAP 156 e MAAP 244
dado_maap = [];
for j = 1:length(data)
    
    %date = strcat(data{j}(:,1),{' '},data{j}(:,2));
    date = data{j}(:,1);
    date_conv = datenum(date,'dd-mmm-yyyy HH:MM');
    
    %BC = cell2mat(data{j}(:,4))*1000; % x1000 p converter para ng/cm3
    BC = str2mat(data{j}(:,3));
    BC = str2num(BC);
    BC(find(BC == 0)) = NaN;
    
    dado_maap{j} = vertcat([date_conv,BC]);
    
end

SSA = data_nef(:,4)./(data_nef(:,4) + dado_maap{1}(:,2)*6.6/1000);
SSA(find(SSA > 0.95)) = NaN;
SSA(find(SSA < 0.5)) = NaN;

figure, plot(data_nef(:,1),SSA,'.k')
ylabel('SSA');
title('Single Scattering Albedo (SSA) at Rio Branco - AC, Brazil, PM 2.5, 2018')
dynamicDateTicks


%%

lev3 = true;
lambda=[450 525 635];
caminho2 = 'D:\Documentos\Google Drive\Doutorado\Campanha Rio Branco - AC\Dados\Processamento_LFA\Nef_AirPhoton\0_Ascii-files\Level2\'
caminho3 = 'D:\Documentos\Google Drive\Doutorado\Campanha Rio Branco - AC\Queimadas em Rio Branco 2018\Dados_processados\AE33\Level2\'
%if(lev3)
% Lê os arquivos Level 2 que acabou de salvar
aux=importdata([caminho2,'T0rb_Nephelometer_AirPhoton_Level2_5min.csv']);
level2_5min=aux.data;
level2_5min(:,1)=level2_5min(:,1)/(24*60*60)+datenum(1970,1,1); % unix time to matlab time
orig=level2_5min;
origsize=sum(~isnan(level2_5min(:,2)));
%
C=importdata([caminho3,'RioBranco_Aethalometer_AE33_Level2_5min.csv']);
maap=C.data;
maap(:,1)=maap(:,1)/(24*60*60)+datenum(1970,1,1); % unix time to matlab time (utc)
maap(:,2)=[]; % Mm-1
maap(maap(:,17)<=0,2)=NaN; %637nm
clear C aux
%
% Specific correction for 2014
aux=find(level2_5min(:,1)>=datenum(2014,8,22) & level2_5min(:,1)<=datenum(2015,2,10,18,0,0));
level2_5min(aux,2)=(level2_5min(aux,2)-1.0)/0.3; % correction based on Mie calculations
level2_5min(aux,3)=(level2_5min(aux,3)-1.0)/0.3;
level2_5min(aux,4)=(level2_5min(aux,4)-1.0)/0.3;
%
% Calculate Angstrom by power law fit
angstrom=nan(length(level2_5min),1);
scat637=nan(length(level2_5min),1);
for i=1:length(level2_5min)
    if (mod(i,100)==0)
        disp(i)
    end
    if(level2_5min(i,4)>0)
        p=polyfit(log10(lambda),log10(level2_5min(i,2:4)),1);
        angstrom(i)=-p(1);
        scat637(i)=10^polyval(p,log10(637));
    end
end
%
ssa=scat637./(scat637+maap(:,17));
%
rem1=find(ssa<0.1);
rem2=find(angstrom<=-1 | angstrom>=10);
level2_5min(rem1,2:size(level2_5min,2))=NaN;
level2_5min(rem2,2:size(level2_5min,2))=NaN;
ssa(rem1)=NaN;
ssa(rem2)=NaN;
angstrom(rem1)=NaN;
angstrom(rem2)=NaN;
clear rem1 rem2
%
%cd 'D:\Documentos\Google Drive\Doutorado\Tratamento\2018\325m\nefelometro\'
%cd 'D:\Documentos\Google Drive\Doutorado\Tratamento\2018\325m\NEPH\2018\formato_txt\'
%cd 'D:\Documentos\Google Drive\Doutorado\Tratamento\2017\60m\Neph\LFA\Processamento\Processamento_LFA\'
figure
subplot(4,1,1),plot(level2_5min(:,1),level2_5min(:,2),'b.'),grid on, %   dynamicDateTicks
subplot(4,1,2),plot(maap(:,1),maap(:,2),'k.'),grid on, %dynamicDateTicks
subplot(4,1,3),plot(level2_5min(:,1),angstrom,'g.'),grid on, %dynamicDateTicks
subplot(4,1,4),plot(level2_5min(:,1),ssa,'r.'),grid on, %dynamicDateTicks
%
% Additional data periods to remove:
% Há períodos com ang>5 e ssa<0.5 que coincidem com períodos em que scat e abs
% estão próximos de zero. Esses períodos não coincidem com as trocas de sílica
% reportadas. Às vezes no início desses períodos há descontinuidades nas séries
% temporais de neph e maap (queda súbita). Não sei o que pode ter
% causado esses efeitos, e por isso decidi manter esses pontos no Level3.
%
%%
% Save Level 3 data
cd 'D:\Documentos\Google Drive\Doutorado\Campanha Rio Branco - AC\Dados\Processamento_LFA\Nef_AirPhoton\'

Instrument='Nephelometer';
Model='AirPhoton';
serial_num='x';
Site='T0a';
Val_level='3';
Proc_vers='1.0';
labels={'scat450';'scat532';'scat632';'bscat450';'bscat532';'bscat632';'RH';'Angstrom';'scat637';'SSA'};
units={'(Mm-1)';'(Mm-1)';'(Mm-1)';'(Mm-1)';'(Mm-1)';'(Mm-1)';'(prct)';'(adim)';'(Mm-1)';'(adim)'};
   
export_LFA_Leveis(labels,units,Instrument,Model,serial_num,Site,Val_level,Proc_vers,origsize,level2_5min(:,1),[level2_5min(:,2:10) ssa])
    