% Append Aethalometer AE33 raw data files. Correct for multiple scattering effects and STP.
% Save files in LFA standard format.
% Written by Luciana Rizzo
% 
% 1-Time(utc) 2-Timebase 3-RefCh1 4-Sen1Ch1 5-Sen2Ch1 6-RefCh2 7-Sen1Ch2 8-Sen2Ch2
% 9-RefCh3 10-Sen1Ch3 11-Sen2Ch3 12-RefCh4 13-Sen1Ch4 14-Sen2Ch4 15-RefCh5 16-Sen1Ch5
% 17-Sen2Ch5 18-RefCh6 19-Sen1Ch6 20-Sen2Ch6 21-RefCh7 22-Sen1Ch7 23-Sen2Ch7 24-Flow1
% 25-Flow2 26-FlowC 27-Pressure(Pa) 28-Temperature(°C) 29-BB(%) 30-ContTemp 31-SupplyTemp 32-Status	
% 33-ContStatus 34-DetectStatus 35-LedStatus 36-ValveStatus 37-LedTemp 
% 38-BC11 39-BC12 40-BC1 41-BC21 42-BC22 43-BC2 44-BC31 45-BC32 46-BC3 47-BC41 45-BC42
% 46-BC4 47-BC51 48-BC52 49-BC5 50-BC61 51-BC62 52-BC6 53-BC71 54-BC72 55-BC7
% 56-K1 57-K2 58-K3 59-K4 60-K5 61-K6 62-K7 63-TapeAdvCount 64-69 ?
%
caminho='C:\A\GOAMAZON\T1\AE33_S0200193sn\rawdata\';
caminho1='C:\A\GOAMAZON\T1\AE33_S0200193sn\';
site='T1_2015-2016_BC';
wl=([370 470 520 590 660 880 950]);  % 7 wavelengths AE33
%
data=zeros(1,78);
files=dir(caminho);
disp('-----Reading raw data files')
for ifile=3:size(files,1)
    nome=getfield(files, {ifile,1}, 'name');
    disp(nome)
    delimiterIn=' ';
    headerlinesIn=8;
    C=importdata(strcat(caminho,nome),delimiterIn,headerlinesIn);
    data_(:,1)=datenum([char(C.textdata(9:length(C.textdata),1)),char(C.textdata(9:length(C.textdata),2))],'yyyy/mm/ddHH:MM:SS');
    data_(:,2:78)=C.data;
    data=cat(1,data,data_);
    clear  data_ C
