% Code used for Huang, Y., et al. (2018). Fine dust emissions from Oceano
% Sand Dunes, California.
%
% 
% Main Input data: ProcessedData_OceanoDust_CalibrationDays.mat
%             ProcessedData_OceanoDust_MeasurementDays.mat
%             ProcessedData_OceanoSaltation.mat
%             OPC_diameter.mat
% Scripts used: 
%             *\scripts\

clearvars; clc; close all;
%% load input data
% % load raw secondly measurements during Oceano 2015
datapath = 'C:\Users\Yue\Documents\MATLAB\Oceano_Dunes_2017\';
dataname1 = 'ProcessedData_OceanoDust_CalibrationDays';  % load measurements of dust during June 5, 6 and 7, 2015
dataname2 = 'ProcessedData_OceanoDust_MeasurementDays'; % load measuremnts of dust during May 15 to June 4, 2015
dataname3 = 'ProcessedData_OceanoSaltation'; % load measurements of saltation during May 26 to June 4, 2015, from R.L.Martin.
dataname4 = 'OPC_diameter'; % load OPC bin size of dust geometric, after using Lorentz-Mie code
dataname5 = 'NonSandySoil_PM10Flux'; % load mass flux within PM10 geometric diameter of non-sany soils, from J.F.Kok 
dataname6 = 'NonSandySoil_PSD'; % load PSD of non-sany soils by in situ field measurements, from J.F.Kok et al. (2017)'s Fig. 1(C) 
dataname7 = 'African_PSD'; % load PSD of African aircraft measurements, from Mahowald et al. (2014)'s Fig. 5(A)
dataname8 = 'Bullard04_R35'; % load PSD of produced dust by clay-coating removal, from Bullard et al. (2004)'s Fig. 17

