make master up to date

config/pollyDefaults/template_polly_defaults.json

@@ -3,14 +3,24 @@
     "depolCaliConstStd532": 0.0,
     "depolCaliConst355": 0.1201,
     "depolCaliConstStd355": 0.0,
+	"polCaliEta355": 33,
+    "polCaliEtaStd355": 1,
+	"polCaliEta532": 11,
+    "polCaliEtaStd532": 1,
+    "polCaliEta1064": 0.7,
+    "polCaliEtaStd1064": 1,
     "LC": [16893959541573.252000, 1, 1, 1, 28466210203696.203000, 1, 1, 67876166905915.734000, 1, 1, 1, 1],
     "LCStd": [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
     "overlapFile532": "template_polly_overlap_532.txt",
     "overlapFile355": "template_polly_overlap_355.txt",
-    "molDepol532": 0.0053,
+    "molDepol532": 0.0044,
     "molDepolStd532": 0.001,
-    "molDepol355": 0.0075,
+    "molDepol355": 0.0077,
     "molDepolStd355": 0.0015,
+	"molDepol1064": 0.0036,
     "wvconst": 11.1,
-    "wvconstStd": 0
+    "wvconstStd": 0,
+	"volDepolerror355": 0.003,
+	"volDepolerror532": 0.004,
+	"volDepolerror1064": 0.005
 }

lib/calibration/depolCaliGHK.m

