diff --git a/lib/interface/Run_picasso_test.m b/lib/interface/Run_picasso_test.m
new file mode 100644
index 00000000..33dc841e
--- /dev/null
+++ b/lib/interface/Run_picasso_test.m
@@ -0,0 +1,201 @@
+% File is originally from Jimenez to start PPC with Inputfiles
+% Modified by Nathan to run in parallel and do a partial run first to get LC values
+% These methods should probably get implemented in "picassoProcTodolist.m" instead of this additional file
+
+
+clc
+clear all
+
+path(path,'E:\Antarctis CCN\picasso_update\Pollynet_Processing_Chain\lib\interface')
+path(path,'E:\Antarctis CCN\picasso_update\Pollynet_Processing_Chain\')
+ 
+javaaddpath('E:\Antarctis CCN\picasso_update\Pollynet_Processing_Chain\include\sqlite-jdbc-3.30.1.jar')
+
+
+initPicassoToolbox
+
+
+Inputfiles={
+'2023_02_27_Mon_ARI_00_00_01.nc'
+'2023_02_27_Mon_ARI_06_00_01.nc'
+'2023_02_27_Mon_ARI_12_00_01.nc'
+'2023_02_27_Mon_ARI_18_00_01.nc'  
+};
+
+
+lidar_path='E:\Antarctis CCN\Datasets\2023 - Rawdata\Cases\';
+
+% Set number of parallel workers (0 and 1 will be ignored -> program will run sequentially)
+par_pool_size = 0;
+
+if par_pool_size > 1
+    pool_obj = gcp("nocreate");
+    if isempty(pool_obj)
+        pool_obj = parpool('local',par_pool_size);
+    elseif pool_obj.NumWorkers ~= par_pool_size
+        delete(gcp("nocreate"));
+        pool_obj = parpool('local',par_pool_size);
+    end
+end
+
+% First run to get the LC values / stops midway
+if par_pool_size > 1
+    % Preallocate output container
+    LidarLCexport_temp = cell(1, length(Inputfiles));
+
+    parfor i = 1:length(Inputfiles)
+        javaaddpath('E:\Antarctis CCN\picasso_update\Pollynet_Processing_Chain\include\sqlite-jdbc-3.30.1.jar');
+
+        InputFiles_lidar = fullfile(lidar_path, Inputfiles{i});  % Construct path
+
+        lidarFileExists = isfile(InputFiles_lidar);
+        borrar = 0;
+
+        if ~lidarFileExists
+            borrar = 1;
+            try
+                unzip(strcat(InputFiles_lidar, '.zip'), lidar_path);
+            catch
+                fprintf('Unable to unzip %s.zip\n', InputFiles_lidar);
+                continue;
+            end
+        end
+
+        try
+            [report, LidarLCexport] = picassoProcV3( ...
+                InputFiles_lidar, ...
+                'arielle', ...
+                'E:\Antarctis CCN\picasso\template_pollynet_processing_chain_overlaptest2_Neumayer.json', ...
+                true, 0);
+
+            LidarLCexport_temp{i} = LidarLCexport;
+
+        catch ME
+            fprintf('Error processing file %s: %s\n', InputFiles_lidar, ME.message);
+        end
+
+        if borrar
+            try
+                delete(InputFiles_lidar);
+            catch
+                fprintf('Could not delete %s\n', InputFiles_lidar);
+            end
+        end
+    end
+    
+    LidarLCexport_all = [LidarLCexport_temp{:}];
+    
+else
+    LidarLCexport_all=[];
+
+    for i=1:length(Inputfiles)
+        InputFiles_lidar=strcat(lidar_path,Inputfiles{i});
+
+        if isfile(InputFiles_lidar)
+            borrar=0;
+        else
+            borrar=1;
+            try
+                unzip(strcat(InputFiles_lidar,'.zip'),lidar_path)
+            catch
+                message='unable to unzip'
+                continue
+            end
+        end
+
+        [report,LidarLCexport]=picassoProcV3(InputFiles_lidar, 'arielle', 'E:\Antarctis CCN\picasso\template_pollynet_processing_chain_overlaptest2_Neumayer.json',true,0);
+
+        LidarLCexport_all=[LidarLCexport_all LidarLCexport];
+
+        if borrar
+            delete(InputFiles_lidar)
+        end
+    end
+end
+
+%%
+LidarLC532=nanmean(LidarLCexport_all(1,:));
+LidarLC607=nanmean(LidarLCexport_all(2,:));
+LidarLC532NR=nanmean(LidarLCexport_all(3,:));
+LidarLC607NR=nanmean(LidarLCexport_all(4,:));
+LidarLC1064=nanmean(LidarLCexport_all(5,:));
+
+LC=[LidarLC532,LidarLC607,LidarLC532NR,LidarLC607NR, LidarLC1064];
+
+%%
+% Second run with the LC values obtained / full runthrough
+if par_pool_size > 1
+    LidarLCexport_temp = cell(1, length(Inputfiles));
+
+    parfor i = 1:length(Inputfiles)
+        javaaddpath('E:\Antarctis CCN\picasso\Pollynet_Processing_Chain\include\sqlite-jdbc-3.30.1.jar');
+
+        InputFiles_lidar = fullfile(lidar_path, Inputfiles{i});
+
+        lidarFileExists = isfile(InputFiles_lidar);
+        borrar = 0;
+
+        if ~lidarFileExists
+            borrar = 1;
+            try
+                unzip(strcat(InputFiles_lidar, '.zip'), lidar_path);
+            catch
+                fprintf('Unable to unzip %s.zip\n', InputFiles_lidar);
+                continue;
+            end
+        end
+
+        try
+            [report, LidarLCexport] = picassoProcV3( ...
+                InputFiles_lidar, ...
+                'arielle', ...
+                'E:\Antarctis CCN\picasso\template_pollynet_processing_chain_overlaptest2_Neumayer.json', ...
+                false, LC);
+
+            LidarLCexport_temp{i} = LidarLCexport;
+
+        catch ME
+            fprintf('Error processing file %s: %s\n', InputFiles_lidar, ME.message);
+        end
+
+        if borrar
+            try
+                delete(InputFiles_lidar);
+            catch
+                fprintf('Could not delete %s\n', InputFiles_lidar);
+            end
+        end
+    end
+
+    % Combine all LidarLCexport structs if needed
+    LidarLCexport_all = [LidarLCexport_temp{:}];
+else
+    for i=1:length(Inputfiles)
+        InputFiles_lidar=strcat(lidar_path,Inputfiles{i});
+
+        if isfile(InputFiles_lidar)
+            borrar=0;
+        else
+            borrar=1;
+            try
+                unzip(strcat(InputFiles_lidar,'.zip'),lidar_path)
+            catch
+                message='unable to unzip'
+                continue
+            end
+        end
+
+        LC=[LidarLC532,LidarLC607,LidarLC532NR,LidarLC607NR, LidarLC1064];
+
+        [report,temp]=picassoProcV3(InputFiles_lidar, 'arielle', 'E:\Antarctis CCN\picasso\template_pollynet_processing_chain_overlaptest2_Neumayer.json',false, LC);
+
+
+        if borrar
+            delete(InputFiles_lidar)
+        end
+    end
+end
+
+
+clear all
+return
diff --git a/lib/interface/picassoProcV3.m b/lib/interface/picassoProcV3.m
index 33c49150..f4a83e76 100644
--- a/lib/interface/picassoProcV3.m
+++ b/lib/interface/picassoProcV3.m
@@ -1,4 +1,4 @@
-function [report] = picassoProcV3(pollyDataFile, pollyType, PicassoConfigFile, varargin)
+function [report,LidarLCexport] = picassoProcV3(pollyDataFile, pollyType, PicassoConfigFile, calibration_on, LC_input, varargin)
 % PICASSOPROCV3 Picasso processing main program (Version 3.0).
 %End changed to UNIX Line changed to LF
 % USAGE:
@@ -1673,6 +1673,11 @@
         continue;
     end
 
+    %here the lower end of the extinction profiles is set to contant values
+    %better be done via config file
+    sig532(1:10)=sig532(10);
+    sig607(1:10)=sig607(10);
+    
     thisAerExt532_raman_tmp = thisAerExt532_raman;
     thisAerExt532_raman(1:hBaseInd532) = thisAerExt532_raman(hBaseInd532);
     [thisAerBsc532_raman, thisAerBscStd532_raman, ~] = pollyRamanBsc_smart_MC(data.distance0, sig532, sig607, thisAerExt532_raman, PollyConfig.angstrexp, mExt532(iGrp,:), mBsc532(iGrp,:),mExt607(iGrp,:), mBsc607(iGrp,:), refH532,      PollyConfig.refBeta532, PollyConfig.smoothWin_raman_532,  true, 532, 607,  bg532, bg607, thisAerExtStd532_raman, sigma_angstroem, MC_count, 'monte-carlo');
@@ -1738,6 +1743,11 @@
     if (SNRRef1064 < PollyConfig.minRamanRefSNR1064) || (SNRRef607 < PollyConfig.minRamanRefSNR607)
         continue;
     end
+    
+    %again lower end of the extinction profiles is set to contant values
+    %better be done via config file
+    sig1064(1:10)=sig1064(10);
+    sig607(1:10)=sig607(10);
 
     thisAerExt1064_raman_tmp = thisAerExt1064_raman;
     thisAerExt1064_raman(1:hBaseInd1064) = thisAerExt1064_raman(hBaseInd1064);
@@ -2558,7 +2568,7 @@
     data.LRStd355_OC_raman(iGrp, :) = thisLRStd355_OC_raman;
 
 end
-clearvars bgOLCor355 bgOLCor355 bgOLCor387 sigOLCor387
+clearvars bgOLCor355 bgOLCor355 bgOLCor387 %sigOLCor387
 %% Raman method (overlap corrected 532 nm)
 data.aerBsc532_OC_raman = NaN(size(clFreGrps, 1), length(data.height));
 data.aerBscStd532_OC_raman = NaN(size(clFreGrps, 1), length(data.height));
@@ -2689,7 +2699,7 @@
     data.LRStd1064_OC_raman(iGrp, :) = thisLRStd1064_OC_raman;
 
 end
-clearvars bgOLCor607 bgOLCor1064 sigOLCor607
+clearvars bgOLCor607 bgOLCor1064 %sigOLCor607
 %% Volume depolarization ratio new implemantation 
 if flagGHK
     print_msg('Calculating volume depolarization ratio with GHK.\n', 'flagTimestamp', true);
@@ -4177,6 +4187,27 @@
 
 print_msg('Finish\n', 'flagTimestamp', true);
 
+%% export lidar constants
+if calibration_on
+    LidarLCexport(1,:)=data.LC.LC_raman_532;
+    LidarLCexport(2,:)=data.LC.LC_raman_607;
+    LidarLCexport(3,:)=data.LC.LC_raman_532_NR;
+    LidarLCexport(4,:)=data.LC.LC_raman_607_NR;
+    LidarLCexport(5,:)=data.LC.LC_raman_1064;
+    return
+else
+    data.LCUsed.LCUsed532=LC_input(1);
+    data.LCUsed.LCUsed607=LC_input(2);
+    data.LCUsed.LCUsed532NR=LC_input(3);
+    data.LCUsed.LCUsed607NR=LC_input(4);
+    data.LCUsed.LCUsed1064=LC_input(5);
+    LidarLCexport(1,:)=data.LCUsed.LCUsed532;
+    LidarLCexport(2,:)=data.LCUsed.LCUsed607;
+    LidarLCexport(3,:)=data.LCUsed.LCUsed532NR;
+    LidarLCexport(4,:)=data.LCUsed.LCUsed607NR;
+    LidarLCexport(5,:)=data.LCUsed.LCUsed1064;
+end
+
 %% attenuated backscatter
 print_msg('Start calculating attenuated backscatter.\n', 'flagTimestamp', true);
 
@@ -4228,17 +4259,17 @@
     data.att_beta_OC_1064(:, data.depCalMask) = NaN;
 end
 
-% att_beta_OC_387 = NaN(length(data.distance0), length(data.mTime));
-% if (sum(flag387FR) == 1)
-%     att_beta_OC_387 = sigOLCor387 .* repmat(transpose(data.distance0), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed387;
-%     att_beta_OC_387(:, data.depCalMask) = NaN;
-% end
+att_beta_OC_387 = NaN(length(data.distance0), length(data.mTime));
+if (sum(flag387FR) == 1) %this might be incorrect flag
+    att_beta_OC_387 = sigOLCor387 .* repmat(transpose(data.distance0), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed387;
+    att_beta_OC_387(:, data.depCalMask) = NaN;
+end
 
-% att_beta_OC_607 = NaN(length(data.distance0), length(data.mTime));
-% if (sum(flag607FR) == 1)
-%     att_beta_OC_607 = sigOLCor607 .* repmat(transpose(data.distance0), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed607;
-%     att_beta_OC_607(:, data.depCalMask) = NaN;
-% end
+att_beta_OC_607 = NaN(length(data.distance0), length(data.mTime));
+if (sum(flag607FR) == 1) %this might be incorrect flag
+    att_beta_OC_607 = sigOLCor607 .* repmat(transpose(data.distance0), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed607;
+    att_beta_OC_607(:, data.depCalMask) = NaN;
+end
 
 data.att_beta_NR_355 = NaN(length(data.distance0), length(data.mTime));
 if (sum(flag355NR) == 1)
@@ -4252,6 +4283,19 @@
     data.att_beta_NR_532(:, data.depCalMask) = NaN;
 end
 
+% added 387 and 607 near-range for NR quasi retrieval V2
+att_beta_NR_387 = NaN(length(data.distance0), length(data.mTime));
+if (sum(flag387NR) == 1)
+    att_beta_NR_387 = squeeze(data.signal(flag387NR, :, :)) .* repmat(transpose(data.distance0), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed387NR;
+    att_beta_NR_387(:, data.depCalMask) = NaN;
+end
+
+att_beta_NR_607 = NaN(length(data.distance0), length(data.mTime));
+if (sum(flag607NR) == 1)
+    att_beta_NR_607 = squeeze(data.signal(flag607NR, :, :)) .* repmat(transpose(data.distance0), 1, length(data.mTime)).^2 / data.LCUsed.LCUsed607NR;
+    att_beta_NR_607(:, data.depCalMask) = NaN;
+end
+
 print_msg('Finish.\n', 'flagTimestamp', true);
 
 %% Volume linear depolarization ratio with high temporal resolution
@@ -4440,6 +4484,136 @@
     [data.qsiBsc1064V1, ~] = quasiRetrieval(data.height, att_beta_1064_qsi, mExt1064, mBsc1064, PollyConfig.LR1064, 'nIters', 6);
 end
 
+% quasi-retrieved backscatter at 355 nm overlap corrected
+data.qsiBscOC355V1 = NaN(length(data.height), length(data.mTime));
+att_beta_OC_355_qsi = data.att_beta_OC_355;
+if (sum(flag355t) == 1)
+    att_beta_OC_355_qsi(data.quality_mask_355 ~= 0) = NaN;
+    att_beta_OC_355_qsi = smooth2(att_beta_OC_355_qsi, PollyConfig.quasi_smooth_h(flag355t), PollyConfig.quasi_smooth_t(flag355t));
+
+    % Rayleigh scattering
+%---------------achtung
+    [mBsc355, mExt355] = rayleigh_scattering(355, pressure, temperature + 273.15, 380, 70);
+    mBsc355 = transpose(interp2(TimeMg, HeightMg, mBsc355, mTimeg, Heightg, 'linear'));
+    mExt355 = transpose(interp2(TimeMg, HeightMg, mExt355, mTimeg, Heightg, 'linear'));
+    data.quasiAttri.flagGDAS1 = strcmpi(thisMeteorAttri.dataSource, 'gdas1');
+    data.quasiAttri.meteorSource = thisMeteorAttri.dataSource;
+    data.quasiAttri.timestamp = thisMeteorAttri.datetime;
+
+    hIndOL = find(data.height >= PollyConfig.heightFullOverlap(flag355t), 1);
+    if ~ isempty(hIndOL)
+        att_beta_OC_355_qsi(1:hIndOL, :) = repmat(att_beta_OC_355_qsi(hIndOL, :), hIndOL, 1);
+    else
+        warning('Full overlap height is too large.');
+    end
+
+    [data.qsiBscOC355V1, ~] = quasiRetrieval(data.height, att_beta_OC_355_qsi, mExt355, mBsc355, PollyConfig.LR355, 'nIters', 6);
+end
+
+% quasi-retrieved backscatter at 532 nm overlap corrected
+data.qsiBscOC532V1 = NaN(length(data.height), length(data.mTime));
+att_beta_OC_532_qsi = data.att_beta_OC_532;
+if (sum(flag532t) == 1)
+    att_beta_OC_532_qsi(data.quality_mask_532 ~= 0) = NaN;
+    att_beta_OC_532_qsi = smooth2(att_beta_OC_532_qsi, PollyConfig.quasi_smooth_h(flag532t), PollyConfig.quasi_smooth_t(flag532t));
+
+    % Rayleigh scattering
+    [mBsc532, mExt532] = rayleigh_scattering(532, pressure, temperature + 273.15, 380, 70);
+  %achtung
+    mBsc532 = transpose(interp2(TimeMg, HeightMg, mBsc532, mTimeg, Heightg, 'linear'));
+    mExt532 = transpose(interp2(TimeMg, HeightMg, mExt532, mTimeg, Heightg, 'linear'));
+    data.quasiAttri.flagGDAS1 = strcmpi(thisMeteorAttri.dataSource, 'gdas1');
+    data.quasiAttri.meteorSource = thisMeteorAttri.dataSource;
+    data.quasiAttri.timestamp = thisMeteorAttri.datetime;
+
+    hIndOL = find(data.height >= PollyConfig.heightFullOverlap(flag532t), 1);
+    if ~ isempty(hIndOL)
+        att_beta_OC_532_qsi(1:hIndOL, :) = repmat(att_beta_OC_532_qsi(hIndOL, :), hIndOL, 1);
+    else
+        warning('Full overlap height is too large.');
+    end
+
+    [data.qsiBscOC532V1, ~] = quasiRetrieval(data.height, att_beta_OC_532_qsi, mExt532, mBsc532, PollyConfig.LR532, 'nIters', 6);
+end
+
+% quasi-retrieved backscatter at 1064 nm overlap corrected
+data.qsiBscOC1064V1 = NaN(length(data.height), length(data.mTime));
+att_beta_OC_1064_qsi = data.att_beta_OC_1064;
+if (sum(flag1064t) == 1)
+    att_beta_OC_1064_qsi(data.quality_mask_1064 ~= 0) = NaN;
+    att_beta_OC_1064_qsi = smooth2(att_beta_OC_1064_qsi, PollyConfig.quasi_smooth_h(flag1064t), PollyConfig.quasi_smooth_t(flag1064t));
+
+    % Rayleigh scattering
+%achtung
+    [mBsc1064, mExt1064] = rayleigh_scattering(1064, pressure, temperature + 273.15, 380, 70);
+    mBsc1064 = transpose(interp2(TimeMg, HeightMg, mBsc1064, mTimeg, Heightg, 'linear'));
+    mExt1064 = transpose(interp2(TimeMg, HeightMg, mExt1064, mTimeg, Heightg, 'linear'));
+    data.quasiAttri.flagGDAS1 = strcmpi(thisMeteorAttri.dataSource, 'gdas1');
+    data.quasiAttri.meteorSource = thisMeteorAttri.dataSource;
+    data.quasiAttri.timestamp = thisMeteorAttri.datetime;
+
+    hIndOL = find(data.height >= PollyConfig.heightFullOverlap(flag1064t), 1);
+    if ~ isempty(hIndOL)
+        att_beta_OC_1064_qsi(1:hIndOL, :) = repmat(att_beta_OC_1064_qsi(hIndOL, :), hIndOL, 1);
+    else
+        warning('Full overlap height is too large.');
+    end
+
+    [data.qsiBscOC1064V1, ~] = quasiRetrieval(data.height, att_beta_OC_1064_qsi, mExt1064, mBsc1064, PollyConfig.LR1064, 'nIters', 6);
+end
+
+% quasi-retrieved backscatter at 355 nm near range
+data.qsiBscNR355V1 = NaN(length(data.height), length(data.mTime));
+att_beta_NR_355_qsi = data.att_beta_NR_355;
+if (sum(flag355NR) == 1)
+    att_beta_NR_355_qsi(data.quality_mask_NR_355 ~= 0) = NaN;
+    att_beta_NR_355_qsi = smooth2(att_beta_NR_355_qsi, PollyConfig.quasi_smooth_h(flag355NR), PollyConfig.quasi_smooth_t(flag355NR));
+
+    % Rayleigh scattering
+%---------------achtung
+    [mBsc355, mExt355] = rayleigh_scattering(355, pressure, temperature + 273.15, 380, 70);
+    mBsc355 = transpose(interp2(TimeMg, HeightMg, mBsc355, mTimeg, Heightg, 'linear'));
+    mExt355 = transpose(interp2(TimeMg, HeightMg, mExt355, mTimeg, Heightg, 'linear'));
+    data.quasiAttri.flagGDAS1 = strcmpi(thisMeteorAttri.dataSource, 'gdas1');
+    data.quasiAttri.meteorSource = thisMeteorAttri.dataSource;
+    data.quasiAttri.timestamp = thisMeteorAttri.datetime;
+
+    hIndOL = find(data.height >= PollyConfig.heightFullOverlap(flag355NR), 1);
+    if ~ isempty(hIndOL)
+        att_beta_NR_355_qsi(1:hIndOL, :) = repmat(att_beta_NR_355_qsi(hIndOL, :), hIndOL, 1);
+    else
+        warning('Full overlap height is too large.');
+    end
+
+    [data.qsiBscNR355V1, ~] = quasiRetrieval(data.height, att_beta_NR_355_qsi, mExt355, mBsc355, PollyConfig.LR355, 'nIters', 6);
+end
+
+% quasi-retrieved backscatter at 532 nm near range
+data.qsiBscNR532V1 = NaN(length(data.height), length(data.mTime));
+att_beta_NR_532_qsi = data.att_beta_NR_532;
+if (sum(flag532NR) == 1)
+    att_beta_NR_532_qsi(data.quality_mask_NR_532 ~= 0) = NaN;
+    att_beta_NR_532_qsi = smooth2(att_beta_NR_532_qsi, PollyConfig.quasi_smooth_h(flag532NR), PollyConfig.quasi_smooth_t(flag532NR));
+
+    % Rayleigh scattering
+    [mBsc532, mExt532] = rayleigh_scattering(532, pressure, temperature + 273.15, 380, 70);
+  %achtung
+    mBsc532 = transpose(interp2(TimeMg, HeightMg, mBsc532, mTimeg, Heightg, 'linear'));
+    mExt532 = transpose(interp2(TimeMg, HeightMg, mExt532, mTimeg, Heightg, 'linear'));
+    data.quasiAttri.flagGDAS1 = strcmpi(thisMeteorAttri.dataSource, 'gdas1');
+    data.quasiAttri.meteorSource = thisMeteorAttri.dataSource;
+    data.quasiAttri.timestamp = thisMeteorAttri.datetime;
+
+    hIndOL = find(data.height >= PollyConfig.heightFullOverlap(flag532NR), 1);
+    if ~ isempty(hIndOL)
+        att_beta_NR_532_qsi(1:hIndOL, :) = repmat(att_beta_NR_532_qsi(hIndOL, :), hIndOL, 1);
+    else
+        warning('Full overlap height is too large.');
+    end
+
+    [data.qsiBscNR532V1, ~] = quasiRetrieval(data.height, att_beta_NR_532_qsi, mExt532, mBsc532, PollyConfig.LR532, 'nIters', 6);
+end
+
 % quasi-retrieved particle depolarization ratio at 532 nm
 if flagGHK
     data.qsiPDR532V1 = NaN(length(data.height), length(data.mTime));
@@ -4618,6 +4792,127 @@
     [data.qsiBsc1064V2, ~] = quasiRetrieval2(data.height, att_beta_1064_qsi, att_beta_607_qsi, 1064, mExt1064, mBsc1064, mExt607, 0.5, PollyConfig.LR1064, 'nIters', 3);
     data.qsiBsc1064V2 = smooth2(data.qsiBsc1064V2, PollyConfig.quasi_smooth_h(flag1064t), PollyConfig.quasi_smooth_t(flag1064t));
 end
+
+% quasi-retrieved backscatter at 355 nm (V2) overlap corrected
+data.qsiBscOC355V2 = NaN(length(data.height), length(data.mTime));
+att_beta_OC_355_qsi = data.att_beta_OC_355;
+att_beta_OC_387_qsi = att_beta_OC_387;
+if (sum(flag355t) == 1) && (sum(flag387FR) == 1)
+    att_beta_OC_355_qsi(data.quality_mask_355 ~= 0) = NaN;
+    att_beta_OC_387_qsi(data.quality_mask_387 ~= 0) = NaN;
+    att_beta_OC_355_qsi = smooth2(att_beta_OC_355_qsi, PollyConfig.quasi_smooth_h(flag355t), PollyConfig.quasi_smooth_t(flag355t));
+    att_beta_OC_387_qsi = smooth2(att_beta_OC_387_qsi, PollyConfig.quasi_smooth_h(flag387FR), PollyConfig.quasi_smooth_t(flag387FR));
+
+    % Rayleigh scattering
+    [mBsc355, mExt355] = rayleigh_scattering(355, pressure, temperature + 273.15, 380, 70);
+    [~, mExt387] = rayleigh_scattering(387, pressure, temperature + 273.15, 380, 70);
+    mBsc355 = transpose(interp2(TimeMg, HeightMg, mBsc355, mTimeg, Heightg, 'linear'));
+    mExt355 = transpose(interp2(TimeMg, HeightMg, mExt355, mTimeg, Heightg, 'linear'));
+    mExt387 = transpose(interp2(TimeMg, HeightMg, mExt387, mTimeg, Heightg, 'linear'));
+    data.quasiAttri.flagGDAS1 = strcmpi(thisMeteorAttri.dataSource, 'gdas1');
+    data.quasiAttri.meteorSource = thisMeteorAttri.dataSource;
+    data.quasiAttri.timestamp = thisMeteorAttri.datetime;
+
+    [data.qsiBscOC355V2, ~] = quasiRetrieval2(data.height, att_beta_OC_355_qsi, att_beta_OC_387_qsi, 355, mExt355, mBsc355, mExt387, 0.5, PollyConfig.LR355, 'nIters', 3);
+    data.qsiBscOC355V2 = smooth2(data.qsiBscOC355V2, PollyConfig.quasi_smooth_h(flag355t), PollyConfig.quasi_smooth_t(flag355t));
+end
+
+% quasi-retrieved backscatter at 532 nm (V2) overlap corrected
+data.qsiBscOC532V2 = NaN(length(data.height), length(data.mTime));
+att_beta_OC_532_qsi = data.att_beta_OC_532;
+att_beta_OC_607_qsi = att_beta_OC_607;
+if (sum(flag532t) == 1) && (sum(flag607FR) == 1)
+    att_beta_OC_532_qsi(data.quality_mask_532 ~= 0) = NaN;
+    att_beta_OC_607_qsi(data.quality_mask_607 ~= 0) = NaN;
+    att_beta_OC_532_qsi = smooth2(att_beta_OC_532_qsi, PollyConfig.quasi_smooth_h(flag532t), PollyConfig.quasi_smooth_t(flag532t));
+    att_beta_OC_607_qsi = smooth2(att_beta_OC_607_qsi, PollyConfig.quasi_smooth_h(flag607FR), PollyConfig.quasi_smooth_t(flag607FR));
+
+    % Rayleigh scattering
+    [mBsc532, mExt532] = rayleigh_scattering(532, pressure, temperature + 273.15, 380, 70);
+    [~, mExt607] = rayleigh_scattering(607, pressure, temperature + 273.15, 380, 70);
+    mBsc532 = transpose(interp2(TimeMg, HeightMg, mBsc532, mTimeg, Heightg, 'linear'));
+    mExt532 = transpose(interp2(TimeMg, HeightMg, mExt532, mTimeg, Heightg, 'linear'));
+    mExt607 = transpose(interp2(TimeMg, HeightMg, mExt607, mTimeg, Heightg, 'linear'));
+    data.quasiAttri.flagGDAS1 = strcmpi(thisMeteorAttri.dataSource, 'gdas1');
+    data.quasiAttri.meteorSource = thisMeteorAttri.dataSource;
+    data.quasiAttri.timestamp = thisMeteorAttri.datetime;
+
+    [data.qsiBscOC532V2, ~] = quasiRetrieval2(data.height, att_beta_OC_532_qsi, att_beta_OC_607_qsi, 532, mExt532, mBsc532, mExt607, 0.5, PollyConfig.LR532, 'nIters', 3);
+    data.qsiBscOC532V2 = smooth2(data.qsiBscOC532V2, PollyConfig.quasi_smooth_h(flag532t), PollyConfig.quasi_smooth_t(flag532t));
+end
+
+% quasi-retrieved backscatter at 1064 nm (V2) overlap corrected
+data.qsiBscOC1064V2 = NaN(length(data.height), length(data.mTime));
+att_beta_OC_1064_qsi = data.att_beta_OC_1064;
+att_beta_OC_607_qsi = att_beta_OC_607;
+if (sum(flag1064t) == 1) && (sum(flag607FR) == 1)
+    att_beta_OC_1064_qsi(data.quality_mask_1064 ~= 0) = NaN;
+    att_beta_OC_607_qsi(data.quality_mask_607 ~= 0) = NaN;
+    att_beta_OC_1064_qsi = smooth2(att_beta_OC_1064_qsi, PollyConfig.quasi_smooth_h(flag1064t), PollyConfig.quasi_smooth_t(flag1064t));
+    att_beta_OC_607_qsi = smooth2(att_beta_OC_607_qsi, PollyConfig.quasi_smooth_h(flag607FR), PollyConfig.quasi_smooth_t(flag607FR));
+
+    % Rayleigh scattering
+    [mBsc1064, mExt1064] = rayleigh_scattering(1064, pressure, temperature + 273.15, 380, 70);
+    [~, mExt607] = rayleigh_scattering(607, pressure, temperature + 273.15, 380, 70);
+    mBsc1064 = transpose(interp2(TimeMg, HeightMg, mBsc1064, mTimeg, Heightg, 'linear'));
+    mExt1064 = transpose(interp2(TimeMg, HeightMg, mExt1064, mTimeg, Heightg, 'linear'));
+    mExt607 = transpose(interp2(TimeMg, HeightMg, mExt607, mTimeg, Heightg, 'linear'));
+    data.quasiAttri.flagGDAS1 = strcmpi(thisMeteorAttri.dataSource, 'gdas1');
+    data.quasiAttri.meteorSource = thisMeteorAttri.dataSource;
+    data.quasiAttri.timestamp = thisMeteorAttri.datetime;
+
+    [data.qsiBscOC1064V2, ~] = quasiRetrieval2(data.height, att_beta_OC_1064_qsi, att_beta_OC_607_qsi, 1064, mExt1064, mBsc1064, mExt607, 0.5, PollyConfig.LR1064, 'nIters', 3);
+    data.qsiBscOC1064V2 = smooth2(data.qsiBscOC1064V2, PollyConfig.quasi_smooth_h(flag1064t), PollyConfig.quasi_smooth_t(flag1064t));
+end
+
+% quasi-retrieved backscatter at 355 nm (V2) near range
+data.qsiBscNR355V2 = NaN(length(data.height), length(data.mTime));
+att_beta_NR_355_qsi = data.att_beta_NR_355;
+att_beta_NR_387_qsi = att_beta_NR_387;
+if (sum(flag355t) == 1) && (sum(flag387FR) == 1)
+    att_beta_NR_355_qsi(data.quality_mask_NR_355 ~= 0) = NaN;
+    %att_beta_NR_387_qsi(data.quality_mask_NR_387 ~= 0) = NaN;
+    att_beta_NR_355_qsi = smooth2(att_beta_NR_355_qsi, PollyConfig.quasi_smooth_h(flag355t), PollyConfig.quasi_smooth_t(flag355t));
+    att_beta_NR_387_qsi = smooth2(att_beta_NR_387_qsi, PollyConfig.quasi_smooth_h(flag387FR), PollyConfig.quasi_smooth_t(flag387FR));
+
+    % Rayleigh scattering
+    [mBsc355, mExt355] = rayleigh_scattering(355, pressure, temperature + 273.15, 380, 70);
+    [~, mExt387] = rayleigh_scattering(387, pressure, temperature + 273.15, 380, 70);
+    mBsc355 = transpose(interp2(TimeMg, HeightMg, mBsc355, mTimeg, Heightg, 'linear'));
+    mExt355 = transpose(interp2(TimeMg, HeightMg, mExt355, mTimeg, Heightg, 'linear'));
+    mExt387 = transpose(interp2(TimeMg, HeightMg, mExt387, mTimeg, Heightg, 'linear'));
+    data.quasiAttri.flagGDAS1 = strcmpi(thisMeteorAttri.dataSource, 'gdas1');
+    data.quasiAttri.meteorSource = thisMeteorAttri.dataSource;
+    data.quasiAttri.timestamp = thisMeteorAttri.datetime;
+
+    [data.qsiBscNR355V2, ~] = quasiRetrieval2(data.height, att_beta_NR_355_qsi, att_beta_NR_387_qsi, 355, mExt355, mBsc355, mExt387, 0.5, PollyConfig.LR355, 'nIters', 3);
+    data.qsiBscNR355V2 = smooth2(data.qsiBscNR355V2, PollyConfig.quasi_smooth_h(flag355t), PollyConfig.quasi_smooth_t(flag355t));
+end
+
+% quasi-retrieved backscatter at 532 nm (V2) near range
+data.qsiBscNR532V2 = NaN(length(data.height), length(data.mTime));
+att_beta_NR_532_qsi = data.att_beta_NR_532;
+att_beta_NR_607_qsi = att_beta_NR_607;
+if (sum(flag532t) == 1) && (sum(flag607FR) == 1)
+    att_beta_NR_532_qsi(data.quality_mask_NR_532 ~= 0) = NaN;
+    %att_beta_NR_607_qsi(data.quality_mask_NR_607 ~= 0) = NaN;
+    att_beta_NR_532_qsi = smooth2(att_beta_NR_532_qsi, PollyConfig.quasi_smooth_h(flag532t), PollyConfig.quasi_smooth_t(flag532t));
+    att_beta_NR_607_qsi = smooth2(att_beta_NR_607_qsi, PollyConfig.quasi_smooth_h(flag607FR), PollyConfig.quasi_smooth_t(flag607FR));
+
+    % Rayleigh scattering
+    [mBsc532, mExt532] = rayleigh_scattering(532, pressure, temperature + 273.15, 380, 70);
+    [~, mExt607] = rayleigh_scattering(607, pressure, temperature + 273.15, 380, 70);
+    mBsc532 = transpose(interp2(TimeMg, HeightMg, mBsc532, mTimeg, Heightg, 'linear'));
+    mExt532 = transpose(interp2(TimeMg, HeightMg, mExt532, mTimeg, Heightg, 'linear'));
+    mExt607 = transpose(interp2(TimeMg, HeightMg, mExt607, mTimeg, Heightg, 'linear'));
+    data.quasiAttri.flagGDAS1 = strcmpi(thisMeteorAttri.dataSource, 'gdas1');
+    data.quasiAttri.meteorSource = thisMeteorAttri.dataSource;
+    data.quasiAttri.timestamp = thisMeteorAttri.datetime;
+
+    [data.qsiBscNR532V2, ~] = quasiRetrieval2(data.height, att_beta_NR_532_qsi, att_beta_NR_607_qsi, 532, mExt532, mBsc532, mExt607, 0.5, PollyConfig.LR532, 'nIters', 3);
+    data.qsiBscNR532V2 = smooth2(data.qsiBscNR532V2, PollyConfig.quasi_smooth_h(flag532t), PollyConfig.quasi_smooth_t(flag532t));
+end
+
 clearvars att_beta_1064_qsi att_beta_355_qsi att_beta_387_qsi att_beta_532_qsi att_beta_607_qsi att_beta_387 att_beta_607;
 clearvars  mBsc355 mExt355 mBsc387 mExt387 mBsc407 mExt407 mBsc532 mExt532 mBsc607 mExt607 mBsc1058 mExt1058 mBsc1064 mExt1064 number_density ;
 % quasi-retrieved particle depolarization ratio at 532 nm (V2)
diff --git a/lib/io/pollySaveQsiV1.m b/lib/io/pollySaveQsiV1.m
index c47b6610..ed24956f 100644
--- a/lib/io/pollySaveQsiV1.m
+++ b/lib/io/pollySaveQsiV1.m
@@ -34,8 +34,14 @@ function pollySaveQsiV1(data)
 varID_latitude = netcdf.defVar(ncID, 'latitude', 'NC_FLOAT', dimID_constant);
 varID_time = netcdf.defVar(ncID, 'time', 'NC_DOUBLE', dimID_time);
 varID_height = netcdf.defVar(ncID, 'height', 'NC_FLOAT', dimID_height);
+varID_quasi_bsc_355 = netcdf.defVar(ncID, 'quasi_bsc_355', 'NC_FLOAT', [dimID_height, dimID_time]);
 varID_quasi_bsc_532 = netcdf.defVar(ncID, 'quasi_bsc_532', 'NC_FLOAT', [dimID_height, dimID_time]);
 varID_quasi_bsc_1064 = netcdf.defVar(ncID, 'quasi_bsc_1064', 'NC_FLOAT', [dimID_height, dimID_time]);
+varID_quasi_bsc_OC_355 = netcdf.defVar(ncID, 'quasi_bsc_OC_355', 'NC_FLOAT', [dimID_height, dimID_time]);
+varID_quasi_bsc_OC_532 = netcdf.defVar(ncID, 'quasi_bsc_OC_532', 'NC_FLOAT', [dimID_height, dimID_time]);
+varID_quasi_bsc_OC_1064 = netcdf.defVar(ncID, 'quasi_bsc_OC_1064', 'NC_FLOAT', [dimID_height, dimID_time]);
+varID_quasi_bsc_NR_355 = netcdf.defVar(ncID, 'quasi_bsc_NR_355', 'NC_FLOAT', [dimID_height, dimID_time]);
+varID_quasi_bsc_NR_532 = netcdf.defVar(ncID, 'quasi_bsc_NR_532', 'NC_FLOAT', [dimID_height, dimID_time]);
 varID_quasi_pardepol_532 = netcdf.defVar(ncID, 'quasi_pardepol_532', 'NC_FLOAT', [dimID_height, dimID_time]);
 varID_quasi_ang_532_1064 = netcdf.defVar(ncID, 'quasi_ang_532_1064', 'NC_FLOAT', [dimID_height, dimID_time]);
 varID_quality_mask_532 = netcdf.defVar(ncID, 'quality_mask_532', 'NC_BYTE', [dimID_height, dimID_time]);
@@ -43,8 +49,14 @@ function pollySaveQsiV1(data)
 varID_quality_mask_voldepol_532 = netcdf.defVar(ncID, 'quality_mask_voldepol_532', 'NC_BYTE', [dimID_height, dimID_time]);
 
 % define the filling value
+netcdf.defVarFill(ncID, varID_quasi_bsc_355, false, -999);
 netcdf.defVarFill(ncID, varID_quasi_bsc_532, false, -999);
 netcdf.defVarFill(ncID, varID_quasi_bsc_1064, false, -999);
+netcdf.defVarFill(ncID, varID_quasi_bsc_OC_355, false, -999);
+netcdf.defVarFill(ncID, varID_quasi_bsc_OC_532, false, -999);
+netcdf.defVarFill(ncID, varID_quasi_bsc_OC_1064, false, -999);
+netcdf.defVarFill(ncID, varID_quasi_bsc_NR_355, false, -999);
+netcdf.defVarFill(ncID, varID_quasi_bsc_NR_532, false, -999);
 netcdf.defVarFill(ncID, varID_quasi_pardepol_532, false, -999);
 netcdf.defVarFill(ncID, varID_quasi_ang_532_1064, false, -999);
 netcdf.defVarFill(ncID, varID_quality_mask_532, false, 1);
@@ -52,8 +64,14 @@ function pollySaveQsiV1(data)
 netcdf.defVarFill(ncID, varID_quality_mask_voldepol_532, false, 1);
 
 % define the data compression
+netcdf.defVarDeflate(ncID, varID_quasi_bsc_355, true, true, 5);
 netcdf.defVarDeflate(ncID, varID_quasi_bsc_532, true, true, 5);
 netcdf.defVarDeflate(ncID, varID_quasi_bsc_1064, true, true, 5);
+netcdf.defVarDeflate(ncID, varID_quasi_bsc_OC_355, true, true, 5);
+netcdf.defVarDeflate(ncID, varID_quasi_bsc_OC_532, true, true, 5);
+netcdf.defVarDeflate(ncID, varID_quasi_bsc_OC_1064, true, true, 5);
+netcdf.defVarDeflate(ncID, varID_quasi_bsc_NR_355, true, true, 5);
+netcdf.defVarDeflate(ncID, varID_quasi_bsc_NR_532, true, true, 5);
 netcdf.defVarDeflate(ncID, varID_quasi_pardepol_532, true, true, 5);
 netcdf.defVarDeflate(ncID, varID_quasi_ang_532_1064, true, true, 5);
 netcdf.defVarDeflate(ncID, varID_quality_mask_532, true, true, 5);
@@ -69,8 +87,14 @@ function pollySaveQsiV1(data)
 netcdf.putVar(ncID, varID_latitude, single(data.lat));
 netcdf.putVar(ncID, varID_time, datenum_2_unix_timestamp(data.mTime));   % do the conversion
 netcdf.putVar(ncID, varID_height, single(data.height));
+netcdf.putVar(ncID, varID_quasi_bsc_355, single(data.qsiBsc355V1));
 netcdf.putVar(ncID, varID_quasi_bsc_532, single(data.qsiBsc532V1));
 netcdf.putVar(ncID, varID_quasi_bsc_1064, single(data.qsiBsc1064V1));
+netcdf.putVar(ncID, varID_quasi_bsc_OC_355, single(data.qsiBscOC355V1));
+netcdf.putVar(ncID, varID_quasi_bsc_OC_532, single(data.qsiBscOC532V1));
+netcdf.putVar(ncID, varID_quasi_bsc_OC_1064, single(data.qsiBscOC1064V1));
+netcdf.putVar(ncID, varID_quasi_bsc_NR_355, single(data.qsiBscNR355V1));
+netcdf.putVar(ncID, varID_quasi_bsc_NR_532, single(data.qsiBscNR532V1));
 netcdf.putVar(ncID, varID_quasi_pardepol_532, single(data.qsiPDR532V1));
 netcdf.putVar(ncID, varID_quasi_ang_532_1064, single(data.qsiAE_532_1064_V1));
 netcdf.putVar(ncID, varID_quality_mask_532, int8(data.quality_mask_532));
@@ -112,12 +136,27 @@ function pollySaveQsiV1(data)
 netcdf.putAtt(ncID, varID_height, 'standard_name', 'height');
 netcdf.putAtt(ncID, varID_height, 'axis', 'Z');
 
+% quasi_bsc_355
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'unit', 'sr^-1 m^-1');
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'unit_html', 'sr<sup>-1</sup> m<sup>-1</sup>');
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'long_name', 'quasi aerosol backscatter coefficients at 355 nm');
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'standard_name', 'quasi_bsc_355');
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'plot_range', PollyConfig.zLim_quasi_beta_355/1e6);
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'plot_scale', 'linear');
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'source', CampaignConfig.name);
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'error_variable', 'quasi_beta_355_error');
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'bias_variable', 'quasi_beta_355_bias');
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'Lidar_calibration_constant_used', data.LCUsed.LCUsed355);
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'Lidar_ratio_used', PollyConfig.LR355);
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'retrieved_info', sprintf('Fixed Lidar ratio: %5.1f[Sr]', PollyConfig.LR355));
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'comment', 'This parameter is retrieved by the method demonstrated in (Holger, ATM, 2017). The retrieved results are dependent on the lidar constants and the AOD below the current bin. If the AOD is greater than 0.2, the relative uncertainty can be as large as 20%. Be careful about that!');
+
 % quasi_bsc_532
 netcdf.putAtt(ncID, varID_quasi_bsc_532, 'unit', 'sr^-1 m^-1');
 netcdf.putAtt(ncID, varID_quasi_bsc_532, 'unit_html', 'sr<sup>-1</sup> m<sup>-1</sup>');
 netcdf.putAtt(ncID, varID_quasi_bsc_532, 'long_name', 'quasi aerosol backscatter coefficients at 532 nm');
 netcdf.putAtt(ncID, varID_quasi_bsc_532, 'standard_name', 'quasi_bsc_532');
