%%% This script loops over each data file extracts absorbance, transmittance, and PSAP instrument data.
%%% Angstrom exponenets are then calculated from this.  Half-hour rolling
%%% average data vectors are also calculated here.
%%% D Bonanno, 2018
%% Initialization
%these must start as empty matricies to allow the vertical matrix
%concatenation to work
PSAPtime = [];

PSAPspotsize = [];
PSAPimpactor = [];
PSAPflowrate = [];

PSAPabsred = [];
qc_absR = [];
PSAPtransred = [];

PSAPabsgreen = [];
qc_absG = [];
PSAPtransgreen = [];

PSAPabsblue = [];
qc_absB = [];
PSAPtransblue = [];

qc_rawblue = [];

%% Looping over each file in BonDir
%this makes sure you can run this program from anywhere, you don't have to
%be in the right folder anymore (as long as you are on UOP computer);
cd('D:\Google Drive\Projects\HI-SCALE SGP\Data\Peripheral Data\PSAP_matthewdownload');

%% Uncomment if you need to use this script on another folder
% datadir = uipickfiles();
% cd(datadir{1});

BonDir=dir;

for i=1:length(BonDir)
    %% Testing if file "i" is a .nc file
    ncdatatest = strfind(BonDir(i).name,'.nc');
    if isempty(ncdatatest)
        continue
    end
    
    %% Finding Datetime Vector
%     periodPos = strfind(BonDir(i).name,'.2017');
	periodPos = strfind(BonDir(i).name,'.2016');  %%%%%%%%%%%%%%%%THIS IS SUPPOSED TO BE TEMPORARY
    startDate = BonDir(i).name(periodPos+1:periodPos+8);
    startDate = datetime(startDate,'InputFormat','yyyyMMdd');
    startDatevec_single = datevec(startDate);
           
    PSAPtimetest = ncread(BonDir(i).name,'time');%% seconds since midnight
    PSAPtimesize = size(PSAPtimetest,1);
    
    armDateVec = repmat(startDatevec_single,PSAPtimesize,1);
    
    armDateVec(:,6) = PSAPtimetest;
    armDateVec = datetime(armDateVec);  
    
    %making large time vector
    PSAPtime = [PSAPtime ; armDateVec];
    
    %% Pulling out Relevant Variables
    %instrument values
    currPSAPimpactor = ncread(BonDir(i).name,'impactor_state');
    currPSAPflowrate = ncread(BonDir(i).name,'sample_flow_rate'); %length(PSAPtimetest)x1 value in liter/min
    currPSAPspotsize = ncread(BonDir(i).name,'spot_size_area'); % 1x1 value in m^2
    
    %red
    currPSAPAbsRed = ncread(BonDir(i).name,'Ba_R_Weiss');
	currAred_QC = ncread(BonDir(i).name,'qc_Ba_R_Weiss');
% 	currPSAPAbsRed = ncread(BonDir(i).name,'absorbance_red');%%%%%%%%%%%%%%%%THIS IS SUPPOSED TO BE TEMPORARY
    currPSAPTransRed = ncread(BonDir(i).name,'transmittance_red');
    
    %green
    currPSAPAbsGreen = ncread(BonDir(i).name,'Ba_G_Weiss'); 
	currAgreen_QC = ncread(BonDir(i).name,'qc_Ba_G_Weiss');
% 	currPSAPAbsGreen = ncread(BonDir(i).name,'absorbance_green'); %%%%%%%%%%%%%%%%THIS IS SUPPOSED TO BE TEMPORARY
    currPSAPTransGreen = ncread(BonDir(i).name,'transmittance_green');
    
    %blue
    currPSAPAbsBlue = ncread(BonDir(i).name,'Ba_B_Weiss');
	currAblue_QC = ncread(BonDir(i).name,'qc_Ba_B_Weiss');
	currAblue_QCraw = ncread(BonDir(i).name,'qc_Ba_B_raw');