load(strcat(datapath,dataname1));
load(strcat(datapath,dataname2));
load(strcat(datapath,dataname3));
load(strcat(datapath,dataname4));
load(strcat(datapath,dataname5));
load(strcat(datapath,dataname6));
load(strcat(datapath,dataname7));
load(strcat(datapath,dataname8));
%% add scripts path
addpath(strcat(datapath,'scripts\'));
%% indicate output figure path
figurepath = strcat(datapath,'figures\');

fprintf(['Load raw data, and obtain two structs "OceanoDust_CalibrationDays" \n'...
    'and "OceanoDust_MeasurementDays" \n' ...
    'Press "Enter" to continue \n']);
% pause;

%% ========== Part 1: Data pre-processing  =========== 
% %====== data preprocessing on Calibration Days (June 5th, 6th and 7th, 2015)
[StartTimes, EndTimes] = TimeSeriesContinuous(OceanoDust_CalibrationDays);
% fprintf('Press "Enter" to continue \n');pause;

TimeInterval = duration(0,29,59); % 30 minutes in seconds using Timedate format
[StartTimeBlocks, EndTimeBlocks] = CreateTimeBlocks(StartTimes, EndTimes, TimeInterval);
% fprintf('Press "Enter" to continue \n');pause;

TimeInterval_in_Seconds = 1799; % TimeInterval in seconds
[OceanoDust_CaliDays_Org] = DistributeIntoBlocks(OceanoDust_CalibrationDays, StartTimeBlocks, EndTimeBlocks, TimeInterval_in_Seconds);
% fprintf('Press "Enter" to continue \n');pause;

[OceanoDust_CaliDays_Avg, OceanoDust_CaliDays_Convr] = Average_N_Convert_inBlocks(OceanoDust_CaliDays_Org);
% fprintf('Press "Enter" to continue \n');pause;

[CaliDays_5avg, CaliDays_5std, CaliDays_5se] = CalculateStatVars(OceanoDust_CaliDays_Convr);
% fprintf('Press "Enter" to continue \n');pause;

% %====== linear regression on Calibration Days (June 5th, 6th and 7th, 2015)
[a_log10, b_log10, SE_a_log10, SE_b_log10, Cov_squared_log10, OceanoDust_CaliDays_Convr_Fit] =...
    Calibrate3Days_Powerfit(OceanoDust_CaliDays_Convr, CaliDays_5avg, CaliDays_5se);
% fprintf('Press "Enter" to continue \n');pause;

% %====== data pre-processing on Measurement Days (May 26th to June 4th, 2015)
[OceanoDust_MeasDays_IDOrder] = DistributeOPC_IDOrder(OceanoDust_MeasurementDays); 
% fprintf('Press "Enter" to continue \n');pause;

[OceanoDust_MeasDays_TimePre] = DistributeIntoQBlocks(OceanoDust_MeasDays_IDOrder, OceanoSaltation.StartTimes, OceanoSaltation.EndTimes);
% fprintf('Press "Enter" to continue \n');pause;

[OceanoDust_MeasDays_Time] = Org_Avg_N_Convert_inQBlocks(OceanoDust_MeasDays_TimePre, TimeInterval_in_Seconds);
% fprintf('Press "Enter" to continue \n');pause;

% %====== linear regression on Measurement Days
[OceanoDust_MeasDays_Cali] = CalibrateMDays_Powerfit(OceanoDust_MeasDays_Time, a_log10, b_log10, SE_a_log10, SE_b_log10, Cov_squared_log10);
% fprintf('Press "Enter" to continue \n');pause;

%% ====== Part 2: Calibration plotting ============
% %====== plot of Calibration Days (June 5th, 6th and 7th, 2015)
for i = 1:7 % lope for 7 bins
    PlotCalibrationDays(OceanoDust_CaliDays_Convr, CaliDays_5avg, CaliDays_5se, OceanoDust_CaliDays_Convr_Fit, i);
    print([figurepath,'CalibrationDays_Bin_',num2str(i),'.png'],'-dpng'); 
%     fprintf('Press "Enter" to continue \n');
%     pause;    
    clf;
end
close all;

% %====== plot of Measurement Days (May 26th to June 4th, 2015)
for i = 1:7 % lope for 7 bins
    PlotMeasurementDays(OceanoDust_CaliDays_Convr, CaliDays_5avg, CaliDays_5se, OceanoDust_CaliDays_Convr_Fit, OceanoDust_MeasDays_Cali, i);
    print([figurepath,'MeasurementDays_Bin_',num2str(i),'.png'],'-dpng'); 
% fprintf('Press "Enter" to continue \n');
% pause;
clf;
end
close all;

%% ====== Part 3: Data-quality control criteria 
% %======#1 pick up data with at least N (4, for example) sensors running
N_sensor = 4;
[OceanoDust_MeasDays_NOPC] = PickupNOPCs(OceanoDust_MeasDays_Cali, OceanoSaltation, N_sensor); %OceanoSaltation struct here for time series
% fprintf('Press "Enter" to continue \n');pause;

% % ======#2 pick up data with theta within -/+45 degrees
[OceanoDust_MeasDays_Dir] = Direction45Degrees(OceanoDust_MeasDays_NOPC); 
% fprintf('Press "Enter" to continue \n');pause;

% % ======#3 pick up data after May 29th, 2015, when lowest sensor z1 at same height
[OceanoDust_MeasDays_Z1Same] = PickupTimeBlocks(OceanoDust_MeasDays_Dir);
% fprintf('Press "Enter" to continue \n');pause;

% % ======#4 pick up data without time overlapping
[OceanoDust_MeasDays_Nonlap] = PickupTimeNonOverlap(OceanoDust_MeasDays_Z1Same);
% fprintf('Press "Enter" to continue \n');pause;

%% ====== Part 4: Mass flux from concentration profile
% % ====== order the struct with increaing of height of OPCs
D_rou = 2.5*10^3; % unit: kg/m^3, Kok et al. (2017) from Nature Geoscience

[OceanoDust_MeasDays_Z] = DistributeOPC_HeightOrder(OceanoDust_MeasDays_Nonlap);
% fprintf('Press "Enter" to continue \n');pause;

% % ====== plot mass concentration profile (kg/m^3) and fit line
for i = 1:7
   PlotMassConcenProfile(OceanoDust_MeasDays_Z, D_dust, D_rou, i, figurepath); 
end
close all;

%% ====== Part 5: Diagnose and remove sea-salt deposition flux
% % ====== divide into (1) when saltation occurs, Q>0 and (2) saltation not occur (Q==0)
[OceanoDust_MeasDays_Salta, OceanoDust_MeasDays_NoSalta] = SaltationTwoTypes(OceanoDust_MeasDays_Z);
% fprintf('Press "Enter" to continue \n');pause;

% % ====== plot norm. volume PSD at each of the 4 heights to diagnoze the existence of sea salt aerosols
for N = 1:4
    PlotPSDeaOPC(OceanoDust_MeasDays_NoSalta, D_dust, D_dust_range, D_rou, N);
end
% in prepare for legend
d1 = plot(-10,-10,'o','MarkerFaceColor',[1,0,0.2],'MarkerEdgeColor','none','Markersize',7,'linewidth',1); hold on;
d2 = plot(-10,-10,'o','MarkerFaceColor',[0.2,0.6,0.2],'MarkerEdgeColor','none','Markersize',7,'linewidth',1); hold on;
d3 = plot(-10,-10,'o','MarkerFaceColor',[0,0.6,0.8],'MarkerEdgeColor','none','Markersize',7,'linewidth',1); hold on;
d4 = plot(-10,-10,'o','MarkerFaceColor',[0.8,0.6,0.4],'MarkerEdgeColor','none','Markersize',7,'linewidth',1); hold on;
the = legend([d1 d2 d3 d4],{'D1 = 0.74 m','D2 = 1.48 m','D3 = 2.29 m','D4 = 6.44 m'});
set(the,'location','Southeast','FontSize',8,'FontWeight','bold');
title('Inactive saltation (Q = 0)','FontSize',12,'FontWeight','bold');
print([figurepath,'FourHeights_PSD_nosalta.png'],'-dpng'); 
clf;
    
for N = 1:4
    PlotPSDeaOPC(OceanoDust_MeasDays_Salta, D_dust, D_dust_range, D_rou, N);
end
% in prepare for legend
d1 = plot(-10,-10,'o','MarkerFaceColor',[1,0,0.2],'MarkerEdgeColor','none','Markersize',7,'linewidth',1); hold on;
d2 = plot(-10,-10,'o','MarkerFaceColor',[0.2,0.6,0.2],'MarkerEdgeColor','none','Markersize',7,'linewidth',1); hold on;
d3 = plot(-10,-10,'o','MarkerFaceColor',[0,0.6,0.8],'MarkerEdgeColor','none','Markersize',7,'linewidth',1); hold on;
d4 = plot(-10,-10,'o','MarkerFaceColor',[0.8,0.6,0.4],'MarkerEdgeColor','none','Markersize',7,'linewidth',1); hold on;
the = legend([d1 d2 d3 d4],{'D1 = 0.74 m','D2 = 1.48 m','D3 = 2.29 m','D4 = 6.44 m'});
set(the,'location','Southeast','FontSize',8,'FontWeight','bold');
title('Active saltation (Q > 0)','FontSize',12,'FontWeight','bold');
print([figurepath,'FourHeights_PSD_salta.png'],'-dpng'); 
clf; 

% % ====== calculate mass flux using concentration profile
N_height = 1:2;
[DustFlux_Salta] = NumConcenProfile2Flux(OceanoDust_MeasDays_Salta, D_dust, D_rou, N_height); %input number concentration(#/m^3), output mass flux(kg/m^2/s)

[DustFlux_NoSalta] = NumConcenProfile2Flux(OceanoDust_MeasDays_NoSalta, D_dust, D_rou, N_height); %input number concentration(#/m^3), output mass flux(kg/m^2/s)

% % ====== plot mass flux within both saltation flux and non-saltation flux 
for i = 1:7 % lope for 7 bins; i=8 means the sum of mass flux of 7 bins
    PlotAllMassFlux(DustFlux_Salta, i, 1);% "1" = 'closed circles' as plot Marker
    PlotAllMassFlux(DustFlux_NoSalta, i, -1); % "-1" = 'open cycles' as plot Marker
    title(['Bin - ',num2str(i),' (Plan: D1 and D2)']);
    print([figurepath,'AllMassFlux_Bin_',num2str(i),'_2dn.png'],'-dpng'); 
    clf;
end
close all;

% % ====== subtract sea-salt effect on the mass flux, uing (1) powerfit and (2) assume constant depostion flux with shear velocity
[DustFlux_rSS_powerfit, DustFlux_Salta_rSS_powerfit_Plot, DustFlux_NoSalta_rSS_powerfit_Plot] = RemoveSeaSalt_powerfit(DustFlux_Salta, DustFlux_NoSalta);% input unit mass flux(kg/m^2/s), output unit mass flux(kg/m^2/s)

[DustFlux_rSS_constfit, DustFlux_Salta_rSS_constfit_Plot, DustFlux_NoSalta_rSS_constfit_Plot] = RemoveSeaSalt_constantfit(DustFlux_Salta, DustFlux_NoSalta);% input unit mass flux(kg/m^2/s), output unit mass flux(kg/m^2/s)

% % ====== plot of sea-salt mass flux removal
for i = 1:7 % lope for 7 bins
    PlotRemoveSeaSalt(DustFlux_Salta_rSS_powerfit_Plot, DustFlux_NoSalta_rSS_powerfit_Plot, i, 1);% 1 indicates powerfit, in preparation for legend
    print([figurepath,'SeaSaltFlux_Bin_',num2str(i),'_powerfit.png'],'-dpng'); 
    clf;
end
for i = 1:7 % lope for 7 bins
    PlotRemoveSeaSalt(DustFlux_Salta_rSS_constfit_Plot, DustFlux_NoSalta_rSS_constfit_Plot, i, -1);% -1 indicates constfit, in preparation for legend    
    print([figurepath,'SeaSaltFlux_Bin_',num2str(i),'_constfit.png'],'-dpng'); 
    clf;
end
close all;

%% ====== Part 6: Compare dust mass flux with other non-sandy soils
% % ===== plot mass flux of each bin after addressing sea salt issue
for i = 1:7 % lope for 7 bins; i=8 means the sum of mass flux of 7 bins
    PlotDustMassFlux(DustFlux_rSS_powerfit, i);
    print([figurepath,'DustFlux_Bin_',num2str(i),'_powerfit.png'],'-dpng'); 
    clf;
end

for i = 1:7% lope for 7 bins; i=8 means the sum of mass flux of 7 bins
    PlotDustMassFlux(DustFlux_rSS_constfit, i);    
    print([figurepath,'DustFlux_Bin_',num2str(i),'_constfit.png'],'-dpng'); 
    clf;
end
% % ===== calculate summed flux over seveal bins
bin_nn = 6; % sum over the first 6 bins
PlotDustBulkMassFlux(DustFlux_rSS_powerfit, bin_nn);
print([figurepath,'DustBulkFlux_first',num2str(bin_nn),'bins_powerfit.png'],'-dpng'); 
clf;

PlotDustBulkMassFlux(DustFlux_rSS_constfit, bin_nn);
print([figurepath,'DustBulkFlux_first',num2str(bin_nn),'bins_constfit.png'],'-dpng'); 
clf;
close all;

% % ===== plot PM_{10,g} mass flux and other measurements
PlotDustFluxPM10Comparison(DustFlux_rSS_powerfit, D_dust_range, NonSandyFlux);
print([figurepath,'MassFluxComparison_PM10_powerfit.png'],'-dpng');   
clf;

PlotDustFluxPM10Comparison(DustFlux_rSS_constfit, D_dust_range, NonSandyFlux);
print([figurepath,'MassFluxComparison_PM10_constfit.png'],'-dpng');   
clf;

% % ===== plot PM_{9,g} mass flux and other measurements (first 6 bins)
PlotDustFluxPM9Comparison(DustFlux_rSS_powerfit, NonSandyFlux);
print([figurepath,'MassFluxComparison_PM9_powerfit.png'],'-dpng');   
clf;

PlotDustFluxPM9Comparison(DustFlux_rSS_constfit, NonSandyFlux);
print([figurepath,'MassFluxComparison_PM9_constfit.png'],'-dpng');   
clf;

% % ===== convert geometric diameter to aerodynamic diameter in preparation for PM_{2.5, a} calculation
[Dae_dust] = GeometricD2AerodynamicD(D_dust, D_rou);
[Dae_dust_range] = GeometricD2AerodynamicD(D_dust_range, D_rou);

% % ===== plot PM_{2.5,a} mass flux and other measurements
PlotDustFluxPM2_5Comparison(DustFlux_rSS_powerfit, Dae_dust_range, NonSandyFlux, 0.052); 
% 0.052 as the weight Factor for PM_{10,g} to PM_{2.5,a} for non-sandy soils, after email to J.F.kok (10/5/17)
print([figurepath,'MassFluxComparison_PM2_5_powerfit.png'],'-dpng');   
clf;

PlotDustFluxPM2_5Comparison(DustFlux_rSS_constfit, Dae_dust_range, NonSandyFlux, 0.052); 
% 0.052 as the weight Factor for PM_{10,g} to PM_{2.5,a} for non-sandy soils, after email to J.F.kok (10/5/17)
print([figurepath,'MassFluxComparison_PM2_5_constfit.png'],'-dpng');   
clf;
close all;

%% ====== Part 7: PSD and alpha = Fd/Q
% %====== calculate sandblasting efficiency for the first 6 bins and the bulk alpha (sum of the 6 bins)
[Sandblasting_alpha_powerfit] = CalculateAlpha(DustFlux_rSS_powerfit);

[Sandblasting_alpha_constfit] = CalculateAlpha(DustFlux_rSS_constfit);

% %====== plot sandblasting efficiency of each individual bin
for i = 1:7 % lope over 7 bins
   PlotAlpha(Sandblasting_alpha_powerfit, i, D_dust_range, D_rou);
   print([figurepath,'Alpha_Bin_',num2str(i),'_powerfit.png'],'-dpng'); 
   clf;

   PlotAlpha(Sandblasting_alpha_constfit, i, D_dust_range, D_rou);
   print([figurepath,'Alpha_Bin_',num2str(i),'_constfit.png'],'-dpng'); 
   clf;
end
close all;
% %====== plot the bulk sandblasting efficiency 
PlotBulkAlpha(DustFlux_rSS_powerfit, bin_nn, D_dust_range, D_rou);
print([figurepath,'BulkAlpha_first',num2str(bin_nn),'bins_powerfit.png'],'-dpng'); 
clf;

PlotBulkAlpha(DustFlux_rSS_constfit, bin_nn, D_dust_range, D_rou);
print([figurepath,'BulkAlpha_first',num2str(bin_nn),'bins_constfit.png'],'-dpng'); 
clf;

% % % % % % % ====== calculate PSD
% % % % % [Dust_dvdlnd_powerfit] = CalculatePSD(DustFlux_rSS_powerfit, D_dust_range, bin_nn);
% % % % % 
% % % % % [Dust_dvdlnd_constfit] = CalculatePSD(DustFlux_rSS_constfit, D_dust_range, bin_nn);
% % % % % 
% % % % % % %====== plot PSD, mean of all samples
% % % % % PlotPSDMean(Dust_dvdlnd_powerfit, D_dust(1:bin_nn,1), [0,0,0]);
% % % % % print([figurepath,'PSD_powerfit.png'],'-dpng'); 
% % % % % clf;
% % % % % 
% % % % % PlotPSDMean(Dust_dvdlnd_constfit, D_dust(1:bin_nn,1), [0,0,0]);
% % % % % print([figurepath,'PSD_constfit.png'],'-dpng'); 
% % % % % clf;
% % % % % close all;
% % % % % 
% % % % % % % ====== plot PSD within 3 ustar bins
% % % % % PlotPSD_in3Groups(Dust_dvdlnd_powerfit, D_dust(1:bin_nn,1));
% % % % % print([figurepath,'PSD_powerfit_3groups.png'],'-dpng'); 
% % % % % clf;
% % % % % 
% % % % % PlotPSD_in3Groups(Dust_dvdlnd_constfit, D_dust(1:bin_nn,1));
% % % % % print([figurepath,'PSD_constfit_3groups.png'],'-dpng'); 
% % % % % clf;
% % % % % close all;
% % % % % 
% % % % % % %====== plot PSD (mean) with other non-sandy soils
% % % % % PlotPSDComparison(Dust_dvdlnd_powerfit, NonSandyPSD, D_dust(1:bin_nn,1));
% % % % % print([figurepath,'PSD_nonsandy_comparison_powerfit.png'],'-dpng'); 
% % % % % clf;
% % % % % 
% % % % % PlotPSDComparison(Dust_dvdlnd_constfit, NonSandyPSD, D_dust(1:bin_nn,1));
% % % % % print([figurepath,'PSD_nonsandy_comparison_constfit.png'],'-dpng'); 
% % % % % clf;
% % % % % close all;
% % % % % 
% % % % % % %====== plot PSD (mean) with other non-sandy soils, and Bullard et al. (2004)
% % % % % PlotPSDComparisonBullard(Dust_dvdlnd_powerfit, NonSandyPSD, D_dust(1:bin_nn,1), D_save72h, dvdlnd_save72h, D_save120h, dvdlnd_save120h);
% % % % % print([figurepath,'PSD_nonsandyNbullard_comparison_powerfit.png'],'-dpng'); 
% % % % % clf;
% % % % % 
% % % % % PlotPSDComparisonBullard(Dust_dvdlnd_constfit, NonSandyPSD, D_dust(1:bin_nn,1), D_save72h, dvdlnd_save72h, D_save120h, dvdlnd_save120h);
% % % % % print([figurepath,'PSD_nonsandyNbullard_comparison_constfit.png'],'-dpng'); 
% % % % % clf;
% % % % % close all;
% % % % % 
% % % % % % %====== plot PSD (mean) with African aircraft measurements 
% % % % % PlotPSDCompareAfrica(Dust_dvdlnd_powerfit, AfricanPSD, D_dust(1:bin_nn,1));
% % % % % print([figurepath,'PSD_Africa_comparison_powerfit.png'],'-dpng'); 
% % % % % clf;
% % % % % 
% % % % % PlotPSDCompareAfrica(Dust_dvdlnd_constfit, AfricanPSD, D_dust(1:bin_nn,1));
% % % % % print([figurepath,'PSD_Africa_comparison_constfit.png'],'-dpng'); 
% % % % % clf;
% % % % % close all;
% % % % % 
% % % % % % %====== plot PSD (mean) with African aircraft measurements , and Bullard et al. (2004)
% % % % % PlotPSDCompareAfricaBullard(Dust_dvdlnd_powerfit, AfricanPSD, D_dust(1:bin_nn,1), D_save72h, dvdlnd_save72h, D_save120h, dvdlnd_save120h);
% % % % % print([figurepath,'PSD_AfricaNbullard_comparison_powerfit.png'],'-dpng'); 
% % % % % clf;
% % % % % 
% % % % % PlotPSDCompareAfricaBullard(Dust_dvdlnd_constfit, AfricanPSD, D_dust(1:bin_nn,1), D_save72h, dvdlnd_save72h, D_save120h, dvdlnd_save120h);
% % % % % print([figurepath,'PSD_AfricaNbullard_comparison_constfit.png'],'-dpng'); 
% % % % % clf;
% % % % % close all;