-netcdf.putAtt(ncID, varID_quasi_bsc_532, 'plot_range', PollyConfig.zLim_quasi_beta_355/1e6);
+netcdf.putAtt(ncID, varID_quasi_bsc_532, 'plot_range', PollyConfig.zLim_quasi_beta_532/1e6);
 netcdf.putAtt(ncID, varID_quasi_bsc_532, 'plot_scale', 'linear');
 netcdf.putAtt(ncID, varID_quasi_bsc_532, 'source', CampaignConfig.name);
 netcdf.putAtt(ncID, varID_quasi_bsc_532, 'error_variable', 'quasi_beta_532_error');
@@ -142,6 +181,81 @@ function pollySaveQsiV1(data)
 netcdf.putAtt(ncID, varID_quasi_bsc_1064, 'retrieved_info', sprintf('Fixed Lidar ratio: %5.1f[Sr]', PollyConfig.LR1064));
 netcdf.putAtt(ncID, varID_quasi_bsc_1064, 'comment', 'This parameter is retrieved by the method demonstrated in (Holger, ATM, 2017). The retrieved results are dependent on the lidar constants and the AOD below the current bin. If the AOD is greater than 0.2, the relative uncertainty can be as large as 20%. Be careful about that!');
 
+% quasi_bsc_OC_355
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'unit', 'sr^-1 m^-1');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'unit_html', 'sr<sup>-1</sup> m<sup>-1</sup>');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'long_name', 'quasi aerosol backscatter coefficients overlap corrected at 355 nm');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'standard_name', 'quasi_bsc_OC_355');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'plot_range', PollyConfig.zLim_quasi_beta_355/1e6);
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'plot_scale', 'linear');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'source', CampaignConfig.name);
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'error_variable', 'quasi_beta_355_error');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'bias_variable', 'quasi_beta_355_bias');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'Lidar_calibration_constant_used', data.LCUsed.LCUsed355);
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'Lidar_ratio_used', PollyConfig.LR355);
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'retrieved_info', sprintf('Fixed Lidar ratio: %5.1f[Sr]', PollyConfig.LR355));
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'comment', 'This parameter is retrieved by the method demonstrated in (Holger, ATM, 2017). The retrieved results are dependent on the lidar constants and the AOD below the current bin. If the AOD is greater than 0.2, the relative uncertainty can be as large as 20%. Be careful about that!');
+
+% quasi_bsc_OC_532
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'unit', 'sr^-1 m^-1');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'unit_html', 'sr<sup>-1</sup> m<sup>-1</sup>');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'long_name', 'quasi aerosol backscatter coefficients overlap corrected at 532 nm');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'standard_name', 'quasi_bsc_OC_532');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'plot_range', PollyConfig.zLim_quasi_beta_532/1e6);
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'plot_scale', 'linear');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'source', CampaignConfig.name);
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'error_variable', 'quasi_beta_532_error');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'bias_variable', 'quasi_beta_532_bias');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'Lidar_calibration_constant_used', data.LCUsed.LCUsed532);
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'Lidar_ratio_used', PollyConfig.LR532);
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'retrieved_info', sprintf('Fixed Lidar ratio: %5.1f[Sr]', PollyConfig.LR532));
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'comment', 'This parameter is retrieved by the method demonstrated in (Holger, ATM, 2017). The retrieved results are dependent on the lidar constants and the AOD below the current bin. If the AOD is greater than 0.2, the relative uncertainty can be as large as 20%. Be careful about that!');
+
+% quasi_bsc_OC_1064
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'unit', 'sr^-1 m^-1');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'unit_html', 'sr<sup>-1</sup> m<sup>-1</sup>');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'long_name', 'quasi aerosol backscatter coefficients overlap corrected at 1064 nm');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'standard_name', 'quasi_bsc_OC_1064');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'plot_range', PollyConfig.zLim_quasi_beta_1064/1e6);
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'plot_scale', 'linear');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'source', CampaignConfig.name);
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'error_variable', 'quasi_beta_1064_error');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'bias_variable', 'quasi_beta_1064_bias');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'Lidar_calibration_constant_used', data.LCUsed.LCUsed1064);
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'Lidar_ratio_used', PollyConfig.LR1064);
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'retrieved_info', sprintf('Fixed Lidar ratio: %5.1f[Sr]', PollyConfig.LR1064));
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'comment', 'This parameter is retrieved by the method demonstrated in (Holger, ATM, 2017). The retrieved results are dependent on the lidar constants and the AOD below the current bin. If the AOD is greater than 0.2, the relative uncertainty can be as large as 20%. Be careful about that!');
+
+% quasi_bsc_NR_355
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'unit', 'sr^-1 m^-1');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'unit_html', 'sr<sup>-1</sup> m<sup>-1</sup>');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'long_name', 'quasi aerosol backscatter coefficients near range at 355 nm');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'standard_name', 'quasi_bsc_NR_355');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'plot_range', PollyConfig.zLim_quasi_beta_355/1e6);
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'plot_scale', 'linear');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'source', CampaignConfig.name);
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'error_variable', 'quasi_beta_355_error');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'bias_variable', 'quasi_beta_355_bias');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'Lidar_calibration_constant_used', data.LCUsed.LCUsed355NR);
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'Lidar_ratio_used', PollyConfig.LR355);
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'retrieved_info', sprintf('Fixed Lidar ratio: %5.1f[Sr]', PollyConfig.LR355));
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'comment', 'This parameter is retrieved by the method demonstrated in (Holger, ATM, 2017). The retrieved results are dependent on the lidar constants and the AOD below the current bin. If the AOD is greater than 0.2, the relative uncertainty can be as large as 20%. Be careful about that!');
+
+% quasi_bsc_NR_532
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'unit', 'sr^-1 m^-1');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'unit_html', 'sr<sup>-1</sup> m<sup>-1</sup>');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'long_name', 'quasi aerosol backscatter coefficients near range at 532 nm');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'standard_name', 'quasi_bsc_NR_532');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'plot_range', PollyConfig.zLim_quasi_beta_532/1e6);
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'plot_scale', 'linear');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'source', CampaignConfig.name);
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'error_variable', 'quasi_beta_532_error');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'bias_variable', 'quasi_beta_532_bias');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'Lidar_calibration_constant_used', data.LCUsed.LCUsed532NR);
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'Lidar_ratio_used', PollyConfig.LR532);
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'retrieved_info', sprintf('Fixed Lidar ratio: %5.1f[Sr]', PollyConfig.LR532));
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'comment', 'This parameter is retrieved by the method demonstrated in (Holger, ATM, 2017). The retrieved results are dependent on the lidar constants and the AOD below the current bin. If the AOD is greater than 0.2, the relative uncertainty can be as large as 20%. Be careful about that!');
+
 % quasi_pardepol_532
 netcdf.putAtt(ncID, varID_quasi_pardepol_532, 'unit', '');
 netcdf.putAtt(ncID, varID_quasi_pardepol_532, 'long_name', 'quasi particle depolarization ratio at 532 nm');