@@ -0,0 +1,264 @@
+function [polCaliEta, polCaliEtaStd, polCaliStartTime, polCaliStopTime,cali_status, globalAttri ] = depolCaliGHK(signal_t, ...
+        bg_t, signal_x, bg_x, time, polCaliPAngStartTime, ...
+        polCaliPAngStopTime, polCaliNAngStartTime, ...
+        polCaliNAngStopTime, K, caliHIndxRange, ...
+        SNRmin, sigMax, rel_std_dplus, rel_std_dminus, segmentLen, smoothWin)
+% DEPOLCALI polarization calibration for PollyXT lidar system.
+%
+% USAGE:
+%    [polCaliEta, polCaliEtaStd, depol_cal_time] = depolCaliGHK(signal_t, 
+%        bg_t, signal_x, bg_x, time, polCaliPAngStartTime, 
+%        polCaliPAngStopTime, polCaliNAngStartTime, 
+%        polCaliNAngStopTime, K, caliHIndxRange, 
+%        SNRmin, sigMax, rel_std_dplus, rel_std_dminus, 
+%        segmentLen, smoothWin, flagShowResults)
+%
+% INPUTS:
+%    signal_t: matrix
+%        background-removed photon count signal at total channel.
+%        (nBins * nProfiles)
+%    bg_t: matrix
+%        background at total channel. (nBins * nProfiles)
+%    signal_x: matrix
+%        background-removed photon count signal at cross channel.
+%        (nBins * nProfiles)
+%    bg_x: matrix
+%        background at cross channel. (nBins * nProfiles)
+%    time: array
+%        datenum array represents the measurement time of each profile.
+%    polCaliPAngStartTime: array
+%        datenum array represents the start time that the polarizer
+%        rotates to the positive angle. 
+%    polCaliPAngStopTime: array
+%        datenum array represents the stop time that the polarizer
+%        rotates to the positive angle.
+%    polCaliNAngStartTime array
+%        datenum array represents the start time that the polarizer
+%        rotates to the negative angle.
+%    polCaliNAngStopTime: array
+%        datenum array represents the end time that the polarizer
+%        rotates to the negative angle.
+%    K: float
+%        K parameter from GHK to correct the calibration.
+%    caliHIndxRange: 2-element array
+%        range of height indexes at which the signal can be used for
+%        polarization calibration.
+%    SNRmin: array
+%        minimum SNR for calibration.
+%    sigMax: array
+%        maximum signal that could be used in the calibration to prevent
+%        pulse pileup effects. (Photon Count)
+%    rel_std_dplus: float
+%        maximum relative uncertainty of dplus that is allowed.
+%    rel_std_dplus: float
+%        maximum relative uncertainty of dminus that is allowed.
+%    segmentLen: integer
+%        segement length for testing the variability of the calibration results
+%        to prevent of cloud contamintaion.
+%    smoothWin: integer
+%        width of the sliding window for smoothing the signal.
+%    flagShowResults: logical
+%        flag to control whether to save the intermediate results.
+%
+% OUTPUTS:
+%    polCaliEta: array
+%        eta from polarization calibration.
+%    polCaliEtaStd: array
+%        uncertainty of eta from polarizatrion calibration.
+%    polCaliStartTime: array
+%        start time for each successful calibration.
+%    polCaliStopTime: array
+%        stop time for each successful calibration.
+%    cali_status: integer
+%        1 if calibration is successfull, 0 otherwise.
+%    globalAttri: struct
+%        all the information about the depol calibration.
+%
+% REFERENCE:
+%    1. Freudenthaler, V.: About the effects of polarising optics on lidar signals and the Delta90��calibration, Atmos. Meas. Tech., 9, 4181-4255, 10.5194/amt-9-4181-2016, 2016.
+%
+% HISTORY:
+%    - 2018-07-25: First edition by Zhenping.
+%    - 2019-06-08: If no depol cali, return empty array.
+%    - 2019-09-06: Remove the part to replace the bins of low SNR with NaN, because it will lead to bias when doing smoothing.
+%    - 2024-08-27: Transfrom to GHK parameters using only eta and not V* (polCaliFac) any more 
+%    - 2024-10-02: Add calibration status variable to output. 
+%
+% .. Authors: - zhenping@tropos.de, haarig@tropos.de
+
+%% parameters initialization
+polCaliEta = [];
+polCaliEtaStd = [];
+mean_dminus = [];
+mean_dplus = [];
+std_dminus = [];
+std_dplus = [];
+polCaliStartTime = [];
+polCaliStopTime = [];
+cali_status = 0;
+globalAttri = struct();
+globalAttri.sig_t_p = cell(0);
+globalAttri.sig_t_m = cell(0);
+globalAttri.sig_x_p = cell(0);
+globalAttri.sig_x_m = cell(0);
+globalAttri.caliHIndxRange = cell(0);
+globalAttri.indx_45m = cell(0);
+globalAttri.indx_45p = cell(0);
+globalAttri.dplus = cell(0);
+globalAttri.dminus = cell(0);
+globalAttri.segmentLen = cell(0);
+globalAttri.indx = cell(0);
+globalAttri.mean_dplus_tmp = cell(0);
+globalAttri.std_dplus_tmp = cell(0);
+globalAttri.mean_dminus_tmp = cell(0);
+globalAttri.std_dminus_tmp = cell(0);
+globalAttri.K = cell(0);
+globalAttri.segIndx = cell(0);
+globalAttri.caliTime = cell(0);
+
+if isempty(signal_t) || isempty(signal_x) 
+    warning('No data for polarization calibration.');
+    return;
+end
+
+days = unique(fix(time));
+nDays = length(days);
+
+for iDay = 1:nDays
+    for iDepolCal = 1:length(polCaliNAngStartTime)
+        indx_45p = find(time >= days(iDay) & time < (days(iDay) + 1) & ...
+                        time >= polCaliPAngStartTime(iDepolCal) & ...
+                        time <= polCaliPAngStopTime(iDepolCal));
+        indx_45m = find(time >= days(iDay) & time < (days(iDay) + 1) & ...
+                        time >= polCaliNAngStartTime(iDepolCal) & ... 
+                        time <= polCaliNAngStopTime(iDepolCal));
+
+        if (length(indx_45p) < 4) || (length(indx_45m) < 4)
+            % if not enough depol cali profiles were found, break the loop
+            break;
+        end
+
+        thisCaliStartTime = min([polCaliPAngStartTime(iDepolCal), ...
+                                 polCaliNAngStartTime(iDepolCal)]);
+        thisCaliStopTime = max([polCaliPAngStopTime(iDepolCal), ...
+                                polCaliNAngStopTime(iDepolCal)]);
+
+        % neglect the first and last profile which could be unstable due to
+        % the rotation of the polarizer
+        indx_45m = indx_45m(2:end-1);
+        indx_45p = indx_45p(2:end-1);
+
+        sig_t_p = nanmean(signal_t(:, indx_45p), 2);
+        bg_t_p = nanmean(bg_t(:, indx_45p), 2);
+        SNR_t_p = pollySNR(sig_t_p, bg_t_p);
+        indxBad_t_p = (SNR_t_p <= SNRmin(1)) | (sig_t_p >= sigMax(1));
+
+        sig_t_m = nanmean(signal_t(:, indx_45m), 2);
+        bg_t_m = nanmean(bg_t(:, indx_45m), 2);
+        SNR_t_m = pollySNR(sig_t_m, bg_t_m);
+        indxBad_t_m = (SNR_t_m <= SNRmin(2)) | (sig_t_m >= sigMax(2));
+
+        sig_x_p = nanmean(signal_x(:, indx_45p), 2);
+        bg_x_p = nanmean(bg_x(:, indx_45p), 2);
+        SNR_x_p = pollySNR(sig_x_p, bg_x_p);
+        indxBad_x_p = (SNR_x_p <= SNRmin(3)) | (sig_x_p >= sigMax(3));
+
+        sig_x_m = nanmean(signal_x(:, indx_45m), 2);
+        bg_x_m = nanmean(bg_x(:, indx_45m), 2);
+        SNR_x_m = pollySNR(sig_x_m, bg_x_m);
+        indxBad_x_m = (SNR_x_m <= SNRmin(4)) | (sig_x_m >= sigMax(4));
+
+
+        dplus = smooth(sig_x_p, 'moving', smoothWin) ./ ...
+                smooth(sig_t_p, 'moving', smoothWin);
+        dminus = smooth(sig_x_m, 'moving', smoothWin) ./ ...
+                 smooth(sig_t_m, 'moving', smoothWin);
+        dplus(isinf(dplus)) = NaN;
+        dminus(isinf(dminus)) = NaN;
+        dplus(indxBad_t_p | indxBad_x_p) = NaN;
+        dminus(indxBad_t_m | indxBad_x_m) = NaN;
+        dplus = dplus(caliHIndxRange(1):caliHIndxRange(2));
+        dminus = dminus(caliHIndxRange(1):caliHIndxRange(2));
+
+        mean_dplus_tmp = [];
+        std_dplus_tmp = [];
+        mean_dminus_tmp = [];
+        std_dminus_tmp = [];
+        segIndx_tmp = [];
+        % find the most stable region where the realtive std of the signal
+        % is less than rel_std_dminus and rel_std_dplus
+        for iReg = 1:(caliHIndxRange(2) - caliHIndxRange(1) - segmentLen)
+
+            if sum(~ isnan(dplus(iReg:(iReg + segmentLen)))) ...
+                <= segmentLen/4 || ...
+               sum(~ isnan(dminus(iReg:(iReg + segmentLen)))) <= segmentLen/4
+                continue;
+            end
+
+            this_mean_dplus = nanmean(dplus(iReg:(iReg + segmentLen)));
+            this_std_dplus = nanstd(dplus(iReg:(iReg + segmentLen)));
+            this_mean_dminus = nanmean(dminus(iReg:(iReg + segmentLen)));
+            this_std_dminus = nanstd(dminus(iReg:(iReg + segmentLen)));
+
+            if abs(this_std_dminus / this_mean_dminus) <= rel_std_dminus && ...
+               abs(this_std_dplus / this_mean_dplus) <= rel_std_dplus
+                segIndx_tmp = cat(2, segIndx_tmp, iReg);
+                mean_dplus_tmp = cat(2, mean_dplus_tmp, this_mean_dplus);
+                mean_dminus_tmp = cat(2, mean_dminus_tmp, this_mean_dminus);
+                std_dplus_tmp = cat(2, std_dplus_tmp, this_std_dplus);
+                std_dminus_tmp = cat(2, std_dminus_tmp, this_std_dminus);
+            end 
+        end
+
+        % if there is no stable calibration segments, start the next 
+        % calibration 
+        if isempty(mean_dplus_tmp)
+            continue;
+        end
+
+        % find the most stable calbiration region
+        [~, segIndx] = min(sqrt((std_dplus_tmp./mean_dplus_tmp).^2 + ...
+                                (std_dminus_tmp./mean_dminus_tmp).^2));
+        indx = segIndx_tmp(segIndx);
+        polCaliStartTime = cat(2, polCaliStartTime, thisCaliStartTime);
+        polCaliStopTime = cat(2, polCaliStopTime, thisCaliStopTime);
+        mean_dplus = cat(2, mean_dplus, mean_dplus_tmp(segIndx));
+        std_dplus = cat(2, std_dplus, std_dplus_tmp(segIndx));
+        mean_dminus = cat(2, mean_dminus, mean_dminus_tmp(segIndx));
+        std_dminus = cat(2, std_dminus, std_dminus_tmp(segIndx));
+
+        % save the intermediate results
+        globalAttri.sig_t_p{end + 1} = sig_t_p;
+        globalAttri.sig_t_m{end + 1} = sig_t_m;
+        globalAttri.sig_x_p{end + 1} = sig_x_p;
+        globalAttri.sig_x_m{end + 1} = sig_x_m;
+        globalAttri.caliHIndxRange{end + 1} = caliHIndxRange;
+        globalAttri.indx_45m{end + 1} = indx_45m;
+        globalAttri.indx_45p{end + 1} = indx_45p;
+        globalAttri.dplus{end + 1} = dplus;
+        globalAttri.dminus{end + 1} = dminus;
+        globalAttri.segmentLen{end + 1} = segmentLen;
+        globalAttri.indx{end + 1} = indx;
+        globalAttri.mean_dplus_tmp{end + 1} = mean_dplus_tmp;
+        globalAttri.std_dplus_tmp{end + 1} = std_dplus_tmp;
+        globalAttri.mean_dminus_tmp{end + 1} = mean_dminus_tmp;
+        globalAttri.std_dminus_tmp{end + 1} = std_dminus_tmp;
+        globalAttri.K{end + 1} = K;
+        globalAttri.segIndx{end + 1} = segIndx;
+        globalAttri.caliTime{end + 1} = mean([thisCaliStartTime, thisCaliStopTime]);
+
+    end
+end
+
+if isempty(mean_dminus) || isempty(mean_dplus)
+    print_msg('Plus or minus 45� calibration is missing.\n', 'flagTimestamp', true);
+    cali_status =0; %Calibration was not successfull.
+    return;
+end
+
+% calculate the depol-calibration factor and unceratinty and correcting it
+% with K
+polCaliEta = 1/K .* sqrt(mean_dplus .* mean_dminus);
+polCaliEtaStd = 0.5 .* (mean_dplus .* std_dminus + mean_dminus .* std_dplus) ./ sqrt(mean_dplus .* mean_dminus);
+cali_status =1; % Calibration successfull.
+end