end
data(1,:)=[];
[~,n]=sort(data(:,1)); 
data=data(n,:);
rawdata=data;
save([caminho1,site,'alldata_level0.mat'],'data','-mat')
%
% Plots for quality control ------------------------------
figure,subplot(3,1,1),plot(data(:,1),data(:,24),'-b.'),grid on,dynamicDateTicks
ylabel('Flow1'), axis tight
subplot(3,1,2),plot(data(:,1),data(:,25),'-c.'),grid on,dynamicDateTicks
ylabel('Flow2'), axis tight
subplot(3,1,3),plot(data(:,1),data(:,26),'-m.'),grid on,dynamicDateTicks
ylabel('Flow3'), axis tight
figure,subplot(2,1,1),plot(data(:,1),data(:,27),'-r.'),grid on,dynamicDateTicks
ylabel('Press(Pa)'), axis tight
subplot(2,1,2),plot(data(:,1),data(:,28),'-g.'),grid on,dynamicDateTicks
ylabel('Temp(oC)'), axis tight
figure,subplot(3,1,1),plot(data(:,1),data(:,30),'-m.'),grid on,dynamicDateTicks
ylabel('ContTemp'), axis tight
subplot(3,1,2),plot(data(:,1),data(:,31),'-r.'),grid on,dynamicDateTicks
ylabel('SupplyTemp'), axis tight
subplot(3,1,3),plot(data(:,1),data(:,37),'-m.'),grid on,dynamicDateTicks
ylabel('LedTemp'), axis tight
figure,subplot(3,1,1),plot(data(:,1),data(:,29),'-r.'),grid on,dynamicDateTicks
ylabel('BB(%)'), axis tight
subplot(3,1,2),plot(data(:,1),data(:,40),'-b.'),grid on,dynamicDateTicks
ylabel('BC370'), axis tight
subplot(3,1,3),plot(data(:,1),data(:,55),'-k.'),grid on,dynamicDateTicks
ylabel('BC880'), axis tight
figure,plot(data(:,1),data(:,36),'-k.'),grid on,dynamicDateTicks
ylabel('ValveStatus'), axis tight
%
% Level 1 criteria ----------------------------------------
% Exclude ValveStatus=0 and the subsequent point
% Exclude Flow3=0
% BC880 data must be >-1 ug.m-3 and <100 ug.m-3
% 
orig=data;
aux=find(data(:,36)>0);
aux=cat(1,aux,aux+1,find(data(:,26)==0),find(data(:,55)<-1000),find(data(:,55)>100000));
data(aux,:)=[];
figure,subplot(3,1,1),plot(orig(:,1),orig(:,29),'-r.'),grid on,dynamicDateTicks, hold on
plot(data(:,1),data(:,29),'k.')
ylabel('BB(%)'), axis tight
subplot(3,1,2),plot(orig(:,1),orig(:,40),'-r.'),grid on,dynamicDateTicks, hold on
plot(data(:,1),data(:,40),'k.')
ylabel('BC370'), axis tight
subplot(3,1,3),plot(orig(:,1),orig(:,55),'-r.'),grid on,dynamicDateTicks, hold on
plot(data(:,1),data(:,55),'k.')
ylabel('BC880'), axis tight
%
% Save Level 1 data 
data1=[data(:,1),data(:,40),data(:,43),data(:,46),data(:,49),data(:,52),data(:,55),data(:,58),data(:,29)]; 
data1(:,2:8)=data1(:,2:8)/1000; % ng/m3 to ug/m3
origsize=length(orig);
Instrument='Aethalometer';
Model='AE33';
serial_num='S0200193';
Site=site;
Val_level='1';
Proc_vers='1.1';
labels={'BC370','BC470','BC520','BC590','BC660','BC880','BC950','BB'};
units={'(ug.m-3)','(ug.m-3)','(ug.m-3)','(ug.m-3)','(ug.m-3)','(ug.m-3)','(ug.m-3)','(percent)'};
export_LFA_GoAmazon_jan2015(labels,units,Instrument,Model,serial_num,Site,Val_level,Proc_vers,origsize,data1(:,1),data1(:,2:9))
%
%
% Level 2 -----------------------------------------
% Lê os arquivos Level 1 
caminho1='C:\A\GOAMAZON\T1\AE33_S0200193sn\';
site='T1_2015-2016_BC';
wl=([370 470 520 590 660 880 950]);  % 7 wavelengths AE33
csa=[39.5 31.1 28.1 24.8 22.2 16.6 15.4]; % aethalometer cross section (m2/g)
aux=importdata('C:\A\GOAMAZON\T1\LFA_standard\T1_2015-2016_BC_Aethalometer_AE33_Level1_5min.csv');
data1=cat(2,aux.data(:,1:2),aux.data(:,9),aux.data(:,16),aux.data(:,23),aux.data(:,30),aux.data(:,37),aux.data(:,44),aux.data(:,51));
data1(:,1)=data1(:,1)/(24*60*60)+datenum(1970,1,1); % unix time to matlab time
for jj=1:7
    data1(data1(:,jj)<=0,2:size(data1,2))=NaN;
end
origsize=sum(~isnan(data1(:,2)));
orig=data1;
% Calculate attenuation coefficient (Mm-1) from BC concentrations
for jj=1:7
    data1(:,jj+8)=data1(:,jj+1).*csa(jj); % 
end
%
aux=importdata('C:\A\GOAMAZON\T1\LFA_standard\T1_2015-2016_Nephelometer_Aurora_Ecotech_2000_Level2_5min.csv');
data2=aux.data; % 
data2(:,1)=data2(:,1)/(24*60*60)+datenum(1970,1,1); % unix time to matlab time
data2(:,7)=data2(:,6); % SAE (Scattering Angstrom exponent)
data2(:,9)=nan; % SSA (single scattering albedo)
%
% Cref correction factor (multiple scattering) --> eq 10 from Schmid 2006
% Since there was not a MAAP running at T1, Cref from T2 measurements will be used.
aux=importdata('C:\A\GOAMAZON\T2\Level3\T2_Aethalometer_AE33_193_Level1_30min.csv');
aetht2=cat(2,aux.data(:,1:2),aux.data(:,9),aux.data(:,16),aux.data(:,23),aux.data(:,30),aux.data(:,37),aux.data(:,44));
aetht2(:,1)=aetht2(:,1)/(24*60*60)+datenum(1970,1,1);
for jj=1:7
    aetht2(:,jj+8)=aetht2(:,jj+1).*csa(jj)/1000; % 