diff --git a/lib/io/pollySaveQsiV2.m b/lib/io/pollySaveQsiV2.m
index 51b04924..63aa649f 100644
--- a/lib/io/pollySaveQsiV2.m
+++ b/lib/io/pollySaveQsiV2.m
@@ -33,8 +33,14 @@ function pollySaveQsiV2(data)
 varID_latitude = netcdf.defVar(ncID, 'latitude', 'NC_FLOAT', dimID_constant);
 varID_time = netcdf.defVar(ncID, 'time', 'NC_DOUBLE', dimID_time);
 varID_height = netcdf.defVar(ncID, 'height', 'NC_FLOAT', dimID_height);
+varID_quasi_bsc_355 = netcdf.defVar(ncID, 'quasi_bsc_355', 'NC_FLOAT', [dimID_height, dimID_time]);
 varID_quasi_bsc_532 = netcdf.defVar(ncID, 'quasi_bsc_532', 'NC_FLOAT', [dimID_height, dimID_time]);
 varID_quasi_bsc_1064 = netcdf.defVar(ncID, 'quasi_bsc_1064', 'NC_FLOAT', [dimID_height, dimID_time]);
+varID_quasi_bsc_OC_355 = netcdf.defVar(ncID, 'quasi_bsc_OC_355', 'NC_FLOAT', [dimID_height, dimID_time]);
+varID_quasi_bsc_OC_532 = netcdf.defVar(ncID, 'quasi_bsc_OC_532', 'NC_FLOAT', [dimID_height, dimID_time]);
+varID_quasi_bsc_OC_1064 = netcdf.defVar(ncID, 'quasi_bsc_OC_1064', 'NC_FLOAT', [dimID_height, dimID_time]);
+varID_quasi_bsc_NR_355 = netcdf.defVar(ncID, 'quasi_bsc_NR_355', 'NC_FLOAT', [dimID_height, dimID_time]);
+varID_quasi_bsc_NR_532 = netcdf.defVar(ncID, 'quasi_bsc_NR_532', 'NC_FLOAT', [dimID_height, dimID_time]);
 varID_quasi_pardepol_532 = netcdf.defVar(ncID, 'quasi_pardepol_532', 'NC_FLOAT', [dimID_height, dimID_time]);
 varID_quasi_ang_532_1064 = netcdf.defVar(ncID, 'quasi_ang_532_1064', 'NC_FLOAT', [dimID_height, dimID_time]);
 varID_quality_mask_532 = netcdf.defVar(ncID, 'quality_mask_532', 'NC_BYTE', [dimID_height, dimID_time]);
