Moved "from __future__ import print_function" 0.9.4

Tue, 24 Jan 2017 02:26:19 +0100

author
Volker Freudenthaler <volker.freudenthaler@lmu.de>
date
Tue, 24 Jan 2017 02:26:19 +0100
changeset 21
857c95060313
parent 20
161490e56a2c
child 22
74ea89fcbeb2

Moved "from __future__ import print_function"
before "# !/usr/bin/env python3"
Now it should work with Python 3.x .

Added some comments and a compact code description at the top.

lidar_correction_ghk.py file | annotate | diff | comparison | revisions
output_files/output_example lidar.dat file | annotate | diff | comparison | revisions
system_settings/optic_input_example_lidar.py file | annotate | diff | comparison | revisions
--- a/lidar_correction_ghk.py	Mon Nov 28 16:45:48 2016 +0200
+++ b/lidar_correction_ghk.py	Tue Jan 24 02:26:19 2017 +0100
@@ -18,11 +18,31 @@
 Equation reference: http://www.atmos-meas-tech-discuss.net/amt-2015-338/amt-2015-338.pdf
 With equations code from Appendix C
 Python 3.4.2
+
+Code description:
+
+From measured lidar signals we cannot directly determine the desired backscatter coefficient (F11) and the linear depolarization ratio (LDR)
+because of the cross talk between the channles and systematic errors of a lidar system.
+http://www.atmos-meas-tech-discuss.net/amt-2015-338/amt-2015-338.pdf provides an analytical model for the description of these errors,
+with which the measured signals can be corrected.
+This code simulates the lidar measurements with "assumed true" model parameters from an input file, and calculates the correction parameters (G,H, and K).
+The "assumed true" system parameters are the ones we think are the right ones, but in reality these parameters probably deviate from the assumed truth due to
+uncertainties. The uncertainties of the "assumed true" parameters can be described in the input file. Then this code calculates the lidar signals and the
+gain ratio eta* with all possible combinations of "errors", which represents the distribution of "possibly real" signals, and "corrects" them with the "assumed true"
+GHK parameters (GT0, GR0, HT0, HR0, and K0) to derive finally the distributions of "possibly real" linear depolarization ratios (LDRcorr),
+which are plotted for five different input linear depolarization ratios (LDRtrue). The red bars in the plots represent the input values of LDRtrue.
+A complication arises from the fact that the correction parameter K = eta*/eta (Eq. 83) can depend on the LDR during the calibration measurement, i.e. LDRcal or aCal
+in the code (see e.g. Eqs. (103), (115), and (141); mind the mistake in Eq. (116)). Therefor values of K for LDRcal = 0.004, 0.2, and 0.45 are calculated for
+"assumed true" system parameters and printed in the output file behind the GH parameters. The full impact of the LDRcal dependent K can be considered in the error
+calculation by specifying a range of possible LDRcal values in the input file. For the real calibration measurements a calibration range with low or no aerosol
+content should be chosen, and the default in the input file is a range of LDRcal between 0.004 and 0.014 (i.e. 0.009 +-0.005).
+
+Tip: In case you run the code with Spyder, all output text and plots can be displayed together in an IPython console, which can be saved as an html file.
 """
-# !/usr/bin/env python3
 
 # Comment:  The code works with Python 2.7  with the help of following line, which enables Python2 to correctly interpret the Python 3 print statements.
 from __future__ import print_function
+# !/usr/bin/env python3
 
 import os
 import sys
@@ -49,7 +69,6 @@
 
 from contextlib import contextmanager
 
-
 @contextmanager
 def redirect_stdout(new_target):
     old_target, sys.stdout = sys.stdout, new_target  # replace sys.stdout
@@ -59,7 +78,6 @@
         sys.stdout.flush()
         sys.stdout = old_target  # restore to the previous value
 
-
 '''
 real_raw_input = vars(__builtins__).get('raw_input',input)
 '''
@@ -93,9 +111,9 @@
 
 sqr05 = 0.5 ** 0.5
 
+# ---- Initial definition of variables; the actual values will be read in with exec(open('./optic_input.py').read()) below
 # Do you want to calculate the errors? If not, just the GHK-parameters are determined.
 Error_Calc = True
-# ---- Initial definition of variables; the actual values will be read in with exec(open('./optic_input.py').read()) below
 LID = "internal"
 EID = "internal"
 # --- IL Laser IL and +-Uncertainty
@@ -184,40 +202,8 @@
 #  end of initial definition of variables
 # *******************************************************************************************************************************
 
-# --- Read actual lidar system parameters from ./optic_input.py  (must be in the same directory)
+# --- Read actual lidar system parameters from optic_input.py  (should be in sub-directory 'system_settings')
 InputFile = 'optic_input_example_lidar.py'
-# InputFile = 'optic_input_ver6e_POLIS_355.py'
-# InputFile = 'optic_input_ver6e_POLIS_355_JA.py'
-# InputFile = 'optic_input_ver6c_POLIS_532.py'
-# InputFile = 'optic_input_ver6e_POLIS_532.py'
-# InputFile = 'optic_input_ver8c_POLIS_532.py'
-# InputFile = 'optic_input_ver6e_MUSA.py'
-# InputFile = 'optic_input_ver6e_MUSA_JA.py'
-# InputFile = 'optic_input_ver6e_PollyXTSea.py'
-# InputFile = 'optic_input_ver6e_PollyXTSea_JA.py'
-# InputFile = 'optic_input_ver6e_PollyXT_RALPH.py'
-# InputFile = 'optic_input_ver8c_PollyXT_RALPH.py'
-# InputFile = 'optic_input_ver8c_PollyXT_RALPH_2.py'
-# InputFile = 'optic_input_ver8c_PollyXT_RALPH_3.py'
-# InputFile = 'optic_input_ver8c_PollyXT_RALPH_4.py'
-# InputFile = 'optic_input_ver8c_PollyXT_RALPH_5.py'
-# InputFile = 'optic_input_ver8c_PollyXT_RALPH_6.py'
-# InputFile = 'optic_input_ver8c_PollyXT_RALPH_7.py'
-# InputFile = 'optic_input_ver8a_MOHP_DPL_355.py'
-# InputFile = 'optic_input_ver9_MOHP_DPL_355.py'
-# InputFile = 'optic_input_ver6e_RALI.py'
-# InputFile = 'optic_input_ver6e_RALI_JA.py'
-# InputFile = 'optic_input_ver6e_RALI_new.py'
-# InputFile = 'optic_input_ver6e_RALI_act.py'
-# InputFile = 'optic_input_ver6e_MULHACEN.py'
-# InputFile = 'optic_input_ver6e_MULHACEN_JA.py'
-# InputFile = 'optic_input_ver6e-IPRAL.py'
-# InputFile = 'optic_input_ver6e-IPRAL_JA.py'
-# InputFile = 'optic_input_ver6e-LB21.py'
-# InputFile = 'optic_input_ver6e-LB21_JA.py'
-# InputFile = 'optic_input_ver6e_Bertha_b_355.py'
-# InputFile = 'optic_input_ver6e_Bertha_b_532.py'
-# InputFile = 'optic_input_ver6e_Bertha_b_1064.py'
 
 '''
 print("From ", dname)
@@ -225,10 +211,11 @@
 print("Reading input file ", InputFile, " for")
 '''
 input_path = os.path.join('.', 'system_settings', InputFile)
-# this works with Python 2 - and 3?
+# this works with Python 2 and 3!
 exec (open(input_path).read(), globals())
 #  end of read actual system parameters
 
+
 # --- Manual Parameter Change ---
 #  (use for quick parameter changes without changing the input file )
 # DiO = 0.
@@ -368,10 +355,10 @@
     # -------------------------
 
     # analyser
-    # For POLLY_XTs
     if (RS_RP_depend_on_TS_TP):
         RS = 1 - TS
         RP = 1 - TP
+
     TiT = 0.5 * (TP + TS)
     DiT = (TP - TS) / (TP + TS)
     ZiT = (1. - DiT ** 2) ** 0.5
@@ -853,14 +840,14 @@
     TRa = TiR * TaR  # *ARP1
 
     F11sim = 1 / (TiO * TiE) * (
-    (HR * Etax / K * It / TTa - HT * Ir / TRa) / (HR * GT - HT * GR))  # IL = 1, Etat = Etar = 1
+    (HR * Etax / K * It / TTa - HT * Ir / TRa) / (HR * GT - HT * GR))  # IL = 1, Etat = Etar = 1  ;  AMT Eq.64; what is Etax/K? => see about 20 lines above: = Eta
 
     return (GT, HT, GR, HR, K, Eta, LDRsimx, LDRCorr, DTa, DRa, TTa, TRa, F11sim)
 
 
 # *******************************************************************************************************************************
 
-# --- CALC truth
+# --- CALC truth with assumed true parameters from the input file
 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0,
                                                                                        RetO0, DiO0, RotC0, RetC0, DiC0,
                                                                                        TP0, TS0, RP0, RS0, ERaT0,
@@ -901,6 +888,8 @@
         print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}, {5:1.0f}".format("DT,TT,DR,TR,Y   :", DiT, TiT, DiR, TiR, Y))
         print("{0:12}".format(" --- Combined PBS + Pol.-filter ---"))
         print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}".format("DT,TT,DR,TR     :", DTa0, TTa0, DRa0, TRa0))
+        print("{0:26}: {1:6.3f}± {2:5.3f}/{3:2d}".format("LDRCal during calibration in calibration range", LDRCal0,
+                                                        dLDRCal, nLDRCal))
         print()
         print("Rotation Error Epsilon For Normal Measurements = ", RotationErrorEpsilonForNormalMeasurements)
         print(Type[TypeC], Loc[LocC])
@@ -917,8 +906,6 @@
         # print("{0:6.3f}±{1:5.3f}/{2:2d}|{3:8.5f}->{4:8.5f}->{5:8.5f}".format(LDRCal0, dLDRCal, nLDRCal, LDRtrue, LDRsimx, LDRCorr))
         # print("{0:8}       |{1:8}->{2:8}->{3:8}".format(" LDRCal"," LDRtrue", " LDRsimx", " LDRcorr"))
         # print(" --- LDRCal during calibration")
-        print("{0:26}: {1:6.3f}±{2:5.3f}/{3:2d}".format("LDRCal during calibration in calibration range", LDRCal0,
-                                                        dLDRCal, nLDRCal))
 
         # print("{0:8}={1:8.5f};{2:8}={3:8.5f}".format(" IinP",IinP," F11sim",F11sim))
         print()
@@ -926,6 +913,8 @@
         K0List = np.zeros(3)
         LDRsimxList = np.zeros(3)
         LDRCalList = 0.004, 0.2, 0.45
+        # The loop over LDRCalList is ony for checking whether and how much the LDR depends on the LDRCal during calibration and whether the corrections work.
+        # Still with assumed true parameters in input file
         for i, LDRCal in enumerate(LDRCalList):
             GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0,
                                                                                                    DiE0, RotO0, RetO0,
@@ -943,10 +932,18 @@
         print("{0:8.5f},{1:8.5f},{2:8.5f},{3:8.5f},{4:9.5f},{5:9.5f},{6:9.5f}".format(GR0, GT0, HR0, HT0, K0List[0],
                                                                                       K0List[1], K0List[2]))
         print('========================================================================')
+        print("{0:9},{1:9},{2:9}".format("  LDRtrue", "  LDRsimx", "  LDRCorr"))
 
-        print("{0:9},{1:9},{2:9}".format("  LDRtrue", "  LDRsimx", "  LDRCorr"))
-        LDRtrueList = 0.004, 0.02, 0.2, 0.45
-        for i, LDRtrue in enumerate(LDRtrueList):
+        #LDRtrueList = 0.004, 0.02, 0.2, 0.45
+        aF11sim0 = np.zeros(5)
+        LDRrange = np.zeros(5)
+        LDRrange = 0.004, 0.02, 0.1, 0.3, 0.45
+
+        # The loop over LDRtrueList is ony for checking how much the uncorrected LDRsimx deviates from LDRtrue ... and whether the corrections work.
+        # LDRsimx = LDRsim = Ir / It    or      1/LDRsim
+        # Still with assumed true parameters in input file
+        for i, LDRtrue in enumerate(LDRrange):
+        #for LDRtrue in LDRrange:
             GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0,
                                                                                                    DiE0, RotO0, RetO0,
                                                                                                    DiO0, RotC0, RetC0,
@@ -955,6 +952,8 @@
                                                                                                    RetT0, ERaR0, RotaR0,
                                                                                                    RetR0, LDRCal0)
             print("{0:9.5f},{1:9.5f},{2:9.5f}".format(LDRtrue, LDRsimx, LDRCorr))
+            aF11sim0[i] = F11sim0
+            # the assumed true aF11sim0 results will be used below to calc the deviation from the real signals
 
 file = open('output_files\output_' + LID + '.dat', 'r')
 print(file.read())
@@ -974,14 +973,13 @@
     sys.stdout = old_target
 '''
 if (Error_Calc):
-    # --- CALC again truth with LDRCal0 to reset all 0-values
+    # --- CALC again assumed truth with LDRCal0 and with assumed true parameters in input file to reset all 0-values
     GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0,
                                                                                            RotO0, RetO0, DiO0, RotC0,
                                                                                            RetC0, DiC0, TP0, TS0, RP0,
                                                                                            RS0, ERaT0, RotaT0, RetT0,
                                                                                            ERaR0, RotaR0, RetR0,
                                                                                            LDRCal0)
-
     # --- Start Errors calculation with variable parameters ------------------------------------------------------------------
 
     iN = -1
@@ -1008,8 +1006,8 @@
     F11min = np.zeros(5)
     F11max = np.zeros(5)
 
-    LDRrange = np.zeros(5)
-    LDRrange = 0.004, 0.02, 0.1, 0.3, 0.45
+    #LDRrange = np.zeros(5)
+    #LDRrange = 0.004, 0.02, 0.1, 0.3, 0.45
     # aLDRsimx = np.zeros(N)
     # aLDRsimx2 = np.zeros(N)
     # aLDRcorr = np.zeros(N)
@@ -1043,21 +1041,22 @@
     dtime = clock()
 
     # --- Calc Error signals
-    # GT, HT, GR, HR, K, Eta, LDRsim = Calc(RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC, TP, TS)
     # ---- Do the calculations of bra-ket vectors
     h = -1. if TypeC == 2 else 1
 
+    '''
     # from input file: measured LDRm and true LDRtrue, LDRtrue2  =>
     ameas = (1. - LDRmeas) / (1 + LDRmeas)
     atrue = (1. - LDRtrue) / (1 + LDRtrue)
     atrue2 = (1. - LDRtrue2) / (1 + LDRtrue2)
+    '''
 
     for iLDRCal in range(-nLDRCal, nLDRCal + 1):
         # from input file:  assumed LDRCal for calibration measurements
         LDRCal = LDRCal0
         if nLDRCal > 0: LDRCal = LDRCal0 + iLDRCal * dLDRCal / nLDRCal
 
-        GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0,
+        GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim = Calc(RotL0, RotE0, RetE0,
                                                                                                DiE0, RotO0, RetO0, DiO0,
                                                                                                RotC0, RetC0, DiC0, TP0,
                                                                                                TS0, RP0, RS0, ERaT0,
@@ -1673,14 +1672,19 @@
                         LDRsimx2 = 1./LDRsim2
                     '''
                     # ----- Backward correction
-                    # Corrected LDRCorr from forward simulated LDRsim (atrue) with assumed true G0,H0,K0
+                    # Corrected LDRCorr  with assumed true G0,H0,K0 from forward simulated (real) LDRsim (atrue)
                     LDRCorr = (LDRsim * K0 / Etax * (GT0 + HT0) - (GR0 + HR0)) / (
                     (GR0 - HR0) - LDRsim * K0 / Etax * (GT0 - HT0))
-
+                    '''
                     # -- F11corr from It and Ir and calibration EtaX
                     Text1 = "!!! EXPERIMENTAL !!!  F11corr from It and Ir with calibration EtaX: x-axis: F11corr(LDRtrue) / F11corr(LDRtrue = 0.004) - 1"
                     F11corr = 1 / (TiO * TiE) * (
-                    (HR0 * Etax / K0 * It / TTa - HT0 * Ir / TRa) / (HR0 * GT0 - HT0 * GR0))  # IL = 1  Eq.(64)
+                    (HR0 * Etax / K0 * It / TTa - HT0 * Ir / TRa) / (HR0 * GT0 - HT0 * GR0))  # IL = 1  Eq.(64); Etax/K0 = Eta0.
+                    '''
+                    # Corrected F11corr  with assumed true G0,H0,K0 from forward simulated (real) It and Ir (atrue)
+                    Text1 = "!!! EXPERIMENTAL !!!  F11corr from real It and Ir with real calibration EtaX: x-axis: F11corr(LDRtrue) / aF11sim0(LDRtrue) - 1"
+                    F11corr = 1 / (TiO * TiE) * (
+                    (HR0 * Etax / K0 * It / TTa - HT0 * Ir / TRa) / (HR0 * GT0 - HT0 * GR0))  # IL = 1  Eq.(64); Etax/K0 = Eta0.
 
                     # Text1 = "F11corr from It and Ir without corrections but with calibration EtaX: x-axis: F11corr(LDRtrue) devided by F11corr(LDRtrue = 0.004)"
                     # F11corr = 0.5/(TiO*TiE)*(Etax*It/TTa+Ir/TRa)    # IL = 1  Eq.(64)
@@ -1806,9 +1810,10 @@
 
     plt.show()
     plt.close
+
     '''
+    # --- Plot F11 histograms
     print()
-    #print("IT(LDRtrue) devided by IT(LDRtrue = 0.004)")
     print(" ############################################################################## ")
     print(Text1)
     print()
@@ -1816,8 +1821,8 @@
     iLDR = 5
     for LDRTrue in LDRrange:
         iLDR = iLDR - 1
-        aF11corr[iLDR,:] = aF11corr[iLDR,:] / aF11corr[0,:] - 1.0
-
+        #aF11corr[iLDR,:] = aF11corr[iLDR,:] / aF11corr[0,:] - 1.0
+        aF11corr[iLDR,:] = aF11corr[iLDR,:] / aF11sim0[iLDR] - 1.0
     # Plot F11
     def PlotSubHistF11(aVar, aX, X0, daX, iaX, naX):
         fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2))
@@ -1882,8 +1887,8 @@
 
     plt.show()
     plt.close
+
     '''
-
     '''
     # only histogram
     #print("******************* " + aVar + " *******************")
@@ -1903,6 +1908,7 @@
         plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2)
     plt.show()
     plt.close
+     # --- End of Plot F11 histograms
     '''
 
     # --- Plot LDRmin, LDRmax
--- a/output_files/output_example lidar.dat	Mon Nov 28 16:45:48 2016 +0200
+++ b/output_files/output_example lidar.dat	Tue Jan 24 02:26:19 2017 +0100
@@ -17,12 +17,12 @@
 DT,TT,DR,TR,Y   : 0.9588, 0.4850, -0.9029, 0.5150, 1
  --- Combined PBS + Pol.-filter ---
 DT,TT,DR,TR     : 1.0000, 0.2427, -0.9999, 0.2578
+LDRCal during calibration in calibration range:  0.009± 0.005/ 1
 
 Rotation Error Epsilon For Normal Measurements =  False
 linear polarizer before receiver
 Parallel signal detected in transmitted channel
 RS_RP_depend_on_TS_TP =  False
-LDRCal during calibration in calibration range:  0.009±0.005/ 1
 
 ========================================================================
  GR     , GT     , HR     , HT     ,  K(0.004), K(0.2) ,  K(0.45)
@@ -31,5 +31,6 @@
   LDRtrue,  LDRsimx,  LDRCorr
   0.00400,  0.00418,  0.00400
   0.02000,  0.02068,  0.02000
-  0.20000,  0.20637,  0.20000
+  0.10000,  0.10321,  0.10000
+  0.30000,  0.30952,  0.30000
   0.45000,  0.46426,  0.45000
--- a/system_settings/optic_input_example_lidar.py	Mon Nov 28 16:45:48 2016 +0200
+++ b/system_settings/optic_input_example_lidar.py	Tue Jan 24 02:26:19 2017 +0100
@@ -110,7 +110,7 @@
 
 # --- LDRCal assumed atmospheric linear depolarization ratio during the calibration measurements in calibration range with almost clean air (first guess)
 LDRCal,dLDRCal,nLDRCal= 0.009, 0.005, 1     # spans the interference filter influence 
-
+ 
 # ====================================================
 # NOTE: there is no need to change anything below.
 

mercurial