Mercurial > public > atmospheric_lidar_ghk / changeset
changeset
Moved "from __future__ import print_function"
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.
--- 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.