@@ -42,8 +48,14 @@ function pollySaveQsiV2(data)
 varID_quality_mask_voldepol_532 = netcdf.defVar(ncID, 'quality_mask_voldepol_532', 'NC_BYTE', [dimID_height, dimID_time]);
 
 % define the filling value
+netcdf.defVarFill(ncID, varID_quasi_bsc_355, false, -999);
 netcdf.defVarFill(ncID, varID_quasi_bsc_532, false, -999);
 netcdf.defVarFill(ncID, varID_quasi_bsc_1064, false, -999);
+netcdf.defVarFill(ncID, varID_quasi_bsc_OC_355, false, -999);
+netcdf.defVarFill(ncID, varID_quasi_bsc_OC_532, false, -999);
+netcdf.defVarFill(ncID, varID_quasi_bsc_OC_1064, false, -999);
+netcdf.defVarFill(ncID, varID_quasi_bsc_NR_355, false, -999);
+netcdf.defVarFill(ncID, varID_quasi_bsc_NR_532, false, -999);
 netcdf.defVarFill(ncID, varID_quasi_pardepol_532, false, -999);
 netcdf.defVarFill(ncID, varID_quasi_ang_532_1064, false, -999);
 netcdf.defVarFill(ncID, varID_quality_mask_532, false, 1);