lib/calibration/pollyMDRGHK.m

@@ -0,0 +1,77 @@
+function [MDR, MDRStd, flagDeft] = pollyMDRGHK(sigT, bgT, sigC, bgC, flagT,flagC, eta, voldepol_error, minSNR, deftMDR, deftMDRStd, PollyConfig)
+% POLLYMDR etimate the molecular depolarization ratio according to the measurements at reference height.
+%
+% USAGE:
+%    [MDR, MDRStd, flagDeft] = pollyMDRGHK(sigT, bgT, ...
+%    sigC, bgC, flagT,flagC, eta, voldep_sys_uncertainty, ...
+%    minSNR, deftMDR, deftMDRStd)
+%
+% INPUTS:
+%    sigT: array
+%        signal strength of the total channel at reference height. 
+%        [photon count]
+%    bgT: array
+%        background of the total channel at reference height. [photon count]
+%    sigCross: array
+%        signal strength of the cross channel at reference height. 
+%        [photon count]
+%    bgCross: array
+%        background of the cross channel at reference height. [photon count]
+%    flagT: 
+%        flag the total channel for the respective wavelength.
+%    flagC: 
+%        flag the cross channel for the respective wavelength.
+%    eta: scalar
+%        depolarzation calibration constant. 
+%    voldep_sys_uncertainty: scalar
+%        systematic uncertainty of the volume depolarization ratio (in
+%        future it should be given in the config file)
+%    minSNR: float
+%        the SNR constrain for the the signal strength at reference height. 
+%        Choose a strong constrain for ensuring a stable result, 
+%        like 50 or 100.
+%    deftMDR: float
+%        default molecular depolarization ratio.
+%    deftMDRStd: float
+%        default std of molecular depolarization ratio.
+%    Polly.Config:
+%        contains the GHK parameters for further calcualtions. 
+%
+% OUTPUTS:
+%    MDR: float
+%        retrieved molecular depolarization ratio. 
+%    MDRStd: float
+%        std of retrieved molecular depolarization ratio. 
+%    flagDeft: logical
+%        flag to show whether using the default values. If true, 
+%        it means default MDR and MDRStd were used.
+%
+% HISTORY:
+%    - 2021-05-31: first edition by Zhenping
+%    - 2024-08-13: MH changed to GHK
+%
+% .. Authors: - zhenping@tropos.de, haarig@tropos.de
+
+MDR = deftMDR;
+MDRStd = deftMDRStd;
+flagDeft = true;
+
+snrTot = pollySNR(sum(sigT), sum(bgT));
+snrCro = pollySNR(sum(sigC), sum(bgC));
+
+flagValidPointTot = (snrTot >= minSNR);
+flagValidPointCro = (snrCro >= minSNR);
+
+if (~ flagValidPointCro) || (~ flagValidPointTot)
+    fprintf(['Too noisy in the reference height to calculate ' ...
+             'the molecular depol ratio.\n']);
+    return;
+end
+
+[MDR, MDRStd] = pollyVDRGHK(sum(sigT), sum(sigC), ...
+                            PollyConfig.G(flagT), PollyConfig.G(flagC), ...
+                            PollyConfig.H(flagT), PollyConfig.H(flagC), ...
+                            eta, voldepol_error(1), voldepol_error(2), voldepol_error(3), 1);
+flagDeft = false;
+
+end