end
aetht2(aetht2(:,13)<=0,13)=NaN;
aetht2(aetht2(:,12)<=0,12)=NaN;
aae=-log10(aetht2(:,13)./aetht2(:,12))/log10(660/590);
aetht2(:,16)=10.^(log10(aetht2(:,13))+aae.*log10(660/637)); % abs637
aux=importdata('C:\A\GOAMAZON\T2\Level3\T2_MAAP_5012_Level2_30min.csv');
maapt2=aux.data; % BC637, Abs637
maapt2(:,1)=maapt2(:,1)/(24*60*60)+datenum(1970,1,1);
aux=find(maapt2(:,1)>=datenum(2014,3,1) & maapt2(:,1)<datenum(2014,4,30));
Cref=aetht2(aux,16)./maapt2(aux,3);
Cref=Cref(~isnan(Cref));
Cref(Cref<0)=[];
outlier=find(Cref>(mean(Cref)+2*std(Cref)));
Cref(outlier)=[];
figure,plot(Cref,'k.'),grid on, ylabel('Cref (637nm)')
disp('Mean Cref will be calculated from the data in the plot.')
Crefmean=mean(Cref)
clear aetht2 maapt2 aae aux
%
% Filter loading correction is not needed, because AE33 is dual spot
%
% SAE (Scattering Angstrom exponent) and SSA (single scattering albedo) from maap and neph data
% Since there was not a MAAP running at T1, mean SSA from T2 measurements will be used.
aux=importdata('C:\A\GOAMAZON\T2\Level3\T2_Nephelometer_Aurora_Ecotech_3000_Level3_30min.csv');
ssat2=[aux.data(:,1) aux.data(:,9)]; % SSA637
ssat2(:,1)=ssat2(:,1)/(24*60*60)+datenum(1970,1,1);
aux=find(ssat2(:,1)>=datenum(2014,3,1) & ssat2(:,1)<datenum(2014,4,30));
SSAmean=nanmedian(ssat2(aux,2))
SAEmean=nanmedian(data2(:,7))
clear ssat2 aux
%
% Parameterization of C(lambda) as a function of AAE (1.0, 1.25 1.5, 1.75 2.0) --> Fig4a from Schmid 2006
Arnott=[370 470 521 590 660 880 950; 1.813 2.073 2.076 2.104 2.182 2.226 2.199; 0.0335 0.0457 0.0523 0.0616 0.0713 0.1038 0.1148];    
cor=['r','g','b','c','m','y','k'];
Cparam=zeros(3,length(wl));
AAEr=[1.0 1.25 1.5 1.75 2.0];  % assumed AAE
A=zeros(1,length(AAEr));
B=zeros(1,length(AAEr));
D=zeros(1,length(AAEr));
opt2=fitoptions('Method','NonlinearLeastSquares');
f2=fittype('A*(x^2) + B*x + D','options',opt2);
%
% Iterations to calculate AAE and corrected abs(lambda)  
maxit=10; % maximum number of iterations to calculate AAE (Angstrom Absorption Exponent)
prec=0.005; % precision for AAE convergence
AAE=nan(size(data1,1),maxit);  
AAEconv=nan(size(data1,1),1);  % final estimative AAE, after iterations converged
C=zeros(size(data1,1),length(wl));
corrected=nan(size(data1,1),length(wl));
meanSSA_SAE=false;
disp('Looping through data files to do iterative calculations until convergence of AAE')
for ii=1:size(data1,1)
    if(mod(ii,50)==0)
        disp(ii)
    end
    if(sum(isnan(data1(ii,:))>0))
        continue
    end
    % Raw AAE (first guess)
    b=glmfit(log(wl),log(data1(ii,9:15))); % linear fit (y=b1+ b2*x)
    AAE(ii,1)=-b(2);
    clear b
    % Obtain parameterization of C(lambda,AAE) according to SSA and SAE, point by point.
    % If SSA or SAE are not available, use SSAmean and SAEmean.
    if(~meanSSA_SAE)
        if(isnan(data2(ii,9)))
            SSA=SSAmean;
        else
            SSA=data2(ii,9);
        end
        if(isnan(data2(ii,7)))
            SAE=SAEmean;
        else
            SAE=data2(ii,7);
        end
        Cparam=zeros(3,length(wl));
        AAEr=[1.0 1.25 1.5 1.75 2.0];  % assumed AAE
        A=zeros(1,length(AAEr));
        B=zeros(1,length(AAEr));
        D=zeros(1,length(AAEr));
        opt2=fitoptions('Method','NonlinearLeastSquares');
        f2=fittype('A*(x^2) + B*x + D','options',opt2);
        for k=1:length(AAEr)
            ratio=wl./637;
            w0=SSA*(ratio.^(-SAE))./(SSA*(ratio.^(-SAE))+(1-SSA)*(ratio.^(-AAEr(k))));
            Cparam(k,:)=Arnott(2,:)+(Arnott(3,:).*(w0./(1-w0)));
            [cc,gof] = fit(log(wl)',log(Cparam(k,:)'),f2,'Startpoint',[-0.2 2.0 -9]);
            A(k)=getfield(cc,'A');   % parameters A,B,D
            B(k)=getfield(cc,'B');
            D(k)=getfield(cc,'D');
        end
        [cc,gof]=fit(AAEr',A',f2,'Startpoint',[0.1 0.1 0.1]);
        AB(1)=getfield(cc,'A');   % parameters A,B,D
        AB(2)=getfield(cc,'B');
        AB(3)=getfield(cc,'D');
        [cc,gof]=fit(AAEr',B',f2,'Startpoint',[-1 2 2]);
        AB(4)=getfield(cc,'A');   % parameters A,B,D
        AB(5)=getfield(cc,'B');
        AB(6)=getfield(cc,'D');
        clear SSA SAE Cparam A B D f2 opt2 curva cc gof AAEr ratio w0
    end
    % Iterations to obtain AAE
    dif=1;
    it=1;  % count iterations
    while(dif>prec && it<=maxit)
        %disp(it)
        it=it+1;
        x=AAE(ii,it-1);
        Acoef=(AB(1)*(x^2)) + (AB(2)*x) + AB(3); %Acoef= (0.102*(x^2)) - (0.187*x) - 0.141;
        Bcoef=(AB(4)*(x^2)) + (AB(5)*x) + AB(6); %Bcoef=(-1.275*(x^2)) + (2.564*x) + 1.827;
        for jj=1:length(wl)
            C(ii,jj)=Crefmean*(wl(jj)^(Acoef*log(wl(jj))+Bcoef))/(637^(Acoef*log(637)+Bcoef));
            corrected(ii,jj)=data1(ii,jj+8)/(C(ii,jj));  % tentative corrected absorption coefficient
        end
        b=glmfit(log(wl),log(corrected(ii,:))); % linear fit (y=b1+ b2*x)
        AAE(ii,it)=-b(2);   % tentative AAE
        clear b
        dif=abs(AAE(ii,it)-AAE(ii,it-1));
    end
    if(dif<=prec); AAEconv(ii)=AAE(ii,it); end