@@ -51,8 +63,14 @@ function pollySaveQsiV2(data)
 netcdf.defVarFill(ncID, varID_quality_mask_voldepol_532, false, 1);
 
 % define the data compression
+netcdf.defVarDeflate(ncID, varID_quasi_bsc_355, true, true, 5);
 netcdf.defVarDeflate(ncID, varID_quasi_bsc_532, true, true, 5);
 netcdf.defVarDeflate(ncID, varID_quasi_bsc_1064, true, true, 5);
+netcdf.defVarDeflate(ncID, varID_quasi_bsc_OC_355, true, true, 5);
+netcdf.defVarDeflate(ncID, varID_quasi_bsc_OC_532, true, true, 5);
+netcdf.defVarDeflate(ncID, varID_quasi_bsc_OC_1064, true, true, 5);
+netcdf.defVarDeflate(ncID, varID_quasi_bsc_NR_355, true, true, 5);
+netcdf.defVarDeflate(ncID, varID_quasi_bsc_NR_532, true, true, 5);
 netcdf.defVarDeflate(ncID, varID_quasi_pardepol_532, true, true, 5);
 netcdf.defVarDeflate(ncID, varID_quasi_ang_532_1064, true, true, 5);
 netcdf.defVarDeflate(ncID, varID_quality_mask_532, true, true, 5);
@@ -68,8 +86,14 @@ function pollySaveQsiV2(data)
 netcdf.putVar(ncID, varID_latitude, single(data.lat));
 netcdf.putVar(ncID, varID_time, datenum_2_unix_timestamp(data.mTime));   % do the conversion
 netcdf.putVar(ncID, varID_height, single(data.height));
+netcdf.putVar(ncID, varID_quasi_bsc_355, single(data.qsiBsc355V2));
 netcdf.putVar(ncID, varID_quasi_bsc_532, single(data.qsiBsc532V2));
 netcdf.putVar(ncID, varID_quasi_bsc_1064, single(data.qsiBsc1064V2));
+netcdf.putVar(ncID, varID_quasi_bsc_OC_355, single(data.qsiBscOC355V2));
+netcdf.putVar(ncID, varID_quasi_bsc_OC_532, single(data.qsiBscOC532V2));
+netcdf.putVar(ncID, varID_quasi_bsc_OC_1064, single(data.qsiBscOC1064V2));
+netcdf.putVar(ncID, varID_quasi_bsc_NR_355, single(data.qsiBscNR355V2));
+netcdf.putVar(ncID, varID_quasi_bsc_NR_532, single(data.qsiBscNR532V2));
 netcdf.putVar(ncID, varID_quasi_pardepol_532, single(data.qsiPDR532V2));
 netcdf.putVar(ncID, varID_quasi_ang_532_1064, single(data.qsiAE_532_1064_V2));
 netcdf.putVar(ncID, varID_quality_mask_532, int8(data.quality_mask_532_V2));
@@ -111,6 +135,21 @@ function pollySaveQsiV2(data)
 netcdf.putAtt(ncID, varID_height, 'standard_name', 'height');
 netcdf.putAtt(ncID, varID_height, 'axis', 'Z');
 
+% quasi_bsc_355
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'unit', 'sr^-1 m^-1');
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'unit_html', 'sr<sup>-1</sup> m<sup>-1</sup>');
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'long_name', 'quasi aerosol backscatter coefficients at 355 nm');
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'standard_name', 'quasi_bsc_355');
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'plot_range', PollyConfig.zLim_quasi_beta_355/1e6);
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'plot_scale', 'linear');
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'source', CampaignConfig.name);
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'error_variable', 'quasi_beta_355_error');
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'bias_variable', 'quasi_beta_355_bias');
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'Lidar_calibration_constant_used', data.LCUsed.LCUsed355);
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'Lidar_ratio_used', PollyConfig.LR355);
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'retrieved_info', sprintf('Fixed Lidar ratio: %5.1f[Sr]', PollyConfig.LR355));
+netcdf.putAtt(ncID, varID_quasi_bsc_355, 'comment', 'This parameter is retrieved by the method demonstrated in (Holger, ATM, 2017). The retrieved results are dependent on the lidar constants and the AOD below the current bin. If the AOD is greater than 0.2, the relative uncertainty can be as large as 20%. Be careful about that!');
+
 % quasi_bsc_532
 netcdf.putAtt(ncID, varID_quasi_bsc_532, 'unit', 'sr^-1 m^-1');
 netcdf.putAtt(ncID, varID_quasi_bsc_532, 'unit_html', 'sr<sup>-1</sup> m<sup>-1</sup>');
@@ -121,6 +160,8 @@ function pollySaveQsiV2(data)
 netcdf.putAtt(ncID, varID_quasi_bsc_532, 'source', CampaignConfig.name);
 netcdf.putAtt(ncID, varID_quasi_bsc_532, 'error_variable', 'quasi_beta_532_error');
 netcdf.putAtt(ncID, varID_quasi_bsc_532, 'bias_variable', 'quasi_beta_532_bias');
