Mercurial > public > atmospheric_lidar_ghk / changeset
changeset
Bug fixes and further options. Read Improvements_of_lidar_correction_ghk_ver.0.9.8_190124.pdf
Thu, 24 Jan 2019 21:24:25 +0100
- author
- Volker Freudenthaler <volker.freudenthaler@lmu.de>
- date
- Thu, 24 Jan 2019 21:24:25 +0100
- changeset 28
- 79fa4a41420f
- parent 27
- e7b3f4bf631d
- child 29
- 785232e69af7
Bug fixes and further options. Read Improvements_of_lidar_correction_ghk_ver.0.9.8_190124.pdf
--- a/.hgignore Fri Nov 24 22:54:15 2017 +0100 +++ b/.hgignore Thu Jan 24 21:24:25 2019 +0100 @@ -2,3 +2,5 @@ *~ re:^\.idea/ *.orig +re:^lidar_correction_ghk_0\.9\.8_published\.py$ +re:^lidar_correction_ghk_0\.9\.6\.py$
--- a/lidar_correction_ghk.py Fri Nov 24 22:54:15 2017 +0100 +++ b/lidar_correction_ghk.py Thu Jan 24 21:24:25 2019 +0100 @@ -1,6 +1,6 @@ # -*- coding: utf-8 -*- """ -Copyright 2016, 2017 Volker Freudenthaler +Copyright 2016, 2019 Volker Freudenthaler Licensed under the EUPL, Version 1.1 only (the "Licence"). @@ -17,7 +17,7 @@ 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 +Python 3.7, seaborn 0.9.0 Code description: @@ -38,9 +38,13 @@ 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. + +Ver. 0.9.8: - for details, see "Improvements_of_lidar_correction_ghk_ver.0.9.8_190124.pdf" + - correct calculation of Eta for cleaned anaylsers considering the combined transmission Eta = (TaT* TiT)(1 + cos2RotaT * DaT * DiT) and (TaR * TiR)(1 + cos2RotaR * DaR * DiR) according to the papers supplement Eqs. (S.10.10.1) ff + - ND-filters can be added for the calibration measurements in the transmitted (TCalT) and the reflected path (TCalR) in order to include their uncertainties in the error calculation. + - includes the simulation of signal noise """ - -# 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. +# Comment: The code might 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 @@ -58,7 +62,8 @@ sns_loaded = False import matplotlib.pyplot as plt -from time import clock +# from time import clock # python 2 +from timeit import default_timer as clock # from matplotlib.backends.backend_pdf import PdfPages # pdffile = '{}.pdf'.format('path') @@ -113,12 +118,14 @@ # ---- 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. +ScriptVersion = "0.9.8d" Error_Calc = True LID = "internal" EID = "internal" # --- IL Laser IL and +-Uncertainty DOLP, dDOLP, nDOLP = 0.995, 0.005, 1 # degree of linear polarization; default 1 RotL, dRotL, nRotL = 0.0, 0.0, 1 # alpha; rotation of laser polarization in degrees; default 0 +# IL = 1e5 #photons in the laser beam, including detection efficiency of the telescope, atmodspheric and r^2 attenuation # --- ME Emitter and +-Uncertainty DiE, dDiE, nDiE = 0., 0.00, 1 # Diattenuation TiE = 1. # Unpolarized transmittance @@ -157,12 +164,33 @@ DaR = (1 - ERaR) / (1 + ERaR) TaR = 0.5 * (1 + ERaR) -# Parellel signal detected in the transmitted channel => Y = 1, or in the reflected channel => Y = -1 -Y = -1. +# +++ Orientation of the PBS with respect to the reference plane (see Polarisation-orientation.png and Polarisation-orientation-2.png in /system_settings) +# Y = +1: polarisation in reference plane is finally transmitted, +# Y = -1: polarisation in reference plane is finally reflected. +Y = 1. # Calibrator = type defined by matrix values LocC = 4 # location of calibrator: behind laser = 1; behind emitter = 2; before receiver = 3; before PBS = 4 +# --- Additional attenuation (transmission of the ND-filter) during the calibration +TCalT, dTCalT, nTCalT = 1, 0, 0 # transmitting path; error calc not working yet +TCalR, dTCalR, nTCalR = 1, 0, 0 # reflecting path; error calc not working yet + +# *** signal noise error calculation +# --- number of photon counts in the signal summed up in the calibration range during the calibration measurements +NCalT = 1e6 # default 1e6, assumed the same in +45° and -45° signals +NCalR = 1e6 # default 1e6, assumed the same in +45° and -45° signals +NILfac = 200 # duration of standard (0°) measurement relative to calibration measurements +nNCal = 0 # error nNCal: one-sigma in steps to left and right for calibration signals +nNI = 0 # error nNI: one-sigma in steps to left and right for 0° signals +IoutTp0, IoutTp, dIoutTp0 = 0.5, 0.5, 0 +IoutTm0, IoutTm, dIoutTm0 = 0.5, 0.5, 0 +IoutRp0, IoutRp, dIoutRp0 = 0.5, 0.5, 0 +IoutRm0, IoutRm, dIoutRm0 = 0.5, 0.5, 0 +It0, It, dIt0 = 1 , 1, 0 +Ir0, Ir, dTr0 = 1 , 1, 0 +CalcFrom0deg = True + TypeC = 3 # linear polarizer calibrator # example with extinction ratio 0.001 DiC, dDiC, nDiC = 1.0, 0., 0 # ideal 1.0 @@ -173,7 +201,12 @@ # Rotation error without calibrator: if False, then epsilon = 0 for normal measurements RotationErrorEpsilonForNormalMeasurements = True - +# BSR backscatter ratio +# BSR, dBSR, nBSR = 10, 0.05, 1 +BSR = np.zeros(5) +BSR = [1.1, 2, 5, 10, 50] +# theoretical molecular LDR LDRm +LDRm, dLDRm, nLDRm = 0.004, 0.001, 1 # LDRCal assumed atmospheric linear depolarization ratio during the calibration measurements (first guess) LDRCal0, dLDRCal, nLDRCal = 0.25, 0.04, 1 LDRCal = LDRCal0 @@ -182,7 +215,7 @@ # LDRtrue for simulation of measurement => LDRsim LDRtrue = 0.5 LDRtrue2 = 0.004 - +LDRunCorr = 1 # Initialize other values to 0 ER, nER, dER = 0.001, 0, 0.001 K = 0. @@ -198,22 +231,18 @@ Type = ['', 'mechanical rotator', 'hwp rotator', 'linear polarizer', 'qwp rotator', 'circular polarizer', 'real HWP +-22.5°'] dY = ['reflected channel', '', 'transmitted channel'] +bPlotEtax = False # end of initial definition of variables # ******************************************************************************************************************************* -# --- Read actual lidar system parameters from optic_input.py (should be in sub-directory 'system_settings') -InputFile = 'optic_input_crossed_mirrors_test.py' -InputFile = 'optic_input_crossed_mirrors_test_combined.py' -InputFile = 'optic_input_ver8c_POLIS_355_Mar2017.py' -# InputFile = 'optic_input_ver8c_POLIS_532_Mar2017.py' -InputFile = 'optic_input_example_lidar_0.9.5.py' -InputFile = 'optic_input_ver8c_PollyXT_532_Lacros.py' -InputFile = 'optic_input_ver8c_CUT_532_May2017.py' -InputFile = 'optic_input_ver8c_MUSA.py' -InputFile = 'optic_input_ver10_RALI_JA.py' -InputFile = 'optic_input_ver10_RALI_act.py' -InputFile = 'optic_input_ver8c-IPRAL-170331.py' +# --- Read actual lidar system parameters from optic_input.py (must be in the programs sub-directory 'system_settings') +#InputFile = 'optic_input_example_lidar_2.py' +#InputFile = 'optic_input_example_lidar_3.py' +#InputFile = 'optic_input_example_lidar_4.py' +#InputFile = 'optic_input_example_lidar_5.py' +InputFile = 'optic_input_example_lidar.py' + ''' print("From ", dname) print("Running ", fname) @@ -266,13 +295,16 @@ RotaR0, dRotaR, nRotaR = RotaR, dRotaR, nRotaR LDRCal0, dLDRCal, nLDRCal = LDRCal, dLDRCal, nLDRCal -# LDRCal0,dLDRCal,nLDRCal=LDRCal,dLDRCal,0 + +# BSR0, dBSR, nBSR = BSR, dBSR, nBSR +LDRm0, dLDRm, nLDRm = LDRm, dLDRm, nLDRm # ---------- End of manual parameter change RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC = RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0 TP, TS, RP, RS, ERaT, RotaT, RetT, ERaR, RotaR, RetR = TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0 LDRCal = LDRCal0 DTa0, TTa0, DRa0, TRa0, LDRsimx, LDRCorr = 0, 0, 0, 0, 0, 0 +TCalT0, TCalR0 = TCalT, TCalR TiT = 0.5 * (TP + TS) DiT = (TP - TS) / (TP + TS) @@ -281,10 +313,21 @@ DiR = (RP - RS) / (RP + RS) ZiR = (1. - DiR ** 2) ** 0.5 - +C2aT = np.cos(np.deg2rad(2 * RotaT)) +C2aR = np.cos(np.deg2rad(2 * RotaR)) +ATPT = (1 + C2aT * DaT * DiT) +ARPT = (1 + C2aR * DaR * DiR) +TTa = TiT * TaT * ATPT # unpolarized transmission +TRa = TiR * TaR * ARPT # unpolarized transmission +Eta0 = TRa / TTa +# --- this subroutine is for the calculation of the PLDR from LDR, BSR, and LDRm ----------------------------------------------------- +def CalcPLDR(LDR, BSR, LDRm): + PLDR = (BSR * (1. + LDRm) * LDR - LDRm * (1. + LDR)) / (BSR * (1. + LDRm) - (1. + LDR)) + return (PLDR) # --- this subroutine is for the calculation with certain fixed parameters ----------------------------------------------------- -def Calc(DOLP, RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC, TP, TS, RP, RS, ERaT, RotaT, RetT, ERaR, RotaR, RetR, - LDRCal): +def Calc(TCalT, TCalR, NCalT, NCalR, DOLP, RotL, RotE, RetE, DiE, RotO, RetO, DiO, + RotC, RetC, DiC, TP, TS, RP, RS, + ERaT, RotaT, RetT, ERaR, RotaR, RetR, LDRCal): # ---- Do the calculations of bra-ket vectors h = -1. if TypeC == 2 else 1 # from input file: assumed LDRCal for calibration measurements @@ -389,21 +432,25 @@ S2aR = np.sin(np.deg2rad(h * 2 * RotaR)) C2aR = np.cos(np.deg2rad(2 * RotaR)) - # Analyzer As before the PBS Eq. D.5 - ATP1 = (1 + C2aT * DaT * DiT) - ATP2 = Y * (DiT + C2aT * DaT) - ATP3 = Y * S2aT * DaT * ZiT * CosT - ATP4 = S2aT * DaT * ZiT * SinT + # Analyzer As before the PBS Eq. D.5; combined PBS and cleaning pol-filter + ATPT = (1 + C2aT * DaT * DiT) # unpolarized transmission correction + TTa = TiT * TaT * ATPT # unpolarized transmission + ATP1 = 1 + ATP2 = Y * (DiT + C2aT * DaT) / ATPT + ATP3 = Y * S2aT * DaT * ZiT * CosT / ATPT + ATP4 = S2aT * DaT * ZiT * SinT / ATPT ATP = np.array([ATP1, ATP2, ATP3, ATP4]) + DTa = ATP2 * Y - ARP1 = (1 + C2aR * DaR * DiR) - ARP2 = Y * (DiR + C2aR * DaR) - ARP3 = Y * S2aR * DaR * ZiR * CosR - ARP4 = S2aR * DaR * ZiR * SinR + ARPT = (1 + C2aR * DaR * DiR) # unpolarized transmission correction + TRa = TiR * TaR * ARPT # unpolarized transmission + ARP1 = 1 + ARP2 = Y * (DiR + C2aR * DaR) / ARPT + ARP3 = Y * S2aR * DaR * ZiR * CosR / ARPT + ARP4 = S2aR * DaR * ZiR * SinR / ARPT ARP = np.array([ARP1, ARP2, ARP3, ARP4]) + DRa = ARP2 * Y - DTa = ATP2 * Y / ATP1 - DRa = ARP2 * Y / ARP1 # ---- Calculate signals and correction parameters for diffeent locations and calibrators if LocC == 4: # Calibrator before the PBS @@ -808,18 +855,17 @@ print("Calibrator location not implemented yet") sys.exit() - # Determination of the correction K of the calibration factor - IoutTp = TaT * TiT * TiO * TiE * (AT + BT) - IoutTm = TaT * TiT * TiO * TiE * (AT - BT) - IoutRp = TaR * TiR * TiO * TiE * (AR + BR) - IoutRm = TaR * TiR * TiO * TiE * (AR - BR) - + # Determination of the correction K of the calibration factor. + IoutTp = TTa * TiC * TiO * TiE * (AT + BT) + IoutTm = TTa * TiC * TiO * TiE * (AT - BT) + IoutRp = TRa * TiC * TiO * TiE * (AR + BR) + IoutRm = TRa * TiC * TiO * TiE * (AR - BR) # --- Results and Corrections; electronic etaR and etaT are assumed to be 1 Etapx = IoutRp / IoutTp Etamx = IoutRm / IoutTm Etax = (Etapx * Etamx) ** 0.5 - Eta = (TaR * TiR) / (TaT * TiT) # Eta = Eta*/K Eq. 84 + Eta = (TRa / TTa) # = TRa / TTa; Eta = Eta*/K Eq. 84 => K = Eta* / Eta; equation corrected according to the papers supplement Eqs. (S.10.10.1) ff K = Etax / Eta # For comparison with Volkers Libreoffice Müller Matrix spreadsheet @@ -827,10 +873,33 @@ # Eta_test_m = (IoutRm/IoutTm) # Eta_test = (Eta_test_p*Eta_test_m)**0.5 - # ----- Forward simulated signals and LDRsim with atrue; from input file - It = TaT * TiT * TiO * TiE * (GT + atrue * HT) - Ir = TaR * TiR * TiO * TiE * (GR + atrue * HR) - # LDRsim = 1/Eta*Ir/It # simulated LDR* with Y from input file + # ----- random error calculation ---------- + # noise must be calculated with the photon counts of measured signals; + # relative standard deviation of calibration signals with LDRcal; assumed to be statisitcally independent + # normalised noise errors + if (CalcFrom0deg): + dIoutTp = (NCalT * IoutTp) ** -0.5 + dIoutTm = (NCalT * IoutTm) ** -0.5 + dIoutRp = (NCalR * IoutRp) ** -0.5 + dIoutRm = (NCalR * IoutRm) ** -0.5 + else: + dIoutTp = (NCalT ** -0.5) + dIoutTm = (NCalT ** -0.5) + dIoutRp = (NCalR ** -0.5) + dIoutRm = (NCalR ** -0.5) + # Forward simulated 0°-signals with LDRCal with atrue; from input file + + It = TTa * TiO * TiE * (GT + atrue * HT) + Ir = TRa * TiO * TiE * (GR + atrue * HR) + # relative standard deviation of standard signals with LDRmeas; assumed to be statisitcally independent + if (CalcFrom0deg): + dIt = ((NCalT * It / IoutTp * NILfac / TCalT) ** -0.5) + dIr = ((NCalR * Ir / IoutRp * NILfac / TCalR) ** -0.5) + else: + dIt = ((NCalT * 2 * NILfac / TCalT ) ** -0.5) * It + dIr = ((NCalR * 2 * NILfac / TCalR) ** -0.5) * Ir + + # ----- Forward simulated LDRsim = 1/Eta*Ir/It # simulated LDR* with Y from input file LDRsim = Ir / It # simulated uncorrected LDR with Y from input file # Corrected LDRsimCorr from forward simulated LDRsim (atrue) # LDRsimCorr = (1./Eta*LDRsim*(GT+HT)-(GR+HR))/((GR-HR)-1./Eta*LDRsim*(GT-HT)) @@ -839,45 +908,48 @@ LDRsimx = 1. / LDRsim / Etax else: LDRsimx = LDRsim / Etax - ''' + ''' LDRsimx = LDRsim # The following is correct without doubt - # LDRCorr = (LDRsim*K/Etax*(GT+HT)-(GR+HR))/((GR-HR)-LDRsim*K/Etax*(GT-HT)) + # LDRCorr = (LDRsim/(Etax/K)*(GT+HT)-(GR+HR))/((GR-HR)-LDRsim/(Etax/K)*(GT-HT)) # The following is a test whether the equations for calibration Etax and normal signal (GHK, LDRsim) are consistent - LDRCorr = (LDRsim / Eta * (GT + HT) - (GR + HR)) / ((GR - HR) - LDRsim * K / Etax * (GT - HT)) - #LDRCorr = LDRsimx # for test only - TTa = TiT * TaT # *ATP1 - TRa = TiR * TaR # *ARP1 + LDRCorr = (LDRsim / (Etax / K) * (GT + HT) - (GR + HR)) / ((GR - HR) - LDRsim / (Etax / K) * (GT - HT)) + # here we could also use Eta instead of Etax / K => how to test whether Etax is correct? => comparison with MüllerMatrix simulation! + # Without any correction: only measured It, Ir, EtaX are used + LDRunCorr = (LDRsim / Etax * (GT / abs(GT) + HT / abs(HT)) - (GR / abs(GR) + HR / abs(HR))) / ((GR / abs(GR) - HR / abs(HR)) - LDRsim / Etax * (GT / abs(GT) - HT / abs(HT))) - F11sim = 1 / (TiO * TiE) * ( - (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 + #LDRCorr = LDRsimx # for test only + + F11sim = 1 / (TiO * TiE) * ((HR * Eta * It - HT * Ir) / (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) + return (IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr, + GT, HT, GR, HR, K, Eta, LDRsimx, LDRCorr, DTa, DRa, TTa, TRa, F11sim, LDRunCorr) + # ******************************************************************************************************************************* -# --- CALC truth with assumed true parameters from the input file -GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(DOLP0, RotL0, RotE0, RetE0, DiE0, RotO0, - RetO0, DiO0, RotC0, RetC0, DiC0, - TP0, TS0, RP0, RS0, ERaT0, - RotaT0, RetT0, ERaR0, RotaR0, - RetR0, LDRCal0) - +# --- CALC with assumed true parameters from the input file +IoutTp0, IoutTm0, IoutRp0, IoutRm0, It0, Ir0, dIoutTp0, dIoutTm0, dIoutRp0, dIoutRm0, dIt0, dIr0, \ +GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \ +Calc(TCalT, TCalR, NCalT, NCalR, DOLP0, RotL0, RotE0, RetE0, DiE0, + RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, + ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0) +Etax0 = K0 * Eta0 # --- Print parameters to console and output file -with open('output_files\output_' + LID + '.dat', 'w') as f: +with open('output_files\\' + LID + '-' + InputFile[0:-3] + '-GHK.dat', 'w') as f: with redirect_stdout(f): - print("From ", dname) - print("Running ", fname) + print("From folder", dname) + print("Running prog", fname) + print("Version", ScriptVersion) print("Reading input file ", InputFile) # , " for Lidar system :", EID, ", ", LID) print("for Lidar system: ", EID, ", ", LID) # --- Print iput information********************************* print(" --- Input parameters: value ±error / ±steps ----------------------") - print("{0:5}{1:5} {2:6.4f}±{3:7.4f}/{4:2d}; {5:8} {6:8.4f}±{7:7.4f}/{8:2d}".format("Laser: ", "DOLP =", DOLP0, dDOLP, nDOLP, - " Rotation alpha = ", RotL0, dRotL, - nRotL)) + print("{0:5}{1:5} {2:6.4f}±{3:7.4f}/{4:2d}; {5:8} {6:8.4f}±{7:7.4f}/{8:2d}".format( + "Laser: ", "DOLP =", DOLP0, dDOLP, nDOLP," Rotation alpha = ", RotL0, dRotL, nRotL)) print(" Diatt., Tunpol, Retard., Rotation (deg)") print("{0:12} {1:7.4f}±{2:7.4f}/{8:2d}, {3:7.4f}, {4:3.0f}±{5:3.0f}/{9:2d}, {6:7.4f}±{7:7.4f}/{10:2d}".format( "Emitter ", DiE0, dDiE, TiE, RetE0, dRetE, RotE0, dRotE, nDiE, nRetE, nRotE)) @@ -886,22 +958,25 @@ print("{0:12} {1:7.4f}±{2:7.4f}/{8:2d}, {3:7.4f}, {4:3.0f}±{5:3.0f}/{9:2d}, {6:7.4f}±{7:7.4f}/{10:2d}".format( "Calibrator ", DiC0, dDiC, TiC, RetC0, dRetC, RotC0, dRotC, nDiC, nRetC, nRotC)) print("{0:12}".format(" --- Pol.-filter ---")) - print( - "{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format("ERT, RotT :", ERaT0, dERaT, nERaT, - RotaT0, dRotaT, nRotaT)) - print( - "{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format("ERR, RotR :", ERaR0, dERaR, nERaR, - RotaR0, dRotaR, nRotaR)) + print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format( + "ERT, RotT :", ERaT0, dERaT, nERaT, RotaT0, dRotaT, nRotaT)) + print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format( + "ERR, RotR :", ERaR0, dERaR, nERaR, RotaR0, dRotaR, nRotaR)) print("{0:12}".format(" --- PBS ---")) - print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format("TP,TS :", TP0, dTP, nTP, TS0, - dTS, nTS)) - print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format("RP,RS :", RP0, dRP, nRP, RS0, - dRS, nRS)) - 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}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format( + "TP,TS :", TP0, dTP, nTP, TS0, dTS, nTS)) + print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format( + "RP,RS :", RP0, dRP, nRP, RS0, dRS, nRS)) + 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("{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("{0:12}".format(" --- Additional ND filter attenuation (transmission) during the calibration ---")) + print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format( + "TCalT,TCalR :", TCalT0, dTCalT, nTCalT, TCalR0, dTCalR, nTCalR)) print() print("Rotation Error Epsilon For Normal Measurements = ", RotationErrorEpsilonForNormalMeasurements) print(Type[TypeC], Loc[LocC]) @@ -922,53 +997,85 @@ # print("{0:8}={1:8.5f};{2:8}={3:8.5f}".format(" IinP",IinP," F11sim",F11sim)) print() - K0List = np.zeros(3) - LDRsimxList = np.zeros(3) - LDRCalList = 0.004, 0.2, 0.45 + K0List = np.zeros(6) + LDRsimxList = np.zeros(6) + LDRCalList = 0.004, 0.02, 0.1, 0.2, 0.3, 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 + + facIt = NCalT / TCalT0 * NILfac + facIr = NCalR / TCalR0 * NILfac + print("IoutTp, IoutTm, IoutRp, IoutRm, It , Ir , dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr") + for i, LDRCal in enumerate(LDRCalList): - GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(DOLP0, RotL0, RotE0, RetE0, - DiE0, RotO0, RetO0, - DiO0, RotC0, RetC0, - DiC0, TP0, TS0, RP0, - RS0, ERaT0, RotaT0, - RetT0, ERaR0, RotaR0, - RetR0, LDRCal) + IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr, \ + GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \ + Calc(TCalT0, TCalR0, NCalT, NCalR, DOLP0, RotL0, RotE0, RetE0, DiE0, + RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, + ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal) K0List[i] = K0 LDRsimxList[i] = LDRsimx - print('========================================================================') - print("{0:8},{1:8},{2:8},{3:8},{4:9},{5:8},{6:9}".format(" GR", " GT", " HR", " HT", " K(0.004)", " K(0.2)", - " K(0.45)")) - 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:10},{1:10},{2:10},{3:10}".format(" LDRtrue", " LDRsimx", " 1/LDRsimx", " LDRCorr")) + # check error signals + # print( "{:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}" + # .format(IoutTp * NCalT, IoutTm * NCalT, IoutRp * NCalR, IoutRm * NCalR, It * facIt, Ir * facIr, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr)) + #print( "{:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}".format(IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr)) + # end check error signals + print('===========================================================================================================') + print("{0:8},{1:8},{2:8},{3:8},{4:9},{5:8},{6:9},{7:9},{8:9},{9:9}".format( + " GR", " GT", " HR", " HT", " K(0.004)", " K(0.02)", " K(0.1)", " K(0.2)", " K(0.3)", " K(0.45)")) + print("{0:8.5f},{1:8.5f},{2:8.5f},{3:8.5f},{4:9.5f},{5:9.5f},{6:9.5f},{7:9.5f},{8:9.5f},{9:9.5f}".format( + GR0, GT0, HR0, HT0, K0List[0], K0List[1], K0List[2], K0List[3], K0List[4], K0List[5])) + print('===========================================================================================================') + print() + print("Errors from neglecting GHK corrections and/or calibration:") + print("{0:>10},{1:>10},{2:>10},{3:>10},{4:>10},{5:>10}".format( + "LDRtrue", "LDRunCorr", "1/LDRunCorr", "LDRsimx", "1/LDRsimx", "LDRCorr")) #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 + LDRrange = [0.004, 0.02, 0.1, 0.3, 0.45] # list - # The loop over LDRtrueList is ony for checking how much the uncorrected LDRsimx deviates from LDRtrue ... and whether the corrections work. + # The loop over LDRtrueList is only 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(DOLP0, RotL0, RotE0, RetE0, - DiE0, RotO0, RetO0, - DiO0, RotC0, RetC0, - DiC0, TP0, TS0, RP0, - RS0, ERaT0, RotaT0, - RetT0, ERaR0, RotaR0, - RetR0, LDRCal0) - print("{0:10.5f},{1:10.5f},{2:10.5f},{3:10.5f}".format(LDRtrue, LDRsimx, 1/LDRsimx, LDRCorr)) + IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr, \ + GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \ + Calc(TCalT0, TCalR0, NCalT, NCalR, DOLP0, RotL0, RotE0, RetE0, DiE0, + RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, + ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0) + print("{0:10.5f},{1:10.5f},{2:10.5f},{3:10.5f},{4:10.5f},{5:10.5f}".format(LDRtrue, LDRunCorr, 1/LDRunCorr, LDRsimx, 1/LDRsimx, LDRCorr)) aF11sim0[i] = F11sim0 # the assumed true aF11sim0 results will be used below to calc the deviation from the real signals - print("Note: LDRsimx = LDR of the nominal system directly from measured signals without GHK-corrections") + print("LDRsimx = LDR of the nominal system directly from measured signals without calibration and GHK-corrections") + print("LDRunCorr = LDR of the nominal system directly from measured signals with calibration but without GHK-corrections; electronic amplifications = 1 assumed") + print("LDRCorr = LDR calibrated and GHK-corrected") + print() + print("Errors from signal noise:") + print("Signal counts: NCalT, NCalR, NILfac, nNCal, nNI = {0:10.0f},{1:10.0f},{2:3.0f},{3:2.0f},{4:2.0f}".format( + NCalT, NCalR, NILfac, nNCal, nNI)) -file = open('output_files\output_' + LID + '.dat', 'r') + '''# das muß wieder weg + print("IoutTp, IoutTm, IoutRp, IoutRm, It , Ir , dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr") + LDRCal = 0.01 + for i, LDRtrue in enumerate(LDRrange): + IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr, \ + GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \ + Calc(TCalT0, TCalR0, NCalT, NCalR, DOLP0, RotL0, RotE0, RetE0, DiE0, + RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, + ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0) + print( "{:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}".format( + IoutTp * NCalT, IoutTm * NCalT, IoutRp * NCalR, IoutRm * NCalR, It * facIt, Ir * facIr, + dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr)) + aF11sim0[i] = F11sim0 + # the assumed true aF11sim0 results will be used below to calc the deviation from the real signals + # bis hierher weg + ''' + +file = open('output_files\\' + LID + '-' + InputFile[0:-3] + '-GHK.dat', 'r') print(file.read()) file.close() @@ -987,22 +1094,26 @@ ''' if (Error_Calc): # --- 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(DOLP0, 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 ------------------------------------------------------------------ + IoutTp0, IoutTm0, IoutRp0, IoutRm0, It0, Ir0, dIoutTp0, dIoutTm0, dIoutRp0, dIoutRm0, dIt0, dIr0, \ + GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \ + Calc(TCalT0, TCalR0, NCalT, NCalR, DOLP0, RotL0, RotE0, RetE0, DiE0, + RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, + ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0) + Etax0 = K0 * Eta0 + # --- Start Error calculation with variable parameters ------------------------------------------------------------------ + # error nNCal: one-sigma in steps to left and right for calibration signals + # error nNI: one-sigma in steps to left and right for 0° signals iN = -1 - N = ((nDOLP * 2 + 1) * (nRotL * 2 + 1) * + N = ((nTCalT * 2 + 1) * (nTCalR * 2 + 1) * + (nNCal * 2 + 1) ** 4 * (nNI * 2 + 1) ** 2 * + (nDOLP * 2 + 1) * (nRotL * 2 + 1) * (nRotE * 2 + 1) * (nRetE * 2 + 1) * (nDiE * 2 + 1) * (nRotO * 2 + 1) * (nRetO * 2 + 1) * (nDiO * 2 + 1) * (nRotC * 2 + 1) * (nRetC * 2 + 1) * (nDiC * 2 + 1) * (nTP * 2 + 1) * (nTS * 2 + 1) * (nRP * 2 + 1) * (nRS * 2 + 1) * (nERaT * 2 + 1) * (nERaR * 2 + 1) * (nRotaT * 2 + 1) * (nRotaR * 2 + 1) * (nRetT * 2 + 1) * (nRetR * 2 + 1) * (nLDRCal * 2 + 1)) - print("N = ", N, " ", end="") + print("number of system variations N = ", N, " ", end="") if N > 1e6: if user_yes_no_query('Warning: processing ' + str( @@ -1018,9 +1129,11 @@ LDRmax = np.zeros(5) F11min = np.zeros(5) F11max = np.zeros(5) + Etaxmin = np.zeros(5) + Etaxmax = 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) @@ -1047,9 +1160,23 @@ aRetO = np.zeros(N) aRotO = np.zeros(N) aLDRCal = np.zeros(N) - aA = np.zeros((5, N)) - aX = np.zeros((5, N)) + aNCalTp = np.zeros(N) + aNCalTm = np.zeros(N) + aNCalRp = np.zeros(N) + aNCalRm = np.zeros(N) + aNIt = np.zeros(N) + aNIr = np.zeros(N) + aTCalT = np.zeros(N) + aTCalR = np.zeros(N) + + # each np.zeros((LDRrange, N)) array has the same N-dependency + aLDRcorr = np.zeros((5, N)) aF11corr = np.zeros((5, N)) + aPLDR = np.zeros((5, N)) + aEtax = np.zeros((5, N)) + + # np.zeros((GHKs, N)) + aGHK = np.zeros((5, N)) atime = clock() dtime = clock() @@ -1066,16 +1193,19 @@ ''' for iLDRCal in range(-nLDRCal, nLDRCal + 1): - # from input file: assumed LDRCal for calibration measurements + # from input file: 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, F11sim = Calc(DOLP0, RotL0, RotE0, RetE0, - DiE0, RotO0, RetO0, DiO0, - RotC0, RetC0, DiC0, TP0, - TS0, RP0, RS0, ERaT0, - RotaT0, RetT0, ERaR0, - RotaR0, RetR0, LDRCal) + if nLDRCal > 0: + LDRCal = LDRCal0 + iLDRCal * dLDRCal / nLDRCal + # provides the intensities of the calibration measurements at various LDRCal for signal noise errors + # IoutTp, IoutTm, IoutRp, IoutRm, dIoutTp, dIoutTm, dIoutRp, dIoutRm + ''' + IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr, \ + GT, HT, GR, HR, K, Eta, LDRsimx, LDRCorr, DTa, DRa, TTa, TRa, F11sim, LDRunCorr = \ + Calc(TCalT, TCalR, NCalT, NCalR, DOLP0, RotL0, RotE0, RetE0, DiE0, + RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, + ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal) + ''' aCal = (1. - LDRCal) / (1 + LDRCal) for iDOLP, iRotL, iRotE, iRetE, iDiE \ in [(iDOLP, iRotL, iRotE, iRetE, iDiE) @@ -1151,16 +1281,6 @@ for iRotaR in range(-nRotaR, nRotaR + 1) for iRetR in range(-nRetR, nRetR + 1)]: - iN = iN + 1 - if (iN == 10001): - ctime = clock() - print(" estimated time ", "{0:4.2f}".format(N / 10000 * (ctime - atime)), "sec ") # , end="") - print("\r elapsed time ", "{0:5.0f}".format((ctime - atime)), "sec ", end="\r") - ctime = clock() - if ((ctime - dtime) > 10): - print("\r elapsed time ", "{0:5.0f}".format((ctime - atime)), "sec ", end="\r") - dtime = ctime - if nRotO > 0: RotO = RotO0 + iRotO * dRotO / nRotO if nRetO > 0: RetO = RetO0 + iRetO * dRetO / nRetO if nDiO > 0: DiO = DiO0 + iDiO * dDiO / nDiO @@ -1209,6 +1329,7 @@ CosR = np.cos(np.deg2rad(RetR)) SinR = np.sin(np.deg2rad(RetR)) + # cleaning pol-filter DaT = (1 - ERaT) / (1 + ERaT) DaR = (1 - ERaR) / (1 + ERaR) TaT = 0.5 * (1 + ERaT) @@ -1219,21 +1340,24 @@ S2aR = np.sin(np.deg2rad(h * 2 * RotaR)) C2aR = np.cos(np.deg2rad(2 * RotaR)) - # Aanalyzer As before the PBS Eq. D.5 - ATP1 = (1 + C2aT * DaT * DiT) - ATP2 = Y * (DiT + C2aT * DaT) - ATP3 = Y * S2aT * DaT * ZiT * CosT - ATP4 = S2aT * DaT * ZiT * SinT + # Analyzer As before the PBS Eq. D.5; combined PBS and cleaning pol-filter + ATPT = (1 + C2aT * DaT * DiT) # unpolarized transmission correction + TTa = TiT * TaT * ATPT # unpolarized transmission + ATP1 = 1 + ATP2 = Y * (DiT + C2aT * DaT) / ATPT + ATP3 = Y * S2aT * DaT * ZiT * CosT / ATPT + ATP4 = S2aT * DaT * ZiT * SinT / ATPT ATP = np.array([ATP1, ATP2, ATP3, ATP4]) + DTa = ATP2 * Y - ARP1 = (1 + C2aR * DaR * DiR) - ARP2 = Y * (DiR + C2aR * DaR) - ARP3 = Y * S2aR * DaR * ZiR * CosR - ARP4 = S2aR * DaR * ZiR * SinR + ARPT = (1 + C2aR * DaR * DiR) # unpolarized transmission correction + TRa = TiR * TaR * ARPT # unpolarized transmission + ARP1 = 1 + ARP2 = Y * (DiR + C2aR * DaR) / ARPT + ARP3 = Y * S2aR * DaR * ZiR * CosR / ARPT + ARP4 = S2aR * DaR * ZiR * SinR / ARPT ARP = np.array([ARP1, ARP2, ARP3, ARP4]) - - TTa = TiT * TaT # *ATP1 - TRa = TiR * TaR # *ARP1 + DRa = ARP2 * Y # ---- Calculate signals and correction parameters for diffeent locations and calibrators if LocC == 4: # Calibrator before the PBS @@ -1644,150 +1768,237 @@ GR = ARE1o * D + ARE4o * C HT = ATE2a * A + ATE3a * B + ATE4a * C HR = ARE2a * A + ARE3a * B + ARE4a * C - else: print('Calibrator not implemented yet') sys.exit() - # Calibration signals with aCal => Determination of the correction K of the real calibration factor - IoutTp = TaT * TiT * TiO * TiE * (AT + BT) - IoutTm = TaT * TiT * TiO * TiE * (AT - BT) - IoutRp = TaR * TiR * TiO * TiE * (AR + BR) - IoutRm = TaR * TiR * TiO * TiE * (AR - BR) - # --- Results and Corrections; electronic etaR and etaT are assumed to be 1 - # Eta = TiR/TiT # Eta = Eta*/K Eq. 84 - Etapx = IoutRp / IoutTp - Etamx = IoutRm / IoutTm - Etax = (Etapx * Etamx) ** 0.5 - K = Etax / Eta0 - # print("{0:6.3f},{1:6.3f},{2:6.3f},{3:6.3f},{4:6.3f},{5:6.3f},{6:6.3f},{7:6.3f},{8:6.3f},{9:6.3f},{10:6.3f}".format(AT, BT, AR, BR, DiC, ZiC, RetO, TP, TS, Kp, Km)) - # print("{0:6.3f},{1:6.3f},{2:6.3f},{3:6.3f}".format(DiC, ZiC, Kp, Km)) + for iTCalT, iTCalR, iNCalTp, iNCalTm, iNCalRp, iNCalRm, iNIt, iNIr \ + in [ + (iTCalT, iTCalR, iNCalTp, iNCalTm, iNCalRp, iNCalRm, iNIt, iNIr) + for iTCalT in range(-nTCalT, nTCalT + 1) # Etax + for iTCalR in range(-nTCalR, nTCalR + 1) # Etax + for iNCalTp in range(-nNCal, nNCal + 1) # noise error of calibration signals => Etax + for iNCalTm in range(-nNCal, nNCal + 1) # noise error of calibration signals => Etax + for iNCalRp in range(-nNCal, nNCal + 1) # noise error of calibration signals => Etax + for iNCalRm in range(-nNCal, nNCal + 1) # noise error of calibration signals => Etax + for iNIt in range(-nNI, nNI + 1) + for iNIr in range(-nNI, nNI + 1)]: - # For comparison with Volkers Libreoffice Müller Matrix spreadsheet - # Eta_test_p = (IoutRp/IoutTp) - # Eta_test_m = (IoutRm/IoutTm) - # Eta_test = (Eta_test_p*Eta_test_m)**0.5 + # Calibration signals with aCal => Determination of the correction K of the real calibration factor + IoutTp = TTa * TiC * TiO * TiE * (AT + BT) + IoutTm = TTa * TiC * TiO * TiE * (AT - BT) + IoutRp = TRa * TiC * TiO * TiE * (AR + BR) + IoutRm = TRa * TiC * TiO * TiE * (AR - BR) - # ************************************************************************* - iLDR = -1 - for LDRTrue in LDRrange: - iLDR = iLDR + 1 - atrue = (1 - LDRTrue) / (1 + LDRTrue) - # ----- Forward simulated signals and LDRsim with atrue; from input file - It = TaT * TiT * TiO * TiE * (GT + atrue * HT) # TaT*TiT*TiC*TiO*IinL*(GT+atrue*HT) - Ir = TaR * TiR * TiO * TiE * (GR + atrue * HR) # TaR*TiR*TiC*TiO*IinL*(GR+atrue*HR) + if nTCalT > 0: TCalT = TCalT0 + iTCalT * dTCalT / nTCalT + if nTCalR > 0: TCalR = TCalR0 + iTCalR * dTCalR / nTCalR + # signal noise errors + # ----- random error calculation ---------- + # noise must be calculated from/with the actually measured signals; influence of TCalT, TCalR errors on nouse are not considered ? + # actually measured signals are in input file and don't change + # relative standard deviation of calibration signals with LDRcal; assumed to be statisitcally independent + # error nNCal: one-sigma in steps to left and right for calibration signals + if nNCal > 0: + if (CalcFrom0deg): + dIoutTp = (NCalT * IoutTp) ** -0.5 + dIoutTm = (NCalT * IoutTm) ** -0.5 + dIoutRp = (NCalR * IoutRp) ** -0.5 + dIoutRm = (NCalR * IoutRm) ** -0.5 + else: + dIoutTp = dIoutTp0 * (IoutTp / IoutTp0) + dIoutTm = dIoutTm0 * (IoutTm / IoutTm0) + dIoutRp = dIoutRp0 * (IoutRp / IoutRp0) + dIoutRm = dIoutRm0 * (IoutRm / IoutRm0) + # print(iTCalT, iTCalR, iNCalTp, iNCalTm, iNCalRp, iNCalRm, iNIt, iNIr, IoutTp, dIoutTp) + IoutTp = IoutTp * (1 + iNCalTp * dIoutTp / nNCal) + IoutTm = IoutTm * (1 + iNCalTm * dIoutTm / nNCal) + IoutRp = IoutRp * (1 + iNCalRp * dIoutRp / nNCal) + IoutRm = IoutRm * (1 + iNCalRm * dIoutRm / nNCal) - # LDRsim = 1/Eta*Ir/It # simulated LDR* with Y from input file - LDRsim = Ir / It # simulated uncorrected LDR with Y from input file + IoutTp = IoutTp * TCalT / TCalT0 + IoutTm = IoutTm * TCalT / TCalT0 + IoutRp = IoutRp * TCalR / TCalR0 + IoutRm = IoutRm * TCalR / TCalR0 + # --- Results and Corrections; electronic etaR and etaT are assumed to be 1 for true and assumed true systems + # calibration factor + Eta = (TRa / TTa) # = TRa / TTa; Eta = Eta*/K Eq. 84; corrected according to the papers supplement Eqs. (S.10.10.1) ff + # possibly real calibration factor + Etapx = IoutRp / IoutTp + Etamx = IoutRm / IoutTm + Etax = (Etapx * Etamx) ** 0.5 + K = Etax / Eta + # print("{0:6.3f},{1:6.3f},{2:6.3f},{3:6.3f},{4:6.3f},{5:6.3f},{6:6.3f},{7:6.3f},{8:6.3f},{9:6.3f},{10:6.3f}".format(AT, BT, AR, BR, DiC, ZiC, RetO, TP, TS, Kp, Km)) + # print("{0:6.3f},{1:6.3f},{2:6.3f},{3:6.3f}".format(DiC, ZiC, Kp, Km)) + + # For comparison with Volkers Libreoffice Müller Matrix spreadsheet + # Eta_test_p = (IoutRp/IoutTp) + # Eta_test_m = (IoutRm/IoutTm) + # Eta_test = (Eta_test_p*Eta_test_m)**0.5 ''' - if Y == 1.: - LDRsimx = LDRsim - LDRsimx2 = LDRsim2 - else: - LDRsimx = 1./LDRsim - LDRsimx2 = 1./LDRsim2 - ''' - # ----- Backward correction - # 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)) + for iIt, iIr \ + in [(iIt, iIr) + for iIt in range(-nNI, nNI + 1) + for iIr in range(-nNI, nNI + 1)]: ''' - # -- 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); 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. + + iN = iN + 1 + if (iN == 10001): + ctime = clock() + print(" estimated time ", "{0:4.2f}".format(N / 10000 * (ctime - atime)), "sec ") # , end="") + print("\r elapsed time ", "{0:5.0f}".format((ctime - atime)), "sec ", end="\r") + ctime = clock() + if ((ctime - dtime) > 10): + print("\r elapsed time ", "{0:5.0f}".format((ctime - atime)), "sec ", end="\r") + dtime = ctime - # 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) + # *** loop for different real LDRs ********************************************************************** + iLDR = -1 + for LDRTrue in LDRrange: + iLDR = iLDR + 1 + atrue = (1 - LDRTrue) / (1 + LDRTrue) + # ----- Forward simulated signals and LDRsim with atrue; from input file; not considering TiC. + It = TTa * TiO * TiE * (GT + atrue * HT) # TaT*TiT*TiC*TiO*IinL*(GT+atrue*HT) + Ir = TRa * TiO * TiE * (GR + atrue * HR) # TaR*TiR*TiC*TiO*IinL*(GR+atrue*HR) + # # signal noise errors; standard deviation of signals; assumed to be statisitcally independent + # because the signals depend on LDRtrue, the errors dIt and dIr must be calculated for each LDRtrue + if (CalcFrom0deg): + dIt = ((NCalT * It / IoutTp * NILfac / TCalT) ** -0.5) + dIr = ((NCalR * Ir / IoutRp * NILfac / TCalR) ** -0.5) + else: + dIt = ((NCalT * 2 * NILfac / TCalT ) ** -0.5) * It + dIr = ((NCalR * 2 * NILfac / TCalR) ** -0.5) * Ir + # error nNI: one-sigma in steps to left and right for 0° signals + if nNI > 0: + It = It * (1 + iNIt * dIt / nNI) + Ir = Ir * (1 + iNIr * dIr / nNI) - # -- It from It only with atrue without corrections - for BERTHA (and PollyXTs) - # Text1 = " x-axis: IT(LDRtrue) / IT(LDRtrue = 0.004) - 1" - # F11corr = It/(TaT*TiT*TiO*TiE) #/(TaT*TiT*TiO*TiE*(GT0+atrue*HT0)) - # !!! see below line 1673ff + # LDRsim = 1/Eta*Ir/It # simulated LDR* with Y from input file + LDRsim = Ir / It # simulated uncorrected LDR with Y from input file - aF11corr[iLDR, iN] = F11corr - aA[iLDR, iN] = LDRCorr # LDRCorr # LDRsim # for test only + # ----- Backward correction + # Corrected LDRCorr with assumed true G0,H0,K0,Eta0 from forward simulated (real) LDRsim(atrue) + LDRCorr = (LDRsim / (Etax / K0) * (GT0 + HT0) - (GR0 + HR0)) / ((GR0 - HR0) - LDRsim / (Etax / K0) * (GT0 - HT0)) + + # The following is a test whether the equations for calibration Etax and normal signal (GHK, LDRsim) are consistent + # LDRCorr = (LDRsim / Eta * (GT + HT) - (GR + HR)) / ((GR - HR) - LDRsim / Eta * (GT - HT)) + # Without any correction + LDRunCorr = (LDRsim / Etax * (GT / abs(GT) + HT / abs(HT)) - (GR / abs(GR) + HR / abs(HR))) / ((GR / abs(GR) - HR / abs(HR)) - LDRsim / Etax * (GT / abs(GT) - HT / abs(HT))) + - aX[0, iN] = GR - aX[1, iN] = GT - aX[2, iN] = HR - aX[3, iN] = HT - aX[4, iN] = K + ''' + # -- 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); 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) - aLDRCal[iN] = iLDRCal - aDOLP[iN] = iDOLP - aERaT[iN] = iERaT - aERaR[iN] = iERaR - aRotaT[iN] = iRotaT - aRotaR[iN] = iRotaR - aRetT[iN] = iRetT - aRetR[iN] = iRetR + # -- It from It only with atrue without corrections - for BERTHA (and PollyXTs) + # Text1 = " x-axis: IT(LDRtrue) / IT(LDRtrue = 0.004) - 1" + # F11corr = It/(TaT*TiT*TiO*TiE) #/(TaT*TiT*TiO*TiE*(GT0+atrue*HT0)) + # ! see below line 1673ff + + aF11corr[iLDR, iN] = F11corr + aLDRcorr[iLDR, iN] = LDRCorr # LDRCorr # LDRsim # for test only + # aPLDR[iLDR, iN] = CalcPLDR(LDRCorr, BSR[iLDR], LDRm0) + aEtax[iLDR, iN] = Etax + + aGHK[0, iN] = GR + aGHK[1, iN] = GT + aGHK[2, iN] = HR + aGHK[3, iN] = HT + aGHK[4, iN] = K - aRotL[iN] = iRotL - aRotE[iN] = iRotE - aRetE[iN] = iRetE - aRotO[iN] = iRotO - aRetO[iN] = iRetO - aRotC[iN] = iRotC - aRetC[iN] = iRetC - aDiO[iN] = iDiO - aDiE[iN] = iDiE - aDiC[iN] = iDiC - aTP[iN] = iTP - aTS[iN] = iTS - aRP[iN] = iRP - aRS[iN] = iRS + aLDRCal[iN] = iLDRCal + aDOLP[iN] = iDOLP + aERaT[iN] = iERaT + aERaR[iN] = iERaR + aRotaT[iN] = iRotaT + aRotaR[iN] = iRotaR + aRetT[iN] = iRetT + aRetR[iN] = iRetR + + aRotL[iN] = iRotL + aRotE[iN] = iRotE + aRetE[iN] = iRetE + aRotO[iN] = iRotO + aRetO[iN] = iRetO + aRotC[iN] = iRotC + aRetC[iN] = iRetC + aDiO[iN] = iDiO + aDiE[iN] = iDiE + aDiC[iN] = iDiC + aTP[iN] = iTP + aTS[iN] = iTS + aRP[iN] = iRP + aRS[iN] = iRS + aTCalT[iN] = iTCalT + aTCalR[iN] = iTCalR + + aNCalTp[iN] = iNCalTp # IoutTp, IoutTm, IoutRp, IoutRm => Etax + aNCalTm[iN] = iNCalTm # IoutTp, IoutTm, IoutRp, IoutRm => Etax + aNCalRp[iN] = iNCalRp # IoutTp, IoutTm, IoutRp, IoutRm => Etax + aNCalRm[iN] = iNCalRm # IoutTp, IoutTm, IoutRp, IoutRm => Etax + aNIt[iN] = iNIt # It, Tr + aNIr[iN] = iNIr # It, Tr # --- END loop btime = clock() - print("\r done in ", "{0:5.0f}".format(btime - atime), "sec") # , end="\r") - + # print("\r done in ", "{0:5.0f}".format(btime - atime), "sec. => producing plots now .... some more seconds ..."), # , end="\r"); + print(" done in ", "{0:5.0f}".format(btime - atime), "sec. => producing plots now .... some more seconds ...") # --- Plot ----------------------------------------------------------------- + print("Errors from GHK correction uncertainties:") if (sns_loaded): sns.set_style("whitegrid") - sns.set_palette("bright", 6) + sns.set_palette("bright6", 6) + # for older seaborn versions: + # sns.set_palette("bright", 6) ''' fig2 = plt.figure() - plt.plot(aA[2,:],'b.') - plt.plot(aA[3,:],'r.') - plt.plot(aA[4,:],'g.') - #plt.plot(aA[6,:],'c.') + plt.plot(aLDRcorr[2,:],'b.') + plt.plot(aLDRcorr[3,:],'r.') + plt.plot(aLDRcorr[4,:],'g.') + #plt.plot(aLDRcorr[6,:],'c.') plt.show ''' - # Plot LDR def PlotSubHist(aVar, aX, X0, daX, iaX, naX): + # aVar is the name of the parameter and aX is the subset of aLDRcorr which is coloured in the plot + # example: PlotSubHist("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP) fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) iLDR = -1 for LDRTrue in LDRrange: iLDR = iLDR + 1 - - LDRmin[iLDR] = np.min(aA[iLDR, :]) - LDRmax[iLDR] = np.max(aA[iLDR, :]) - Rmin = LDRmin[iLDR] * 0.995 # np.min(aA[iLDR,:]) * 0.995 - Rmax = LDRmax[iLDR] * 1.005 # np.max(aA[iLDR,:]) * 1.005 + LDRmin[iLDR] = np.amin(aLDRcorr[iLDR, :]) + LDRmax[iLDR] = np.amax(aLDRcorr[iLDR, :]) + Rmin = LDRmin[iLDR] * 0.995 # np.min(aLDRcorr[iLDR,:]) * 0.995 + Rmax = LDRmax[iLDR] * 1.005 # np.max(aLDRcorr[iLDR,:]) * 1.005 # plt.subplot(5,2,iLDR+1) plt.subplot(1, 5, iLDR + 1) - (n, bins, patches) = plt.hist(aA[iLDR, :], + (n, bins, patches) = plt.hist(aLDRcorr[iLDR, :], bins=100, log=False, range=[Rmin, Rmax], - alpha=0.5, normed=False, color='0.5', histtype='stepfilled') + alpha=0.5, density=False, color='0.5', histtype='stepfilled') for iaX in range(-naX, naX + 1): - plt.hist(aA[iLDR, aX == iaX], + plt.hist(aLDRcorr[iLDR, aX == iaX], range=[Rmin, Rmax], - bins=100, log=False, alpha=0.3, normed=False, histtype='stepfilled', + bins=100, log=False, alpha=0.3, density=False, histtype='stepfilled', label=str(round(X0 + iaX * daX / naX, 5))) - if (iLDR == 2): plt.legend() + if (iLDR == 2): + leg = plt.legend() + leg.get_frame().set_alpha(0.1) + plt.tick_params(axis='both', labelsize=9) plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2) @@ -1802,6 +2013,37 @@ # plt.close return + # Plot Etax + def PlotEtax(aVar, aX, X0, daX, iaX, naX): + # aVar is the name of the parameter and aX is the subset of aLDRcorr which is coloured in the plot + # example: PlotSubHist("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP) + fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) + iLDR = -1 + for LDRTrue in LDRrange: + iLDR = iLDR + 1 + Etaxmin[iLDR] = np.amin(aEtax[iLDR, :]) + Etaxmax[iLDR] = np.amax(aEtax[iLDR, :]) + Rmin = Etaxmin[iLDR] * 0.995 # np.min(aLDRcorr[iLDR,:]) * 0.995 + Rmax = Etaxmax[iLDR] * 1.005 # np.max(aLDRcorr[iLDR,:]) * 1.005 + + # plt.subplot(5,2,iLDR+1) + plt.subplot(1, 5, iLDR + 1) + (n, bins, patches) = plt.hist(aEtax[iLDR, :], + bins=100, log=False, + range=[Rmin, Rmax], + alpha=0.5, density=False, color='0.5', histtype='stepfilled') + for iaX in range(-naX, naX + 1): + plt.hist(aEtax[iLDR, aX == iaX], + range=[Rmin, Rmax], + bins=100, log=False, alpha=0.3, density=False, histtype='stepfilled', + label=str(round(X0 + iaX * daX / naX, 5))) + if (iLDR == 2): + leg = plt.legend() + leg.get_frame().set_alpha(0.1) + plt.tick_params(axis='both', labelsize=9) + plt.plot([Etax0, Etax0], [0, np.max(n)], 'r-', lw=2) + fig.suptitle('Etax - ' + LID + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar, fontsize=14, y=1.05) + return if (nDOLP > 0): PlotSubHist("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP) if (nRotL > 0): PlotSubHist("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) @@ -1825,7 +2067,14 @@ if (nRotaT > 0): PlotSubHist("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) if (nRotaR > 0): PlotSubHist("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) if (nLDRCal > 0): PlotSubHist("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) - + if (nTCalT > 0): PlotSubHist("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT) + if (nTCalR > 0): PlotSubHist("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR) + if (nNCal > 0): PlotSubHist("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal) + if (nNCal > 0): PlotSubHist("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal) + if (nNCal > 0): PlotSubHist("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal) + if (nNCal > 0): PlotSubHist("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal) + if (nNI > 0): PlotSubHist("SigNoiseIt", aNIt, 0, 1, iNIt, nNI) + if (nNI > 0): PlotSubHist("SigNoiseIr", aNIr, 0, 1, iNIr, nNI) plt.show() plt.close @@ -1850,8 +2099,8 @@ #F11min[iLDR] = np.min(aF11corr[iLDR,:]) #F11max[iLDR] = np.max(aF11corr[iLDR,:]) - #Rmin = F11min[iLDR] * 0.995 # np.min(aA[iLDR,:]) * 0.995 - #Rmax = F11max[iLDR] * 1.005 # np.max(aA[iLDR,:]) * 1.005 + #Rmin = F11min[iLDR] * 0.995 # np.min(aLDRcorr[iLDR,:]) * 0.995 + #Rmax = F11max[iLDR] * 1.005 # np.max(aLDRcorr[iLDR,:]) * 1.005 #Rmin = 0.8 #Rmax = 1.2 @@ -1860,11 +2109,11 @@ plt.subplot(1,5,iLDR+1) (n, bins, patches) = plt.hist(aF11corr[iLDR,:], bins=100, log=False, - alpha=0.5, normed=False, color = '0.5', histtype='stepfilled') + alpha=0.5, density=False, color = '0.5', histtype='stepfilled') for iaX in range(-naX,naX+1): plt.hist(aF11corr[iLDR,aX == iaX], - bins=100, log=False, alpha=0.3, normed=False, histtype='stepfilled', label = str(round(X0 + iaX*daX/naX,5))) + bins=100, log=False, alpha=0.3, density=False, histtype='stepfilled', label = str(round(X0 + iaX*daX/naX,5))) if (iLDR == 2): plt.legend() @@ -1903,6 +2152,11 @@ if (nRotaT > 0): PlotSubHistF11("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) if (nRotaR > 0): PlotSubHistF11("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) if (nLDRCal > 0): PlotSubHistF11("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) + if (nTCalT > 0): PlotSubHistF11("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT) + if (nTCalR > 0): PlotSubHistF11("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR) + if (nNCal > 0): PlotSubHistF11("CalNoise", aNCal, 0, 1/nNCal, iNCal, nNCal) + if (nNI > 0): PlotSubHistF11("SigNoise", aNI, 0, 1/nNI, iNI, nNI) + plt.show() plt.close @@ -1915,14 +2169,14 @@ iLDR = -1 for LDRTrue in LDRrange: iLDR = iLDR + 1 - LDRmin[iLDR] = np.min(aA[iLDR,:]) - LDRmax[iLDR] = np.max(aA[iLDR,:]) - Rmin = np.min(aA[iLDR,:]) * 0.999 - Rmax = np.max(aA[iLDR,:]) * 1.001 + LDRmin[iLDR] = np.min(aLDRcorr[iLDR,:]) + LDRmax[iLDR] = np.max(aLDRcorr[iLDR,:]) + Rmin = np.min(aLDRcorr[iLDR,:]) * 0.999 + Rmax = np.max(aLDRcorr[iLDR,:]) * 1.001 plt.subplot(5,2,iLDR+1) - (n, bins, patches) = plt.hist(aA[iLDR,:], + (n, bins, patches) = plt.hist(aLDRcorr[iLDR,:], range=[Rmin, Rmax], - bins=200, log=False, alpha=0.2, normed=False, color = '0.5', histtype='stepfilled') + bins=200, log=False, alpha=0.2, density=False, color = '0.5', histtype='stepfilled') plt.tick_params(axis='both', labelsize=9) plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2) plt.show() @@ -1931,9 +2185,18 @@ ''' # --- Plot LDRmin, LDRmax + iLDR = -1 + for LDRTrue in LDRrange: + iLDR = iLDR + 1 + LDRmin[iLDR] = np.amin(aLDRcorr[iLDR, :]) + LDRmax[iLDR] = np.amax(aLDRcorr[iLDR, :]) + fig2 = plt.figure() - plt.plot(LDRrange, LDRmax - LDRrange, linewidth=2.0, color='b') - plt.plot(LDRrange, LDRmin - LDRrange, linewidth=2.0, color='g') + LDRrangeA = np.array(LDRrange) + if((np.amax(LDRmax - LDRrangeA)-np.amin(LDRmin - LDRrangeA)) < 0.001): + plt.ylim(-0.001,0.001) + plt.plot(LDRrangeA, LDRmax - LDRrangeA, linewidth=2.0, color='b') + plt.plot(LDRrangeA, LDRmin - LDRrangeA, linewidth=2.0, color='g') plt.xlabel('LDRtrue', fontsize=18) plt.ylabel('LDRTrue-LDRmin, LDRTrue-LDRmax', fontsize=14) @@ -1944,17 +2207,62 @@ # --- Save LDRmin, LDRmax to file # http://stackoverflow.com/questions/4675728/redirect-stdout-to-a-file-in-python - with open('output_files\LDR_min_max_' + LID + '.dat', 'w') as f: + with open('output_files\\' + LID + '-' + InputFile[0:-3] + '-LDR_min_max.dat', 'w') as f: with redirect_stdout(f): print(LID) print("LDRtrue, LDRmin, LDRmax") - for i in range(len(LDRrange)): - print("{0:7.4f},{1:7.4f},{2:7.4f}".format(LDRrange[i], LDRmin[i], LDRmax[i])) + for i in range(len(LDRrangeA)): + print("{0:7.4f},{1:7.4f},{2:7.4f}".format(LDRrangeA[i], LDRmin[i], LDRmax[i])) + + if (bPlotEtax): + if (nDOLP > 0): PlotEtax("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP) + if (nRotL > 0): PlotEtax("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) + if (nRetE > 0): PlotEtax("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) + if (nRotE > 0): PlotEtax("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) + if (nDiE > 0): PlotEtax("DiE", aDiE, DiE0, dDiE, iDiE, nDiE) + if (nRetO > 0): PlotEtax("RetO", aRetO, RetO0, dRetO, iRetO, nRetO) + if (nRotO > 0): PlotEtax("RotO", aRotO, RotO0, dRotO, iRotO, nRotO) + if (nDiO > 0): PlotEtax("DiO", aDiO, DiO0, dDiO, iDiO, nDiO) + if (nDiC > 0): PlotEtax("DiC", aDiC, DiC0, dDiC, iDiC, nDiC) + if (nRotC > 0): PlotEtax("RotC", aRotC, RotC0, dRotC, iRotC, nRotC) + if (nRetC > 0): PlotEtax("RetC", aRetC, RetC0, dRetC, iRetC, nRetC) + if (nTP > 0): PlotEtax("TP", aTP, TP0, dTP, iTP, nTP) + if (nTS > 0): PlotEtax("TS", aTS, TS0, dTS, iTS, nTS) + if (nRP > 0): PlotEtax("RP", aRP, RP0, dRP, iRP, nRP) + if (nRS > 0): PlotEtax("RS", aRS, RS0, dRS, iRS, nRS) + if (nRetT > 0): PlotEtax("RetT", aRetT, RetT0, dRetT, iRetT, nRetT) + if (nRetR > 0): PlotEtax("RetR", aRetR, RetR0, dRetR, iRetR, nRetR) + if (nERaT > 0): PlotEtax("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) + if (nERaR > 0): PlotEtax("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) + if (nRotaT > 0): PlotEtax("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) + if (nRotaR > 0): PlotEtax("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) + if (nLDRCal > 0): PlotEtax("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) + if (nTCalT > 0): PlotEtax("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT) + if (nTCalR > 0): PlotEtax("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR) + if (nNCal > 0): PlotEtax("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal) + if (nNCal > 0): PlotEtax("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal) + if (nNCal > 0): PlotEtax("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal) + if (nNCal > 0): PlotEtax("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal) + if (nNI > 0): PlotEtax("SigNoiseIt", aNIt, 0, 1, iNIt, nNI) + if (nNI > 0): PlotEtax("SigNoiseIr", aNIr, 0, 1, iNIr, nNI) + plt.show() + plt.close + + #Etaxmin = np.amin(aEtax[1, :]) + Etaxmin = np.amin(aEtax[1, :]) + Etaxmax = np.amax(aEtax[1, :]) + Etaxstd = np.std(aEtax[1, :]) + Etaxmean = np.mean(aEtax[1, :]) + Etaxmedian = np.mean(aEtax[1, :]) + + print("Etax: mean±std, median, max-mean, mean-min") + print("{0:7.4f}±{1:7.4f},{2:7.4f},+{3:7.4f},-{4:7.4f}".format(Etaxmean, Etaxstd, Etaxmedian, Etaxmax-Etaxmean, Etaxmean-Etaxmin, )) + ''' # --- Plot K over LDRCal fig3 = plt.figure() - plt.plot(LDRCal0+aLDRCal*dLDRCal/nLDRCal,aX[4,:], linewidth=2.0, color='b') + plt.plot(LDRCal0+aLDRCal*dLDRCal/nLDRCal,aGHK[4,:], linewidth=2.0, color='b') plt.xlabel('LDRCal', fontsize=18) plt.ylabel('K', fontsize=14)
--- a/output_files/LDR_min_max_example lidar.dat Fri Nov 24 22:54:15 2017 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,7 +0,0 @@ -example lidar -LDRtrue, LDRmin, LDRmax - 0.0040, 0.0039, 0.0057 - 0.0200, 0.0197, 0.0219 - 0.1000, 0.0988, 0.1026 - 0.3000, 0.2965, 0.3042 - 0.4500, 0.4448, 0.4554
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/output_files/example lidar-optic_input_example_lidar-GHK.dat Thu Jan 24 21:24:25 2019 +0100 @@ -0,0 +1,48 @@ +From folder C:\Users\volker\Documents\atmospheric_lidar_ghk +Running prog lidar_correction_ghk_0.9.8d_Py3.7.py +Version 0.9.8d +Reading input file optic_input_example_lidar.py +for Lidar system: xx , example lidar + --- Input parameters: value ±error / ±steps ---------------------- +Laser: DOLP = 0.9950± 0.0050/ 1; Rotation alpha = 0.0000± 2.0000/ 1 + Diatt., Tunpol, Retard., Rotation (deg) +Emitter 0.0000± 0.1000/ 0, 1.0000, 0±180/ 0, 0.0000± 1.0000/ 0 +Receiver 0.0000± 0.2000/ 1, 1.0000, 0±180/ 0, 0.0000± 0.1000/ 0 +Calibrator 0.9998± 0.0002/ 1, 0.4000, 0±180/ 3, 0.0000± 0.1000/ 0 + --- Pol.-filter --- +ERT, RotT : 0.0010± 0.0010/ 0, 0.0000± 1.0000/ 0 +ERR, RotR : 0.0010± 0.0010/ 1, 90.0000± 1.0000/ 1 + --- PBS --- +TP,TS : 0.9500± 0.0100/ 1, 0.0050± 0.0010/ 1 +RP,RS : 0.0500± 0.0000/ 0, 0.9950± 0.0000/ 0 +DT,TT,DR,TR,Y : 0.9895, 0.4775, -0.9043, 0.5225, 1 + --- Combined PBS + Pol.-filter --- +DT,TT,DR,TR : 1.0000, 0.4750, -0.9999, 0.4975 +LDRCal during calibration in calibration range: 0.200± 0.150/ 1 + --- Additional ND filter attenuation (transmission) during the calibration --- +TCalT,TCalR : 1.0000± 0.0100/ 0, 0.1000± 0.0010/ 1 + +Rotation Error Epsilon For Normal Measurements = False +linear polarizer before receiver +Parallel signal detected in transmitted channel +RS_RP_depend_on_TS_TP = True + +IoutTp, IoutTm, IoutRp, IoutRm, It , Ir , dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr +=========================================================================================================== + GR , GT , HR , HT , K(0.004), K(0.02), K(0.1) , K(0.2) , K(0.3) , K(0.45) + 1.00000, 1.00000,-0.99490, 0.99499, 0.96129, 0.96248, 0.96796, 0.97382, 0.97880, 0.98502 +=========================================================================================================== + +Errors from neglecting GHK corrections and/or calibration: + LDRtrue, LDRunCorr,1/LDRunCorr, LDRsimx, 1/LDRsimx, LDRCorr + 0.00400, 0.00673, 148.53471, 0.00687, 145.62357, 0.00400 + 0.02000, 0.02316, 43.17641, 0.02362, 42.33019, 0.02000 + 0.10000, 0.10528, 9.49823, 0.10739, 9.31207, 0.10000 + 0.30000, 0.31044, 3.22120, 0.31665, 3.15806, 0.30000 + 0.45000, 0.46418, 2.15434, 0.47346, 2.11212, 0.45000 +LDRsimx = LDR of the nominal system directly from measured signals without calibration and GHK-corrections +LDRunCorr = LDR of the nominal system directly from measured signals with calibration but without GHK-corrections; electronic amplifications = 1 assumed +LDRCorr = LDR calibrated and GHK-corrected + +Errors from signal noise: +Signal counts: NCalT, NCalR, NILfac, nNCal, nNI = 50000, 50000, 2, 0, 0
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/output_files/example lidar-optic_input_example_lidar-LDR_min_max.dat Thu Jan 24 21:24:25 2019 +0100 @@ -0,0 +1,7 @@ +example lidar +LDRtrue, LDRmin, LDRmax + 0.0040, 0.0011, 0.0086 + 0.0200, 0.0158, 0.0252 + 0.1000, 0.0891, 0.1078 + 0.3000, 0.2726, 0.3143 + 0.4500, 0.4104, 0.4689
--- a/output_files/output_example lidar.dat Fri Nov 24 22:54:15 2017 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,36 +0,0 @@ -From C:\Users\volker\Documents\atmospheric_lidar_ghk -Running lidar_correction_ghk.py -Reading input file optic_input_example_lidar.py -for Lidar system: xx , example lidar - --- Input parameters: value ±error / ±steps ---------------------- -Laser: DOLP = 1.00000; rotation alpha = 0.0000± 2.0000/ 1 - Diatt., Tunpol, Retard., Rotation (deg) -Emitter 0.0000± 0.1000/ 0, 1.0000, 0±180/ 0, 0.0000± 1.0000/ 0 -Receiver 0.0000± 0.1000/ 1, 1.0000, 0±180/ 0, 0.0000± 0.5000/ 0 -Calibrator 0.9998± 0.0001/ 1, 0.4000, 0± 0/ 0, 0.0000± 0.1000/ 0 - --- Pol.-filter --- -ERT, RotT : 0.0010± 0.0010/ 1, 0.0000± 1.0000/ 1 -ERR, RotR : 0.0010± 0.0010/ 1, 90.0000± 1.0000/ 1 - --- PBS --- -TP,TS : 0.9500± 0.0100/ 1, 0.0200± 0.0100/ 1 -RP,RS : 0.0500± 0.0100/ 1, 0.9800± 0.0100/ 1 -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 - -======================================================================== - GR , GT , HR , HT , K(0.004), K(0.2) , K(0.45) - 1.90111, 1.95685,-1.90091, 1.95676, 0.93372, 0.94595, 0.95689 -======================================================================== - LDRtrue, LDRsimx, LDRCorr - 0.00400, 0.00418, 0.00400 - 0.02000, 0.02068, 0.02000 - 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 Fri Nov 24 22:54:15 2017 +0100 +++ b/system_settings/optic_input_example_lidar.py Thu Jan 24 21:24:25 2019 +0100 @@ -3,50 +3,49 @@ # which might improve the portability of the code within an executable. # Due to problems I had with some two letter variables, most variables are now with at least # three letters mixed small and capital. -# To be used with lidar_correction_ghk.py ver. 0.9.5 and larger - -# Do you want to calculate the errors? If not, just the GHK-parameters are determined. +# Error calculation? +# False: only calculate the GHK-parameters. True: calculate also errors. Can take a long time. Error_Calc = True # Header to identify the lidar system # -EID = "xx" # Earlinet station ID +EID = "xx" # Earlinet station ID LID = "example lidar" # Additional lidar ID (short descriptive text) print(" Lidar system :", EID, ", ", LID) -# +++ IL Laser and +-Uncertainty +# +++ IL Laser beam polarisation and +-Uncertainty DOLP, dDOLP, nDOLP = 0.995, 0.005, 1 #degree of linear polarization; default 1 RotL, dRotL, nRotL = 0., 2., 1 #alpha; rotation of laser polarization in degrees; default 0 # +++ ME Emitter optics and +-Uncertainty; default = no emitter optics -DiE, dDiE, nDiE = 0.0, 0.1, 0 # Diattenuation -TiE = 1.0 # Unpolarized transmittance -RetE, dRetE, nRetE = 0., 180., 0 # Retardance in degrees -RotE, dRotE, nRotE = 0., 1.0, 0 # beta: Rotation of optical element in degrees +DiE, dDiE, nDiE = 0.0, 0.1, 0 # Diattenuation; default 0 +TiE = 1.0 # Unpolarized transmittance; default 1 +RetE, dRetE, nRetE = 0., 180., 0 # Retardance in degrees; default 0 +RotE, dRotE, nRotE = 0., 1., 0 # beta: Rotation of optical element in degrees; default 0 # +++ MO Receiver optics including telescope -DiO, dDiO, nDiO = 0.0, 0.1, 1 # Diattenuation -TiO = 1.0 # Unpolarized transmittance -RetO, dRetO, nRetO = 0., 180., 0 # Retardance in degrees -RotO, dRotO, nRotO = 0., 0.5, 0 #gamma: Rotation of the optical element in degrees +DiO, dDiO, nDiO = 0.0, 0.2, 1 # Diattenuation; default 0 +TiO = 1.0 # Unpolarized transmittance; default 1 +RetO, dRetO, nRetO = 0., 180., 0 # Retardance in degrees; default 0 +RotO, dRotO, nRotO = 0., 0.1, 0 #gamma: Rotation of the optical element in degrees; default 0 # +++++ PBS MT Transmitting path defined with TS, TP, PolFilter extinction ratio ERaT, and +-Uncertainty # --- Polarizing beam splitter transmitting path TP, dTP, nTP = 0.95, 0.01, 1 # transmittance of the PBS for parallel polarized light -TS, dTS, nTS = 0.02, 0.01, 1 # transmittance of the PBS for cross polarized light +TS, dTS, nTS = 0.005, 0.001, 1 # transmittance of the PBS for cross polarized light RetT, dRetT, nRetT = 0.0, 180., 0 # Retardance in degrees # --- Pol.Filter behind transmitted path of PBS -ERaT, dERaT, nERaT = 0.001, 0.001, 1 # Extinction ratio -RotaT, dRotaT, nRotaT = 0., 1., 1 # Rotation of the Pol.-filter in degrees; usually close to 0° because TP >> TS, but for PollyXTs it can also be close to 90° +ERaT, dERaT, nERaT = 0.001, 0.001, 0 # Extinction ratio +RotaT, dRotaT, nRotaT = 0., 1., 0 # Rotation of the Pol.-filter in degrees; usually close to 0° because TP >> TS, but for PollyXTs it can also be close to 90° # -- -TiT = 0.5 * (TP + TS) -DiT = (TP-TS)/(TP+TS) -DaT = (1-ERaT)/(1+ERaT) -TaT = 0.5*(1+ERaT) +TiT = 0.5 * (TP + TS) # do not change this +DiT = (TP-TS)/(TP+TS) # do not change this +DaT = (1-ERaT)/(1+ERaT) # do not change this +TaT = 0.5*(1+ERaT) # do not change this -# +++++ PBS MR Reflecting path defined with RS, RP, PolFilter extinction ratio ERaR and +-Uncertainty +# +++++ PBS MR Reflecting path defined with RS, RP, and cleaning PolFilter extinction ratio ERaR and +-Uncertainty # ---- for PBS without absorption the change of RS and RP must depend on the change of TP and TS. Hence the values and uncertainties are not independent. -RS_RP_depend_on_TS_TP = False +RS_RP_depend_on_TS_TP = True # --- Polarizing beam splitter reflecting path if(RS_RP_depend_on_TS_TP): RP, dRP, nRP = 1-TP, 0.00, 0 # do not change this @@ -59,20 +58,26 @@ ERaR, dERaR, nERaR = 0.001, 0.001, 1 # Extinction ratio RotaR, dRotaR, nRotaR = 90., 1., 1 # Rotation of the Pol.-filter in degrees; usually 90° because RS >> RP, but for PollyXTs it can also be 0° # -- -TiR = 0.5 * (RP + RS) -DiR = (RP-RS)/(RP+RS) -DaR = (1-ERaR)/(1+ERaR) -TaR = 0.5*(1+ERaR) +TiR = 0.5 * (RP + RS) # do not change this +DiR = (RP-RS)/(RP+RS) # do not change this +DaR = (1-ERaR)/(1+ERaR) # do not change this +TaR = 0.5*(1+ERaR) # do not change this +# NEW --- Additional ND filter transmission (attenuation) during the calibration +TCalT, dTCalT, nTCalT = 1, 0.01, 0 # transmitting path, default 1, 0, 0 +TCalR, dTCalR, nTCalR = 0.1, 0.001, 1 # reflecting path, default 1, 0, 0 -# +++ Orientation of the PBS with respect to the reference plane (see Polarisation-orientation.png and Polarisation-orientation-2.png in /system_settings) -# Y = +1: polarisation in reference plane is transmitted, -# Y = -1: polarisation in reference plane is reflected. +# +++ Orientation of the PBS with respect to the reference plane (see Improvements_of_lidar_correction_ghk_ver.0.9.8_190124.pdf) +# Y = +1: polarisation in reference plane is finally transmitted, +# Y = -1: polarisation in reference plane is finally reflected. Y = +1. -# +++ Calibrator Location -LocC = 3 #location of calibrator: 1 = behind laser; 2 = behind emitter; 3 = before receiver; 4 = before PBS +# +++ Calibrator + # --- Calibrator Type used; defined by matrix values below TypeC = 3 #Type of calibrator: 1 = mechanical rotator; 2 = hwp rotator (fixed retardation); 3 = linear polarizer; 4 = qwp; 5 = circular polarizer; 6 = real HWP calibration +-22.5° + +# --- Calibrator Location +LocC = 3 #location of calibrator: 1 = behind laser; 2 = behind emitter; 3 = before receiver; 4 = before PBS # --- MC Calibrator parameters if TypeC == 1: #mechanical rotator DiC, dDiC, nDiC = 0., 0., 0 # Diattenuation @@ -81,7 +86,7 @@ RotC, dRotC, nRotC = 0., 0.1, 1 #constant calibrator rotation offset epsilon # Rotation error without calibrator: if False, then epsilon = 0 for normal measurements RotationErrorEpsilonForNormalMeasurements = True # is in general True for TypeC == 1 calibrator -elif TypeC == 2: # HWP rotator without retardance! +elif TypeC == 2: # HWP simulated by rotator without retardance! DiC, dDiC, nDiC = 0., 0., 0 # Diattenuation; ideal 0.0 TiC = 1. RetC, dRetC, nRetC = 180., 0., 0 # Retardance in degrees @@ -89,9 +94,9 @@ RotC, dRotC, nRotC = 0.0, 0.1, 1 #constant calibrator rotation offset epsilon RotationErrorEpsilonForNormalMeasurements = True # is in general True for TypeC == 2 calibrator elif TypeC == 3: # linear polarizer calibrator. Diattenuation DiC = (1-ERC)/(1+ERC); ERC = extinction ratio of calibrator - DiC, dDiC, nDiC = 0.9998, 0.0001, 1 # Diattenuation; ideal 1.0 + DiC, dDiC, nDiC = 0.9998, 0.00019, 1 # Diattenuation; ideal 1.0 TiC = 0.4 # ideal 0.5 - RetC, dRetC, nRetC = 0., 0., 0 # Retardance in degrees + RetC, dRetC, nRetC = 0., 180., 3 # Retardance in degrees RotC, dRotC, nRotC = 0.0, 0.1, 0 #constant calibrator rotation offset epsilon RotationErrorEpsilonForNormalMeasurements = False # is in general False for TypeC == 3 calibrator elif TypeC == 4: # QWP calibrator @@ -112,10 +117,40 @@ sys.exit() # --- 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 +LDRCal,dLDRCal,nLDRCal= 0.2, 0.15, 1 # spans most of the atmospheric depolarisation variability +# LDRCal,dLDRCal,nLDRCal= 0.009, 0.005, 1 # spans the interference filter influence # ==================================================== # NOTE: there is no need to change anything below. +# ==================================================== +# !!! don't change anything in this section !!! +# NEW *** +bPlotEtax = False # plot error histogramms for Etax + +# *** Only for signal noise errors *** +nNCal = 0 # error nNCal, calibration signals: one-sigma in steps to left and right +nNI = 0 # error nNI, 0° signals: one-sigma in steps to left and right; NI signals are calculated from NCalT and NCalR in main programm, but noise is assumed to be independent. + +# --- number of photon counts in the signal summed up in the calibration range during the calibration measurements +NCalT = 28184 # default 1e6, assumed the same in +45° and -45° signals +NCalR = 28184 # default 1e6, assumed the same in +45° and -45° signals +NILfac = 2 # (relative duration (laser shots) of standard (0°) measurement to calibration measurements) * (range of std. meas. smoothing / calibration range); example: 100000#/5000# * 100/1000 = 2 + # LDRmeas below will be used to calculate IR and IT of 0° signals. +# calculate signal counts only from parallel 0° signal assuming the same electronic amplification in both channels; overwrites above values +CalcFrom0deg = True +NI = 1e5 #number of photon counts in the parallel 0°-signal + +if(CalcFrom0deg): + # either eFactT or eFacR is = 1 => rel. amplification + eFacT = 1 # rel. amplification of transmitted channel, approximate values are sufficient; def. = 1 + eFacR = 10 # rel. amplification of reflected channel, approximate values are sufficient; def. = 1 + NILfac = 2 # (relative duration (laser shots) of standard (0°) measurement to calibration measurements) * (range of std. meas. smoothing / calibration range); example: 100000#/5000# * 100/1000 = 2 + + NCalT = NI / NILfac * TCalT * eFacT # photon counts in transmitted signal during calibration + NCalR = NI / NILfac * TCalR * eFacR # photon counts in reflected signal during calibration + # LDRmeas below will be used to calculate IR and IT of 0° signals. +# NEW *** End of signal noise error parameters *** + # --- LDRtrue for simulation of measurement => LDRsim LDRtrue = 0.4