end
% date+time, bc370-950, %BB, abs370-950, abs637, AAE
level2=cat(2,data1(:,1:8),orig(:,9),corrected);
level2(:,17)=10.^(log10(level2(:,13))-AAEconv.*log10(637/wl(4)));
level2(:,18)=AAEconv;
%
figure(1)
subplot(2,1,1),plot(data1(:,9),level2(:,10),'b.'),grid on, hold on
[a,b,rsquare]=flinreg(data1(:,9),level2(:,10));
x=[min(data1(:,9)) max(data1(:,9))];
plot(x,a*x+b,'-k')
legend('data',[num2str(a,'%1.2f'),'*x+',num2str(b,'%3.0f'),' R^2=',num2str(rsquare,'%1.2f')],'Location','SouthEast')  
title('370nm'), xlabel('raw abs (Mm^-^1)'),ylabel('corrected abs (Mm^-^1)')
subplot(2,1,2),plot(data1(:,14),level2(:,15),'r.'),grid on, hold on
[a,b,rsquare]=flinreg(data1(:,14),level2(:,15));
x=[min(data1(:,14)) max(data1(:,14))];
plot(x,a*x+b,'-k')
legend('data',[num2str(a,'%1.2f'),'*x+',num2str(b,'%3.0f'),' R^2=',num2str(rsquare,'%1.2f')],'Location','SouthEast')   
title('880nm'), xlabel('raw abs (Mm^-^1)'),ylabel('corrected abs (Mm^-^1)')
saveas(1,strcat(caminho1,'rawXcorrected'),'jpg')
%
figure(2)
subplot(3,1,1),plot(level2(:,1),level2(:,15),'b.'),grid on,ylabel('Corrected abs 880nm'),dynamicDateTicks
set(gca,'XLim',[datenum(2015,12,08) datenum(2016,04,26)])
subplot(3,1,2),plot(level2(:,1),AAEconv,'m.'),grid on,ylabel('Corrected AAE'),dynamicDateTicks
set(gca,'XLim',[datenum(2015,12,08) datenum(2016,04,26)])
subplot(3,1,3),plot(level2(:,1),level2(:,9),'k.'),grid on,ylabel('%BB'),dynamicDateTicks
set(gca,'XLim',[datenum(2015,12,08) datenum(2016,04,26)])
saveas(2,strcat(caminho1,'timeseries'),'jpg')
%
figure(3)
hist(AAEconv,0.5:0.1:3.5), grid on
xlabel('Corrected AAE'),ylabel('Frequency')
legend(['Median AAE = ',num2str(nanmedian(AAEconv))])
saveas(3,strcat(caminho1,'hist_AAE'),'jpg')
%
% STP correction
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(level2))
    disp('check length of meteorological data file')
    stop