+netcdf.putAtt(ncID, varID_quasi_bsc_532, 'Lidar_calibration_constant_used', data.LCUsed.LCUsed532);
+netcdf.putAtt(ncID, varID_quasi_bsc_532, 'Lidar_ratio_used', PollyConfig.LR532);
 netcdf.putAtt(ncID, varID_quasi_bsc_532, 'retrieved_info', sprintf('Fixed Lidar ratio: %5.1f[Sr]', PollyConfig.LR532));
 netcdf.putAtt(ncID, varID_quasi_bsc_532, 'comment', 'This parameter is retrieved by the method demonstrated in (Holger, ATM, 2017). The retrieved results are dependent on the lidar constants and the AOD below the current bin. If the AOD is greater than 0.2, the relative uncertainty can be as large as 20%. Be careful about that!');
 
@@ -134,9 +175,86 @@ function pollySaveQsiV2(data)
 netcdf.putAtt(ncID, varID_quasi_bsc_1064, 'source', CampaignConfig.name);
 netcdf.putAtt(ncID, varID_quasi_bsc_1064, 'error_variable', 'quasi_beta_1064_error');
 netcdf.putAtt(ncID, varID_quasi_bsc_1064, 'bias_variable', 'quasi_beta_1064_bias');
+netcdf.putAtt(ncID, varID_quasi_bsc_1064, 'Lidar_calibration_constant_used', data.LCUsed.LCUsed1064);
+netcdf.putAtt(ncID, varID_quasi_bsc_1064, 'Lidar_ratio_used', PollyConfig.LR1064);
 netcdf.putAtt(ncID, varID_quasi_bsc_1064, 'retrieved_info', sprintf('Fixed Lidar ratio: %5.1f[Sr]', PollyConfig.LR1064));
 netcdf.putAtt(ncID, varID_quasi_bsc_1064, 'comment', 'This parameter is retrieved by the method demonstrated in (Holger, ATM, 2017). The retrieved results are dependent on the lidar constants and the AOD below the current bin. If the AOD is greater than 0.2, the relative uncertainty can be as large as 20%. Be careful about that!');
 
+% quasi_bsc_OC_355
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'unit', 'sr^-1 m^-1');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'unit_html', 'sr<sup>-1</sup> m<sup>-1</sup>');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'long_name', 'quasi aerosol backscatter coefficients overlap corrected at 355 nm');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'standard_name', 'quasi_bsc_OC_355');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'plot_range', PollyConfig.zLim_quasi_beta_355/1e6);
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'plot_scale', 'linear');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'source', CampaignConfig.name);
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'error_variable', 'quasi_beta_355_error');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'bias_variable', 'quasi_beta_355_bias');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'Lidar_calibration_constant_used', data.LCUsed.LCUsed355);
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'Lidar_ratio_used', PollyConfig.LR355);
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'retrieved_info', sprintf('Fixed Lidar ratio: %5.1f[Sr]', PollyConfig.LR355));
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_355, 'comment', 'This parameter is retrieved by the method demonstrated in (Holger, ATM, 2017). The retrieved results are dependent on the lidar constants and the AOD below the current bin. If the AOD is greater than 0.2, the relative uncertainty can be as large as 20%. Be careful about that!');
+
+% quasi_bsc_OC_532
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'unit', 'sr^-1 m^-1');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'unit_html', 'sr<sup>-1</sup> m<sup>-1</sup>');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'long_name', 'quasi aerosol backscatter coefficients overlap corrected at 532 nm');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'standard_name', 'quasi_bsc_OC_532');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'plot_range', PollyConfig.zLim_quasi_beta_532/1e6);
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'plot_scale', 'linear');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'source', CampaignConfig.name);
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'error_variable', 'quasi_beta_532_error');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'bias_variable', 'quasi_beta_532_bias');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'Lidar_calibration_constant_used', data.LCUsed.LCUsed532);
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'Lidar_ratio_used', PollyConfig.LR532);
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'retrieved_info', sprintf('Fixed Lidar ratio: %5.1f[Sr]', PollyConfig.LR532));
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_532, 'comment', 'This parameter is retrieved by the method demonstrated in (Holger, ATM, 2017). The retrieved results are dependent on the lidar constants and the AOD below the current bin. If the AOD is greater than 0.2, the relative uncertainty can be as large as 20%. Be careful about that!');
+
+% quasi_bsc_OC_1064
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'unit', 'sr^-1 m^-1');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'unit_html', 'sr<sup>-1</sup> m<sup>-1</sup>');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'long_name', 'quasi aerosol backscatter coefficients overlap corrected at 1064 nm');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'standard_name', 'quasi_bsc_OC_1064');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'plot_range', PollyConfig.zLim_quasi_beta_1064/1e6);
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'plot_scale', 'linear');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'source', CampaignConfig.name);
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'error_variable', 'quasi_beta_1064_error');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'bias_variable', 'quasi_beta_1064_bias');
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'Lidar_calibration_constant_used', data.LCUsed.LCUsed1064);
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'Lidar_ratio_used', PollyConfig.LR1064);
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'retrieved_info', sprintf('Fixed Lidar ratio: %5.1f[Sr]', PollyConfig.LR1064));
+netcdf.putAtt(ncID, varID_quasi_bsc_OC_1064, 'comment', 'This parameter is retrieved by the method demonstrated in (Holger, ATM, 2017). The retrieved results are dependent on the lidar constants and the AOD below the current bin. If the AOD is greater than 0.2, the relative uncertainty can be as large as 20%. Be careful about that!');
+
+% quasi_bsc_NR_355
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'unit', 'sr^-1 m^-1');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'unit_html', 'sr<sup>-1</sup> m<sup>-1</sup>');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'long_name', 'quasi aerosol backscatter coefficients near range at 355 nm');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'standard_name', 'quasi_bsc_NR_355');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'plot_range', PollyConfig.zLim_quasi_beta_355/1e6);
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'plot_scale', 'linear');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'source', CampaignConfig.name);
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'error_variable', 'quasi_beta_355_error');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'bias_variable', 'quasi_beta_355_bias');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'Lidar_calibration_constant_used', data.LCUsed.LCUsed355NR);
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'Lidar_ratio_used', PollyConfig.LR355);
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'retrieved_info', sprintf('Fixed Lidar ratio: %5.1f[Sr]', PollyConfig.LR355));
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_355, 'comment', 'This parameter is retrieved by the method demonstrated in (Holger, ATM, 2017). The retrieved results are dependent on the lidar constants and the AOD below the current bin. If the AOD is greater than 0.2, the relative uncertainty can be as large as 20%. Be careful about that!');
+
+% quasi_bsc_NR_532
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'unit', 'sr^-1 m^-1');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'unit_html', 'sr<sup>-1</sup> m<sup>-1</sup>');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'long_name', 'quasi aerosol backscatter coefficients near range at 532 nm');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'standard_name', 'quasi_bsc_NR_532');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'plot_range', PollyConfig.zLim_quasi_beta_532/1e6);
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'plot_scale', 'linear');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'source', CampaignConfig.name);
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'error_variable', 'quasi_beta_532_error');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'bias_variable', 'quasi_beta_532_bias');
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'Lidar_calibration_constant_used', data.LCUsed.LCUsed532NR);
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'Lidar_ratio_used', PollyConfig.LR532);
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'retrieved_info', sprintf('Fixed Lidar ratio: %5.1f[Sr]', PollyConfig.LR532));
+netcdf.putAtt(ncID, varID_quasi_bsc_NR_532, 'comment', 'This parameter is retrieved by the method demonstrated in (Holger, ATM, 2017). The retrieved results are dependent on the lidar constants and the AOD below the current bin. If the AOD is greater than 0.2, the relative uncertainty can be as large as 20%. Be careful about that!');
+
 % quasi_pardepol_532
 netcdf.putAtt(ncID, varID_quasi_pardepol_532, 'unit', '');
 netcdf.putAtt(ncID, varID_quasi_pardepol_532, 'long_name', 'quasi particle depolarization ratio at 532 nm');