% 	currPSAPAbsBlue = ncread(BonDir(i).name,'absorbance_blue');%%%%%%%%%%%%%%%%THIS IS SUPPOSED TO BE TEMPORARY
    currPSAPTransBlue = ncread(BonDir(i).name,'transmittance_blue');
    
    %making large vector over all files
    PSAPimpactor = [PSAPimpactor ; currPSAPimpactor];
    PSAPflowrate = [PSAPflowrate ; currPSAPflowrate];
    PSAPspotsize = [PSAPspotsize ; currPSAPspotsize];
    
    PSAPabsred = [PSAPabsred ; currPSAPAbsRed];
	qc_absR = [qc_absR ; currAred_QC];
    PSAPtransred = [PSAPtransred ; currPSAPTransRed];
    
    PSAPabsgreen = [PSAPabsgreen ; currPSAPAbsGreen];
	qc_absG = [qc_absG ; currAgreen_QC];
    PSAPtransgreen = [PSAPtransgreen ; currPSAPTransGreen];
    
    PSAPabsblue = [PSAPabsblue ; currPSAPAbsBlue];
	qc_absB = [qc_absB ; currAblue_QC];
	qc_rawblue = [qc_rawblue ; currAblue_QCraw];
    PSAPtransblue = [PSAPtransblue ; currPSAPTransBlue];
    
end
%% Raw Data Corrections
% correctionfactor_red = ones(size(PSAPabsred));
% correctionfactor_green = ones(size(PSAPabsgreen));
% correctionfactor_blue = ones(size(PSAPabsblue));
% 
% correctionfactor_red(qc_absR >0) = NaN;
% correctionfactor_green(qc_absG >0) = NaN;
% correctionfactor_blue(qc_absB >0) = NaN;
% 
% PSAPabsred = PSAPabsred .* correctionfactor_red;
% PSAPabsgreen = PSAPabsgreen .* correctionfactor_green;
% PSAPabsblue = PSAPabsblue .* correctionfactor_blue;



%making -9999 data points (data errors) NaNs
% PSAPimpactor(PSAPimpactor == -9999) = NaN;
% PSAPflowrate(PSAPflowrate == -9999) = NaN;
% PSAPspotsize(PSAPspotsize == -9999) = NaN;

PSAPabsred(PSAPabsred == -9999) = NaN;
% PSAPabsred(PSAPabsred < 0) = NaN;
PSAPabsred(PSAPabsred > 35) = NaN;
% PSAPtransred(PSAPtransred == -9999) = NaN;

PSAPabsgreen(PSAPabsgreen == -9999) = NaN;
% PSAPabsgreen(PSAPabsgreen < 0) = NaN;
PSAPabsgreen(PSAPabsgreen > 35) = NaN;
% PSAPtransgreen(PSAPtransgreen == -9999) = NaN;

PSAPabsblue(PSAPabsblue == -9999) = NaN;
% PSAPabsblue(PSAPabsblue < 0) = NaN;
PSAPabsblue(PSAPabsblue > 35) = NaN;
% PSAPtransblue(PSAPtransblue == -9999) = NaN;

%% Calculations
%correcting PSAP absorbance data
% SigmaRPSAP = PSAPcorrection(PSAPabsred,PSAPtransred,PSAPflowrate,PSAPspotsize);
% SigmaGPSAP = PSAPcorrection(PSAPabsgreen,PSAPtransgreen,PSAPflowrate,PSAPspotsize);
% SigmaBPSAP = PSAPcorrection(PSAPabsblue,PSAPtransblue,PSAPflowrate,PSAPspotsize);

%rolling average of corrected absorbance data (half an hour) 
[avgPSAPtime,avgred] = RollingAverage(PSAPtime,PSAPabsred);
[~,avggreen] = RollingAverage(PSAPtime,PSAPabsgreen);
[~,avgblue] = RollingAverage(PSAPtime,PSAPabsblue);

%calculating Angstrom Exponent
%keeping longer wavelength (lower energy) data on top of ratios
RBangstrom = -((log((PSAPabsred)./(PSAPabsblue)))/(log(660/467))); % Red over Blue
GBangstrom = -((log((PSAPabsgreen)./(PSAPabsblue)))/(log(530/467))); % Green over Blue
RGangstrom = -((log((PSAPabsred)./(PSAPabsgreen)))/(log(660/530))); % Red over Green

%rolling average of angstrom exponents
[~,avgangrb] = RollingAverage(PSAPtime,RBangstrom);
[~,avganggb] = RollingAverage(PSAPtime,GBangstrom);
[~,avgangrg] = RollingAverage(PSAPtime,RGangstrom);