end
for j=2:size(level2,2)
    if(j~=9 && j~=18) % skip %BB, AAE
        level2(:,j)=level2(:,j).*stp(:,2);
    end
end
%
% Save Level 2 data 
origsize=sum(~isnan(orig(:,2)));
Instrument='Aethalometer';
Model='AE33';
serial_num='S0200193';
Site=site;
Val_level='2';
Proc_vers='1.1';
labels={'BC370','BC470','BC520','BC590','BC660','BC880','BC950','BB','abs370','abs470','abs520','abs590','abs660','abs880','abs950','abs637','AAE'};
units={'(ug.m-3)','(ug.m-3)','(ug.m-3)','(ug.m-3)','(ug.m-3)','(ug.m-3)','(ug.m-3)','(percent)','(Mm-1)','(Mm-1)','(Mm-1)','(Mm-1)','(Mm-1)','(Mm-1)','(Mm-1)','(Mm-1)','(adim)'};
export_LFA_GoAmazon_jan2015(labels,units,Instrument,Model,serial_num,Site,Val_level,Proc_vers,origsize,level2(:,1),level2(:,2:size(level2,2)))
%
%
%
% Level 3 -----------------------------------------
% Lê os arquivos Level 2
site='T1_2015-2016_BC';
aux=importdata('C:\A\GOAMAZON\T1\LFA_standard\T1_2015-2016_BC_Aethalometer_AE33_Level2_5min.csv');
level2=aux.data;
level2(:,1)=level2(:,1)/(24*60*60)+datenum(1970,1,1); % unix time to matlab time
origsize=sum(~isnan(level2(:,2)));
%
aux=importdata('C:\A\GOAMAZON\T1\LFA_standard\T1_2015-2016_Nephelometer_Aurora_Ecotech_2000_Level3_5min.csv');
neph=aux.data;
neph(:,1)=neph(:,1)/(24*60*60)+datenum(1970,1,1);
ssa=neph(:,7); % SSA calculado com corrected abs637 e scat637
%
rem1=find(ssa<0.1); 
level2(rem,2:size(level2,2))=NaN;
%
% Save Level 3 data 
Instrument='Aethalometer';
Model='AE33';
serial_num='S0200193, 1/4" Cu line, cutoff PM2.5, Magee dryer, in series with NephAurora, 4LPM';
Site=site;
Val_level='3';
Proc_vers='1.1';
labels={'BC370','BC470','BC520','BC590','BC660','BC880','BC950','BB','abs370','abs470','abs520','abs590','abs660','abs880','abs950','abs637','AAE'};
units={'(ug.m-3)','(ug.m-3)','(ug.m-3)','(ug.m-3)','(ug.m-3)','(ug.m-3)','(ug.m-3)','(percent)','(Mm-1)','(Mm-1)','(Mm-1)','(Mm-1)','(Mm-1)','(Mm-1)','(Mm-1)','(Mm-1)','(adim)'};
export_LFA_GoAmazon_jan2015(labels,units,Instrument,Model,serial_num,Site,Val_level,Proc_vers,origsize,level2(:,1),level2(:,2:size(level2,2)))
%