diff --git a/lib/misc/dynamic_smooth.m b/lib/misc/dynamic_smooth.m
new file mode 100644
index 00000000..e63c4063
--- /dev/null
+++ b/lib/misc/dynamic_smooth.m
@@ -0,0 +1,789 @@
+function [c,ww] = dynamic_smooth(varargin)
+%SMOOTH  Smooth data.
+%   Z = SMOOTH(Y) smooths data Y using a 5-point moving average.
+%
+%   Z = SMOOTH(Y,SPAN) smooths data Y using SPAN as the number of points used
+%   to compute each element of Z.
+%
+%   Z = SMOOTH(Y,SPAN,METHOD) smooths data Y with specified METHOD. The
+%   available methods are:
+%
+%           'moving'   - Moving average (default)
+%           'dynamic_window' - Linearly increasing window, full bins
+%           'dynamic_window2' - Linearly increasing window, partial bins
+%           'lowess'   - Lowess (linear fit)
+%           'loess'    - Loess (quadratic fit)
+%           'sgolay'   - Savitzky-Golay
+%           'rlowess'  - Robust Lowess (linear fit)
+%           'rloess'   - Robust Loess (quadratic fit)
+%
+%   Z = SMOOTH(Y,METHOD) uses the default SPAN 5.
+%
+%   Z = SMOOTH(Y,SPAN,'sgolay',DEGREE) and Z = SMOOTH(Y,'sgolay',DEGREE)
+%   additionally specify the degree of the polynomial to be used in the
+%   Savitzky-Golay method. The default DEGREE is 2. DEGREE must be smaller
+%   than SPAN.
+%
+%   Z = SMOOTH(X,Y,...) additionally specifies the X coordinates.  If X is
+%   not provided, methods that require X coordinates assume X = 1:N, where
+%   N is the length of Y.
+%
+%   Notes:
+%   1. When X is given and X is not uniformly distributed, the default method
+%   is 'lowess'.  The 'moving' method is not recommended.
+%
+%   2. For the 'moving' and 'sgolay' methods, SPAN must be odd.
+%   If an even SPAN is specified, it is reduced by 1.
+%
+%   3. If SPAN is greater than the length of Y, it is reduced to the
+%   length of Y.
+%
+%   4. In the case of (robust) lowess and (robust) loess, it is also
+%   possible to specify the SPAN as a percentage of the total number
+%   of data points. When SPAN is less than or equal to 1, it is
+%   treated as a percentage.
+%
+%   For example:
+%
+%   Z = SMOOTH(Y) uses the moving average method with span 5 and
+%   X=1:length(Y).
+%
+%   Z = SMOOTH(Y,7) uses the moving average method with span 7 and
+%   X=1:length(Y).
+%
+%   Z = SMOOTH(Y,'sgolay') uses the Savitzky-Golay method with DEGREE=2,
+%   SPAN = 5, X = 1:length(Y).
+%
+%   Z = SMOOTH(X,Y,'lowess') uses the lowess method with SPAN=5.
+%
+%   Z = SMOOTH(X,Y,SPAN,'rloess') uses the robust loess method.
+%
+%   Z = SMOOTH(X,Y) where X is unevenly distributed uses the
+%   'lowess' method with span 5.
+%
+%   Z = SMOOTH(X,Y,8,'sgolay') uses the Savitzky-Golay method with
+%   span 7 (8 is reduced by 1 to make it odd).
+%
+%   Z = SMOOTH(X,Y,0.3,'loess') uses the loess method where span is
+%   30% of the data, i.e. span = ceil(0.3*length(Y)).
+%
+%   See also SPLINE.
+
+%   Copyright 2001-2016 The MathWorks, Inc.
+
+if nargin < 1
+    error(message('curvefit:smooth:needMoreArgs'));
+end
+
+if nargout > 1 % Called from the GUI cftool
+    ws = warning('off', 'all'); % turn warning off and record the previous warning state.
+    [lw,lwid] = lastwarn;
+    lastwarn('');
+else
+    ws = warning('query','all'); % Leave warning state alone but save it so resets are no-ops.
+end
+
+% is x given as the first argument?
+if nargin==1 || ( nargin > 1 && (length(varargin{2})==1 || ischar(varargin{2})) )
+    % smooth(Y) | smooth(Y,span,...) | smooth(Y,method,...)
+    is_x = 0; % x is not given
+    y = varargin{1};
+    y = y(:);
+    x = (1:length(y))';
+else % smooth(X,Y,...)
+    is_x = 1;
+    y = varargin{2};
+    x = varargin{1};
+    y = y(:);
+    x = x(:);
+end
+
+% is span given?
+span = [];
+if nargin == 1+is_x || ischar(varargin{2+is_x})
+    % smooth(Y), smooth(X,Y) || smooth(X,Y,method,..), smooth(Y,method)
+    is_span = 0;
+else
+    % smooth(...,SPAN,...)
+    is_span = 1;
+    span = varargin{2+is_x};
+end
+
+% is method given?
+method = [];
+if nargin >= 2+is_x+is_span
+    % smooth(...,Y,method,...) | smooth(...,Y,span,method,...)
+    method = varargin{2+is_x+is_span};
+end
+
+t = length(y);
+if t == 0
+    c = y;
+    ww = '';
+    if nargout > 1
+        ww = lastwarn;
+        lastwarn(lw,lwid);
+        warning(ws);  % turn warning back to the previous state.
+    end
+    return
+elseif length(x) ~= t
+    warning(ws); % reset warn state before erroring
+    error(message('curvefit:smooth:XYmustBeSameLength'));
+end
+
+% sets the rate of increase of the window size
+slope = 0.01;
+
+if isempty(method)
+    diffx = diff(x);
+    if uniformx(diffx,x,y)
+        if 0<slope && slope<=1
+		    method = 'dynamic_window';
+		else
+	        method = 'moving'; % uniformly distributed X.
+        end
+    else
+        method = 'lowess';
+    end
+end
+
+% realize span
+if span <= 0
+    warning(ws); % reset warn state before erroring
+    error(message('curvefit:smooth:spanMustBePositive'));
+end
+if span < 1, span = ceil(span*t); end % percent convention
+if isempty(span), span = 5; end % smooth(Y,[],method)
+
+idx = 1:t;
+
+sortx = any(diff(isnan(x))<0);   % if NaNs not all at end
+if sortx || any(diff(x)<0) % sort x
+    [x,idx] = sort(x);
+    y = y(idx);
+end
+
+if islogical(y)
+    y = double(y);
+end
+
+c = NaN(size(y),'like',y);
+
+ok = ~isnan(x);
+switch method
+    case 'moving'
+        c(ok) = moving(x(ok),y(ok),span);
+    case 'dynamic_window'
+        c(ok) = dynamic_window(x(ok),y(ok),span,slope);
+    case 'dynamic_window2'
+        c(ok) = dynamic_window2(x(ok),y(ok),span,slope);
+    case {'lowess','loess','rlowess','rloess'}
+        robust = 0;
+        iter = 5;
+        if method(1)=='r'
+            robust = 1;
+            method = method(2:end);
+        end
+        c(ok) = lowess(x(ok),y(ok),span, method,robust,iter);
+    case 'sgolay'
+        if nargin >= 3+is_x+is_span
+            degree = varargin{3+is_x+is_span};
+        else
+            degree = 2;
+        end
+        if degree < 0 || degree ~= floor(degree) || degree >= span
+            warning(ws); % reset warn state before erroring
+            error(message('curvefit:smooth:invalidDegree'));
+        end
+        c(ok) = sgolay(x(ok),y(ok),span,degree);
+    otherwise
+        warning(ws); % reset warn state before erroring
+        error(message('curvefit:smooth:unrecognizedMethod'));
+end
+
+c(idx) = c;
+
+if nargout > 1
+    ww = lastwarn;
+    lastwarn(lw,lwid);
+    warning(ws);  % turn warning back to the previous state.
+end
+%--------------------------------------------------------------------
+function c = moving(x,y, span)
+% moving average of the data.
+
+ynan = isnan(y);
+span = floor(span);
+n = length(y);
+span = min(span,n);
+width = span-1+mod(span,2); % force it to be odd
+xreps = any(diff(x)==0);
+if width==1 && ~xreps && ~any(ynan), c = y; return; end
+if ~xreps && ~any(ynan)
+    % simplest method for most common case
+    c = filter(ones(width,1)/width,1,y);
+    cbegin = cumsum(y(1:width-2));
+    cbegin = cbegin(1:2:end)./(1:2:(width-2))';
+    cend = cumsum(y(n:-1:n-width+3));
+    cend = cend(end:-2:1)./(width-2:-2:1)';
+    c = [cbegin;c(width:end);cend];
+elseif ~xreps
+    % with no x repeats, can take ratio of two smoothed sequences
+    yy = y;
+    yy(ynan) = 0;
+    nn = double(~ynan);
+    ynum = moving(x,yy,span);
+    yden = moving(x,nn,span);
+    c = ynum ./ yden;
+else
+    % with some x repeats, loop
+    notnan = ~ynan;
+    yy = y;
+    yy(ynan) = 0;
+    c = zeros(n,1,'like',y);
+    for i=1:n
+        if i>1 && x(i)==x(i-1)
+            c(i) = c(i-1);
+            continue;
+        end
+        R = i;                                 % find rightmost value with same x
+        while(R<n && x(R+1)==x(R))
+            R = R+1;
+        end
+        hf = ceil(max(0,(span - (R-i+1))/2));  % need this many more on each side
+        hf = min(min(hf,(i-1)), (n-R));
+        L = i-hf;                              % find leftmost point needed
+        while(L>1 && x(L)==x(L-1))
+            L = L-1;
+        end
+        R = R+hf;                              % find rightmost point needed
+        while(R<n && x(R)==x(R+1))
+            R = R+1;
+        end
+        c(i) = sum(yy(L:R)) / sum(notnan(L:R));
+    end
+end
+
+%%
+function c = dynamic_window(y,span,slope)
+%disp(size(y));
+if length(y)<10
+    c=y;
+else
+%Linearly increasing averaging window, only full bins
+
+%bins parameter determines the initial target avergaing bins to right and
+%left, window=2*size+1
+bins=(span-1)/2;
+%slope parameter determines the bins increment has to be between 0 and 1
+%slope = slope;
+
+%calculate height-bin (xmax) where max window size (bmax) is reached
+%length(y)-bmax=size+floor(slope*(bmax-size-1))
+xmax=ceil((length(y)+slope+bins*(slope-1))/(slope+1));
+bmax=length(y)-xmax;
+
+%region 1: increase bins up to target value
+reg1y = cumsum(y(1:2*bins+1));
+reg1y = reg1y(1:2:end);
+reg1x = (1:2:2*bins+1);
+
+%region 2: increase bins by slope up to xmax/bmax
+%filter determines when to increase the bins
+filter = floor(slope*(1:xmax-bins-2))-floor(slope*(0:xmax-bins-3));
+
+reg2x = 2*floor(bins+slope*(1:xmax-bins-2))+1;
+reg2y = NaN(1,length(reg2x));
+
+cbins = bins;
+csum = reg1y(end);
+for i = bins+2:xmax-1
+   csum = csum + y(i+cbins);
+   if filter(i-bins-1)
+       cbins = cbins+1;
+       csum = csum + y(i+cbins);
+   else
+       csum = csum - y(i-cbins-1);
+   end
+   reg2y(i-bins-1)=csum;
+end
+
+%region 3: decrease
+reg3y = cumsum(y(length(y):-1:xmax-bmax));
+reg3y = reg3y(end:-2:1);
+reg3x = 2*(bmax:-1:0)+1;
+
+%mean and combine
+c=[(reg1y(:)./reg1x(:)).', (reg2y(:)./reg2x(:)).', (reg3y(:)./reg3x(:)).'].';
+end
+
+%%
+function c = dynamic_window2(y,span,slope)
+%disp(size(y));
+if length(y)<10
+    c=y;
+else
+%Linearly increasing averaging window, also partial bins
+
+%bins parameter determines the initial target avergaing bins to right and
+%left, window=2*size+1
+bins=(span-1)/2;
+%slope parameter determines the bins increment has to be between 0 and 1
+%slope = slope;
+
+%calculate max bin point
+%length(y)-bmax=size+floor(slope*(bmax-size-1))
+xmax=ceil((length(y)+slope+bins*(slope-1))/(slope+1));
+bmax=length(y)-xmax;
+
+%region 1: increase bins up to target value
+reg1y = cumsum(y(1:2*bins+1));
+reg1y = reg1y(1:2:end);
+reg1x = (1:2:2*bins+1);
+
+%region 2: increase bins by slope up to xmax
+reg2x = 2*(bins+slope*(1:xmax-bins-2))+1;
+reg2y = NaN(1,length(reg2x));
+
+cbins = bins;
+csum = reg1y(end);
+factor = 0;
+for i = bins+2:xmax-1
+   csum = csum - y(i-cbins-1) + y(i+cbins);
+   factor=factor+slope;
+   % bins are increased when factor exceeds 1
+   if floor(factor)==1
+       factor = factor - 1;
+       cbins = cbins+1;
+       csum = csum + y(i-cbins) + y(i+cbins);
+   end
+   reg2y(i-bins-1)=csum + factor * (y(i-cbins-1) + y(i+cbins+1));
+end
+
+%region 3: decrease
+reg3y = cumsum(y(length(y):-1:xmax-bmax));
+reg3y = reg3y(end:-2:1);
+reg3x = 2*(bmax:-1:0)+1;
+
+%mean and combine
+c=[(reg1y(:)./reg1x(:)).', (reg2y(:)./reg2x(:)).', (reg3y(:)./reg3x(:)).'].';
+end
+
+%--------------------------------------------------------------------
+function c = lowess(x,y, span, method, robust, iter)
+% LOWESS  Smooth data using Lowess or Loess method.
+%
+% The difference between LOWESS and LOESS is that LOWESS uses a
+% linear model to do the local fitting whereas LOESS uses a
+% quadratic model to do the local fitting. Some other software
+% may not have LOWESS, instead, they use LOESS with order 1 or 2 to
+% represent these two smoothing methods.
+%
+% Reference: 
+% [C79] W.S.Cleveland, "Robust Locally Weighted Regression and Smoothing
+%    Scatterplots", _J. of the American Statistical Ass._, Vol 74, No. 368 
+%    (Dec.,1979), pp. 829-836.
+%    http://www.math.tau.ac.il/~yekutiel/MA%20seminar/Cleveland%201979.pdf
+
+
+n = length(y);
+span = floor(span);
+span = min(span,n);
+c = y;
+if span == 1
+    return;
+end
+
+useLoess = false;
+if isequal(method,'loess')
+    useLoess = true;
+end
+
+diffx = diff(x);
+
+% For problems where x is uniform, there's a faster way
+isuniform = uniformx(diffx,x,y);
+if isuniform
+    % For uniform data, an even span actually covers an odd number of
+    % points.  For example, the four closest points to 5 in the
+    % sequence 1:10 are {3,4,5,6}, but 7 is as close as 3.
+    % Therefore force an odd span.
+    span = 2*floor(span/2) + 1;
+
+    c = unifloess(y,span,useLoess);
+    if ~robust || span<=2
+        return;
+    end
+end
+
+% Turn off warnings when called from command line (already off if called from
+% cftool).
+ws = warning( 'off', 'MATLAB:rankDeficientMatrix' );
+cleanup = onCleanup( @() warning( ws ) );
+
+
+ynan = isnan(y);
+anyNans = any(ynan(:));
+seps = sqrt(eps);
+theDiffs = [1; diffx; 1];
+
+if isuniform
+    % We've already computed the non-robust smooth, so in preparation for
+    % the robust smooth, compute the following arrays directly
+    halfw = floor(span/2);
+    
+    % Each local interval is from |halfw| below the current index to |halfw|
+    % above
+    lbound = (1:n)-halfw;
+    rbound = (1:n)+halfw;
+    % However, there always has to be at least |span| points to the right of the
+    % left bound
+    lbound = min( n+1-span, lbound );
+    % ... and at least |span| points to the left of the right bound
+    rbound = max( span, rbound );
+    % Furthermore, because these bounds index into vectors of length n, they
+    % must contain valid indices
+    lbound = max( 1, lbound );
+    rbound = min( n, rbound );
+    
+    % Since the input is uniform we can use natural numbers for the input when
+    % we need them.
+    x = (1:numel(x))';
+else
+    if robust
+        % pre-allocate space for lower and upper indices for each fit,
+        % to avoid re-computing this information in robust iterations
+        lbound = zeros(n,1,'like',y);
+        rbound = zeros(n,1,'like',y);
+    end
+
+    % Compute the non-robust smooth for non-uniform x
+    for i=1:n
+        % if x(i) and x(i-1) are equal we just use the old value.
+        if theDiffs(i) == 0
+            c(i) = c(i-1);
+            if robust
+                lbound(i) = lbound(i-1);
+                rbound(i) = rbound(i-1);
+            end
+            continue;
+        end
+        
+        % Find nearest neighbours
+        idx = iKNearestNeighbours( span, i, x, ~ynan );
+        if robust
+            % Need to store neighborhoods for robust loop
+            lbound(i) = min(idx);
+            rbound(i) = max(idx);
+        end
+        
+        if isempty(idx)
+            c(i) = NaN;
+            continue
+        end
+
+        x1 = x(idx)-x(i); % center around current point to improve conditioning
+        d1 = abs(x1);
+        y1 = y(idx);
+
+        weight = iTricubeWeights( d1 );
+        if all(weight<seps)
+            weight(:) = 1;    % if all weights are 0, just skip weighting
+        end
+
+        v = [ones(size(x1)) x1];
+        if useLoess
+            v = [v x1.*x1]; %#ok<AGROW> There is no significant growth here
+        end
+        
+        v = weight(:,ones(1,size(v,2))).*v;
+        y1 = weight.*y1;
+        if size(v,1)==size(v,2)
+            % Square v may give infs in the \ solution, so force least squares
+            b = [v;zeros(1,size(v,2))]\[y1;0];
+        else
+            b = v\y1;
+        end
+        c(i) = b(1);
+    end
+end
+
+% now that we have a non-robust fit, we can compute the residual and do
+% the robust fit if required
+maxabsyXeps = max(abs(y))*eps;
+if robust
+    for k = 1:iter
+        r = y-c;
+        
+        % Compute robust weights
+        rweight = iBisquareWeights( r, maxabsyXeps ); 
+        
+        % Find new value for each point.
+        for i=1:n
+            if i>1 && x(i)==x(i-1)
+                c(i) = c(i-1);
+                continue;
+            end
+            if isnan(c(i))
+                continue; 
+            end
+            
+            idx = lbound(i):rbound(i);
+            if anyNans
+                idx = idx(~ynan(idx));
+            end
+            % check robust weights for removed points
+            if any( rweight(idx) <= 0 )
+                idx = iKNearestNeighbours( span, i, x, (rweight > 0) );
+            end
+            
+            x1 = x(idx) - x(i);
+            d1 = abs(x1);
+            y1 = y(idx);
+
+            weight = iTricubeWeights( d1 );
+            if all(weight<seps)
+                weight(:) = 1;    % if all weights 0, just skip weighting
+            end
+
+            v = [ones(size(x1)) x1];
+            if useLoess
+                v = [v x1.*x1]; %#ok<AGROW> There is no significant growth here
+            end
+            
+            % Modify the weights based on x values by multiplying them by
+            % robust weights.
+            weight = weight.*rweight(idx);
+            
+            v = weight(:,ones(1,size(v,2))).*v;
+            y1 = weight.*y1;
+            if size(v,1)==size(v,2)
+                % Square v may give infs in the \ solution, so force least squares
+                b = [v;zeros(1,size(v,2))]\[y1;0];
+            else
+                b = v\y1;
+            end
+            c(i) = b(1);
+        end
+    end
+end
+%--------------------------------------------------------------------
+function c=sgolay(x,y,f,k)
+% savitziki-golay smooth
+% (x,y) are given data. f is the frame length to be taken, should
+% be an odd number. k is the degree of polynomial filter. It should
+% be less than f.
+
+% Reference: Orfanidis, S.J., Introduction to Signal Processing,
+% Prentice-Hall, Englewood Cliffs, NJ, 1996.
+
+n = length(x);
+f = floor(f);
+f = min(f,n);
+f = f-mod(f-1,2); % will subtract 1 if frame is even.
+diffx = diff(x);
+notnan = ~isnan(y);
+nomissing = all(notnan);
+if f <= k && all(diffx>0) && nomissing, c = y; return; end
+hf = (f-1)/2; % half frame length
+
+idx = 1:n;
+if any(diffx<0) % make sure x is monotonically increasing
+    [x,idx]=sort(x);
+    y = y(idx);
+    notnan = notnan(idx);
+    diffx = diff(x);
+end
+% note that x is sorted so max(abs(x)) must be abs(x(1)) or abs(x(end));
+% already calculated diffx for monotonic case, so use it again. Only
+% recalculate if we sort x.
+if nomissing && uniformx(diffx,x,y)
+    v = ones(f,k+1);
+    t=(-hf:hf)';
+    for i=1:k
+        v(:,i+1)=t.^i;
+    end
+    [q,~]=qr(v,0);
+    ymid = filter(q*q(hf+1,:)',1,y);
+    ybegin = q(1:hf,:)*q'*y(1:f);
+    yend = q((hf+2):end,:)*q'*y(n-f+1:n);
+    c = [ybegin;ymid(f:end);yend];
+    return;
+end
+
+% non-uniformly distributed data
+c = y;
+
+% Turn off warnings when called from command line (already off if called from
+% cftool).
+ws = warning('off', 'all');
+[lastwarnmsg,lastwarnid]=lastwarn;
+
+for i = 1:n
+    if i>1 && x(i)==x(i-1)
+        c(i) = c(i-1);
+        continue
+    end
+    L = i; R = i;                          % find leftmost and rightmost values
+    while(R<n && x(R+1)==x(i))
+        R = R+1;
+    end
+    while(L>1 && x(L-1)==x(i))
+        L = L-1;
+    end
+    HF = ceil(max(0,(f - (R-L+1))/2));     % need this many more on each side
+
+    L = min(n-f+1,max(1,L-HF));            % find leftmost point needed
+    while(L>1 && x(L)==x(L-1))
+        L = L-1;
+    end
+    R = min(n,max(R+HF,L+f-1));            % find rightmost point needed
+    while(R<n && x(R)==x(R+1))
+        R = R+1;
+    end
+
+    tidx = L:R;
+    tidx = tidx(notnan(tidx));
+    if isempty(tidx)
+        c(i) = NaN;
+        continue;
+    end
+    q = x(tidx) - x(i);   % center to improve conditioning
+    vrank = 1 + sum(diff(q)>0);
+    ncols = min(k+1, vrank);
+    v = ones(length(q),ncols,'like',q);
+    for j = 1:ncols-1
+        v(:,j+1) = q.^j;
+    end
+    if size(v,1)==size(v,2)
+        % Square v may give infs in the \ solution, so force least squares
+        d = [v;zeros(1,size(v,2))]\[y(tidx);0];
+    else
+        d = v\y(tidx);
+    end
+    c(i) = d(1);
+end
+c(idx) = c;
+
+lastwarn(lastwarnmsg,lastwarnid);
+warning(ws);
+%--------------------------------------------------------------------
+function ys = unifloess(y,span,useLoess)
+%UNIFLOESS Apply loess on uniformly spaced X values
+
+y = y(:);
+
+% Omit points at the extremes, which have zero weight
+halfw = (span-1)/2;              % halfwidth of entire span
+d = abs((1-halfw:halfw-1));      % distances to pts with nonzero weight
+dmax = halfw;                    % max distance for tri-cubic weight
+
+% Set up weighted Vandermonde matrix using equally spaced X values
+x1 = (2:span-1)-(halfw+1);
+weight = (1 - (d/dmax).^3).^1.5; % tri-cubic weight
+v = [ones(length(x1),1) x1(:)];
+if useLoess
+    v = [v x1(:).^2];
+end
+V = v .* repmat(weight',1,size(v,2));
+
+% Do QR decomposition
+[Q,~] = qr(V,0);
+
+% The projection matrix is Q*Q'.  We want to project onto the middle
+% point, so we can take just one row of the first factor.
+alpha = Q(halfw,:)*Q';
+
+% This alpha defines the linear combination of the weighted y values that
+% yields the desired smooth values.  Incorporate the weights into the
+% coefficients of the linear combination, then apply filter.
+alpha = alpha .* weight;
+ys = filter(alpha,1,y);
+
+% We need to slide the values into the center of the array.
+ys(halfw+1:end-halfw) = ys(span-1:end-1);
+
+% Now we have taken care of everything except the end effects.  Loop over
+% the points where we don't have a complete span.  Now the Vandermonde
+% matrix has span-1 points, because only 1 has zero weight.
+x1 = 1:span-1;
+v = [ones(length(x1),1) x1(:)];
+if useLoess
+    v = [v x1(:).^2];
+end
+for j=1:halfw
+    % Compute weights based on deviations from the jth point,
+    % then compute weights and apply them as above.
+    d = abs((1:span-1) - j);
+    weight = (1 - (d/(span-j)).^3).^1.5;
+    V = v .* repmat(weight(:),1,size(v,2));
+    [Q,~] = qr(V,0);
+    alpha = Q(j,:)*Q';
+    alpha = alpha .* weight;
+    ys(j) = alpha * y(1:span-1);
+
+    % These coefficients can be applied to the other end as well
+    ys(end+1-j) = alpha * y(end:-1:end-span+2);
+end
+%--------------------------------------------------------------------
+function isuniform = uniformx(diffx,x,y)
+%ISUNIFORM True if x is of the form a:b:c
+
+if any(isnan(y)) || any(isnan(x))
+    isuniform = false;
+else
+    isuniform = all(abs(diff(diffx)) <= eps*max(abs([x(1),x(end)])));
+end
+%--------------------------------------------------------------------
+function idx = iKNearestNeighbours( k, i, x, in )
+% Find the k points from x(in) closest to x(i)
+
+if nnz( in ) <= k
+    % If we have k points or fewer, then return them all
+    idx = find( in );
+else
+    % Find the distance to the k closest point
+    d = abs( x - x(i) );
+    ds = sort( d(in) );
+    dk = ds(k);
+    
+    % Find all points that are as close as or closer than the k closest point
+    close = (d <= dk);
+    
+    % The required indices are those points that are both close and "in"
+    idx = find( close & in );
+end
+%--------------------------------------------------------------------
+% Bi-square (robust) weight function
+function delta = iBisquareWeights( r, myeps )
+% Convert residuals to weights using the bi-square weight function.
+% NOTE that this function returns the square root of the weights
+
+% Only use non-NaN residuals to compute median
+idx = ~isnan( r );
+% And bound the median away from zero
+s = max( 1e8 * myeps, median( abs( r(idx) ) ) );
+% Covert the residuals to weights
+delta = iBisquare( r/(6*s) );
+% Everything with NaN residual should have zero weight
+delta(~idx) = 0;
+function b = iBisquare( x )
+% This is this bi-square function defined at the top of the left hand
+% column of page 831 in [C79]
+% NOTE that this function returns the square root of the weights
+b = zeros( size( x ) , 'like', x);
+idx = abs( x ) < 1;
+b(idx) = abs( 1 - x(idx).^2 );
+%--------------------------------------------------------------------
+% Tri-cubic weight function
+function w = iTricubeWeights( d )
+% Convert distances into weights using tri-cubic weight function.
+% NOTE that this function returns the square-root of the weights.
+%
+% Protect against divide-by-zero. This can happen if more points than the span
+% are coincident.
+maxD = max( d );
+if maxD > 0
+    d = d/max( d );
+end
+w = (1 - d.^3).^1.5;
diff --git a/lib/retrievals/pollyFernald.m b/lib/retrievals/pollyFernald.m
index 98282505..cb9a504e 100644
--- a/lib/retrievals/pollyFernald.m
+++ b/lib/retrievals/pollyFernald.m
@@ -91,8 +91,13 @@
 RCS = reshape(signal, 1, numel(signal)) .* reshape(alt, 1, numel(alt)).^2;
 
 indRefMid = int32(mean(indRefAlt));
+
 % smooth the signal at the reference height region
-RCS = smooth(RCS, window_size, 'moving');
+% RCS = smooth(RCS, window_size, 'moving'); %old function
+%RCS = dynamic_smooth(RCS, window_size, 'moving'); %new function
+RCS = dynamic_smooth(RCS, window_size, 'dynamic_window'); %new method (full bins)
+%RCS = dynamic_smooth(RCS, window_size, 'dynamic_window2'); %new method (partial bins)
+
 RCS(indRefMid) = mean(RCS(indRefAlt(1):indRefAlt(2)));
 
 % intialize some parameters and set the value at the reference altitude.