Mercurial > public > atmospheric_lidar_ghk / file comparison

comparison: lidar_correction_ghk.py

lidar_correction_ghk.py

changeset 28
79fa4a41420f
parent 26
28b5510492ba
equal deleted inserted replaced
27:e7b3f4bf631d 28:79fa4a41420f
1 # -*- coding: utf-8 -*- 1 # -*- coding: utf-8 -*-
2 """ 2 """
3 Copyright 2016, 2017 Volker Freudenthaler 3 Copyright 2016, 2019 Volker Freudenthaler
4 4
5 Licensed under the EUPL, Version 1.1 only (the "Licence"). 5 Licensed under the EUPL, Version 1.1 only (the "Licence").
6 6
7 You may not use this work except in compliance with the Licence. 7 You may not use this work except in compliance with the Licence.
8 A copy of the licence is distributed with the code. Alternatively, you may obtain 8 A copy of the licence is distributed with the code. Alternatively, you may obtain
15 OF ANY KIND, either express or implied. See the Licence for the specific language governing 15 OF ANY KIND, either express or implied. See the Licence for the specific language governing
16 permissions and limitations under the Licence. 16 permissions and limitations under the Licence.
17 17
18 Equation reference: http://www.atmos-meas-tech-discuss.net/amt-2015-338/amt-2015-338.pdf 18 Equation reference: http://www.atmos-meas-tech-discuss.net/amt-2015-338/amt-2015-338.pdf
19 With equations code from Appendix C 19 With equations code from Appendix C
20 Python 3.4.2 20 Python 3.7, seaborn 0.9.0
21 21
22 Code description: 22 Code description:
23 23
24 From measured lidar signals we cannot directly determine the desired backscatter coefficient (F11) and the linear depolarization ratio (LDR) 24 From measured lidar signals we cannot directly determine the desired backscatter coefficient (F11) and the linear depolarization ratio (LDR)
25 because of the cross talk between the channles and systematic errors of a lidar system. 25 because of the cross talk between the channles and systematic errors of a lidar system.
36 "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 36 "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
37 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 37 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
38 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). 38 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).
39 39
40 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. 40 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.
41
42 Ver. 0.9.8: - for details, see "Improvements_of_lidar_correction_ghk_ver.0.9.8_190124.pdf"
43 - 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
44 - 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.
45 - includes the simulation of signal noise
41 """ 46 """
42 47 # 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.
43 # 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.
44 from __future__ import print_function 48 from __future__ import print_function
45 # !/usr/bin/env python3 49 # !/usr/bin/env python3
46 50
47 import os 51 import os
48 import sys 52 import sys
56 sns_loaded = True 60 sns_loaded = True
57 except ImportError: 61 except ImportError:
58 sns_loaded = False 62 sns_loaded = False
59 63
60 import matplotlib.pyplot as plt 64 import matplotlib.pyplot as plt
61 from time import clock 65 # from time import clock # python 2
66 from timeit import default_timer as clock
62 67
63 # from matplotlib.backends.backend_pdf import PdfPages 68 # from matplotlib.backends.backend_pdf import PdfPages
64 # pdffile = '{}.pdf'.format('path') 69 # pdffile = '{}.pdf'.format('path')
65 # pp = PdfPages(pdffile) 70 # pp = PdfPages(pdffile)
66 ## pp.savefig can be called multiple times to save to multiple pages 71 ## pp.savefig can be called multiple times to save to multiple pages
111 116
112 sqr05 = 0.5 ** 0.5 117 sqr05 = 0.5 ** 0.5
113 118
114 # ---- Initial definition of variables; the actual values will be read in with exec(open('./optic_input.py').read()) below 119 # ---- Initial definition of variables; the actual values will be read in with exec(open('./optic_input.py').read()) below
115 # Do you want to calculate the errors? If not, just the GHK-parameters are determined. 120 # Do you want to calculate the errors? If not, just the GHK-parameters are determined.
121 ScriptVersion = "0.9.8d"
116 Error_Calc = True 122 Error_Calc = True
117 LID = "internal" 123 LID = "internal"
118 EID = "internal" 124 EID = "internal"
119 # --- IL Laser IL and +-Uncertainty 125 # --- IL Laser IL and +-Uncertainty
120 DOLP, dDOLP, nDOLP = 0.995, 0.005, 1 # degree of linear polarization; default 1 126 DOLP, dDOLP, nDOLP = 0.995, 0.005, 1 # degree of linear polarization; default 1
121 RotL, dRotL, nRotL = 0.0, 0.0, 1 # alpha; rotation of laser polarization in degrees; default 0 127 RotL, dRotL, nRotL = 0.0, 0.0, 1 # alpha; rotation of laser polarization in degrees; default 0
128 # IL = 1e5 #photons in the laser beam, including detection efficiency of the telescope, atmodspheric and r^2 attenuation
122 # --- ME Emitter and +-Uncertainty 129 # --- ME Emitter and +-Uncertainty
123 DiE, dDiE, nDiE = 0., 0.00, 1 # Diattenuation 130 DiE, dDiE, nDiE = 0., 0.00, 1 # Diattenuation
124 TiE = 1. # Unpolarized transmittance 131 TiE = 1. # Unpolarized transmittance
125 RetE, dRetE, nRetE = 0., 180.0, 0 # Retardance in degrees 132 RetE, dRetE, nRetE = 0., 180.0, 0 # Retardance in degrees
126 RotE, dRotE, nRotE = 0., 0.0, 0 # beta: Rotation of optical element in degrees 133 RotE, dRotE, nRotE = 0., 0.0, 0 # beta: Rotation of optical element in degrees
155 ERaR, dERaR, nERaR = 0.001, 0.001, 1 162 ERaR, dERaR, nERaR = 0.001, 0.001, 1
156 RotaR, dRotaR, nRotaR = 90., 3., 1 163 RotaR, dRotaR, nRotaR = 90., 3., 1
157 DaR = (1 - ERaR) / (1 + ERaR) 164 DaR = (1 - ERaR) / (1 + ERaR)
158 TaR = 0.5 * (1 + ERaR) 165 TaR = 0.5 * (1 + ERaR)
159 166
160 # Parellel signal detected in the transmitted channel => Y = 1, or in the reflected channel => Y = -1 167 # +++ Orientation of the PBS with respect to the reference plane (see Polarisation-orientation.png and Polarisation-orientation-2.png in /system_settings)
161 Y = -1. 168 # Y = +1: polarisation in reference plane is finally transmitted,
169 # Y = -1: polarisation in reference plane is finally reflected.
170 Y = 1.
162 171
163 # Calibrator = type defined by matrix values 172 # Calibrator = type defined by matrix values
164 LocC = 4 # location of calibrator: behind laser = 1; behind emitter = 2; before receiver = 3; before PBS = 4 173 LocC = 4 # location of calibrator: behind laser = 1; behind emitter = 2; before receiver = 3; before PBS = 4
174
175 # --- Additional attenuation (transmission of the ND-filter) during the calibration
176 TCalT, dTCalT, nTCalT = 1, 0, 0 # transmitting path; error calc not working yet
177 TCalR, dTCalR, nTCalR = 1, 0, 0 # reflecting path; error calc not working yet
178
179 # *** signal noise error calculation
180 # --- number of photon counts in the signal summed up in the calibration range during the calibration measurements
181 NCalT = 1e6 # default 1e6, assumed the same in +45° and -45° signals
182 NCalR = 1e6 # default 1e6, assumed the same in +45° and -45° signals
183 NILfac = 200 # duration of standard (0°) measurement relative to calibration measurements
184 nNCal = 0 # error nNCal: one-sigma in steps to left and right for calibration signals
185 nNI = 0 # error nNI: one-sigma in steps to left and right for 0° signals
186 IoutTp0, IoutTp, dIoutTp0 = 0.5, 0.5, 0
187 IoutTm0, IoutTm, dIoutTm0 = 0.5, 0.5, 0
188 IoutRp0, IoutRp, dIoutRp0 = 0.5, 0.5, 0
189 IoutRm0, IoutRm, dIoutRm0 = 0.5, 0.5, 0
190 It0, It, dIt0 = 1 , 1, 0
191 Ir0, Ir, dTr0 = 1 , 1, 0
192 CalcFrom0deg = True
165 193
166 TypeC = 3 # linear polarizer calibrator 194 TypeC = 3 # linear polarizer calibrator
167 # example with extinction ratio 0.001 195 # example with extinction ratio 0.001
168 DiC, dDiC, nDiC = 1.0, 0., 0 # ideal 1.0 196 DiC, dDiC, nDiC = 1.0, 0., 0 # ideal 1.0
169 TiC = 0.5 # ideal 0.5 197 TiC = 0.5 # ideal 0.5
171 RotC, dRotC, nRotC = 0.0, 0.1, 0 # constant calibrator offset epsilon 199 RotC, dRotC, nRotC = 0.0, 0.1, 0 # constant calibrator offset epsilon
172 RotationErrorEpsilonForNormalMeasurements = False # is in general False for TypeC == 3 calibrator 200 RotationErrorEpsilonForNormalMeasurements = False # is in general False for TypeC == 3 calibrator
173 201
174 # Rotation error without calibrator: if False, then epsilon = 0 for normal measurements 202 # Rotation error without calibrator: if False, then epsilon = 0 for normal measurements
175 RotationErrorEpsilonForNormalMeasurements = True 203 RotationErrorEpsilonForNormalMeasurements = True
176 204 # BSR backscatter ratio
205 # BSR, dBSR, nBSR = 10, 0.05, 1
206 BSR = np.zeros(5)
207 BSR = [1.1, 2, 5, 10, 50]
208 # theoretical molecular LDR LDRm
209 LDRm, dLDRm, nLDRm = 0.004, 0.001, 1
177 # LDRCal assumed atmospheric linear depolarization ratio during the calibration measurements (first guess) 210 # LDRCal assumed atmospheric linear depolarization ratio during the calibration measurements (first guess)
178 LDRCal0, dLDRCal, nLDRCal = 0.25, 0.04, 1 211 LDRCal0, dLDRCal, nLDRCal = 0.25, 0.04, 1
179 LDRCal = LDRCal0 212 LDRCal = LDRCal0
180 # measured LDRm will be corrected with calculated parameters 213 # measured LDRm will be corrected with calculated parameters
181 LDRmeas = 0.015 214 LDRmeas = 0.015
182 # LDRtrue for simulation of measurement => LDRsim 215 # LDRtrue for simulation of measurement => LDRsim
183 LDRtrue = 0.5 216 LDRtrue = 0.5
184 LDRtrue2 = 0.004 217 LDRtrue2 = 0.004
185 218 LDRunCorr = 1
186 # Initialize other values to 0 219 # Initialize other values to 0
187 ER, nER, dER = 0.001, 0, 0.001 220 ER, nER, dER = 0.001, 0, 0.001
188 K = 0. 221 K = 0.
189 Km = 0. 222 Km = 0.
190 Kp = 0. 223 Kp = 0.
196 229
197 Loc = ['', 'behind laser', 'behind emitter', 'before receiver', 'before PBS'] 230 Loc = ['', 'behind laser', 'behind emitter', 'before receiver', 'before PBS']
198 Type = ['', 'mechanical rotator', 'hwp rotator', 'linear polarizer', 'qwp rotator', 'circular polarizer', 231 Type = ['', 'mechanical rotator', 'hwp rotator', 'linear polarizer', 'qwp rotator', 'circular polarizer',
199 'real HWP +-22.5°'] 232 'real HWP +-22.5°']
200 dY = ['reflected channel', '', 'transmitted channel'] 233 dY = ['reflected channel', '', 'transmitted channel']
234 bPlotEtax = False
201 235
202 # end of initial definition of variables 236 # end of initial definition of variables
203 # ******************************************************************************************************************************* 237 # *******************************************************************************************************************************
204 238
205 # --- Read actual lidar system parameters from optic_input.py (should be in sub-directory 'system_settings') 239 # --- Read actual lidar system parameters from optic_input.py (must be in the programs sub-directory 'system_settings')
206 InputFile = 'optic_input_crossed_mirrors_test.py' 240 #InputFile = 'optic_input_example_lidar_2.py'
207 InputFile = 'optic_input_crossed_mirrors_test_combined.py' 241 #InputFile = 'optic_input_example_lidar_3.py'
208 InputFile = 'optic_input_ver8c_POLIS_355_Mar2017.py' 242 #InputFile = 'optic_input_example_lidar_4.py'
209 # InputFile = 'optic_input_ver8c_POLIS_532_Mar2017.py' 243 #InputFile = 'optic_input_example_lidar_5.py'
210 InputFile = 'optic_input_example_lidar_0.9.5.py' 244 InputFile = 'optic_input_example_lidar.py'
211 InputFile = 'optic_input_ver8c_PollyXT_532_Lacros.py' 245
212 InputFile = 'optic_input_ver8c_CUT_532_May2017.py'
213 InputFile = 'optic_input_ver8c_MUSA.py'
214 InputFile = 'optic_input_ver10_RALI_JA.py'
215 InputFile = 'optic_input_ver10_RALI_act.py'
216 InputFile = 'optic_input_ver8c-IPRAL-170331.py'
217 ''' 246 '''
218 print("From ", dname) 247 print("From ", dname)
219 print("Running ", fname) 248 print("Running ", fname)
220 print("Reading input file ", InputFile, " for") 249 print("Reading input file ", InputFile, " for")
221 ''' 250 '''
264 293
265 ERaR0, dERaR, nERaR = ERaR, dERaR, nERaR 294 ERaR0, dERaR, nERaR = ERaR, dERaR, nERaR
266 RotaR0, dRotaR, nRotaR = RotaR, dRotaR, nRotaR 295 RotaR0, dRotaR, nRotaR = RotaR, dRotaR, nRotaR
267 296
268 LDRCal0, dLDRCal, nLDRCal = LDRCal, dLDRCal, nLDRCal 297 LDRCal0, dLDRCal, nLDRCal = LDRCal, dLDRCal, nLDRCal
269 # LDRCal0,dLDRCal,nLDRCal=LDRCal,dLDRCal,0 298
299 # BSR0, dBSR, nBSR = BSR, dBSR, nBSR
300 LDRm0, dLDRm, nLDRm = LDRm, dLDRm, nLDRm
270 # ---------- End of manual parameter change 301 # ---------- End of manual parameter change
271 302
272 RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC = RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0 303 RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC = RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0
273 TP, TS, RP, RS, ERaT, RotaT, RetT, ERaR, RotaR, RetR = TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0 304 TP, TS, RP, RS, ERaT, RotaT, RetT, ERaR, RotaR, RetR = TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0
274 LDRCal = LDRCal0 305 LDRCal = LDRCal0
275 DTa0, TTa0, DRa0, TRa0, LDRsimx, LDRCorr = 0, 0, 0, 0, 0, 0 306 DTa0, TTa0, DRa0, TRa0, LDRsimx, LDRCorr = 0, 0, 0, 0, 0, 0
307 TCalT0, TCalR0 = TCalT, TCalR
276 308
277 TiT = 0.5 * (TP + TS) 309 TiT = 0.5 * (TP + TS)
278 DiT = (TP - TS) / (TP + TS) 310 DiT = (TP - TS) / (TP + TS)
279 ZiT = (1. - DiT ** 2) ** 0.5 311 ZiT = (1. - DiT ** 2) ** 0.5
280 TiR = 0.5 * (RP + RS) 312 TiR = 0.5 * (RP + RS)
281 DiR = (RP - RS) / (RP + RS) 313 DiR = (RP - RS) / (RP + RS)
282 ZiR = (1. - DiR ** 2) ** 0.5 314 ZiR = (1. - DiR ** 2) ** 0.5
283 315
284 316 C2aT = np.cos(np.deg2rad(2 * RotaT))
317 C2aR = np.cos(np.deg2rad(2 * RotaR))
318 ATPT = (1 + C2aT * DaT * DiT)
319 ARPT = (1 + C2aR * DaR * DiR)
320 TTa = TiT * TaT * ATPT # unpolarized transmission
321 TRa = TiR * TaR * ARPT # unpolarized transmission
322 Eta0 = TRa / TTa
323 # --- this subroutine is for the calculation of the PLDR from LDR, BSR, and LDRm -----------------------------------------------------
324 def CalcPLDR(LDR, BSR, LDRm):
325 PLDR = (BSR * (1. + LDRm) * LDR - LDRm * (1. + LDR)) / (BSR * (1. + LDRm) - (1. + LDR))
326 return (PLDR)
285 # --- this subroutine is for the calculation with certain fixed parameters ----------------------------------------------------- 327 # --- this subroutine is for the calculation with certain fixed parameters -----------------------------------------------------
286 def Calc(DOLP, RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC, TP, TS, RP, RS, ERaT, RotaT, RetT, ERaR, RotaR, RetR, 328 def Calc(TCalT, TCalR, NCalT, NCalR, DOLP, RotL, RotE, RetE, DiE, RotO, RetO, DiO,
287 LDRCal): 329 RotC, RetC, DiC, TP, TS, RP, RS,
330 ERaT, RotaT, RetT, ERaR, RotaR, RetR, LDRCal):
288 # ---- Do the calculations of bra-ket vectors 331 # ---- Do the calculations of bra-ket vectors
289 h = -1. if TypeC == 2 else 1 332 h = -1. if TypeC == 2 else 1
290 # from input file: assumed LDRCal for calibration measurements 333 # from input file: assumed LDRCal for calibration measurements
291 aCal = (1. - LDRCal) / (1 + LDRCal) 334 aCal = (1. - LDRCal) / (1 + LDRCal)
292 # from input file: measured LDRm and true LDRtrue, LDRtrue2 => 335 # from input file: measured LDRm and true LDRtrue, LDRtrue2 =>
387 S2aT = np.sin(np.deg2rad(h * 2 * RotaT)) 430 S2aT = np.sin(np.deg2rad(h * 2 * RotaT))
388 C2aT = np.cos(np.deg2rad(2 * RotaT)) 431 C2aT = np.cos(np.deg2rad(2 * RotaT))
389 S2aR = np.sin(np.deg2rad(h * 2 * RotaR)) 432 S2aR = np.sin(np.deg2rad(h * 2 * RotaR))
390 C2aR = np.cos(np.deg2rad(2 * RotaR)) 433 C2aR = np.cos(np.deg2rad(2 * RotaR))
391 434
392 # Analyzer As before the PBS Eq. D.5 435 # Analyzer As before the PBS Eq. D.5; combined PBS and cleaning pol-filter
393 ATP1 = (1 + C2aT * DaT * DiT) 436 ATPT = (1 + C2aT * DaT * DiT) # unpolarized transmission correction
394 ATP2 = Y * (DiT + C2aT * DaT) 437 TTa = TiT * TaT * ATPT # unpolarized transmission
395 ATP3 = Y * S2aT * DaT * ZiT * CosT 438 ATP1 = 1
396 ATP4 = S2aT * DaT * ZiT * SinT 439 ATP2 = Y * (DiT + C2aT * DaT) / ATPT
440 ATP3 = Y * S2aT * DaT * ZiT * CosT / ATPT
441 ATP4 = S2aT * DaT * ZiT * SinT / ATPT
397 ATP = np.array([ATP1, ATP2, ATP3, ATP4]) 442 ATP = np.array([ATP1, ATP2, ATP3, ATP4])
398 443 DTa = ATP2 * Y
399 ARP1 = (1 + C2aR * DaR * DiR) 444
400 ARP2 = Y * (DiR + C2aR * DaR) 445 ARPT = (1 + C2aR * DaR * DiR) # unpolarized transmission correction
401 ARP3 = Y * S2aR * DaR * ZiR * CosR 446 TRa = TiR * TaR * ARPT # unpolarized transmission
402 ARP4 = S2aR * DaR * ZiR * SinR 447 ARP1 = 1
448 ARP2 = Y * (DiR + C2aR * DaR) / ARPT
449 ARP3 = Y * S2aR * DaR * ZiR * CosR / ARPT
450 ARP4 = S2aR * DaR * ZiR * SinR / ARPT
403 ARP = np.array([ARP1, ARP2, ARP3, ARP4]) 451 ARP = np.array([ARP1, ARP2, ARP3, ARP4])
404 452 DRa = ARP2 * Y
405 DTa = ATP2 * Y / ATP1 453
406 DRa = ARP2 * Y / ARP1
407 454
408 # ---- Calculate signals and correction parameters for diffeent locations and calibrators 455 # ---- Calculate signals and correction parameters for diffeent locations and calibrators
409 if LocC == 4: # Calibrator before the PBS 456 if LocC == 4: # Calibrator before the PBS
410 # print("Calibrator location not implemented yet") 457 # print("Calibrator location not implemented yet")
411 458
806 853
807 else: 854 else:
808 print("Calibrator location not implemented yet") 855 print("Calibrator location not implemented yet")
809 sys.exit() 856 sys.exit()
810 857
811 # Determination of the correction K of the calibration factor 858 # Determination of the correction K of the calibration factor.
812 IoutTp = TaT * TiT * TiO * TiE * (AT + BT) 859 IoutTp = TTa * TiC * TiO * TiE * (AT + BT)
813 IoutTm = TaT * TiT * TiO * TiE * (AT - BT) 860 IoutTm = TTa * TiC * TiO * TiE * (AT - BT)
814 IoutRp = TaR * TiR * TiO * TiE * (AR + BR) 861 IoutRp = TRa * TiC * TiO * TiE * (AR + BR)
815 IoutRm = TaR * TiR * TiO * TiE * (AR - BR) 862 IoutRm = TRa * TiC * TiO * TiE * (AR - BR)
816
817 # --- Results and Corrections; electronic etaR and etaT are assumed to be 1 863 # --- Results and Corrections; electronic etaR and etaT are assumed to be 1
818 Etapx = IoutRp / IoutTp 864 Etapx = IoutRp / IoutTp
819 Etamx = IoutRm / IoutTm 865 Etamx = IoutRm / IoutTm
820 Etax = (Etapx * Etamx) ** 0.5 866 Etax = (Etapx * Etamx) ** 0.5
821 867
822 Eta = (TaR * TiR) / (TaT * TiT) # Eta = Eta*/K Eq. 84 868 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
823 K = Etax / Eta 869 K = Etax / Eta
824 870
825 # For comparison with Volkers Libreoffice Müller Matrix spreadsheet 871 # For comparison with Volkers Libreoffice Müller Matrix spreadsheet
826 # Eta_test_p = (IoutRp/IoutTp) 872 # Eta_test_p = (IoutRp/IoutTp)
827 # Eta_test_m = (IoutRm/IoutTm) 873 # Eta_test_m = (IoutRm/IoutTm)
828 # Eta_test = (Eta_test_p*Eta_test_m)**0.5 874 # Eta_test = (Eta_test_p*Eta_test_m)**0.5
829 875
830 # ----- Forward simulated signals and LDRsim with atrue; from input file 876 # ----- random error calculation ----------
831 It = TaT * TiT * TiO * TiE * (GT + atrue * HT) 877 # noise must be calculated with the photon counts of measured signals;
832 Ir = TaR * TiR * TiO * TiE * (GR + atrue * HR) 878 # relative standard deviation of calibration signals with LDRcal; assumed to be statisitcally independent
833 # LDRsim = 1/Eta*Ir/It # simulated LDR* with Y from input file 879 # normalised noise errors
880 if (CalcFrom0deg):
881 dIoutTp = (NCalT * IoutTp) ** -0.5
882 dIoutTm = (NCalT * IoutTm) ** -0.5
883 dIoutRp = (NCalR * IoutRp) ** -0.5
884 dIoutRm = (NCalR * IoutRm) ** -0.5
885 else:
886 dIoutTp = (NCalT ** -0.5)
887 dIoutTm = (NCalT ** -0.5)
888 dIoutRp = (NCalR ** -0.5)
889 dIoutRm = (NCalR ** -0.5)
890 # Forward simulated 0°-signals with LDRCal with atrue; from input file
891
892 It = TTa * TiO * TiE * (GT + atrue * HT)
893 Ir = TRa * TiO * TiE * (GR + atrue * HR)
894 # relative standard deviation of standard signals with LDRmeas; assumed to be statisitcally independent
895 if (CalcFrom0deg):
896 dIt = ((NCalT * It / IoutTp * NILfac / TCalT) ** -0.5)
897 dIr = ((NCalR * Ir / IoutRp * NILfac / TCalR) ** -0.5)
898 else:
899 dIt = ((NCalT * 2 * NILfac / TCalT ) ** -0.5) * It
900 dIr = ((NCalR * 2 * NILfac / TCalR) ** -0.5) * Ir
901
902 # ----- Forward simulated LDRsim = 1/Eta*Ir/It # simulated LDR* with Y from input file
834 LDRsim = Ir / It # simulated uncorrected LDR with Y from input file 903 LDRsim = Ir / It # simulated uncorrected LDR with Y from input file
835 # Corrected LDRsimCorr from forward simulated LDRsim (atrue) 904 # Corrected LDRsimCorr from forward simulated LDRsim (atrue)
836 # LDRsimCorr = (1./Eta*LDRsim*(GT+HT)-(GR+HR))/((GR-HR)-1./Eta*LDRsim*(GT-HT)) 905 # LDRsimCorr = (1./Eta*LDRsim*(GT+HT)-(GR+HR))/((GR-HR)-1./Eta*LDRsim*(GT-HT))
837 ''' 906 '''
838 if ((Y == -1.) and (abs(RotL0) < 45)) or ((Y == +1.) and (abs(RotL0) > 45)): 907 if ((Y == -1.) and (abs(RotL0) < 45)) or ((Y == +1.) and (abs(RotL0) > 45)):
839 LDRsimx = 1. / LDRsim / Etax 908 LDRsimx = 1. / LDRsim / Etax
840 else: 909 else:
841 LDRsimx = LDRsim / Etax 910 LDRsimx = LDRsim / Etax
842 ''' 911 '''
843 LDRsimx = LDRsim 912 LDRsimx = LDRsim
844 913
845 # The following is correct without doubt 914 # The following is correct without doubt
846 # LDRCorr = (LDRsim*K/Etax*(GT+HT)-(GR+HR))/((GR-HR)-LDRsim*K/Etax*(GT-HT)) 915 # LDRCorr = (LDRsim/(Etax/K)*(GT+HT)-(GR+HR))/((GR-HR)-LDRsim/(Etax/K)*(GT-HT))
847 916
848 # The following is a test whether the equations for calibration Etax and normal signal (GHK, LDRsim) are consistent 917 # The following is a test whether the equations for calibration Etax and normal signal (GHK, LDRsim) are consistent
849 LDRCorr = (LDRsim / Eta * (GT + HT) - (GR + HR)) / ((GR - HR) - LDRsim * K / Etax * (GT - HT)) 918 LDRCorr = (LDRsim / (Etax / K) * (GT + HT) - (GR + HR)) / ((GR - HR) - LDRsim / (Etax / K) * (GT - HT))
919 # here we could also use Eta instead of Etax / K => how to test whether Etax is correct? => comparison with MüllerMatrix simulation!
920 # Without any correction: only measured It, Ir, EtaX are used
921 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)))
922
850 #LDRCorr = LDRsimx # for test only 923 #LDRCorr = LDRsimx # for test only
851 TTa = TiT * TaT # *ATP1 924
852 TRa = TiR * TaR # *ARP1 925 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
853 926
854 F11sim = 1 / (TiO * TiE) * ( 927 return (IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr,
855 (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 928 GT, HT, GR, HR, K, Eta, LDRsimx, LDRCorr, DTa, DRa, TTa, TRa, F11sim, LDRunCorr)
856 929
857 return (GT, HT, GR, HR, K, Eta, LDRsimx, LDRCorr, DTa, DRa, TTa, TRa, F11sim)
858 930
859 931
860 # ******************************************************************************************************************************* 932 # *******************************************************************************************************************************
861 933
862 # --- CALC truth with assumed true parameters from the input file 934 # --- CALC with assumed true parameters from the input file
863 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(DOLP0, RotL0, RotE0, RetE0, DiE0, RotO0, 935 IoutTp0, IoutTm0, IoutRp0, IoutRm0, It0, Ir0, dIoutTp0, dIoutTm0, dIoutRp0, dIoutRm0, dIt0, dIr0, \
864 RetO0, DiO0, RotC0, RetC0, DiC0, 936 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \
865 TP0, TS0, RP0, RS0, ERaT0, 937 Calc(TCalT, TCalR, NCalT, NCalR, DOLP0, RotL0, RotE0, RetE0, DiE0,
866 RotaT0, RetT0, ERaR0, RotaR0, 938 RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0,
867 RetR0, LDRCal0) 939 ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0)
868 940 Etax0 = K0 * Eta0
869 # --- Print parameters to console and output file 941 # --- Print parameters to console and output file
870 with open('output_files\output_' + LID + '.dat', 'w') as f: 942 with open('output_files\\' + LID + '-' + InputFile[0:-3] + '-GHK.dat', 'w') as f:
871 with redirect_stdout(f): 943 with redirect_stdout(f):
872 print("From ", dname) 944 print("From folder", dname)
873 print("Running ", fname) 945 print("Running prog", fname)
946 print("Version", ScriptVersion)
874 print("Reading input file ", InputFile) # , " for Lidar system :", EID, ", ", LID) 947 print("Reading input file ", InputFile) # , " for Lidar system :", EID, ", ", LID)
875 print("for Lidar system: ", EID, ", ", LID) 948 print("for Lidar system: ", EID, ", ", LID)
876 # --- Print iput information********************************* 949 # --- Print iput information*********************************
877 print(" --- Input parameters: value ±error / ±steps ----------------------") 950 print(" --- Input parameters: value ±error / ±steps ----------------------")
878 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, 951 print("{0:5}{1:5} {2:6.4f}±{3:7.4f}/{4:2d}; {5:8} {6:8.4f}±{7:7.4f}/{8:2d}".format(
879 " Rotation alpha = ", RotL0, dRotL, 952 "Laser: ", "DOLP =", DOLP0, dDOLP, nDOLP," Rotation alpha = ", RotL0, dRotL, nRotL))
880 nRotL))
881 print(" Diatt., Tunpol, Retard., Rotation (deg)") 953 print(" Diatt., Tunpol, Retard., Rotation (deg)")
882 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( 954 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(
883 "Emitter ", DiE0, dDiE, TiE, RetE0, dRetE, RotE0, dRotE, nDiE, nRetE, nRotE)) 955 "Emitter ", DiE0, dDiE, TiE, RetE0, dRetE, RotE0, dRotE, nDiE, nRetE, nRotE))
884 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( 956 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(
885 "Receiver ", DiO0, dDiO, TiO, RetO0, dRetO, RotO0, dRotO, nDiO, nRetO, nRotO)) 957 "Receiver ", DiO0, dDiO, TiO, RetO0, dRetO, RotO0, dRotO, nDiO, nRetO, nRotO))
886 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( 958 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(
887 "Calibrator ", DiC0, dDiC, TiC, RetC0, dRetC, RotC0, dRotC, nDiC, nRetC, nRotC)) 959 "Calibrator ", DiC0, dDiC, TiC, RetC0, dRetC, RotC0, dRotC, nDiC, nRetC, nRotC))
888 print("{0:12}".format(" --- Pol.-filter ---")) 960 print("{0:12}".format(" --- Pol.-filter ---"))
889 print( 961 print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format(
890 "{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format("ERT, RotT :", ERaT0, dERaT, nERaT, 962 "ERT, RotT :", ERaT0, dERaT, nERaT, RotaT0, dRotaT, nRotaT))
891 RotaT0, dRotaT, nRotaT)) 963 print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format(
892 print( 964 "ERR, RotR :", ERaR0, dERaR, nERaR, RotaR0, dRotaR, nRotaR))
893 "{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format("ERR, RotR :", ERaR0, dERaR, nERaR,
894 RotaR0, dRotaR, nRotaR))
895 print("{0:12}".format(" --- PBS ---")) 965 print("{0:12}".format(" --- PBS ---"))
896 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, 966 print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format(
897 dTS, nTS)) 967 "TP,TS :", TP0, dTP, nTP, TS0, dTS, nTS))
898 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, 968 print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format(
899 dRS, nRS)) 969 "RP,RS :", RP0, dRP, nRP, RS0, dRS, nRS))
900 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)) 970 print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}, {5:1.0f}".format(
971 "DT,TT,DR,TR,Y :", DiT, TiT, DiR, TiR, Y))
901 print("{0:12}".format(" --- Combined PBS + Pol.-filter ---")) 972 print("{0:12}".format(" --- Combined PBS + Pol.-filter ---"))
902 print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}".format("DT,TT,DR,TR :", DTa0, TTa0, DRa0, TRa0)) 973 print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}".format(
903 print("{0:26}: {1:6.3f}± {2:5.3f}/{3:2d}".format("LDRCal during calibration in calibration range", LDRCal0, 974 "DT,TT,DR,TR :", DTa0, TTa0, DRa0, TRa0))
904 dLDRCal, nLDRCal)) 975 print("{0:26}: {1:6.3f}± {2:5.3f}/{3:2d}".format(
976 "LDRCal during calibration in calibration range", LDRCal0, dLDRCal, nLDRCal))
977 print("{0:12}".format(" --- Additional ND filter attenuation (transmission) during the calibration ---"))
978 print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format(
979 "TCalT,TCalR :", TCalT0, dTCalT, nTCalT, TCalR0, dTCalR, nTCalR))
905 print() 980 print()
906 print("Rotation Error Epsilon For Normal Measurements = ", RotationErrorEpsilonForNormalMeasurements) 981 print("Rotation Error Epsilon For Normal Measurements = ", RotationErrorEpsilonForNormalMeasurements)
907 print(Type[TypeC], Loc[LocC]) 982 print(Type[TypeC], Loc[LocC])
908 print("Parallel signal detected in", dY[int(Y + 1)]) 983 print("Parallel signal detected in", dY[int(Y + 1)])
909 print("RS_RP_depend_on_TS_TP = ", RS_RP_depend_on_TS_TP) 984 print("RS_RP_depend_on_TS_TP = ", RS_RP_depend_on_TS_TP)
920 # print(" --- LDRCal during calibration") 995 # print(" --- LDRCal during calibration")
921 996
922 # print("{0:8}={1:8.5f};{2:8}={3:8.5f}".format(" IinP",IinP," F11sim",F11sim)) 997 # print("{0:8}={1:8.5f};{2:8}={3:8.5f}".format(" IinP",IinP," F11sim",F11sim))
923 print() 998 print()
924 999
925 K0List = np.zeros(3) 1000 K0List = np.zeros(6)
926 LDRsimxList = np.zeros(3) 1001 LDRsimxList = np.zeros(6)
927 LDRCalList = 0.004, 0.2, 0.45 1002 LDRCalList = 0.004, 0.02, 0.1, 0.2, 0.3, 0.45
928 # 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. 1003 # 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.
929 # Still with assumed true parameters in input file 1004 # Still with assumed true parameters in input file
1005
1006 facIt = NCalT / TCalT0 * NILfac
1007 facIr = NCalR / TCalR0 * NILfac
1008 print("IoutTp, IoutTm, IoutRp, IoutRm, It , Ir , dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr")
1009
930 for i, LDRCal in enumerate(LDRCalList): 1010 for i, LDRCal in enumerate(LDRCalList):
931 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(DOLP0, RotL0, RotE0, RetE0, 1011 IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr, \
932 DiE0, RotO0, RetO0, 1012 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \
933 DiO0, RotC0, RetC0, 1013 Calc(TCalT0, TCalR0, NCalT, NCalR, DOLP0, RotL0, RotE0, RetE0, DiE0,
934 DiC0, TP0, TS0, RP0, 1014 RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0,
935 RS0, ERaT0, RotaT0, 1015 ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal)
936 RetT0, ERaR0, RotaR0,
937 RetR0, LDRCal)
938 K0List[i] = K0 1016 K0List[i] = K0
939 LDRsimxList[i] = LDRsimx 1017 LDRsimxList[i] = LDRsimx
940 1018
941 print('========================================================================') 1019 # check error signals
942 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)", 1020 # 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}"
943 " K(0.45)")) 1021 # .format(IoutTp * NCalT, IoutTm * NCalT, IoutRp * NCalR, IoutRm * NCalR, It * facIt, Ir * facIr, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr))
944 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], 1022 #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))
945 K0List[1], K0List[2])) 1023 # end check error signals
946 print('========================================================================') 1024 print('===========================================================================================================')
947 print("{0:10},{1:10},{2:10},{3:10}".format(" LDRtrue", " LDRsimx", " 1/LDRsimx", " LDRCorr")) 1025 print("{0:8},{1:8},{2:8},{3:8},{4:9},{5:8},{6:9},{7:9},{8:9},{9:9}".format(
1026 " GR", " GT", " HR", " HT", " K(0.004)", " K(0.02)", " K(0.1)", " K(0.2)", " K(0.3)", " K(0.45)"))
1027 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(
1028 GR0, GT0, HR0, HT0, K0List[0], K0List[1], K0List[2], K0List[3], K0List[4], K0List[5]))
1029 print('===========================================================================================================')
1030 print()
1031 print("Errors from neglecting GHK corrections and/or calibration:")
1032 print("{0:>10},{1:>10},{2:>10},{3:>10},{4:>10},{5:>10}".format(
1033 "LDRtrue", "LDRunCorr", "1/LDRunCorr", "LDRsimx", "1/LDRsimx", "LDRCorr"))
948 1034
949 #LDRtrueList = 0.004, 0.02, 0.2, 0.45 1035 #LDRtrueList = 0.004, 0.02, 0.2, 0.45
950 aF11sim0 = np.zeros(5) 1036 aF11sim0 = np.zeros(5)
951 LDRrange = np.zeros(5) 1037 LDRrange = np.zeros(5)
952 LDRrange = 0.004, 0.02, 0.1, 0.3, 0.45 1038 LDRrange = [0.004, 0.02, 0.1, 0.3, 0.45] # list
953 1039
954 # The loop over LDRtrueList is ony for checking how much the uncorrected LDRsimx deviates from LDRtrue ... and whether the corrections work. 1040 # The loop over LDRtrueList is only for checking how much the uncorrected LDRsimx deviates from LDRtrue ... and whether the corrections work.
955 # LDRsimx = LDRsim = Ir / It or 1/LDRsim 1041 # LDRsimx = LDRsim = Ir / It or 1/LDRsim
956 # Still with assumed true parameters in input file 1042 # Still with assumed true parameters in input file
957 for i, LDRtrue in enumerate(LDRrange): 1043 for i, LDRtrue in enumerate(LDRrange):
958 #for LDRtrue in LDRrange: 1044 #for LDRtrue in LDRrange:
959 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(DOLP0, RotL0, RotE0, RetE0, 1045 IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr, \
960 DiE0, RotO0, RetO0, 1046 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \
961 DiO0, RotC0, RetC0, 1047 Calc(TCalT0, TCalR0, NCalT, NCalR, DOLP0, RotL0, RotE0, RetE0, DiE0,
962 DiC0, TP0, TS0, RP0, 1048 RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0,
963 RS0, ERaT0, RotaT0, 1049 ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0)
964 RetT0, ERaR0, RotaR0, 1050 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))
965 RetR0, LDRCal0)
966 print("{0:10.5f},{1:10.5f},{2:10.5f},{3:10.5f}".format(LDRtrue, LDRsimx, 1/LDRsimx, LDRCorr))
967 aF11sim0[i] = F11sim0 1051 aF11sim0[i] = F11sim0
968 # the assumed true aF11sim0 results will be used below to calc the deviation from the real signals 1052 # the assumed true aF11sim0 results will be used below to calc the deviation from the real signals
969 print("Note: LDRsimx = LDR of the nominal system directly from measured signals without GHK-corrections") 1053 print("LDRsimx = LDR of the nominal system directly from measured signals without calibration and GHK-corrections")
970 1054 print("LDRunCorr = LDR of the nominal system directly from measured signals with calibration but without GHK-corrections; electronic amplifications = 1 assumed")
971 file = open('output_files\output_' + LID + '.dat', 'r') 1055 print("LDRCorr = LDR calibrated and GHK-corrected")
1056 print()
1057 print("Errors from signal noise:")
1058 print("Signal counts: NCalT, NCalR, NILfac, nNCal, nNI = {0:10.0f},{1:10.0f},{2:3.0f},{3:2.0f},{4:2.0f}".format(
1059 NCalT, NCalR, NILfac, nNCal, nNI))
1060
1061 '''# das muß wieder weg
1062 print("IoutTp, IoutTm, IoutRp, IoutRm, It , Ir , dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr")
1063 LDRCal = 0.01
1064 for i, LDRtrue in enumerate(LDRrange):
1065 IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr, \
1066 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \
1067 Calc(TCalT0, TCalR0, NCalT, NCalR, DOLP0, RotL0, RotE0, RetE0, DiE0,
1068 RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0,
1069 ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0)
1070 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(
1071 IoutTp * NCalT, IoutTm * NCalT, IoutRp * NCalR, IoutRm * NCalR, It * facIt, Ir * facIr,
1072 dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr))
1073 aF11sim0[i] = F11sim0
1074 # the assumed true aF11sim0 results will be used below to calc the deviation from the real signals
1075 # bis hierher weg
1076 '''
1077
1078 file = open('output_files\\' + LID + '-' + InputFile[0:-3] + '-GHK.dat', 'r')
972 print(file.read()) 1079 print(file.read())
973 file.close() 1080 file.close()
974 1081
975 ''' 1082 '''
976 if(PrintToOutputFile): 1083 if(PrintToOutputFile):
985 f.close 1092 f.close
986 sys.stdout = old_target 1093 sys.stdout = old_target
987 ''' 1094 '''
988 if (Error_Calc): 1095 if (Error_Calc):
989 # --- CALC again assumed truth with LDRCal0 and with assumed true parameters in input file to reset all 0-values 1096 # --- CALC again assumed truth with LDRCal0 and with assumed true parameters in input file to reset all 0-values
990 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(DOLP0, RotL0, RotE0, RetE0, DiE0, 1097 IoutTp0, IoutTm0, IoutRp0, IoutRm0, It0, Ir0, dIoutTp0, dIoutTm0, dIoutRp0, dIoutRm0, dIt0, dIr0, \
991 RotO0, RetO0, DiO0, RotC0, 1098 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \
992 RetC0, DiC0, TP0, TS0, RP0, 1099 Calc(TCalT0, TCalR0, NCalT, NCalR, DOLP0, RotL0, RotE0, RetE0, DiE0,
993 RS0, ERaT0, RotaT0, RetT0, 1100 RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0,
994 ERaR0, RotaR0, RetR0, 1101 ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0)
995 LDRCal0) 1102 Etax0 = K0 * Eta0
996 # --- Start Errors calculation with variable parameters ------------------------------------------------------------------ 1103 # --- Start Error calculation with variable parameters ------------------------------------------------------------------
1104 # error nNCal: one-sigma in steps to left and right for calibration signals
1105 # error nNI: one-sigma in steps to left and right for 0° signals
997 1106
998 iN = -1 1107 iN = -1
999 N = ((nDOLP * 2 + 1) * (nRotL * 2 + 1) * 1108 N = ((nTCalT * 2 + 1) * (nTCalR * 2 + 1) *
1109 (nNCal * 2 + 1) ** 4 * (nNI * 2 + 1) ** 2 *
1110 (nDOLP * 2 + 1) * (nRotL * 2 + 1) *
1000 (nRotE * 2 + 1) * (nRetE * 2 + 1) * (nDiE * 2 + 1) * 1111 (nRotE * 2 + 1) * (nRetE * 2 + 1) * (nDiE * 2 + 1) *
1001 (nRotO * 2 + 1) * (nRetO * 2 + 1) * (nDiO * 2 + 1) * 1112 (nRotO * 2 + 1) * (nRetO * 2 + 1) * (nDiO * 2 + 1) *
1002 (nRotC * 2 + 1) * (nRetC * 2 + 1) * (nDiC * 2 + 1) * 1113 (nRotC * 2 + 1) * (nRetC * 2 + 1) * (nDiC * 2 + 1) *
1003 (nTP * 2 + 1) * (nTS * 2 + 1) * (nRP * 2 + 1) * (nRS * 2 + 1) * (nERaT * 2 + 1) * (nERaR * 2 + 1) * 1114 (nTP * 2 + 1) * (nTS * 2 + 1) * (nRP * 2 + 1) * (nRS * 2 + 1) * (nERaT * 2 + 1) * (nERaR * 2 + 1) *
1004 (nRotaT * 2 + 1) * (nRotaR * 2 + 1) * (nRetT * 2 + 1) * (nRetR * 2 + 1) * (nLDRCal * 2 + 1)) 1115 (nRotaT * 2 + 1) * (nRotaR * 2 + 1) * (nRetT * 2 + 1) * (nRetR * 2 + 1) * (nLDRCal * 2 + 1))
1005 print("N = ", N, " ", end="") 1116 print("number of system variations N = ", N, " ", end="")
1006 1117
1007 if N > 1e6: 1118 if N > 1e6:
1008 if user_yes_no_query('Warning: processing ' + str( 1119 if user_yes_no_query('Warning: processing ' + str(
1009 N) + ' samples will take very long. Do you want to proceed?') == 0: sys.exit() 1120 N) + ' samples will take very long. Do you want to proceed?') == 0: sys.exit()
1010 if N > 5e6: 1121 if N > 5e6:
1016 # --- Arrays for plotting ------ 1127 # --- Arrays for plotting ------
1017 LDRmin = np.zeros(5) 1128 LDRmin = np.zeros(5)
1018 LDRmax = np.zeros(5) 1129 LDRmax = np.zeros(5)
1019 F11min = np.zeros(5) 1130 F11min = np.zeros(5)
1020 F11max = np.zeros(5) 1131 F11max = np.zeros(5)
1021 1132 Etaxmin = np.zeros(5)
1022 #LDRrange = np.zeros(5) 1133 Etaxmax = np.zeros(5)
1023 #LDRrange = 0.004, 0.02, 0.1, 0.3, 0.45 1134
1135 # LDRrange = np.zeros(5)
1136 # LDRrange = 0.004, 0.02, 0.1, 0.3, 0.45
1024 # aLDRsimx = np.zeros(N) 1137 # aLDRsimx = np.zeros(N)
1025 # aLDRsimx2 = np.zeros(N) 1138 # aLDRsimx2 = np.zeros(N)
1026 # aLDRcorr = np.zeros(N) 1139 # aLDRcorr = np.zeros(N)
1027 # aLDRcorr2 = np.zeros(N) 1140 # aLDRcorr2 = np.zeros(N)
1028 aDOLP = np.zeros(N) 1141 aDOLP = np.zeros(N)
1045 aRetE = np.zeros(N) 1158 aRetE = np.zeros(N)
1046 aRotE = np.zeros(N) 1159 aRotE = np.zeros(N)
1047 aRetO = np.zeros(N) 1160 aRetO = np.zeros(N)
1048 aRotO = np.zeros(N) 1161 aRotO = np.zeros(N)
1049 aLDRCal = np.zeros(N) 1162 aLDRCal = np.zeros(N)
1050 aA = np.zeros((5, N)) 1163 aNCalTp = np.zeros(N)
1051 aX = np.zeros((5, N)) 1164 aNCalTm = np.zeros(N)
1165 aNCalRp = np.zeros(N)
1166 aNCalRm = np.zeros(N)
1167 aNIt = np.zeros(N)
1168 aNIr = np.zeros(N)
1169 aTCalT = np.zeros(N)
1170 aTCalR = np.zeros(N)
1171
1172 # each np.zeros((LDRrange, N)) array has the same N-dependency
1173 aLDRcorr = np.zeros((5, N))
1052 aF11corr = np.zeros((5, N)) 1174 aF11corr = np.zeros((5, N))
1175 aPLDR = np.zeros((5, N))
1176 aEtax = np.zeros((5, N))
1177
1178 # np.zeros((GHKs, N))
1179 aGHK = np.zeros((5, N))
1053 1180
1054 atime = clock() 1181 atime = clock()
1055 dtime = clock() 1182 dtime = clock()
1056 1183
1057 # --- Calc Error signals 1184 # --- Calc Error signals
1064 atrue = (1. - LDRtrue) / (1 + LDRtrue) 1191 atrue = (1. - LDRtrue) / (1 + LDRtrue)
1065 atrue2 = (1. - LDRtrue2) / (1 + LDRtrue2) 1192 atrue2 = (1. - LDRtrue2) / (1 + LDRtrue2)
1066 ''' 1193 '''
1067 1194
1068 for iLDRCal in range(-nLDRCal, nLDRCal + 1): 1195 for iLDRCal in range(-nLDRCal, nLDRCal + 1):
1069 # from input file: assumed LDRCal for calibration measurements 1196 # from input file: LDRCal for calibration measurements
1070 LDRCal = LDRCal0 1197 LDRCal = LDRCal0
1071 if nLDRCal > 0: LDRCal = LDRCal0 + iLDRCal * dLDRCal / nLDRCal 1198 if nLDRCal > 0:
1072 1199 LDRCal = LDRCal0 + iLDRCal * dLDRCal / nLDRCal
1073 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim = Calc(DOLP0, RotL0, RotE0, RetE0, 1200 # provides the intensities of the calibration measurements at various LDRCal for signal noise errors
1074 DiE0, RotO0, RetO0, DiO0, 1201 # IoutTp, IoutTm, IoutRp, IoutRm, dIoutTp, dIoutTm, dIoutRp, dIoutRm
1075 RotC0, RetC0, DiC0, TP0, 1202 '''
1076 TS0, RP0, RS0, ERaT0, 1203 IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr, \
1077 RotaT0, RetT0, ERaR0, 1204 GT, HT, GR, HR, K, Eta, LDRsimx, LDRCorr, DTa, DRa, TTa, TRa, F11sim, LDRunCorr = \
1078 RotaR0, RetR0, LDRCal) 1205 Calc(TCalT, TCalR, NCalT, NCalR, DOLP0, RotL0, RotE0, RetE0, DiE0,
1206 RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0,
1207 ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal)
1208 '''
1079 aCal = (1. - LDRCal) / (1 + LDRCal) 1209 aCal = (1. - LDRCal) / (1 + LDRCal)
1080 for iDOLP, iRotL, iRotE, iRetE, iDiE \ 1210 for iDOLP, iRotL, iRotE, iRetE, iDiE \
1081 in [(iDOLP, iRotL, iRotE, iRetE, iDiE) 1211 in [(iDOLP, iRotL, iRotE, iRetE, iDiE)
1082 for iDOLP in range(-nDOLP, nDOLP + 1) 1212 for iDOLP in range(-nDOLP, nDOLP + 1)
1083 for iRotL in range(-nRotL, nRotL + 1) 1213 for iRotL in range(-nRotL, nRotL + 1)
1149 for iRetT in range(-nRetT, nRetT + 1) 1279 for iRetT in range(-nRetT, nRetT + 1)
1150 for iERaR in range(-nERaR, nERaR + 1) 1280 for iERaR in range(-nERaR, nERaR + 1)
1151 for iRotaR in range(-nRotaR, nRotaR + 1) 1281 for iRotaR in range(-nRotaR, nRotaR + 1)
1152 for iRetR in range(-nRetR, nRetR + 1)]: 1282 for iRetR in range(-nRetR, nRetR + 1)]:
1153 1283
1154 iN = iN + 1
1155 if (iN == 10001):
1156 ctime = clock()
1157 print(" estimated time ", "{0:4.2f}".format(N / 10000 * (ctime - atime)), "sec ") # , end="")
1158 print("\r elapsed time ", "{0:5.0f}".format((ctime - atime)), "sec ", end="\r")
1159 ctime = clock()
1160 if ((ctime - dtime) > 10):
1161 print("\r elapsed time ", "{0:5.0f}".format((ctime - atime)), "sec ", end="\r")
1162 dtime = ctime
1163
1164 if nRotO > 0: RotO = RotO0 + iRotO * dRotO / nRotO 1284 if nRotO > 0: RotO = RotO0 + iRotO * dRotO / nRotO
1165 if nRetO > 0: RetO = RetO0 + iRetO * dRetO / nRetO 1285 if nRetO > 0: RetO = RetO0 + iRetO * dRetO / nRetO
1166 if nDiO > 0: DiO = DiO0 + iDiO * dDiO / nDiO 1286 if nDiO > 0: DiO = DiO0 + iDiO * dDiO / nDiO
1167 if nRotC > 0: RotC = RotC0 + iRotC * dRotC / nRotC 1287 if nRotC > 0: RotC = RotC0 + iRotC * dRotC / nRotC
1168 if nRetC > 0: RetC = RetC0 + iRetC * dRetC / nRetC 1288 if nRetC > 0: RetC = RetC0 + iRetC * dRetC / nRetC
1207 CosT = np.cos(np.deg2rad(RetT)) 1327 CosT = np.cos(np.deg2rad(RetT))
1208 SinT = np.sin(np.deg2rad(RetT)) 1328 SinT = np.sin(np.deg2rad(RetT))
1209 CosR = np.cos(np.deg2rad(RetR)) 1329 CosR = np.cos(np.deg2rad(RetR))
1210 SinR = np.sin(np.deg2rad(RetR)) 1330 SinR = np.sin(np.deg2rad(RetR))
1211 1331
1332 # cleaning pol-filter
1212 DaT = (1 - ERaT) / (1 + ERaT) 1333 DaT = (1 - ERaT) / (1 + ERaT)
1213 DaR = (1 - ERaR) / (1 + ERaR) 1334 DaR = (1 - ERaR) / (1 + ERaR)
1214 TaT = 0.5 * (1 + ERaT) 1335 TaT = 0.5 * (1 + ERaT)
1215 TaR = 0.5 * (1 + ERaR) 1336 TaR = 0.5 * (1 + ERaR)
1216 1337
1217 S2aT = np.sin(np.deg2rad(h * 2 * RotaT)) 1338 S2aT = np.sin(np.deg2rad(h * 2 * RotaT))
1218 C2aT = np.cos(np.deg2rad(2 * RotaT)) 1339 C2aT = np.cos(np.deg2rad(2 * RotaT))
1219 S2aR = np.sin(np.deg2rad(h * 2 * RotaR)) 1340 S2aR = np.sin(np.deg2rad(h * 2 * RotaR))
1220 C2aR = np.cos(np.deg2rad(2 * RotaR)) 1341 C2aR = np.cos(np.deg2rad(2 * RotaR))
1221 1342
1222 # Aanalyzer As before the PBS Eq. D.5 1343 # Analyzer As before the PBS Eq. D.5; combined PBS and cleaning pol-filter
1223 ATP1 = (1 + C2aT * DaT * DiT) 1344 ATPT = (1 + C2aT * DaT * DiT) # unpolarized transmission correction
1224 ATP2 = Y * (DiT + C2aT * DaT) 1345 TTa = TiT * TaT * ATPT # unpolarized transmission
1225 ATP3 = Y * S2aT * DaT * ZiT * CosT 1346 ATP1 = 1
1226 ATP4 = S2aT * DaT * ZiT * SinT 1347 ATP2 = Y * (DiT + C2aT * DaT) / ATPT
1348 ATP3 = Y * S2aT * DaT * ZiT * CosT / ATPT
1349 ATP4 = S2aT * DaT * ZiT * SinT / ATPT
1227 ATP = np.array([ATP1, ATP2, ATP3, ATP4]) 1350 ATP = np.array([ATP1, ATP2, ATP3, ATP4])
1228 1351 DTa = ATP2 * Y
1229 ARP1 = (1 + C2aR * DaR * DiR) 1352
1230 ARP2 = Y * (DiR + C2aR * DaR) 1353 ARPT = (1 + C2aR * DaR * DiR) # unpolarized transmission correction
1231 ARP3 = Y * S2aR * DaR * ZiR * CosR 1354 TRa = TiR * TaR * ARPT # unpolarized transmission
1232 ARP4 = S2aR * DaR * ZiR * SinR 1355 ARP1 = 1
1356 ARP2 = Y * (DiR + C2aR * DaR) / ARPT
1357 ARP3 = Y * S2aR * DaR * ZiR * CosR / ARPT
1358 ARP4 = S2aR * DaR * ZiR * SinR / ARPT
1233 ARP = np.array([ARP1, ARP2, ARP3, ARP4]) 1359 ARP = np.array([ARP1, ARP2, ARP3, ARP4])
1234 1360 DRa = ARP2 * Y
1235 TTa = TiT * TaT # *ATP1
1236 TRa = TiR * TaR # *ARP1
1237 1361
1238 # ---- Calculate signals and correction parameters for diffeent locations and calibrators 1362 # ---- Calculate signals and correction parameters for diffeent locations and calibrators
1239 if LocC == 4: # Calibrator before the PBS 1363 if LocC == 4: # Calibrator before the PBS
1240 # print("Calibrator location not implemented yet") 1364 # print("Calibrator location not implemented yet")
1241 1365
1642 C = -ZiC * (SinC * UinEe - CosC * VinE) 1766 C = -ZiC * (SinC * UinEe - CosC * VinE)
1643 GT = ATE1o * D + ATE4o * C 1767 GT = ATE1o * D + ATE4o * C
1644 GR = ARE1o * D + ARE4o * C 1768 GR = ARE1o * D + ARE4o * C
1645 HT = ATE2a * A + ATE3a * B + ATE4a * C 1769 HT = ATE2a * A + ATE3a * B + ATE4a * C
1646 HR = ARE2a * A + ARE3a * B + ARE4a * C 1770 HR = ARE2a * A + ARE3a * B + ARE4a * C
1647
1648 else: 1771 else:
1649 print('Calibrator not implemented yet') 1772 print('Calibrator not implemented yet')
1650 sys.exit() 1773 sys.exit()
1651 1774
1652 # Calibration signals with aCal => Determination of the correction K of the real calibration factor 1775 for iTCalT, iTCalR, iNCalTp, iNCalTm, iNCalRp, iNCalRm, iNIt, iNIr \
1653 IoutTp = TaT * TiT * TiO * TiE * (AT + BT) 1776 in [
1654 IoutTm = TaT * TiT * TiO * TiE * (AT - BT) 1777 (iTCalT, iTCalR, iNCalTp, iNCalTm, iNCalRp, iNCalRm, iNIt, iNIr)
1655 IoutRp = TaR * TiR * TiO * TiE * (AR + BR) 1778 for iTCalT in range(-nTCalT, nTCalT + 1) # Etax
1656 IoutRm = TaR * TiR * TiO * TiE * (AR - BR) 1779 for iTCalR in range(-nTCalR, nTCalR + 1) # Etax
1657 # --- Results and Corrections; electronic etaR and etaT are assumed to be 1 1780 for iNCalTp in range(-nNCal, nNCal + 1) # noise error of calibration signals => Etax
1658 # Eta = TiR/TiT # Eta = Eta*/K Eq. 84 1781 for iNCalTm in range(-nNCal, nNCal + 1) # noise error of calibration signals => Etax
1659 Etapx = IoutRp / IoutTp 1782 for iNCalRp in range(-nNCal, nNCal + 1) # noise error of calibration signals => Etax
1660 Etamx = IoutRm / IoutTm 1783 for iNCalRm in range(-nNCal, nNCal + 1) # noise error of calibration signals => Etax
1661 Etax = (Etapx * Etamx) ** 0.5 1784 for iNIt in range(-nNI, nNI + 1)
1662 K = Etax / Eta0 1785 for iNIr in range(-nNI, nNI + 1)]:
1663 # 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)) 1786
1664 # print("{0:6.3f},{1:6.3f},{2:6.3f},{3:6.3f}".format(DiC, ZiC, Kp, Km)) 1787 # Calibration signals with aCal => Determination of the correction K of the real calibration factor
1665 1788 IoutTp = TTa * TiC * TiO * TiE * (AT + BT)
1666 # For comparison with Volkers Libreoffice Müller Matrix spreadsheet 1789 IoutTm = TTa * TiC * TiO * TiE * (AT - BT)
1667 # Eta_test_p = (IoutRp/IoutTp) 1790 IoutRp = TRa * TiC * TiO * TiE * (AR + BR)
1668 # Eta_test_m = (IoutRm/IoutTm) 1791 IoutRm = TRa * TiC * TiO * TiE * (AR - BR)
1669 # Eta_test = (Eta_test_p*Eta_test_m)**0.5 1792
1670 1793 if nTCalT > 0: TCalT = TCalT0 + iTCalT * dTCalT / nTCalT
1671 # ************************************************************************* 1794 if nTCalR > 0: TCalR = TCalR0 + iTCalR * dTCalR / nTCalR
1672 iLDR = -1 1795 # signal noise errors
1673 for LDRTrue in LDRrange: 1796 # ----- random error calculation ----------
1674 iLDR = iLDR + 1 1797 # noise must be calculated from/with the actually measured signals; influence of TCalT, TCalR errors on nouse are not considered ?
1675 atrue = (1 - LDRTrue) / (1 + LDRTrue) 1798 # actually measured signals are in input file and don't change
1676 # ----- Forward simulated signals and LDRsim with atrue; from input file 1799 # relative standard deviation of calibration signals with LDRcal; assumed to be statisitcally independent
1677 It = TaT * TiT * TiO * TiE * (GT + atrue * HT) # TaT*TiT*TiC*TiO*IinL*(GT+atrue*HT) 1800 # error nNCal: one-sigma in steps to left and right for calibration signals
1678 Ir = TaR * TiR * TiO * TiE * (GR + atrue * HR) # TaR*TiR*TiC*TiO*IinL*(GR+atrue*HR) 1801 if nNCal > 0:
1679 1802 if (CalcFrom0deg):
1680 # LDRsim = 1/Eta*Ir/It # simulated LDR* with Y from input file 1803 dIoutTp = (NCalT * IoutTp) ** -0.5
1681 LDRsim = Ir / It # simulated uncorrected LDR with Y from input file 1804 dIoutTm = (NCalT * IoutTm) ** -0.5
1805 dIoutRp = (NCalR * IoutRp) ** -0.5
1806 dIoutRm = (NCalR * IoutRm) ** -0.5
1807 else:
1808 dIoutTp = dIoutTp0 * (IoutTp / IoutTp0)
1809 dIoutTm = dIoutTm0 * (IoutTm / IoutTm0)
1810 dIoutRp = dIoutRp0 * (IoutRp / IoutRp0)
1811 dIoutRm = dIoutRm0 * (IoutRm / IoutRm0)
1812 # print(iTCalT, iTCalR, iNCalTp, iNCalTm, iNCalRp, iNCalRm, iNIt, iNIr, IoutTp, dIoutTp)
1813 IoutTp = IoutTp * (1 + iNCalTp * dIoutTp / nNCal)
1814 IoutTm = IoutTm * (1 + iNCalTm * dIoutTm / nNCal)
1815 IoutRp = IoutRp * (1 + iNCalRp * dIoutRp / nNCal)
1816 IoutRm = IoutRm * (1 + iNCalRm * dIoutRm / nNCal)
1817
1818 IoutTp = IoutTp * TCalT / TCalT0
1819 IoutTm = IoutTm * TCalT / TCalT0
1820 IoutRp = IoutRp * TCalR / TCalR0
1821 IoutRm = IoutRm * TCalR / TCalR0
1822 # --- Results and Corrections; electronic etaR and etaT are assumed to be 1 for true and assumed true systems
1823 # calibration factor
1824 Eta = (TRa / TTa) # = TRa / TTa; Eta = Eta*/K Eq. 84; corrected according to the papers supplement Eqs. (S.10.10.1) ff
1825 # possibly real calibration factor
1826 Etapx = IoutRp / IoutTp
1827 Etamx = IoutRm / IoutTm
1828 Etax = (Etapx * Etamx) ** 0.5
1829 K = Etax / Eta
1830 # 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))
1831 # print("{0:6.3f},{1:6.3f},{2:6.3f},{3:6.3f}".format(DiC, ZiC, Kp, Km))
1832
1833 # For comparison with Volkers Libreoffice Müller Matrix spreadsheet
1834 # Eta_test_p = (IoutRp/IoutTp)
1835 # Eta_test_m = (IoutRm/IoutTm)
1836 # Eta_test = (Eta_test_p*Eta_test_m)**0.5
1682 ''' 1837 '''
1683 if Y == 1.: 1838 for iIt, iIr \
1684 LDRsimx = LDRsim 1839 in [(iIt, iIr)
1685 LDRsimx2 = LDRsim2 1840 for iIt in range(-nNI, nNI + 1)
1686 else: 1841 for iIr in range(-nNI, nNI + 1)]:
1687 LDRsimx = 1./LDRsim
1688 LDRsimx2 = 1./LDRsim2
1689 ''' 1842 '''
1690 # ----- Backward correction 1843
1691 # Corrected LDRCorr with assumed true G0,H0,K0 from forward simulated (real) LDRsim (atrue) 1844 iN = iN + 1
1692 LDRCorr = (LDRsim * K0 / Etax * (GT0 + HT0) - (GR0 + HR0)) / ( 1845 if (iN == 10001):
1693 (GR0 - HR0) - LDRsim * K0 / Etax * (GT0 - HT0)) 1846 ctime = clock()
1694 ''' 1847 print(" estimated time ", "{0:4.2f}".format(N / 10000 * (ctime - atime)), "sec ") # , end="")
1695 # -- F11corr from It and Ir and calibration EtaX 1848 print("\r elapsed time ", "{0:5.0f}".format((ctime - atime)), "sec ", end="\r")
1696 Text1 = "!!! EXPERIMENTAL !!! F11corr from It and Ir with calibration EtaX: x-axis: F11corr(LDRtrue) / F11corr(LDRtrue = 0.004) - 1" 1849 ctime = clock()
1697 F11corr = 1 / (TiO * TiE) * ( 1850 if ((ctime - dtime) > 10):
1698 (HR0 * Etax / K0 * It / TTa - HT0 * Ir / TRa) / (HR0 * GT0 - HT0 * GR0)) # IL = 1 Eq.(64); Etax/K0 = Eta0. 1851 print("\r elapsed time ", "{0:5.0f}".format((ctime - atime)), "sec ", end="\r")
1699 ''' 1852 dtime = ctime
1700 # Corrected F11corr with assumed true G0,H0,K0 from forward simulated (real) It and Ir (atrue) 1853
1701 Text1 = "!!! EXPERIMENTAL !!! F11corr from real It and Ir with real calibration EtaX: x-axis: F11corr(LDRtrue) / aF11sim0(LDRtrue) - 1" 1854 # *** loop for different real LDRs **********************************************************************
1702 F11corr = 1 / (TiO * TiE) * ( 1855 iLDR = -1
1703 (HR0 * Etax / K0 * It / TTa - HT0 * Ir / TRa) / (HR0 * GT0 - HT0 * GR0)) # IL = 1 Eq.(64); Etax/K0 = Eta0. 1856 for LDRTrue in LDRrange:
1704 1857 iLDR = iLDR + 1
1705 # Text1 = "F11corr from It and Ir without corrections but with calibration EtaX: x-axis: F11corr(LDRtrue) devided by F11corr(LDRtrue = 0.004)" 1858 atrue = (1 - LDRTrue) / (1 + LDRTrue)
1706 # F11corr = 0.5/(TiO*TiE)*(Etax*It/TTa+Ir/TRa) # IL = 1 Eq.(64) 1859 # ----- Forward simulated signals and LDRsim with atrue; from input file; not considering TiC.
1707 1860 It = TTa * TiO * TiE * (GT + atrue * HT) # TaT*TiT*TiC*TiO*IinL*(GT+atrue*HT)
1708 # -- It from It only with atrue without corrections - for BERTHA (and PollyXTs) 1861 Ir = TRa * TiO * TiE * (GR + atrue * HR) # TaR*TiR*TiC*TiO*IinL*(GR+atrue*HR)
1709 # Text1 = " x-axis: IT(LDRtrue) / IT(LDRtrue = 0.004) - 1" 1862 # # signal noise errors; standard deviation of signals; assumed to be statisitcally independent
1710 # F11corr = It/(TaT*TiT*TiO*TiE) #/(TaT*TiT*TiO*TiE*(GT0+atrue*HT0)) 1863 # because the signals depend on LDRtrue, the errors dIt and dIr must be calculated for each LDRtrue
1711 # !!! see below line 1673ff 1864 if (CalcFrom0deg):
1712 1865 dIt = ((NCalT * It / IoutTp * NILfac / TCalT) ** -0.5)
1713 aF11corr[iLDR, iN] = F11corr 1866 dIr = ((NCalR * Ir / IoutRp * NILfac / TCalR) ** -0.5)
1714 aA[iLDR, iN] = LDRCorr # LDRCorr # LDRsim # for test only 1867 else:
1715 1868 dIt = ((NCalT * 2 * NILfac / TCalT ) ** -0.5) * It
1716 aX[0, iN] = GR 1869 dIr = ((NCalR * 2 * NILfac / TCalR) ** -0.5) * Ir
1717 aX[1, iN] = GT 1870 # error nNI: one-sigma in steps to left and right for 0° signals
1718 aX[2, iN] = HR 1871 if nNI > 0:
1719 aX[3, iN] = HT 1872 It = It * (1 + iNIt * dIt / nNI)
1720 aX[4, iN] = K 1873 Ir = Ir * (1 + iNIr * dIr / nNI)
1721 1874
1722 aLDRCal[iN] = iLDRCal 1875 # LDRsim = 1/Eta*Ir/It # simulated LDR* with Y from input file
1723 aDOLP[iN] = iDOLP 1876 LDRsim = Ir / It # simulated uncorrected LDR with Y from input file
1724 aERaT[iN] = iERaT 1877
1725 aERaR[iN] = iERaR 1878 # ----- Backward correction
1726 aRotaT[iN] = iRotaT 1879 # Corrected LDRCorr with assumed true G0,H0,K0,Eta0 from forward simulated (real) LDRsim(atrue)
1727 aRotaR[iN] = iRotaR 1880 LDRCorr = (LDRsim / (Etax / K0) * (GT0 + HT0) - (GR0 + HR0)) / ((GR0 - HR0) - LDRsim / (Etax / K0) * (GT0 - HT0))
1728 aRetT[iN] = iRetT 1881
1729 aRetR[iN] = iRetR 1882 # The following is a test whether the equations for calibration Etax and normal signal (GHK, LDRsim) are consistent
1730 1883 # LDRCorr = (LDRsim / Eta * (GT + HT) - (GR + HR)) / ((GR - HR) - LDRsim / Eta * (GT - HT))
1731 aRotL[iN] = iRotL 1884 # Without any correction
1732 aRotE[iN] = iRotE 1885 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)))
1733 aRetE[iN] = iRetE 1886
1734 aRotO[iN] = iRotO 1887
1735 aRetO[iN] = iRetO 1888 '''
1736 aRotC[iN] = iRotC 1889 # -- F11corr from It and Ir and calibration EtaX
1737 aRetC[iN] = iRetC 1890 Text1 = "!!! EXPERIMENTAL !!! F11corr from It and Ir with calibration EtaX: x-axis: F11corr(LDRtrue) / F11corr(LDRtrue = 0.004) - 1"
1738 aDiO[iN] = iDiO 1891 F11corr = 1 / (TiO * TiE) * (
1739 aDiE[iN] = iDiE 1892 (HR0 * Etax / K0 * It / TTa - HT0 * Ir / TRa) / (HR0 * GT0 - HT0 * GR0)) # IL = 1 Eq.(64); Etax/K0 = Eta0.
1740 aDiC[iN] = iDiC 1893 '''
1741 aTP[iN] = iTP 1894 # Corrected F11corr with assumed true G0,H0,K0 from forward simulated (real) It and Ir (atrue)
1742 aTS[iN] = iTS 1895 Text1 = "!!! EXPERIMENTAL !!! F11corr from real It and Ir with real calibration EtaX: x-axis: F11corr(LDRtrue) / aF11sim0(LDRtrue) - 1"
1743 aRP[iN] = iRP 1896 F11corr = 1 / (TiO * TiE) * (
1744 aRS[iN] = iRS 1897 (HR0 * Etax / K0 * It / TTa - HT0 * Ir / TRa) / (HR0 * GT0 - HT0 * GR0)) # IL = 1 Eq.(64); Etax/K0 = Eta0.
1898
1899 # Text1 = "F11corr from It and Ir without corrections but with calibration EtaX: x-axis: F11corr(LDRtrue) devided by F11corr(LDRtrue = 0.004)"
1900 # F11corr = 0.5/(TiO*TiE)*(Etax*It/TTa+Ir/TRa) # IL = 1 Eq.(64)
1901
1902 # -- It from It only with atrue without corrections - for BERTHA (and PollyXTs)
1903 # Text1 = " x-axis: IT(LDRtrue) / IT(LDRtrue = 0.004) - 1"
1904 # F11corr = It/(TaT*TiT*TiO*TiE) #/(TaT*TiT*TiO*TiE*(GT0+atrue*HT0))
1905 # ! see below line 1673ff
1906
1907 aF11corr[iLDR, iN] = F11corr
1908 aLDRcorr[iLDR, iN] = LDRCorr # LDRCorr # LDRsim # for test only
1909 # aPLDR[iLDR, iN] = CalcPLDR(LDRCorr, BSR[iLDR], LDRm0)
1910 aEtax[iLDR, iN] = Etax
1911
1912 aGHK[0, iN] = GR
1913 aGHK[1, iN] = GT
1914 aGHK[2, iN] = HR
1915 aGHK[3, iN] = HT
1916 aGHK[4, iN] = K
1917
1918 aLDRCal[iN] = iLDRCal
1919 aDOLP[iN] = iDOLP
1920 aERaT[iN] = iERaT
1921 aERaR[iN] = iERaR
1922 aRotaT[iN] = iRotaT
1923 aRotaR[iN] = iRotaR
1924 aRetT[iN] = iRetT
1925 aRetR[iN] = iRetR
1926
1927 aRotL[iN] = iRotL
1928 aRotE[iN] = iRotE
1929 aRetE[iN] = iRetE
1930 aRotO[iN] = iRotO
1931 aRetO[iN] = iRetO
1932 aRotC[iN] = iRotC
1933 aRetC[iN] = iRetC
1934 aDiO[iN] = iDiO
1935 aDiE[iN] = iDiE
1936 aDiC[iN] = iDiC
1937 aTP[iN] = iTP
1938 aTS[iN] = iTS
1939 aRP[iN] = iRP
1940 aRS[iN] = iRS
1941 aTCalT[iN] = iTCalT
1942 aTCalR[iN] = iTCalR
1943
1944 aNCalTp[iN] = iNCalTp # IoutTp, IoutTm, IoutRp, IoutRm => Etax
1945 aNCalTm[iN] = iNCalTm # IoutTp, IoutTm, IoutRp, IoutRm => Etax
1946 aNCalRp[iN] = iNCalRp # IoutTp, IoutTm, IoutRp, IoutRm => Etax
1947 aNCalRm[iN] = iNCalRm # IoutTp, IoutTm, IoutRp, IoutRm => Etax
1948 aNIt[iN] = iNIt # It, Tr
1949 aNIr[iN] = iNIr # It, Tr
1745 1950
1746 # --- END loop 1951 # --- END loop
1747 btime = clock() 1952 btime = clock()
1748 print("\r done in ", "{0:5.0f}".format(btime - atime), "sec") # , end="\r") 1953 # print("\r done in ", "{0:5.0f}".format(btime - atime), "sec. => producing plots now .... some more seconds ..."), # , end="\r");
1749 1954 print(" done in ", "{0:5.0f}".format(btime - atime), "sec. => producing plots now .... some more seconds ...")
1750 # --- Plot ----------------------------------------------------------------- 1955 # --- Plot -----------------------------------------------------------------
1956 print("Errors from GHK correction uncertainties:")
1751 if (sns_loaded): 1957 if (sns_loaded):
1752 sns.set_style("whitegrid") 1958 sns.set_style("whitegrid")
1753 sns.set_palette("bright", 6) 1959 sns.set_palette("bright6", 6)
1960 # for older seaborn versions:
1961 # sns.set_palette("bright", 6)
1754 1962
1755 ''' 1963 '''
1756 fig2 = plt.figure() 1964 fig2 = plt.figure()
1757 plt.plot(aA[2,:],'b.') 1965 plt.plot(aLDRcorr[2,:],'b.')
1758 plt.plot(aA[3,:],'r.') 1966 plt.plot(aLDRcorr[3,:],'r.')
1759 plt.plot(aA[4,:],'g.') 1967 plt.plot(aLDRcorr[4,:],'g.')
1760 #plt.plot(aA[6,:],'c.') 1968 #plt.plot(aLDRcorr[6,:],'c.')
1761 plt.show 1969 plt.show
1762 ''' 1970 '''
1763 1971
1764
1765 # Plot LDR 1972 # Plot LDR
1766 def PlotSubHist(aVar, aX, X0, daX, iaX, naX): 1973 def PlotSubHist(aVar, aX, X0, daX, iaX, naX):
1974 # aVar is the name of the parameter and aX is the subset of aLDRcorr which is coloured in the plot
1975 # example: PlotSubHist("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP)
1767 fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) 1976 fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2))
1768 iLDR = -1 1977 iLDR = -1
1769 for LDRTrue in LDRrange: 1978 for LDRTrue in LDRrange:
1770 iLDR = iLDR + 1 1979 iLDR = iLDR + 1
1771 1980 LDRmin[iLDR] = np.amin(aLDRcorr[iLDR, :])
1772 LDRmin[iLDR] = np.min(aA[iLDR, :]) 1981 LDRmax[iLDR] = np.amax(aLDRcorr[iLDR, :])
1773 LDRmax[iLDR] = np.max(aA[iLDR, :]) 1982 Rmin = LDRmin[iLDR] * 0.995 # np.min(aLDRcorr[iLDR,:]) * 0.995
1774 Rmin = LDRmin[iLDR] * 0.995 # np.min(aA[iLDR,:]) * 0.995 1983 Rmax = LDRmax[iLDR] * 1.005 # np.max(aLDRcorr[iLDR,:]) * 1.005
1775 Rmax = LDRmax[iLDR] * 1.005 # np.max(aA[iLDR,:]) * 1.005
1776 1984
1777 # plt.subplot(5,2,iLDR+1) 1985 # plt.subplot(5,2,iLDR+1)
1778 plt.subplot(1, 5, iLDR + 1) 1986 plt.subplot(1, 5, iLDR + 1)
1779 (n, bins, patches) = plt.hist(aA[iLDR, :], 1987 (n, bins, patches) = plt.hist(aLDRcorr[iLDR, :],
1780 bins=100, log=False, 1988 bins=100, log=False,
1781 range=[Rmin, Rmax], 1989 range=[Rmin, Rmax],
1782 alpha=0.5, normed=False, color='0.5', histtype='stepfilled') 1990 alpha=0.5, density=False, color='0.5', histtype='stepfilled')
1783 1991
1784 for iaX in range(-naX, naX + 1): 1992 for iaX in range(-naX, naX + 1):
1785 plt.hist(aA[iLDR, aX == iaX], 1993 plt.hist(aLDRcorr[iLDR, aX == iaX],
1786 range=[Rmin, Rmax], 1994 range=[Rmin, Rmax],
1787 bins=100, log=False, alpha=0.3, normed=False, histtype='stepfilled', 1995 bins=100, log=False, alpha=0.3, density=False, histtype='stepfilled',
1788 label=str(round(X0 + iaX * daX / naX, 5))) 1996 label=str(round(X0 + iaX * daX / naX, 5)))
1789 1997
1790 if (iLDR == 2): plt.legend() 1998 if (iLDR == 2):
1999 leg = plt.legend()
2000 leg.get_frame().set_alpha(0.1)
2001
1791 2002
1792 plt.tick_params(axis='both', labelsize=9) 2003 plt.tick_params(axis='both', labelsize=9)
1793 plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2) 2004 plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2)
1794 2005
1795 # plt.title(LID + ' ' + aVar, fontsize=18) 2006 # plt.title(LID + ' ' + aVar, fontsize=18)
1800 # plt.show() 2011 # plt.show()
1801 # fig.savefig(LID + '_' + aVar + '.png', dpi=150, bbox_inches='tight', pad_inches=0) 2012 # fig.savefig(LID + '_' + aVar + '.png', dpi=150, bbox_inches='tight', pad_inches=0)
1802 # plt.close 2013 # plt.close
1803 return 2014 return
1804 2015
2016 # Plot Etax
2017 def PlotEtax(aVar, aX, X0, daX, iaX, naX):
2018 # aVar is the name of the parameter and aX is the subset of aLDRcorr which is coloured in the plot
2019 # example: PlotSubHist("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP)
2020 fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2))
2021 iLDR = -1
2022 for LDRTrue in LDRrange:
2023 iLDR = iLDR + 1
2024 Etaxmin[iLDR] = np.amin(aEtax[iLDR, :])
2025 Etaxmax[iLDR] = np.amax(aEtax[iLDR, :])
2026 Rmin = Etaxmin[iLDR] * 0.995 # np.min(aLDRcorr[iLDR,:]) * 0.995
2027 Rmax = Etaxmax[iLDR] * 1.005 # np.max(aLDRcorr[iLDR,:]) * 1.005
2028
2029 # plt.subplot(5,2,iLDR+1)
2030 plt.subplot(1, 5, iLDR + 1)
2031 (n, bins, patches) = plt.hist(aEtax[iLDR, :],
2032 bins=100, log=False,
2033 range=[Rmin, Rmax],
2034 alpha=0.5, density=False, color='0.5', histtype='stepfilled')
2035 for iaX in range(-naX, naX + 1):
2036 plt.hist(aEtax[iLDR, aX == iaX],
2037 range=[Rmin, Rmax],
2038 bins=100, log=False, alpha=0.3, density=False, histtype='stepfilled',
2039 label=str(round(X0 + iaX * daX / naX, 5)))
2040 if (iLDR == 2):
2041 leg = plt.legend()
2042 leg.get_frame().set_alpha(0.1)
2043 plt.tick_params(axis='both', labelsize=9)
2044 plt.plot([Etax0, Etax0], [0, np.max(n)], 'r-', lw=2)
2045 fig.suptitle('Etax - ' + LID + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar, fontsize=14, y=1.05)
2046 return
1805 2047
1806 if (nDOLP > 0): PlotSubHist("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP) 2048 if (nDOLP > 0): PlotSubHist("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP)
1807 if (nRotL > 0): PlotSubHist("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) 2049 if (nRotL > 0): PlotSubHist("RotL", aRotL, RotL0, dRotL, iRotL, nRotL)
1808 if (nRetE > 0): PlotSubHist("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) 2050 if (nRetE > 0): PlotSubHist("RetE", aRetE, RetE0, dRetE, iRetE, nRetE)
1809 if (nRotE > 0): PlotSubHist("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) 2051 if (nRotE > 0): PlotSubHist("RotE", aRotE, RotE0, dRotE, iRotE, nRotE)
1823 if (nERaT > 0): PlotSubHist("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) 2065 if (nERaT > 0): PlotSubHist("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT)
1824 if (nERaR > 0): PlotSubHist("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) 2066 if (nERaR > 0): PlotSubHist("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR)
1825 if (nRotaT > 0): PlotSubHist("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) 2067 if (nRotaT > 0): PlotSubHist("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT)
1826 if (nRotaR > 0): PlotSubHist("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) 2068 if (nRotaR > 0): PlotSubHist("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR)
1827 if (nLDRCal > 0): PlotSubHist("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) 2069 if (nLDRCal > 0): PlotSubHist("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal)
1828 2070 if (nTCalT > 0): PlotSubHist("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT)
2071 if (nTCalR > 0): PlotSubHist("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR)
2072 if (nNCal > 0): PlotSubHist("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal)
2073 if (nNCal > 0): PlotSubHist("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal)
2074 if (nNCal > 0): PlotSubHist("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal)
2075 if (nNCal > 0): PlotSubHist("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal)
2076 if (nNI > 0): PlotSubHist("SigNoiseIt", aNIt, 0, 1, iNIt, nNI)
2077 if (nNI > 0): PlotSubHist("SigNoiseIr", aNIr, 0, 1, iNIr, nNI)
1829 plt.show() 2078 plt.show()
1830 plt.close 2079 plt.close
1831 2080
1832 ''' 2081 '''
1833 # --- Plot F11 histograms 2082 # --- Plot F11 histograms
1848 for LDRTrue in LDRrange: 2097 for LDRTrue in LDRrange:
1849 iLDR = iLDR + 1 2098 iLDR = iLDR + 1
1850 2099
1851 #F11min[iLDR] = np.min(aF11corr[iLDR,:]) 2100 #F11min[iLDR] = np.min(aF11corr[iLDR,:])
1852 #F11max[iLDR] = np.max(aF11corr[iLDR,:]) 2101 #F11max[iLDR] = np.max(aF11corr[iLDR,:])
1853 #Rmin = F11min[iLDR] * 0.995 # np.min(aA[iLDR,:]) * 0.995 2102 #Rmin = F11min[iLDR] * 0.995 # np.min(aLDRcorr[iLDR,:]) * 0.995
1854 #Rmax = F11max[iLDR] * 1.005 # np.max(aA[iLDR,:]) * 1.005 2103 #Rmax = F11max[iLDR] * 1.005 # np.max(aLDRcorr[iLDR,:]) * 1.005
1855 2104
1856 #Rmin = 0.8 2105 #Rmin = 0.8
1857 #Rmax = 1.2 2106 #Rmax = 1.2
1858 2107
1859 #plt.subplot(5,2,iLDR+1) 2108 #plt.subplot(5,2,iLDR+1)
1860 plt.subplot(1,5,iLDR+1) 2109 plt.subplot(1,5,iLDR+1)
1861 (n, bins, patches) = plt.hist(aF11corr[iLDR,:], 2110 (n, bins, patches) = plt.hist(aF11corr[iLDR,:],
1862 bins=100, log=False, 2111 bins=100, log=False,
1863 alpha=0.5, normed=False, color = '0.5', histtype='stepfilled') 2112 alpha=0.5, density=False, color = '0.5', histtype='stepfilled')
1864 2113
1865 for iaX in range(-naX,naX+1): 2114 for iaX in range(-naX,naX+1):
1866 plt.hist(aF11corr[iLDR,aX == iaX], 2115 plt.hist(aF11corr[iLDR,aX == iaX],
1867 bins=100, log=False, alpha=0.3, normed=False, histtype='stepfilled', label = str(round(X0 + iaX*daX/naX,5))) 2116 bins=100, log=False, alpha=0.3, density=False, histtype='stepfilled', label = str(round(X0 + iaX*daX/naX,5)))
1868 2117
1869 if (iLDR == 2): plt.legend() 2118 if (iLDR == 2): plt.legend()
1870 2119
1871 plt.tick_params(axis='both', labelsize=9) 2120 plt.tick_params(axis='both', labelsize=9)
1872 #plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2) 2121 #plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2)
1901 if (nERaT > 0): PlotSubHistF11("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) 2150 if (nERaT > 0): PlotSubHistF11("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT)
1902 if (nERaR > 0): PlotSubHistF11("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) 2151 if (nERaR > 0): PlotSubHistF11("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR)
1903 if (nRotaT > 0): PlotSubHistF11("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) 2152 if (nRotaT > 0): PlotSubHistF11("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT)
1904 if (nRotaR > 0): PlotSubHistF11("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) 2153 if (nRotaR > 0): PlotSubHistF11("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR)
1905 if (nLDRCal > 0): PlotSubHistF11("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) 2154 if (nLDRCal > 0): PlotSubHistF11("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal)
2155 if (nTCalT > 0): PlotSubHistF11("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT)
2156 if (nTCalR > 0): PlotSubHistF11("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR)
2157 if (nNCal > 0): PlotSubHistF11("CalNoise", aNCal, 0, 1/nNCal, iNCal, nNCal)
2158 if (nNI > 0): PlotSubHistF11("SigNoise", aNI, 0, 1/nNI, iNI, nNI)
2159
1906 2160
1907 plt.show() 2161 plt.show()
1908 plt.close 2162 plt.close
1909 2163
1910 ''' 2164 '''
1913 #print("******************* " + aVar + " *******************") 2167 #print("******************* " + aVar + " *******************")
1914 fig, ax = plt.subplots(nrows=5, ncols=2, sharex=True, sharey=True, figsize=(10, 10)) 2168 fig, ax = plt.subplots(nrows=5, ncols=2, sharex=True, sharey=True, figsize=(10, 10))
1915 iLDR = -1 2169 iLDR = -1
1916 for LDRTrue in LDRrange: 2170 for LDRTrue in LDRrange:
1917 iLDR = iLDR + 1 2171 iLDR = iLDR + 1
1918 LDRmin[iLDR] = np.min(aA[iLDR,:]) 2172 LDRmin[iLDR] = np.min(aLDRcorr[iLDR,:])
1919 LDRmax[iLDR] = np.max(aA[iLDR,:]) 2173 LDRmax[iLDR] = np.max(aLDRcorr[iLDR,:])
1920 Rmin = np.min(aA[iLDR,:]) * 0.999 2174 Rmin = np.min(aLDRcorr[iLDR,:]) * 0.999
1921 Rmax = np.max(aA[iLDR,:]) * 1.001 2175 Rmax = np.max(aLDRcorr[iLDR,:]) * 1.001
1922 plt.subplot(5,2,iLDR+1) 2176 plt.subplot(5,2,iLDR+1)
1923 (n, bins, patches) = plt.hist(aA[iLDR,:], 2177 (n, bins, patches) = plt.hist(aLDRcorr[iLDR,:],
1924 range=[Rmin, Rmax], 2178 range=[Rmin, Rmax],
1925 bins=200, log=False, alpha=0.2, normed=False, color = '0.5', histtype='stepfilled') 2179 bins=200, log=False, alpha=0.2, density=False, color = '0.5', histtype='stepfilled')
1926 plt.tick_params(axis='both', labelsize=9) 2180 plt.tick_params(axis='both', labelsize=9)
1927 plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2) 2181 plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2)
1928 plt.show() 2182 plt.show()
1929 plt.close 2183 plt.close
1930 # --- End of Plot F11 histograms 2184 # --- End of Plot F11 histograms
1931 ''' 2185 '''
1932 2186
1933 # --- Plot LDRmin, LDRmax 2187 # --- Plot LDRmin, LDRmax
2188 iLDR = -1
2189 for LDRTrue in LDRrange:
2190 iLDR = iLDR + 1
2191 LDRmin[iLDR] = np.amin(aLDRcorr[iLDR, :])
2192 LDRmax[iLDR] = np.amax(aLDRcorr[iLDR, :])
2193
1934 fig2 = plt.figure() 2194 fig2 = plt.figure()
1935 plt.plot(LDRrange, LDRmax - LDRrange, linewidth=2.0, color='b') 2195 LDRrangeA = np.array(LDRrange)
1936 plt.plot(LDRrange, LDRmin - LDRrange, linewidth=2.0, color='g') 2196 if((np.amax(LDRmax - LDRrangeA)-np.amin(LDRmin - LDRrangeA)) < 0.001):
2197 plt.ylim(-0.001,0.001)
2198 plt.plot(LDRrangeA, LDRmax - LDRrangeA, linewidth=2.0, color='b')
2199 plt.plot(LDRrangeA, LDRmin - LDRrangeA, linewidth=2.0, color='g')
1937 2200
1938 plt.xlabel('LDRtrue', fontsize=18) 2201 plt.xlabel('LDRtrue', fontsize=18)
1939 plt.ylabel('LDRTrue-LDRmin, LDRTrue-LDRmax', fontsize=14) 2202 plt.ylabel('LDRTrue-LDRmin, LDRTrue-LDRmax', fontsize=14)
1940 plt.title(LID + ' ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]), fontsize=18) 2203 plt.title(LID + ' ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]), fontsize=18)
1941 # plt.ylimit(-0.07, 0.07) 2204 # plt.ylimit(-0.07, 0.07)
1942 plt.show() 2205 plt.show()
1943 plt.close 2206 plt.close
1944 2207
1945 # --- Save LDRmin, LDRmax to file 2208 # --- Save LDRmin, LDRmax to file
1946 # http://stackoverflow.com/questions/4675728/redirect-stdout-to-a-file-in-python 2209 # http://stackoverflow.com/questions/4675728/redirect-stdout-to-a-file-in-python
1947 with open('output_files\LDR_min_max_' + LID + '.dat', 'w') as f: 2210 with open('output_files\\' + LID + '-' + InputFile[0:-3] + '-LDR_min_max.dat', 'w') as f:
1948 with redirect_stdout(f): 2211 with redirect_stdout(f):
1949 print(LID) 2212 print(LID)
1950 print("LDRtrue, LDRmin, LDRmax") 2213 print("LDRtrue, LDRmin, LDRmax")
1951 for i in range(len(LDRrange)): 2214 for i in range(len(LDRrangeA)):
1952 print("{0:7.4f},{1:7.4f},{2:7.4f}".format(LDRrange[i], LDRmin[i], LDRmax[i])) 2215 print("{0:7.4f},{1:7.4f},{2:7.4f}".format(LDRrangeA[i], LDRmin[i], LDRmax[i]))
1953 2216
2217
2218 if (bPlotEtax):
2219 if (nDOLP > 0): PlotEtax("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP)
2220 if (nRotL > 0): PlotEtax("RotL", aRotL, RotL0, dRotL, iRotL, nRotL)
2221 if (nRetE > 0): PlotEtax("RetE", aRetE, RetE0, dRetE, iRetE, nRetE)
2222 if (nRotE > 0): PlotEtax("RotE", aRotE, RotE0, dRotE, iRotE, nRotE)
2223 if (nDiE > 0): PlotEtax("DiE", aDiE, DiE0, dDiE, iDiE, nDiE)
2224 if (nRetO > 0): PlotEtax("RetO", aRetO, RetO0, dRetO, iRetO, nRetO)
2225 if (nRotO > 0): PlotEtax("RotO", aRotO, RotO0, dRotO, iRotO, nRotO)
2226 if (nDiO > 0): PlotEtax("DiO", aDiO, DiO0, dDiO, iDiO, nDiO)
2227 if (nDiC > 0): PlotEtax("DiC", aDiC, DiC0, dDiC, iDiC, nDiC)
2228 if (nRotC > 0): PlotEtax("RotC", aRotC, RotC0, dRotC, iRotC, nRotC)
2229 if (nRetC > 0): PlotEtax("RetC", aRetC, RetC0, dRetC, iRetC, nRetC)
2230 if (nTP > 0): PlotEtax("TP", aTP, TP0, dTP, iTP, nTP)
2231 if (nTS > 0): PlotEtax("TS", aTS, TS0, dTS, iTS, nTS)
2232 if (nRP > 0): PlotEtax("RP", aRP, RP0, dRP, iRP, nRP)
2233 if (nRS > 0): PlotEtax("RS", aRS, RS0, dRS, iRS, nRS)
2234 if (nRetT > 0): PlotEtax("RetT", aRetT, RetT0, dRetT, iRetT, nRetT)
2235 if (nRetR > 0): PlotEtax("RetR", aRetR, RetR0, dRetR, iRetR, nRetR)
2236 if (nERaT > 0): PlotEtax("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT)
2237 if (nERaR > 0): PlotEtax("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR)
2238 if (nRotaT > 0): PlotEtax("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT)
2239 if (nRotaR > 0): PlotEtax("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR)
2240 if (nLDRCal > 0): PlotEtax("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal)
2241 if (nTCalT > 0): PlotEtax("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT)
2242 if (nTCalR > 0): PlotEtax("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR)
2243 if (nNCal > 0): PlotEtax("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal)
2244 if (nNCal > 0): PlotEtax("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal)
2245 if (nNCal > 0): PlotEtax("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal)
2246 if (nNCal > 0): PlotEtax("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal)
2247 if (nNI > 0): PlotEtax("SigNoiseIt", aNIt, 0, 1, iNIt, nNI)
2248 if (nNI > 0): PlotEtax("SigNoiseIr", aNIr, 0, 1, iNIr, nNI)
2249 plt.show()
2250 plt.close
2251
2252 #Etaxmin = np.amin(aEtax[1, :])
2253 Etaxmin = np.amin(aEtax[1, :])
2254 Etaxmax = np.amax(aEtax[1, :])
2255 Etaxstd = np.std(aEtax[1, :])
2256 Etaxmean = np.mean(aEtax[1, :])
2257 Etaxmedian = np.mean(aEtax[1, :])
2258
2259 print("Etax: mean±std, median, max-mean, mean-min")
2260 print("{0:7.4f}±{1:7.4f},{2:7.4f},+{3:7.4f},-{4:7.4f}".format(Etaxmean, Etaxstd, Etaxmedian, Etaxmax-Etaxmean, Etaxmean-Etaxmin, ))
2261
1954 ''' 2262 '''
1955 # --- Plot K over LDRCal 2263 # --- Plot K over LDRCal
1956 fig3 = plt.figure() 2264 fig3 = plt.figure()
1957 plt.plot(LDRCal0+aLDRCal*dLDRCal/nLDRCal,aX[4,:], linewidth=2.0, color='b') 2265 plt.plot(LDRCal0+aLDRCal*dLDRCal/nLDRCal,aGHK[4,:], linewidth=2.0, color='b')
1958 2266
1959 plt.xlabel('LDRCal', fontsize=18) 2267 plt.xlabel('LDRCal', fontsize=18)
1960 plt.ylabel('K', fontsize=14) 2268 plt.ylabel('K', fontsize=14)
1961 plt.title(LID, fontsize=18) 2269 plt.title(LID, fontsize=18)
1962 plt.show() 2270 plt.show()

EARLINET Repository: Lidar polarisation GHK-parameter and error calculation

mercurial