lib/calibration/pollyPolCali.m

@@ -1,5 +1,5 @@
 function [polCaliEta, polCaliEtaStd, polCaliFac, polCaliFacStd, polCaliTime, polCaliAttri] = pollyPolCali(data, transRatio, varargin)
-% POLLYPOLCALI calibrate the PollyXT cross channels for 355 and 532 nm with ±45° method.
+% POLLYPOLCALI calibrate the PollyXT cross channels for 355, 532 and 1064 nm with Delta90� method.
 %
 % USAGE:
 %    [polCaliEta, polCaliEtaStd,polCaliFac, polCaliFacStd,  polCaliTime, polCaliAttri] = pollyPolCali(data, transRatio)
@@ -12,7 +12,7 @@
 %
 % KEYWORDS:
 %    wavelength: char
-%        '355nm' or '532nm'.
+%        '355nm' or '532nm' or '1064nm'
 %    depolCaliMinBin: numeric
 %        minimum search index for stable polarization calibration constants.
 %    depolCaliMaxBin
@@ -266,7 +266,6 @@
         polCaliTime = [polCaliStartTime1064, polCaliStopTime1064];
         polCaliAttri = polCalAttri1064;
     end
-
 otherwise
     error('Unknown wavelength %s for polarization calibration.', p.Results.wavelength);
 end