1 # -*- coding: utf-8 -*- |
|
2 """ |
|
3 Copyright 2016, 2019 Volker Freudenthaler |
|
4 |
|
5 Licensed under the EUPL, Version 1.1 only (the "Licence"). |
|
6 |
|
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 |
|
9 a copy of the Licence at: |
|
10 |
|
11 https://joinup.ec.europa.eu/community/eupl/og_page/eupl |
|
12 |
|
13 Unless required by applicable law or agreed to in writing, software distributed |
|
14 under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR CONDITIONS |
|
15 OF ANY KIND, either express or implied. See the Licence for the specific language governing |
|
16 permissions and limitations under the Licence. |
|
17 |
|
18 Equation reference: http://www.atmos-meas-tech-discuss.net/amt-2015-338/amt-2015-338.pdf |
|
19 With equations code from Appendix C |
|
20 Python 3.7, seaborn 0.9.0 |
|
21 |
|
22 Code description: |
|
23 |
|
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. |
|
26 http://www.atmos-meas-tech-discuss.net/amt-2015-338/amt-2015-338.pdf provides an analytical model for the description of these errors, |
|
27 with which the measured signals can be corrected. |
|
28 This code simulates the lidar measurements with "assumed true" model parameters from an input file, and calculates the correction parameters (G,H, and K). |
|
29 The "assumed true" system parameters are the ones we think are the right ones, but in reality these parameters probably deviate from the assumed truth due to |
|
30 uncertainties. The uncertainties of the "assumed true" parameters can be described in the input file. Then this code calculates the lidar signals and the |
|
31 gain ratio eta* with all possible combinations of "errors", which represents the distribution of "possibly real" signals, and "corrects" them with the "assumed true" |
|
32 GHK parameters (GT0, GR0, HT0, HR0, and K0) to derive finally the distributions of "possibly real" linear depolarization ratios (LDRCorr), |
|
33 which are plotted for five different input linear depolarization ratios (LDRtrue). The red bars in the plots represent the input values of LDRtrue. |
|
34 A complication arises from the fact that the correction parameter K = eta*/eta (Eq. 83) can depend on the LDR during the calibration measurement, i.e. LDRcal or aCal |
|
35 in the code (see e.g. Eqs. (103), (115), and (141); mind the mistake in Eq. (116)). Therefor values of K for LDRcal = 0.004, 0.2, and 0.45 are calculated for |
|
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 |
|
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 |
|
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.7: includes the random error (signal noise) of the calibration and standard measurements |
|
43 Changes: |
|
44 Line 1687 Eta = (TaR * TiR) / (TaT * TiT) |
|
45 Line 1691 K = Etax / Eta # K of the real system; but correction in Line 1721 with K0 / Etax |
|
46 should work with nTCalT = nTCalR = 0 |
|
47 Ver. 0.9.7b: |
|
48 ToDo: include error due to TCalT und TCalR => determination of NCalT and NCalR etc. in error calculation line 1741ff |
|
49 combined error loops iNI and INCal for signals |
|
50 Ver. 0.9.7c: individual error loops for each of the six signals |
|
51 Ver. 0.9.7c2: different calculation of the signal noise errors |
|
52 Ver. 0.9.7c3: n.a.different calculation of the signal noise errors |
|
53 Ver. 0.9.7c4: test to speed up the loops for error calculation by moving them just before the actual calculation: still some code errors |
|
54 Ver. 0.9.8: |
|
55 - 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 |
|
56 - calculation of the PLDR from LDR and BSR, BSR, and LDRm |
|
57 - 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. |
|
58 Ver. 0.9.8b: change from "TTa = TiT * TaT" to "TTa = TiT * TaT * ATPT" etc. (compare ver 0.9.8 with 0.9.8b) removes |
|
59 - the strong Tp dependence of the errors |
|
60 - the factor 2 in the GH parameters |
|
61 - see c:\technik\Optik\Polarizers\DepCal\ApplOpt\GH-parameters-190114.odt |
|
62 Ver. 0.9.8c: includes error of Etax |
|
63 Ver. 0.9.8d: Eta0, K0 etc in error loop replaced by Eta0y, K0y etc. Changes in signal noise calculations |
|
64 Ver. 0.9.8e: ambiguous laser spec. DOLP (no discrimination between left and right circular polarisation) replaced by Stokes parameters Qin, Uin |
|
65 Ver. 0.9.8e2: Added plot of LDRsim, Etax, Etapx, Etamx; LDRCorr and aLDRcorr consistently named |
|
66 Ver. 0.9.8e3: Change of OutputFile name; Change of Ir and It noise if (CalcFrom0deg) = False; (Different calculation of error contributions tested but not implemented) |
|
67 Ver. 0.9.8e4: text changed for y=+-1 (see line 274 ff and line 1044 ff |
|
68 |
|
69 ======================================================== |
|
70 simulation: LDRsim = Ir / It with variable parameters (possible truths) |
|
71 G,H,Eta,Etax,K |
|
72 It = TaT * TiT * ATP1 * TiO * TiE * (GT + atrue * HT) |
|
73 LDRsim = Ir / It |
|
74 consistency test: is forward simulation and correction consistent? |
|
75 LDRCorr = (LDRsim / Eta * (GT + HT) - (GR + HR)) / ((GR - HR) - LDRsim / Eta * (GT - HT)) => atrue? |
|
76 assumed true: G0,H0,Eta0,Etax0,K0 => actual retrievals of LDRCorr |
|
77 => correct possible truths with assumed true G0,H0,Eta0 |
|
78 measure: It, Ir, EtaX |
|
79 coorect it with: G0,H0,K0 |
|
80 LDRCorr = (LDRsim / (Etax / K0) * (GT0 + HT0) - (GR0 + HR0)) / ((GR0 - HR0) - LDRsim0 / (Etax / K0) * (GT0 - HT0)) |
|
81 """ |
|
82 # 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. |
|
83 from __future__ import print_function |
|
84 # !/usr/bin/env python3 |
|
85 |
|
86 import os |
|
87 import sys |
|
88 |
|
89 from scipy.stats import kurtosis |
|
90 from scipy.stats import skew |
|
91 # use: kurtosis(data, fisher=True,bias=False) => 0; skew(data,bias=False) => 0 |
|
92 # Comment: the seaborn library makes nicer plots, but the code works also without it. |
|
93 import numpy as np |
|
94 import matplotlib.pyplot as plt |
|
95 |
|
96 try: |
|
97 import seaborn as sns |
|
98 |
|
99 sns_loaded = True |
|
100 except ImportError: |
|
101 sns_loaded = False |
|
102 |
|
103 # from time import clock # python 2 |
|
104 from timeit import default_timer as clock |
|
105 |
|
106 # from matplotlib.backends.backend_pdf import PdfPages |
|
107 # pdffile = '{}.pdf'.format('path') |
|
108 # pp = PdfPages(pdffile) |
|
109 ## pp.savefig can be called multiple times to save to multiple pages |
|
110 # pp.savefig() |
|
111 # pp.close() |
|
112 |
|
113 from contextlib import contextmanager |
|
114 |
|
115 @contextmanager |
|
116 def redirect_stdout(new_target): |
|
117 old_target, sys.stdout = sys.stdout, new_target # replace sys.stdout |
|
118 try: |
|
119 yield new_target # run some code with the replaced stdout |
|
120 finally: |
|
121 sys.stdout.flush() |
|
122 sys.stdout = old_target # restore to the previous value |
|
123 |
|
124 ''' |
|
125 real_raw_input = vars(__builtins__).get('raw_input',input) |
|
126 ''' |
|
127 try: |
|
128 import __builtin__ |
|
129 |
|
130 input = getattr(__builtin__, 'raw_input') |
|
131 except (ImportError, AttributeError): |
|
132 pass |
|
133 |
|
134 from distutils.util import strtobool |
|
135 |
|
136 |
|
137 def user_yes_no_query(question): |
|
138 sys.stdout.write('%s [y/n]\n' % question) |
|
139 while True: |
|
140 try: |
|
141 return strtobool(input().lower()) |
|
142 except ValueError: |
|
143 sys.stdout.write('Please respond with \'y\' or \'n\'.\n') |
|
144 |
|
145 |
|
146 # if user_yes_no_query('want to exit?') == 1: sys.exit() |
|
147 |
|
148 abspath = os.path.abspath(__file__) |
|
149 dname = os.path.dirname(abspath) |
|
150 fname = os.path.basename(abspath) |
|
151 os.chdir(dname) |
|
152 |
|
153 # PrintToOutputFile = True |
|
154 |
|
155 sqr05 = 0.5 ** 0.5 |
|
156 |
|
157 # ---- Initial definition of variables; the actual values will be read in with exec(open('./optic_input.py').read()) below |
|
158 # Do you want to calculate the errors? If not, just the GHK-parameters are determined. |
|
159 Error_Calc = True |
|
160 LID = "internal" |
|
161 EID = "internal" |
|
162 # --- IL Laser IL and +-Uncertainty |
|
163 Qin, dQin, nQin = 1., 0.0, 0 # second Stokes vector parameter; default 1 => linear polarization |
|
164 Vin, dVin, nVin = 0., 0.0, 0 # fourth Stokes vector parameter |
|
165 RotL, dRotL, nRotL = 0.0, 0.0, 1 # alpha; rotation of laser polarization in degrees; default 0 |
|
166 # IL = 1e5 #photons in the laser beam, including detection efficiency of the telescope, atmodspheric and r^2 attenuation |
|
167 # --- ME Emitter and +-Uncertainty |
|
168 DiE, dDiE, nDiE = 0., 0.00, 1 # Diattenuation |
|
169 TiE = 1. # Unpolarized transmittance |
|
170 RetE, dRetE, nRetE = 0., 180.0, 0 # Retardance in degrees |
|
171 RotE, dRotE, nRotE = 0., 0.0, 0 # beta: Rotation of optical element in degrees |
|
172 # --- MO Receiver Optics including telescope |
|
173 DiO, dDiO, nDiO = -0.055, 0.003, 1 |
|
174 TiO = 0.9 |
|
175 RetO, dRetO, nRetO = 0., 180.0, 2 |
|
176 RotO, dRotO, nRotO = 0., 0.1, 1 # gamma |
|
177 # --- PBS MT transmitting path defined with (TS,TP); and +-Uncertainty |
|
178 TP, dTP, nTP = 0.98, 0.02, 1 |
|
179 TS, dTS, nTS = 0.001, 0.001, 1 |
|
180 TiT = 0.5 * (TP + TS) |
|
181 DiT = (TP - TS) / (TP + TS) |
|
182 # PolFilter |
|
183 RetT, dRetT, nRetT = 0., 180., 0 |
|
184 ERaT, dERaT, nERaT = 0.001, 0.001, 1 |
|
185 RotaT, dRotaT, nRotaT = 0., 3., 1 |
|
186 DaT = (1 - ERaT) / (1 + ERaT) |
|
187 TaT = 0.5 * (1 + ERaT) |
|
188 # --- PBS MR reflecting path defined with (RS,RP); and +-Uncertainty |
|
189 RS_RP_depend_on_TS_TP = False |
|
190 if (RS_RP_depend_on_TS_TP): |
|
191 RP, dRP, nRP = 1 - TP, 0.0, 0 |
|
192 RS, dRS, nRS = 1 - TS, 0.0, 0 |
|
193 else: |
|
194 RP, dRP, nRP = 0.05, 0.01, 1 |
|
195 RS, dRS, nRS = 0.98, 0.01, 1 |
|
196 TiR = 0.5 * (RP + RS) |
|
197 DiR = (RP - RS) / (RP + RS) |
|
198 # PolFilter |
|
199 RetR, dRetR, nRetR = 0., 180., 0 |
|
200 ERaR, dERaR, nERaR = 0.001, 0.001, 1 |
|
201 RotaR, dRotaR, nRotaR = 90., 3., 1 |
|
202 DaR = (1 - ERaR) / (1 + ERaR) |
|
203 TaR = 0.5 * (1 + ERaR) |
|
204 |
|
205 # +++ Orientation of the PBS with respect to the reference plane (see Polarisation-orientation.png and Polarisation-orientation-2.png in /system_settings) |
|
206 # Y = +1: PBS incidence plane is parallel to reference plane and polarisation in reference plane is finally transmitted. |
|
207 # Y = -1: PBS incidence plane is perpendicular to reference plane and polarisation in reference plane is finally reflected. |
|
208 Y = 1. |
|
209 |
|
210 # Calibrator = type defined by matrix values |
|
211 LocC = 4 # location of calibrator: behind laser = 1; behind emitter = 2; before receiver = 3; before PBS = 4 |
|
212 |
|
213 # --- Additional attenuation (transmission of the ND-filter) during the calibration |
|
214 TCalT, dTCalT, nTCalT = 1, 0., 0 # transmitting path; error calc not working yet |
|
215 TCalR, dTCalR, nTCalR = 1, 0., 0 # reflecting path; error calc not working yet |
|
216 |
|
217 # *** signal noise error calculation |
|
218 # --- number of photon counts in the signal summed up in the calibration range during the calibration measurements |
|
219 NCalT = 1e6 # default 1e6, assumed the same in +45° and -45° signals |
|
220 NCalR = 1e6 # default 1e6, assumed the same in +45° and -45° signals |
|
221 NILfac = 1.0 # duration of standard (0°) measurement relative to calibration measurements |
|
222 nNCal = 0 # error nNCal: one-sigma in steps to left and right for calibration signals |
|
223 nNI = 0 # error nNI: one-sigma in steps to left and right for 0° signals |
|
224 NI = 50000 #number of photon counts in the parallel 0°-signal |
|
225 eFacT = 1.0 # rel. amplification of transmitted channel, approximate values are sufficient; def. = 1 |
|
226 eFacR = 10.0 |
|
227 IoutTp0, IoutTp, dIoutTp0 = 0.5, 0.5, 0.0 |
|
228 IoutTm0, IoutTm, dIoutTm0 = 0.5, 0.5, 0.0 |
|
229 IoutRp0, IoutRp, dIoutRp0 = 0.5, 0.5, 0.0 |
|
230 IoutRm0, IoutRm, dIoutRm0 = 0.5, 0.5, 0.0 |
|
231 It0, It, dIt0 = 1 , 1, 0 |
|
232 Ir0, Ir, dTr0 = 1 , 1, 0 |
|
233 CalcFrom0deg = True |
|
234 |
|
235 TypeC = 3 # linear polarizer calibrator |
|
236 # example with extinction ratio 0.001 |
|
237 DiC, dDiC, nDiC = 1.0, 0., 0 # ideal 1.0 |
|
238 TiC = 0.5 # ideal 0.5 |
|
239 RetC, dRetC, nRetC = 0.0, 0.0, 0 |
|
240 RotC, dRotC, nRotC = 0.0, 0.1, 0 # constant calibrator offset epsilon |
|
241 RotationErrorEpsilonForNormalMeasurements = False # is in general False for TypeC == 3 calibrator |
|
242 |
|
243 # Rotation error without calibrator: if False, then epsilon = 0 for normal measurements |
|
244 RotationErrorEpsilonForNormalMeasurements = True |
|
245 # BSR backscatter ratio |
|
246 # BSR, dBSR, nBSR = 10, 0.05, 1 |
|
247 BSR = np.zeros(5) |
|
248 BSR = [1.1, 2, 5, 10., 50.] |
|
249 # theoretical molecular LDR LDRm |
|
250 LDRm, dLDRm, nLDRm = 0.004, 0.001, 1 |
|
251 # LDRCal assumed atmospheric linear depolarization ratio during the calibration measurements (first guess) |
|
252 LDRCal0, dLDRCal, nLDRCal = 0.25, 0.04, 1 |
|
253 LDRCal = LDRCal0 |
|
254 # measured LDRm will be corrected with calculated parameters |
|
255 LDRmeas = 0.015 |
|
256 # LDRtrue for simulation of measurement => LDRsim |
|
257 LDRtrue = 0.004 |
|
258 LDRtrue2 = 0.004 |
|
259 LDRunCorr = 1. |
|
260 # Initialize other values to 0 |
|
261 ER, nER, dER = 0.001, 0, 0.001 |
|
262 K = 0. |
|
263 Km = 0. |
|
264 Kp = 0. |
|
265 LDRCorr = 0. |
|
266 Eta = 0. |
|
267 Ir = 0. |
|
268 It = 0. |
|
269 h = 1. |
|
270 |
|
271 Loc = ['', 'behind laser', 'behind emitter', 'before receiver', 'before PBS'] |
|
272 Type = ['', 'mechanical rotator', 'hwp rotator', 'linear polarizer', 'qwp rotator', 'circular polarizer', |
|
273 'real HWP +-22.5°'] |
|
274 |
|
275 bPlotEtax = False |
|
276 |
|
277 # end of initial definition of variables |
|
278 # ******************************************************************************************************************************* |
|
279 # --- Read actual lidar system parameters from optic_input.py (must be in the programs sub-directory 'system_settings') |
|
280 # ******************************************************************************************************************************* |
|
281 |
|
282 # InputFile = 'optic_input_example_2_1.py' |
|
283 # InputFile = 'ALidar-355-F-3-3c2-0.9.8d.py' |
|
284 # InputFile = 'Polarimeter-4C3-ver0.98e.py' |
|
285 # InputFile = 'Polarimeter-4A-ver0.98e.py' |
|
286 InputFile = 'optic_input_raym-200-02-18-ver0.9.8e.py' |
|
287 InputFile = 'optic_input_raym-200-04-17-ver0.9.8e.py' |
|
288 InputFile = 'optic_input_raym-200-04-17-ver0.9.8e-extended.py' |
|
289 InputFile = 'Adam_ver0.98.py' |
|
290 InputFile = 'MUSA-B3A-ver0.98e.py' |
|
291 InputFile = 'MUSA-B4A-ver0.98e.py' |
|
292 # InputFile = 'MUSA-A3C-ver0.98e.py' |
|
293 InputFile = 'optic_input_ver0.98e_LILI_532_May2020.py' |
|
294 InputFile = 'optic_input_ver0.98e_LILI_532_May2020_RotL=90.py' |
|
295 InputFile = 'optic_input_0.9.8e4-PollyXT_Lacros.py' |
|
296 InputFile = 'optic_input_UPC-lidar_0.9.8e4.py' |
|
297 InputFile = 'optic_input_UV-Pot-ver0.9.8e.py' |
|
298 InputFile = 'optic_input_0.9.8e4-PollyXT_Lacros.py' |
|
299 InputFile = 'optic_input_example_lidar_ver0.9.8e.py' |
|
300 |
|
301 # ******************************************************************************************************************************* |
|
302 |
|
303 ''' |
|
304 print("From ", dname) |
|
305 print("Running ", fname) |
|
306 print("Reading input file ", InputFile, " for") |
|
307 ''' |
|
308 input_path = os.path.join('.', 'system_settings', InputFile) |
|
309 # this works with Python 2 and 3! |
|
310 exec(open(input_path).read(), globals()) |
|
311 # end of read actual system parameters |
|
312 |
|
313 |
|
314 # --- Manual Parameter Change --- |
|
315 # (use for quick parameter changes without changing the input file ) |
|
316 # DiO = 0. |
|
317 # LDRtrue = 0.45 |
|
318 # LDRtrue2 = 0.004 |
|
319 # Y = -1 |
|
320 # LocC = 4 #location of calibrator: 1 = behind laser; 2 = behind emitter; 3 = before receiver; 4 = before PBS |
|
321 # #TypeC = 6 Don't change the TypeC here |
|
322 # RotationErrorEpsilonForNormalMeasurements = True |
|
323 # LDRCal = 0.25 |
|
324 # # --- Errors |
|
325 Qin0, dQin, nQin = Qin, dQin, nQin |
|
326 Vin0, dVin, nVin = Vin, dVin, nVin |
|
327 RotL0, dRotL, nRotL = RotL, dRotL, nRotL |
|
328 |
|
329 DiE0, dDiE, nDiE = DiE, dDiE, nDiE |
|
330 RetE0, dRetE, nRetE = RetE, dRetE, nRetE |
|
331 RotE0, dRotE, nRotE = RotE, dRotE, nRotE |
|
332 |
|
333 DiO0, dDiO, nDiO = DiO, dDiO, nDiO |
|
334 RetO0, dRetO, nRetO = RetO, dRetO, nRetO |
|
335 RotO0, dRotO, nRotO = RotO, dRotO, nRotO |
|
336 |
|
337 DiC0, dDiC, nDiC = DiC, dDiC, nDiC |
|
338 RetC0, dRetC, nRetC = RetC, dRetC, nRetC |
|
339 RotC0, dRotC, nRotC = RotC, dRotC, nRotC |
|
340 |
|
341 TP0, dTP, nTP = TP, dTP, nTP |
|
342 TS0, dTS, nTS = TS, dTS, nTS |
|
343 RetT0, dRetT, nRetT = RetT, dRetT, nRetT |
|
344 |
|
345 ERaT0, dERaT, nERaT = ERaT, dERaT, nERaT |
|
346 RotaT0, dRotaT, nRotaT = RotaT, dRotaT, nRotaT |
|
347 |
|
348 RP0, dRP, nRP = RP, dRP, nRP |
|
349 RS0, dRS, nRS = RS, dRS, nRS |
|
350 RetR0, dRetR, nRetR = RetR, dRetR, nRetR |
|
351 |
|
352 ERaR0, dERaR, nERaR = ERaR, dERaR, nERaR |
|
353 RotaR0, dRotaR, nRotaR = RotaR, dRotaR, nRotaR |
|
354 |
|
355 LDRCal0, dLDRCal, nLDRCal = LDRCal, dLDRCal, nLDRCal |
|
356 |
|
357 # BSR0, dBSR, nBSR = BSR, dBSR, nBSR |
|
358 LDRm0, dLDRm, nLDRm = LDRm, dLDRm, nLDRm |
|
359 # ---------- End of manual parameter change |
|
360 |
|
361 RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC = RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0 |
|
362 TP, TS, RP, RS, ERaT, RotaT, RetT, ERaR, RotaR, RetR = TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0 |
|
363 LDRCal = LDRCal0 |
|
364 DTa0, TTa0, DRa0, TRa0, LDRsimx, LDRCorr = 0., 0., 0., 0., 0., 0. |
|
365 TCalT0, TCalR0 = TCalT, TCalR |
|
366 |
|
367 TiT = 0.5 * (TP + TS) |
|
368 DiT = (TP - TS) / (TP + TS) |
|
369 ZiT = (1. - DiT ** 2) ** 0.5 |
|
370 TiR = 0.5 * (RP + RS) |
|
371 DiR = (RP - RS) / (RP + RS) |
|
372 ZiR = (1. - DiR ** 2) ** 0.5 |
|
373 |
|
374 C2aT = np.cos(np.deg2rad(2. * RotaT)) |
|
375 C2aR = np.cos(np.deg2rad(2. * RotaR)) |
|
376 ATPT = float(1. + C2aT * DaT * DiT) |
|
377 ARPT = float(1. + C2aR * DaR * DiR) |
|
378 TTa = TiT * TaT * ATPT # unpolarized transmission |
|
379 TRa = TiR * TaR * ARPT # unpolarized transmission |
|
380 Eta0 = TRa / TTa |
|
381 |
|
382 # --- alternative texts for output |
|
383 dY = ['perpendicular', '', 'parallel'] |
|
384 dY2 = ['reflected', '', 'transmitted'] |
|
385 if ((abs(RotL) < 45 and Y == 1) or (abs(RotL) >= 45 and Y == -1)): |
|
386 dY3 = "Parallel laser polarisation is detected in transmitted channel" |
|
387 else: |
|
388 dY3 = "Parallel laser polarisation is detected in reflected channel" |
|
389 |
|
390 # --- check input errors |
|
391 if ((Qin ** 2 + Vin ** 2) ** 0.5) > 1: |
|
392 print("Error: degree of polarisation of laser > 1. Check Qin and Vin! ") |
|
393 sys.exit() |
|
394 |
|
395 # --- this subroutine is for the calculation of the PLDR from LDR, BSR, and LDRm ------------------- |
|
396 def CalcPLDR(LDR, BSR, LDRm): |
|
397 PLDR = (BSR * (1. + LDRm) * LDR - LDRm * (1. + LDR)) / (BSR * (1. + LDRm) - (1. + LDR)) |
|
398 return (PLDR) |
|
399 # --- this subroutine is for the calculation with certain fixed parameters ------------------------ |
|
400 def Calc(TCalT, TCalR, NCalT, NCalR, Qin, Vin, RotL, RotE, RetE, DiE, RotO, RetO, DiO, |
|
401 RotC, RetC, DiC, TP, TS, RP, RS, |
|
402 ERaT, RotaT, RetT, ERaR, RotaR, RetR, LDRCal): |
|
403 # ---- Do the calculations of bra-ket vectors |
|
404 h = -1. if TypeC == 2 else 1 |
|
405 # from input file: assumed LDRCal for calibration measurements |
|
406 aCal = (1. - LDRCal) / (1. + LDRCal) |
|
407 atrue = (1. - LDRtrue) / (1. + LDRtrue) |
|
408 |
|
409 # angles of emitter and laser and calibrator and receiver optics |
|
410 # RotL = alpha, RotE = beta, RotO = gamma, RotC = epsilon |
|
411 S2a = np.sin(2 * np.deg2rad(RotL)) |
|
412 C2a = np.cos(2 * np.deg2rad(RotL)) |
|
413 S2b = np.sin(2 * np.deg2rad(RotE)) |
|
414 C2b = np.cos(2 * np.deg2rad(RotE)) |
|
415 S2ab = np.sin(np.deg2rad(2 * RotL - 2 * RotE)) |
|
416 C2ab = np.cos(np.deg2rad(2 * RotL - 2 * RotE)) |
|
417 S2g = np.sin(np.deg2rad(2 * RotO)) |
|
418 C2g = np.cos(np.deg2rad(2 * RotO)) |
|
419 |
|
420 # Laser with Degree of linear polarization DOLP |
|
421 IinL = 1. |
|
422 QinL = Qin |
|
423 UinL = 0. |
|
424 VinL = Vin |
|
425 # VinL = (1. - DOLP ** 2) ** 0.5 |
|
426 |
|
427 # Stokes Input Vector rotation Eq. E.4 |
|
428 A = C2a * QinL - S2a * UinL |
|
429 B = S2a * QinL + C2a * UinL |
|
430 # Stokes Input Vector rotation Eq. E.9 |
|
431 C = C2ab * QinL - S2ab * UinL |
|
432 D = S2ab * QinL + C2ab * UinL |
|
433 |
|
434 # emitter optics |
|
435 CosE = np.cos(np.deg2rad(RetE)) |
|
436 SinE = np.sin(np.deg2rad(RetE)) |
|
437 ZiE = (1. - DiE ** 2) ** 0.5 |
|
438 WiE = (1. - ZiE * CosE) |
|
439 |
|
440 # Stokes Input Vector after emitter optics equivalent to Eq. E.9 with already rotated input vector from Eq. E.4 |
|
441 # b = beta |
|
442 IinE = (IinL + DiE * C) |
|
443 QinE = (C2b * DiE * IinL + A + S2b * (WiE * D - ZiE * SinE * VinL)) |
|
444 UinE = (S2b * DiE * IinL + B - C2b * (WiE * D - ZiE * SinE * VinL)) |
|
445 VinE = (-ZiE * SinE * D + ZiE * CosE * VinL) |
|
446 |
|
447 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) |
|
448 IinF = IinE |
|
449 QinF = aCal * QinE |
|
450 UinF = -aCal * UinE |
|
451 VinF = (1. - 2. * aCal) * VinE |
|
452 |
|
453 # receiver optics |
|
454 CosO = np.cos(np.deg2rad(RetO)) |
|
455 SinO = np.sin(np.deg2rad(RetO)) |
|
456 ZiO = (1. - DiO ** 2) ** 0.5 |
|
457 WiO = (1. - ZiO * CosO) |
|
458 |
|
459 # calibrator |
|
460 CosC = np.cos(np.deg2rad(RetC)) |
|
461 SinC = np.sin(np.deg2rad(RetC)) |
|
462 ZiC = (1. - DiC ** 2) ** 0.5 |
|
463 WiC = (1. - ZiC * CosC) |
|
464 |
|
465 # Stokes Input Vector before the polarising beam splitter Eq. E.31 |
|
466 A = C2g * QinE - S2g * UinE |
|
467 B = S2g * QinE + C2g * UinE |
|
468 |
|
469 IinP = (IinE + DiO * aCal * A) |
|
470 QinP = (C2g * DiO * IinE + aCal * QinE - S2g * (WiO * aCal * B + ZiO * SinO * (1. - 2. * aCal) * VinE)) |
|
471 UinP = (S2g * DiO * IinE - aCal * UinE + C2g * (WiO * aCal * B + ZiO * SinO * (1. - 2. * aCal) * VinE)) |
|
472 VinP = (ZiO * SinO * aCal * B + ZiO * CosO * (1. - 2. * aCal) * VinE) |
|
473 |
|
474 # ------------------------- |
|
475 # F11 assuemd to be = 1 => measured: F11m = IinP / IinE with atrue |
|
476 # F11sim = TiO*(IinE + DiO*atrue*A)/IinE |
|
477 # ------------------------- |
|
478 |
|
479 # analyser |
|
480 if (RS_RP_depend_on_TS_TP): |
|
481 RS = 1. - TS |
|
482 RP = 1. - TP |
|
483 |
|
484 TiT = 0.5 * (TP + TS) |
|
485 DiT = (TP - TS) / (TP + TS) |
|
486 ZiT = (1. - DiT ** 2) ** 0.5 |
|
487 TiR = 0.5 * (RP + RS) |
|
488 DiR = (RP - RS) / (RP + RS) |
|
489 ZiR = (1. - DiR ** 2) ** 0.5 |
|
490 CosT = np.cos(np.deg2rad(RetT)) |
|
491 SinT = np.sin(np.deg2rad(RetT)) |
|
492 CosR = np.cos(np.deg2rad(RetR)) |
|
493 SinR = np.sin(np.deg2rad(RetR)) |
|
494 |
|
495 DaT = (1. - ERaT) / (1. + ERaT) |
|
496 DaR = (1. - ERaR) / (1. + ERaR) |
|
497 TaT = 0.5 * (1. + ERaT) |
|
498 TaR = 0.5 * (1. + ERaR) |
|
499 |
|
500 S2aT = np.sin(np.deg2rad(h * 2 * RotaT)) |
|
501 C2aT = np.cos(np.deg2rad(2 * RotaT)) |
|
502 S2aR = np.sin(np.deg2rad(h * 2 * RotaR)) |
|
503 C2aR = np.cos(np.deg2rad(2 * RotaR)) |
|
504 |
|
505 # Analyzer As before the PBS Eq. D.5; combined PBS and cleaning pol-filter |
|
506 ATPT = (1. + C2aT * DaT * DiT) # unpolarized transmission correction |
|
507 TTa = TiT * TaT * ATPT # unpolarized transmission |
|
508 ATP1 = 1. |
|
509 ATP2 = Y * (DiT + C2aT * DaT) / ATPT |
|
510 ATP3 = Y * S2aT * DaT * ZiT * CosT / ATPT |
|
511 ATP4 = S2aT * DaT * ZiT * SinT / ATPT |
|
512 ATP = np.array([ATP1, ATP2, ATP3, ATP4]) |
|
513 DTa = ATP2 * Y |
|
514 |
|
515 ARPT = (1 + C2aR * DaR * DiR) # unpolarized transmission correction |
|
516 TRa = TiR * TaR * ARPT # unpolarized transmission |
|
517 ARP1 = 1 |
|
518 ARP2 = Y * (DiR + C2aR * DaR) / ARPT |
|
519 ARP3 = Y * S2aR * DaR * ZiR * CosR / ARPT |
|
520 ARP4 = S2aR * DaR * ZiR * SinR / ARPT |
|
521 ARP = np.array([ARP1, ARP2, ARP3, ARP4]) |
|
522 DRa = ARP2 * Y |
|
523 |
|
524 |
|
525 # ---- Calculate signals and correction parameters for diffeent locations and calibrators |
|
526 if LocC == 4: # Calibrator before the PBS |
|
527 # print("Calibrator location not implemented yet") |
|
528 |
|
529 # S2ge = np.sin(np.deg2rad(2*RotO + h*2*RotC)) |
|
530 # C2ge = np.cos(np.deg2rad(2*RotO + h*2*RotC)) |
|
531 S2e = np.sin(np.deg2rad(h * 2 * RotC)) |
|
532 C2e = np.cos(np.deg2rad(2 * RotC)) |
|
533 # rotated AinP by epsilon Eq. C.3 |
|
534 ATP2e = C2e * ATP2 + S2e * ATP3 |
|
535 ATP3e = C2e * ATP3 - S2e * ATP2 |
|
536 ARP2e = C2e * ARP2 + S2e * ARP3 |
|
537 ARP3e = C2e * ARP3 - S2e * ARP2 |
|
538 ATPe = np.array([ATP1, ATP2e, ATP3e, ATP4]) |
|
539 ARPe = np.array([ARP1, ARP2e, ARP3e, ARP4]) |
|
540 # Stokes Input Vector before the polarising beam splitter Eq. E.31 |
|
541 A = C2g * QinE - S2g * UinE |
|
542 B = S2g * QinE + C2g * UinE |
|
543 # C = (WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE) |
|
544 Co = ZiO * SinO * VinE |
|
545 Ca = (WiO * B - 2 * ZiO * SinO * VinE) |
|
546 # C = Co + aCal*Ca |
|
547 # IinP = (IinE + DiO*aCal*A) |
|
548 # QinP = (C2g*DiO*IinE + aCal*QinE - S2g*C) |
|
549 # UinP = (S2g*DiO*IinE - aCal*UinE + C2g*C) |
|
550 # VinP = (ZiO*SinO*aCal*B + ZiO*CosO*(1-2*aCal)*VinE) |
|
551 IinPo = IinE |
|
552 QinPo = (C2g * DiO * IinE - S2g * Co) |
|
553 UinPo = (S2g * DiO * IinE + C2g * Co) |
|
554 VinPo = ZiO * CosO * VinE |
|
555 |
|
556 IinPa = DiO * A |
|
557 QinPa = QinE - S2g * Ca |
|
558 UinPa = -UinE + C2g * Ca |
|
559 VinPa = ZiO * (SinO * B - 2 * CosO * VinE) |
|
560 |
|
561 IinP = IinPo + aCal * IinPa |
|
562 QinP = QinPo + aCal * QinPa |
|
563 UinP = UinPo + aCal * UinPa |
|
564 VinP = VinPo + aCal * VinPa |
|
565 # Stokes Input Vector before the polarising beam splitter rotated by epsilon Eq. C.3 |
|
566 # QinPe = C2e*QinP + S2e*UinP |
|
567 # UinPe = C2e*UinP - S2e*QinP |
|
568 QinPoe = C2e * QinPo + S2e * UinPo |
|
569 UinPoe = C2e * UinPo - S2e * QinPo |
|
570 QinPae = C2e * QinPa + S2e * UinPa |
|
571 UinPae = C2e * UinPa - S2e * QinPa |
|
572 QinPe = C2e * QinP + S2e * UinP |
|
573 UinPe = C2e * UinP - S2e * QinP |
|
574 |
|
575 # Calibration signals and Calibration correction K from measurements with LDRCal / aCal |
|
576 if (TypeC == 2) or (TypeC == 1): # rotator calibration Eq. C.4 |
|
577 # parameters for calibration with aCal |
|
578 AT = ATP1 * IinP + h * ATP4 * VinP |
|
579 BT = ATP3e * QinP - h * ATP2e * UinP |
|
580 AR = ARP1 * IinP + h * ARP4 * VinP |
|
581 BR = ARP3e * QinP - h * ARP2e * UinP |
|
582 # Correction parameters for normal measurements; they are independent of LDR |
|
583 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out |
|
584 IS1 = np.array([IinPo, QinPo, UinPo, VinPo]) |
|
585 IS2 = np.array([IinPa, QinPa, UinPa, VinPa]) |
|
586 GT = np.dot(ATP, IS1) |
|
587 GR = np.dot(ARP, IS1) |
|
588 HT = np.dot(ATP, IS2) |
|
589 HR = np.dot(ARP, IS2) |
|
590 else: |
|
591 IS1 = np.array([IinPo, QinPo, UinPo, VinPo]) |
|
592 IS2 = np.array([IinPa, QinPa, UinPa, VinPa]) |
|
593 GT = np.dot(ATPe, IS1) |
|
594 GR = np.dot(ARPe, IS1) |
|
595 HT = np.dot(ATPe, IS2) |
|
596 HR = np.dot(ARPe, IS2) |
|
597 elif (TypeC == 3) or (TypeC == 4): # linear polariser calibration Eq. C.5 |
|
598 # parameters for calibration with aCal |
|
599 AT = ATP1 * IinP + ATP3e * UinPe + ZiC * CosC * (ATP2e * QinPe + ATP4 * VinP) |
|
600 BT = DiC * (ATP1 * UinPe + ATP3e * IinP) - ZiC * SinC * (ATP2e * VinP - ATP4 * QinPe) |
|
601 AR = ARP1 * IinP + ARP3e * UinPe + ZiC * CosC * (ARP2e * QinPe + ARP4 * VinP) |
|
602 BR = DiC * (ARP1 * UinPe + ARP3e * IinP) - ZiC * SinC * (ARP2e * VinP - ARP4 * QinPe) |
|
603 # Correction parameters for normal measurements; they are independent of LDR |
|
604 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out |
|
605 IS1 = np.array([IinPo, QinPo, UinPo, VinPo]) |
|
606 IS2 = np.array([IinPa, QinPa, UinPa, VinPa]) |
|
607 GT = np.dot(ATP, IS1) |
|
608 GR = np.dot(ARP, IS1) |
|
609 HT = np.dot(ATP, IS2) |
|
610 HR = np.dot(ARP, IS2) |
|
611 else: |
|
612 IS1e = np.array([IinPo + DiC * QinPoe, DiC * IinPo + QinPoe, ZiC * (CosC * UinPoe + SinC * VinPo), |
|
613 -ZiC * (SinC * UinPoe - CosC * VinPo)]) |
|
614 IS2e = np.array([IinPa + DiC * QinPae, DiC * IinPa + QinPae, ZiC * (CosC * UinPae + SinC * VinPa), |
|
615 -ZiC * (SinC * UinPae - CosC * VinPa)]) |
|
616 GT = np.dot(ATPe, IS1e) |
|
617 GR = np.dot(ARPe, IS1e) |
|
618 HT = np.dot(ATPe, IS2e) |
|
619 HR = np.dot(ARPe, IS2e) |
|
620 elif (TypeC == 6): # diattenuator calibration +-22.5° rotated_diattenuator_X22x5deg.odt |
|
621 # parameters for calibration with aCal |
|
622 AT = ATP1 * IinP + sqr05 * DiC * (ATP1 * QinPe + ATP2e * IinP) + (1. - 0.5 * WiC) * ( |
|
623 ATP2e * QinPe + ATP3e * UinPe) + ZiC * (sqr05 * SinC * (ATP3e * VinP - ATP4 * UinPe) + ATP4 * CosC * VinP) |
|
624 BT = sqr05 * DiC * (ATP1 * UinPe + ATP3e * IinP) + 0.5 * WiC * ( |
|
625 ATP2e * UinPe + ATP3e * QinPe) - sqr05 * ZiC * SinC * (ATP2e * VinP - ATP4 * QinPe) |
|
626 AR = ARP1 * IinP + sqr05 * DiC * (ARP1 * QinPe + ARP2e * IinP) + (1. - 0.5 * WiC) * ( |
|
627 ARP2e * QinPe + ARP3e * UinPe) + ZiC * (sqr05 * SinC * (ARP3e * VinP - ARP4 * UinPe) + ARP4 * CosC * VinP) |
|
628 BR = sqr05 * DiC * (ARP1 * UinPe + ARP3e * IinP) + 0.5 * WiC * ( |
|
629 ARP2e * UinPe + ARP3e * QinPe) - sqr05 * ZiC * SinC * (ARP2e * VinP - ARP4 * QinPe) |
|
630 # Correction parameters for normal measurements; they are independent of LDR |
|
631 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out |
|
632 IS1 = np.array([IinPo, QinPo, UinPo, VinPo]) |
|
633 IS2 = np.array([IinPa, QinPa, UinPa, VinPa]) |
|
634 GT = np.dot(ATP, IS1) |
|
635 GR = np.dot(ARP, IS1) |
|
636 HT = np.dot(ATP, IS2) |
|
637 HR = np.dot(ARP, IS2) |
|
638 else: |
|
639 IS1e = np.array([IinPo + DiC * QinPoe, DiC * IinPo + QinPoe, ZiC * (CosC * UinPoe + SinC * VinPo), |
|
640 -ZiC * (SinC * UinPoe - CosC * VinPo)]) |
|
641 IS2e = np.array([IinPa + DiC * QinPae, DiC * IinPa + QinPae, ZiC * (CosC * UinPae + SinC * VinPa), |
|
642 -ZiC * (SinC * UinPae - CosC * VinPa)]) |
|
643 GT = np.dot(ATPe, IS1e) |
|
644 GR = np.dot(ARPe, IS1e) |
|
645 HT = np.dot(ATPe, IS2e) |
|
646 HR = np.dot(ARPe, IS2e) |
|
647 else: |
|
648 print("Calibrator not implemented yet") |
|
649 sys.exit() |
|
650 |
|
651 elif LocC == 3: # C before receiver optics Eq.57 |
|
652 |
|
653 # S2ge = np.sin(np.deg2rad(2*RotO - 2*RotC)) |
|
654 # C2ge = np.cos(np.deg2rad(2*RotO - 2*RotC)) |
|
655 S2e = np.sin(np.deg2rad(2. * RotC)) |
|
656 C2e = np.cos(np.deg2rad(2. * RotC)) |
|
657 |
|
658 # As with C before the receiver optics (rotated_diattenuator_X22x5deg.odt) |
|
659 AF1 = np.array([1., C2g * DiO, S2g * DiO, 0.]) |
|
660 AF2 = np.array([C2g * DiO, 1. - S2g ** 2 * WiO, S2g * C2g * WiO, -S2g * ZiO * SinO]) |
|
661 AF3 = np.array([S2g * DiO, S2g * C2g * WiO, 1. - C2g ** 2 * WiO, C2g * ZiO * SinO]) |
|
662 AF4 = np.array([0., S2g * SinO, -C2g * SinO, CosO]) |
|
663 |
|
664 ATF = (ATP1 * AF1 + ATP2 * AF2 + ATP3 * AF3 + ATP4 * AF4) |
|
665 ARF = (ARP1 * AF1 + ARP2 * AF2 + ARP3 * AF3 + ARP4 * AF4) |
|
666 ATF2 = ATF[1] |
|
667 ATF3 = ATF[2] |
|
668 ARF2 = ARF[1] |
|
669 ARF3 = ARF[2] |
|
670 |
|
671 # rotated AinF by epsilon |
|
672 ATF1 = ATF[0] |
|
673 ATF4 = ATF[3] |
|
674 ATF2e = C2e * ATF[1] + S2e * ATF[2] |
|
675 ATF3e = C2e * ATF[2] - S2e * ATF[1] |
|
676 ARF1 = ARF[0] |
|
677 ARF4 = ARF[3] |
|
678 ARF2e = C2e * ARF[1] + S2e * ARF[2] |
|
679 ARF3e = C2e * ARF[2] - S2e * ARF[1] |
|
680 |
|
681 ATFe = np.array([ATF1, ATF2e, ATF3e, ATF4]) |
|
682 ARFe = np.array([ARF1, ARF2e, ARF3e, ARF4]) |
|
683 |
|
684 QinEe = C2e * QinE + S2e * UinE |
|
685 UinEe = C2e * UinE - S2e * QinE |
|
686 |
|
687 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) |
|
688 IinF = IinE |
|
689 QinF = aCal * QinE |
|
690 UinF = -aCal * UinE |
|
691 VinF = (1. - 2. * aCal) * VinE |
|
692 |
|
693 IinFo = IinE |
|
694 QinFo = 0. |
|
695 UinFo = 0. |
|
696 VinFo = VinE |
|
697 |
|
698 IinFa = 0. |
|
699 QinFa = QinE |
|
700 UinFa = -UinE |
|
701 VinFa = -2. * VinE |
|
702 |
|
703 # Stokes Input Vector before receiver optics rotated by epsilon Eq. C.3 |
|
704 QinFe = C2e * QinF + S2e * UinF |
|
705 UinFe = C2e * UinF - S2e * QinF |
|
706 QinFoe = C2e * QinFo + S2e * UinFo |
|
707 UinFoe = C2e * UinFo - S2e * QinFo |
|
708 QinFae = C2e * QinFa + S2e * UinFa |
|
709 UinFae = C2e * UinFa - S2e * QinFa |
|
710 |
|
711 # Calibration signals and Calibration correction K from measurements with LDRCal / aCal |
|
712 if (TypeC == 2) or (TypeC == 1): # rotator calibration Eq. C.4 |
|
713 # parameters for calibration with aCal |
|
714 AT = ATF1 * IinF + ATF4 * h * VinF |
|
715 BT = ATF3e * QinF - ATF2e * h * UinF |
|
716 AR = ARF1 * IinF + ARF4 * h * VinF |
|
717 BR = ARF3e * QinF - ARF2e * h * UinF |
|
718 # Correction parameters for normal measurements; they are independent of LDR |
|
719 if (not RotationErrorEpsilonForNormalMeasurements): |
|
720 GT = ATF1 * IinE + ATF4 * VinE |
|
721 GR = ARF1 * IinE + ARF4 * VinE |
|
722 HT = ATF2 * QinE - ATF3 * UinE - ATF4 * 2 * VinE |
|
723 HR = ARF2 * QinE - ARF3 * UinE - ARF4 * 2 * VinE |
|
724 else: |
|
725 GT = ATF1 * IinE + ATF4 * h * VinE |
|
726 GR = ARF1 * IinE + ARF4 * h * VinE |
|
727 HT = ATF2e * QinE - ATF3e * h * UinE - ATF4 * h * 2 * VinE |
|
728 HR = ARF2e * QinE - ARF3e * h * UinE - ARF4 * h * 2 * VinE |
|
729 elif (TypeC == 3) or (TypeC == 4): # linear polariser calibration Eq. C.5 |
|
730 # p = +45°, m = -45° |
|
731 IF1e = np.array([IinF, ZiC * CosC * QinFe, UinFe, ZiC * CosC * VinF]) |
|
732 IF2e = np.array([DiC * UinFe, -ZiC * SinC * VinF, DiC * IinF, ZiC * SinC * QinFe]) |
|
733 AT = np.dot(ATFe, IF1e) |
|
734 AR = np.dot(ARFe, IF1e) |
|
735 BT = np.dot(ATFe, IF2e) |
|
736 BR = np.dot(ARFe, IF2e) |
|
737 |
|
738 # Correction parameters for normal measurements; they are independent of LDR --- the same as for TypeC = 6 |
|
739 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out |
|
740 IS1 = np.array([IinE, 0., 0., VinE]) |
|
741 IS2 = np.array([0., QinE, -UinE, -2. * VinE]) |
|
742 GT = np.dot(ATF, IS1) |
|
743 GR = np.dot(ARF, IS1) |
|
744 HT = np.dot(ATF, IS2) |
|
745 HR = np.dot(ARF, IS2) |
|
746 else: |
|
747 IS1e = np.array([IinFo + DiC * QinFoe, DiC * IinFo + QinFoe, ZiC * (CosC * UinFoe + SinC * VinFo), |
|
748 -ZiC * (SinC * UinFoe - CosC * VinFo)]) |
|
749 IS2e = np.array([IinFa + DiC * QinFae, DiC * IinFa + QinFae, ZiC * (CosC * UinFae + SinC * VinFa), |
|
750 -ZiC * (SinC * UinFae - CosC * VinFa)]) |
|
751 GT = np.dot(ATFe, IS1e) |
|
752 GR = np.dot(ARFe, IS1e) |
|
753 HT = np.dot(ATFe, IS2e) |
|
754 HR = np.dot(ARFe, IS2e) |
|
755 |
|
756 elif (TypeC == 6): # diattenuator calibration +-22.5° rotated_diattenuator_X22x5deg.odt |
|
757 # parameters for calibration with aCal |
|
758 IF1e = np.array([IinF + sqr05 * DiC * QinFe, sqr05 * DiC * IinF + (1. - 0.5 * WiC) * QinFe, |
|
759 (1. - 0.5 * WiC) * UinFe + sqr05 * ZiC * SinC * VinF, |
|
760 -sqr05 * ZiC * SinC * UinFe + ZiC * CosC * VinF]) |
|
761 IF2e = np.array([sqr05 * DiC * UinFe, 0.5 * WiC * UinFe - sqr05 * ZiC * SinC * VinF, |
|
762 sqr05 * DiC * IinF + 0.5 * WiC * QinFe, sqr05 * ZiC * SinC * QinFe]) |
|
763 AT = np.dot(ATFe, IF1e) |
|
764 AR = np.dot(ARFe, IF1e) |
|
765 BT = np.dot(ATFe, IF2e) |
|
766 BR = np.dot(ARFe, IF2e) |
|
767 |
|
768 # Correction parameters for normal measurements; they are independent of LDR |
|
769 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out |
|
770 # IS1 = np.array([IinE,0,0,VinE]) |
|
771 # IS2 = np.array([0,QinE,-UinE,-2*VinE]) |
|
772 IS1 = np.array([IinFo, 0., 0., VinFo]) |
|
773 IS2 = np.array([0., QinFa, UinFa, VinFa]) |
|
774 GT = np.dot(ATF, IS1) |
|
775 GR = np.dot(ARF, IS1) |
|
776 HT = np.dot(ATF, IS2) |
|
777 HR = np.dot(ARF, IS2) |
|
778 else: |
|
779 IS1e = np.array([IinFo + DiC * QinFoe, DiC * IinFo + QinFoe, ZiC * (CosC * UinFoe + SinC * VinFo), |
|
780 -ZiC * (SinC * UinFoe - CosC * VinFo)]) |
|
781 IS2e = np.array([IinFa + DiC * QinFae, DiC * IinFa + QinFae, ZiC * (CosC * UinFae + SinC * VinFa), |
|
782 -ZiC * (SinC * UinFae - CosC * VinFa)]) |
|
783 # IS1e = np.array([IinFo,0,0,VinFo]) |
|
784 # IS2e = np.array([0,QinFae,UinFae,VinFa]) |
|
785 GT = np.dot(ATFe, IS1e) |
|
786 GR = np.dot(ARFe, IS1e) |
|
787 HT = np.dot(ATFe, IS2e) |
|
788 HR = np.dot(ARFe, IS2e) |
|
789 |
|
790 else: |
|
791 print('Calibrator not implemented yet') |
|
792 sys.exit() |
|
793 |
|
794 elif LocC == 2: # C behind emitter optics Eq.57 ------------------------------------------------------- |
|
795 # print("Calibrator location not implemented yet") |
|
796 S2e = np.sin(np.deg2rad(2. * RotC)) |
|
797 C2e = np.cos(np.deg2rad(2. * RotC)) |
|
798 |
|
799 # AS with C before the receiver optics (see document rotated_diattenuator_X22x5deg.odt) |
|
800 AF1 = np.array([1, C2g * DiO, S2g * DiO, 0.]) |
|
801 AF2 = np.array([C2g * DiO, 1. - S2g ** 2 * WiO, S2g * C2g * WiO, -S2g * ZiO * SinO]) |
|
802 AF3 = np.array([S2g * DiO, S2g * C2g * WiO, 1. - C2g ** 2 * WiO, C2g * ZiO * SinO]) |
|
803 AF4 = np.array([0., S2g * SinO, -C2g * SinO, CosO]) |
|
804 |
|
805 ATF = (ATP1 * AF1 + ATP2 * AF2 + ATP3 * AF3 + ATP4 * AF4) |
|
806 ARF = (ARP1 * AF1 + ARP2 * AF2 + ARP3 * AF3 + ARP4 * AF4) |
|
807 ATF1 = ATF[0] |
|
808 ATF2 = ATF[1] |
|
809 ATF3 = ATF[2] |
|
810 ATF4 = ATF[3] |
|
811 ARF1 = ARF[0] |
|
812 ARF2 = ARF[1] |
|
813 ARF3 = ARF[2] |
|
814 ARF4 = ARF[3] |
|
815 |
|
816 # AS with C behind the emitter |
|
817 # terms without aCal |
|
818 ATE1o, ARE1o = ATF1, ARF1 |
|
819 ATE2o, ARE2o = 0., 0. |
|
820 ATE3o, ARE3o = 0., 0. |
|
821 ATE4o, ARE4o = ATF4, ARF4 |
|
822 # terms with aCal |
|
823 ATE1a, ARE1a = 0., 0. |
|
824 ATE2a, ARE2a = ATF2, ARF2 |
|
825 ATE3a, ARE3a = -ATF3, -ARF3 |
|
826 ATE4a, ARE4a = -2. * ATF4, -2. * ARF4 |
|
827 # rotated AinEa by epsilon |
|
828 ATE2ae = C2e * ATF2 + S2e * ATF3 |
|
829 ATE3ae = -S2e * ATF2 - C2e * ATF3 |
|
830 ARE2ae = C2e * ARF2 + S2e * ARF3 |
|
831 ARE3ae = -S2e * ARF2 - C2e * ARF3 |
|
832 |
|
833 ATE1 = ATE1o |
|
834 ATE2e = aCal * ATE2ae |
|
835 ATE3e = aCal * ATE3ae |
|
836 ATE4 = (1 - 2 * aCal) * ATF4 |
|
837 ARE1 = ARE1o |
|
838 ARE2e = aCal * ARE2ae |
|
839 ARE3e = aCal * ARE3ae |
|
840 ARE4 = (1 - 2 * aCal) * ARF4 |
|
841 |
|
842 # rotated IinE |
|
843 QinEe = C2e * QinE + S2e * UinE |
|
844 UinEe = C2e * UinE - S2e * QinE |
|
845 |
|
846 # Calibration signals and Calibration correction K from measurements with LDRCal / aCal |
|
847 if (TypeC == 2) or (TypeC == 1): # +++++++++ rotator calibration Eq. C.4 |
|
848 AT = ATE1o * IinE + (ATE4o + aCal * ATE4a) * h * VinE |
|
849 BT = aCal * (ATE3ae * QinEe - ATE2ae * h * UinEe) |
|
850 AR = ARE1o * IinE + (ARE4o + aCal * ARE4a) * h * VinE |
|
851 BR = aCal * (ARE3ae * QinEe - ARE2ae * h * UinEe) |
|
852 |
|
853 # Correction parameters for normal measurements; they are independent of LDR |
|
854 if (not RotationErrorEpsilonForNormalMeasurements): |
|
855 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) |
|
856 GT = ATE1o * IinE + ATE4o * h * VinE |
|
857 GR = ARE1o * IinE + ARE4o * h * VinE |
|
858 HT = ATE2a * QinE + ATE3a * h * UinEe + ATE4a * h * VinE |
|
859 HR = ARE2a * QinE + ARE3a * h * UinEe + ARE4a * h * VinE |
|
860 else: |
|
861 GT = ATE1o * IinE + ATE4o * h * VinE |
|
862 GR = ARE1o * IinE + ARE4o * h * VinE |
|
863 HT = ATE2ae * QinE + ATE3ae * h * UinEe + ATE4a * h * VinE |
|
864 HR = ARE2ae * QinE + ARE3ae * h * UinEe + ARE4a * h * VinE |
|
865 |
|
866 elif (TypeC == 3) or (TypeC == 4): # +++++++++ linear polariser calibration Eq. C.5 |
|
867 # p = +45°, m = -45° |
|
868 AT = ATE1 * IinE + ZiC * CosC * (ATE2e * QinEe + ATE4 * VinE) + ATE3e * UinEe |
|
869 BT = DiC * (ATE1 * UinEe + ATE3e * IinE) + ZiC * SinC * (ATE4 * QinEe - ATE2e * VinE) |
|
870 AR = ARE1 * IinE + ZiC * CosC * (ARE2e * QinEe + ARE4 * VinE) + ARE3e * UinEe |
|
871 BR = DiC * (ARE1 * UinEe + ARE3e * IinE) + ZiC * SinC * (ARE4 * QinEe - ARE2e * VinE) |
|
872 |
|
873 # Correction parameters for normal measurements; they are independent of LDR |
|
874 if (not RotationErrorEpsilonForNormalMeasurements): |
|
875 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) |
|
876 GT = ATE1o * IinE + ATE4o * VinE |
|
877 GR = ARE1o * IinE + ARE4o * VinE |
|
878 HT = ATE2a * QinE + ATE3a * UinE + ATE4a * VinE |
|
879 HR = ARE2a * QinE + ARE3a * UinE + ARE4a * VinE |
|
880 else: |
|
881 D = IinE + DiC * QinEe |
|
882 A = DiC * IinE + QinEe |
|
883 B = ZiC * (CosC * UinEe + SinC * VinE) |
|
884 C = -ZiC * (SinC * UinEe - CosC * VinE) |
|
885 GT = ATE1o * D + ATE4o * C |
|
886 GR = ARE1o * D + ARE4o * C |
|
887 HT = ATE2a * A + ATE3a * B + ATE4a * C |
|
888 HR = ARE2a * A + ARE3a * B + ARE4a * C |
|
889 |
|
890 elif (TypeC == 6): # real HWP calibration +-22.5° rotated_diattenuator_X22x5deg.odt |
|
891 # p = +22.5°, m = -22.5° |
|
892 IE1e = np.array([IinE + sqr05 * DiC * QinEe, sqr05 * DiC * IinE + (1 - 0.5 * WiC) * QinEe, |
|
893 (1 - 0.5 * WiC) * UinEe + sqr05 * ZiC * SinC * VinE, |
|
894 -sqr05 * ZiC * SinC * UinEe + ZiC * CosC * VinE]) |
|
895 IE2e = np.array([sqr05 * DiC * UinEe, 0.5 * WiC * UinEe - sqr05 * ZiC * SinC * VinE, |
|
896 sqr05 * DiC * IinE + 0.5 * WiC * QinEe, sqr05 * ZiC * SinC * QinEe]) |
|
897 ATEe = np.array([ATE1, ATE2e, ATE3e, ATE4]) |
|
898 AREe = np.array([ARE1, ARE2e, ARE3e, ARE4]) |
|
899 AT = np.dot(ATEe, IE1e) |
|
900 AR = np.dot(AREe, IE1e) |
|
901 BT = np.dot(ATEe, IE2e) |
|
902 BR = np.dot(AREe, IE2e) |
|
903 |
|
904 # Correction parameters for normal measurements; they are independent of LDR |
|
905 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out |
|
906 GT = ATE1o * IinE + ATE4o * VinE |
|
907 GR = ARE1o * IinE + ARE4o * VinE |
|
908 HT = ATE2a * QinE + ATE3a * UinE + ATE4a * VinE |
|
909 HR = ARE2a * QinE + ARE3a * UinE + ARE4a * VinE |
|
910 else: |
|
911 D = IinE + DiC * QinEe |
|
912 A = DiC * IinE + QinEe |
|
913 B = ZiC * (CosC * UinEe + SinC * VinE) |
|
914 C = -ZiC * (SinC * UinEe - CosC * VinE) |
|
915 GT = ATE1o * D + ATE4o * C |
|
916 GR = ARE1o * D + ARE4o * C |
|
917 HT = ATE2a * A + ATE3a * B + ATE4a * C |
|
918 HR = ARE2a * A + ARE3a * B + ARE4a * C |
|
919 |
|
920 else: |
|
921 print('Calibrator not implemented yet') |
|
922 sys.exit() |
|
923 |
|
924 else: |
|
925 print("Calibrator location not implemented yet") |
|
926 sys.exit() |
|
927 |
|
928 # Determination of the correction K of the calibration factor. |
|
929 IoutTp = TTa * TiC * TiO * TiE * (AT + BT) |
|
930 IoutTm = TTa * TiC * TiO * TiE * (AT - BT) |
|
931 IoutRp = TRa * TiC * TiO * TiE * (AR + BR) |
|
932 IoutRm = TRa * TiC * TiO * TiE * (AR - BR) |
|
933 # --- Results and Corrections; electronic etaR and etaT are assumed to be 1 |
|
934 Etapx = IoutRp / IoutTp |
|
935 Etamx = IoutRm / IoutTm |
|
936 Etax = (Etapx * Etamx) ** 0.5 |
|
937 |
|
938 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 |
|
939 K = Etax / Eta |
|
940 |
|
941 # For comparison with Volkers Libreoffice Müller Matrix spreadsheet |
|
942 # Eta_test_p = (IoutRp/IoutTp) |
|
943 # Eta_test_m = (IoutRm/IoutTm) |
|
944 # Eta_test = (Eta_test_p*Eta_test_m)**0.5 |
|
945 |
|
946 # ----- random error calculation ---------- |
|
947 # noise must be calculated with the photon counts of measured signals; |
|
948 # relative standard deviation of calibration signals with LDRcal; assumed to be statisitcally independent |
|
949 # normalised noise errors |
|
950 if (CalcFrom0deg): |
|
951 dIoutTp = (NCalT * IoutTp) ** -0.5 |
|
952 dIoutTm = (NCalT * IoutTm) ** -0.5 |
|
953 dIoutRp = (NCalR * IoutRp) ** -0.5 |
|
954 dIoutRm = (NCalR * IoutRm) ** -0.5 |
|
955 else: |
|
956 dIoutTp = (NCalT ** -0.5) |
|
957 dIoutTm = (NCalT ** -0.5) |
|
958 dIoutRp = (NCalR ** -0.5) |
|
959 dIoutRm = (NCalR ** -0.5) |
|
960 # Forward simulated 0°-signals with LDRCal with atrue; from input file |
|
961 |
|
962 It = TTa * TiO * TiE * (GT + atrue * HT) |
|
963 Ir = TRa * TiO * TiE * (GR + atrue * HR) |
|
964 # relative standard deviation of standard signals with LDRmeas; assumed to be statisitcally independent |
|
965 if (CalcFrom0deg): # this works! |
|
966 dIt = ((It * NI * eFacT) ** -0.5) |
|
967 dIr = ((Ir * NI * eFacR) ** -0.5) |
|
968 ''' |
|
969 dIt = ((NCalT * It / IoutTp * NILfac / TCalT) ** -0.5) |
|
970 dIr = ((NCalR * Ir / IoutRp * NILfac / TCalR) ** -0.5) |
|
971 ''' |
|
972 else: # does this work? Why not as above? |
|
973 dIt = ((NCalT * 2 * NILfac / TCalT ) ** -0.5) |
|
974 dIr = ((NCalR * 2 * NILfac / TCalR) ** -0.5) |
|
975 |
|
976 # ----- Forward simulated LDRsim = 1/Eta*Ir/It # simulated LDR* with Y from input file |
|
977 LDRsim = Ir / It # simulated uncorrected LDR with Y from input file |
|
978 # Corrected LDRsimCorr from forward simulated LDRsim (atrue) |
|
979 # LDRsimCorr = (1./Eta*LDRsim*(GT+HT)-(GR+HR))/((GR-HR)-1./Eta*LDRsim*(GT-HT)) |
|
980 ''' |
|
981 if ((Y == -1.) and (abs(RotL0) < 45)) or ((Y == +1.) and (abs(RotL0) > 45)): |
|
982 LDRsimx = 1. / LDRsim / Etax |
|
983 else: |
|
984 LDRsimx = LDRsim / Etax |
|
985 ''' |
|
986 LDRsimx = LDRsim |
|
987 |
|
988 # The following is correct without doubt |
|
989 # LDRCorr = (LDRsim/(Etax/K)*(GT+HT)-(GR+HR))/((GR-HR)-LDRsim/(Etax/K)*(GT-HT)) |
|
990 |
|
991 # The following is a test whether the equations for calibration Etax and normal signal (GHK, LDRsim) are consistent |
|
992 LDRCorr = (LDRsim / (Etax / K) * (GT + HT) - (GR + HR)) / ((GR - HR) - LDRsim / (Etax / K) * (GT - HT)) |
|
993 # here we could also use Eta instead of Etax / K => how to test whether Etax is correct? => comparison with MüllerMatrix simulation! |
|
994 # Without any correction: only measured It, Ir, EtaX are used |
|
995 LDRunCorr = LDRsim / Etax |
|
996 # 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))) |
|
997 |
|
998 #LDRCorr = LDRsimx # for test only |
|
999 |
|
1000 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 |
|
1001 |
|
1002 return (IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr, |
|
1003 GT, HT, GR, HR, K, Eta, LDRsimx, LDRCorr, DTa, DRa, TTa, TRa, F11sim, LDRunCorr) |
|
1004 |
|
1005 |
|
1006 |
|
1007 # ******************************************************************************************************************************* |
|
1008 |
|
1009 # --- CALC with assumed true parameters from the input file |
|
1010 LDRtrue = LDRtrue2 |
|
1011 IoutTp0, IoutTm0, IoutRp0, IoutRm0, It0, Ir0, dIoutTp0, dIoutTm0, dIoutRp0, dIoutRm0, dIt0, dIr0, \ |
|
1012 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \ |
|
1013 Calc(TCalT, TCalR, NCalT, NCalR, Qin0, Vin0, RotL0, RotE0, RetE0, DiE0, |
|
1014 RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, |
|
1015 ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0) |
|
1016 Etax0 = K0 * Eta0 |
|
1017 Etapx0 = IoutRp0 / IoutTp0 |
|
1018 Etamx0 = IoutRm0 / IoutTm0 |
|
1019 # --- Print parameters to console and output file |
|
1020 OutputFile = 'output_' + InputFile[0:-3] + '_' + fname[0:-3] +'.dat' |
|
1021 with open('output_files\\' + OutputFile, 'w') as f: |
|
1022 with redirect_stdout(f): |
|
1023 print("From ", dname) |
|
1024 print("Running ", fname) |
|
1025 print("Reading input file ", InputFile) # , " for Lidar system :", EID, ", ", LID) |
|
1026 print("for Lidar system: ", EID, ", ", LID) |
|
1027 # --- Print iput information********************************* |
|
1028 print(" --- Input parameters: value ±error / ±steps ----------------------") |
|
1029 print("{0:7}{1:17} {2:6.4f}±{3:7.4f}/{4:2d}".format("Laser: ", "Qin =", Qin0, dQin, nQin)) |
|
1030 print("{0:7}{1:17} {2:6.4f}±{3:7.4f}/{4:2d}".format("", "Vin =", Vin0, dVin, nVin)) |
|
1031 print("{0:7}{1:17} {2:6.4f}±{3:7.4f}/{4:2d}".format("", "Rotation alpha = ", RotL0, dRotL, nRotL)) |
|
1032 print("{0:7}{1:15} {2:8.4f} {3:17}".format("", "=> DOP", ((Qin ** 2 + Vin ** 2) ** 0.5), " (degree of polarisation)")) |
|
1033 |
|
1034 print("Optic: Diatt., Tunpol, Retard., Rotation (deg)") |
|
1035 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( |
|
1036 "Emitter ", DiE0, dDiE, TiE, RetE0, dRetE, RotE0, dRotE, nDiE, nRetE, nRotE)) |
|
1037 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( |
|
1038 "Receiver ", DiO0, dDiO, TiO, RetO0, dRetO, RotO0, dRotO, nDiO, nRetO, nRotO)) |
|
1039 print("{0:12} {1:9.6f}±{2:9.6f}/{8:2d}, {3:7.4f}, {4:3.0f}±{5:3.0f}/{9:2d}, {6:7.4f}±{7:7.4f}/{10:2d}".format( |
|
1040 "Calibrator ", DiC0, dDiC, TiC, RetC0, dRetC, RotC0, dRotC, nDiC, nRetC, nRotC)) |
|
1041 print("{0:12}".format(" Pol.-filter ------ ")) |
|
1042 print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format( |
|
1043 "ERT, RotT :", ERaT0, dERaT, nERaT, RotaT0, dRotaT, nRotaT)) |
|
1044 print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format( |
|
1045 "ERR, RotR :", ERaR0, dERaR, nERaR, RotaR0, dRotaR, nRotaR)) |
|
1046 print("{0:12}".format(" PBS ------ ")) |
|
1047 print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format( |
|
1048 "TP,TS :", TP0, dTP, nTP, TS0, dTS, nTS)) |
|
1049 print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format( |
|
1050 "RP,RS :", RP0, dRP, nRP, RS0, dRS, nRS)) |
|
1051 print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}, {5:1.0f}".format( |
|
1052 "DT,TT,DR,TR,Y :", DiT, TiT, DiR, TiR, Y)) |
|
1053 print("{0:12}".format(" Combined PBS + Pol.-filter ------ ")) |
|
1054 print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}".format( |
|
1055 "DT,TT,DR,TR :", DTa0, TTa0, DRa0, TRa0)) |
|
1056 print("{0:26}: {1:6.3f}± {2:5.3f}/{3:2d}".format( |
|
1057 "LDRCal during calibration in calibration range", LDRCal0, dLDRCal, nLDRCal)) |
|
1058 print("{0:12}".format(" --- Additional ND filter attenuation (transmission) during the calibration ---")) |
|
1059 print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format( |
|
1060 "TCalT,TCalR :", TCalT0, dTCalT, nTCalT, TCalR0, dTCalR, nTCalR)) |
|
1061 print() |
|
1062 print("Rotation Error Epsilon For Normal Measurements = ", RotationErrorEpsilonForNormalMeasurements) |
|
1063 print(Type[TypeC], Loc[LocC]) |
|
1064 print("PBS incidence plane is ", dY[int(Y + 1)], "to reference plane and polarisation in reference plane is finally", dY2[int(Y + 1)]) |
|
1065 print(dY3) |
|
1066 print("RS_RP_depend_on_TS_TP = ", RS_RP_depend_on_TS_TP) |
|
1067 # end of print actual system parameters |
|
1068 # ****************************************************************************** |
|
1069 |
|
1070 |
|
1071 print() |
|
1072 |
|
1073 K0List = np.zeros(7) |
|
1074 LDRsimxList = np.zeros(7) |
|
1075 LDRCalList = 0.0, 0.004, 0.02, 0.1, 0.2, 0.3, 0.45 |
|
1076 # 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. |
|
1077 # Still with assumed true parameters in input file |
|
1078 |
|
1079 ''' |
|
1080 facIt = NCalT / TCalT0 * NILfac |
|
1081 facIr = NCalR / TCalR0 * NILfac |
|
1082 ''' |
|
1083 facIt = NI * eFacT |
|
1084 facIr = NI * eFacR |
|
1085 if (bPlotEtax): |
|
1086 # check error signals |
|
1087 # dIs are relative stdevs |
|
1088 print("LDRCal, IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp,dIoutTm,dIoutRp,dIoutRm,dIt, dIr") |
|
1089 |
|
1090 for i, LDRCal in enumerate(LDRCalList): |
|
1091 IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr, \ |
|
1092 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \ |
|
1093 Calc(TCalT0, TCalR0, NCalT, NCalR, Qin0, Vin0, RotL0, RotE0, RetE0, DiE0, |
|
1094 RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, |
|
1095 ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal) |
|
1096 K0List[i] = K0 |
|
1097 LDRsimxList[i] = LDRsimx |
|
1098 |
|
1099 if (bPlotEtax): |
|
1100 # check error signals |
|
1101 print( "{:0.2f}, {: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(LDRCal, IoutTp * NCalT, IoutTm * NCalT, IoutRp * NCalR, IoutRm * NCalR, It * facIt, Ir * facIr, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr)) |
|
1102 #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)) |
|
1103 # end check error signals |
|
1104 print('===========================================================================================================') |
|
1105 print("{0:8},{1:8},{2:8},{3:8},{4:9},{5:8},{6:9},{7:9},{8:9},{9:9},{10:9}".format( |
|
1106 " GR", " GT", " HR", " HT", " K(0.000)", " K(0.004)", " K(0.02)", " K(0.1)", " K(0.2)", " K(0.3)", " K(0.45)")) |
|
1107 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},{10:9.5f}".format( |
|
1108 GR0, GT0, HR0, HT0, K0List[0], K0List[1], K0List[2], K0List[3], K0List[4], K0List[5], K0List[6])) |
|
1109 print('===========================================================================================================') |
|
1110 print() |
|
1111 print("Errors from neglecting GHK corrections and/or calibration:") |
|
1112 print("{0:>10},{1:>10},{2:>10},{3:>10},{4:>10},{5:>10}".format( |
|
1113 "LDRtrue", "LDRunCorr", "1/LDRunCorr", "LDRsimx", "1/LDRsimx", "LDRCorr")) |
|
1114 |
|
1115 aF11sim0 = np.zeros(5) |
|
1116 LDRrange = np.zeros(5) |
|
1117 LDRsim0 = np.zeros(5) |
|
1118 LDRrange = [0.004, 0.02, 0.1, 0.3, 0.45] # list |
|
1119 LDRrange[0] = LDRtrue2 # value in the input file; default 0.004 |
|
1120 |
|
1121 # The loop over LDRtrueList is only for checking how much the uncorrected LDRsimx deviates from LDRtrue ... and whether the corrections work. |
|
1122 # LDRsimx = LDRsim = Ir / It or 1/LDRsim |
|
1123 # Still with assumed true parameters in input file |
|
1124 for i, LDRtrue in enumerate(LDRrange): |
|
1125 #for LDRtrue in LDRrange: |
|
1126 IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr, \ |
|
1127 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \ |
|
1128 Calc(TCalT0, TCalR0, NCalT, NCalR, Qin0, Vin0, RotL0, RotE0, RetE0, DiE0, |
|
1129 RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, |
|
1130 ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0) |
|
1131 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)) |
|
1132 aF11sim0[i] = F11sim0 |
|
1133 LDRsim0[i] = Ir / It |
|
1134 # the assumed true aF11sim0 results will be used below to calc the deviation from the real signals |
|
1135 print("LDRsimx = LDR of the nominal system directly from measured signals without calibration and GHK-corrections") |
|
1136 print("LDRunCorr = LDR of the nominal system directly from measured signals with calibration but without GHK-corrections; electronic amplifications = 1 assumed") |
|
1137 print("LDRCorr = LDR calibrated and GHK-corrected") |
|
1138 print() |
|
1139 print("Errors from signal noise:") |
|
1140 print("Signal counts: NI, NCalT, NCalR, NILfac, nNCal, nNI, stdev(NI)/NI = {0:10.0f},{1:10.0f},{2:10.0f},{3:3.0f},{4:2.0f},{5:2.0f},{6:8.5f}".format( |
|
1141 NI, NCalT, NCalR, NILfac, nNCal, nNI, 1.0 / NI ** 0.5)) |
|
1142 print() |
|
1143 print() |
|
1144 '''# das muß wieder weg |
|
1145 print("IoutTp, IoutTm, IoutRp, IoutRm, It , Ir , dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr") |
|
1146 LDRCal = 0.01 |
|
1147 for i, LDRtrue in enumerate(LDRrange): |
|
1148 IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr, \ |
|
1149 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \ |
|
1150 Calc(TCalT0, TCalR0, NCalT, NCalR, DOLP0, RotL0, RotE0, RetE0, DiE0, |
|
1151 RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, |
|
1152 ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0) |
|
1153 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( |
|
1154 IoutTp * NCalT, IoutTm * NCalT, IoutRp * NCalR, IoutRm * NCalR, It * facIt, Ir * facIr, |
|
1155 dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr)) |
|
1156 aF11sim0[i] = F11sim0 |
|
1157 # the assumed true aF11sim0 results will be used below to calc the deviation from the real signals |
|
1158 # bis hierher weg |
|
1159 ''' |
|
1160 |
|
1161 file = open('output_files\\' + OutputFile, 'r') |
|
1162 print(file.read()) |
|
1163 file.close() |
|
1164 |
|
1165 # --- CALC again assumed truth with LDRCal0 and with assumed true parameters in input file to reset all 0-values |
|
1166 LDRtrue = LDRtrue2 |
|
1167 IoutTp0, IoutTm0, IoutRp0, IoutRm0, It0, Ir0, dIoutTp0, dIoutTm0, dIoutRp0, dIoutRm0, dIt0, dIr0, \ |
|
1168 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \ |
|
1169 Calc(TCalT0, TCalR0, NCalT, NCalR, Qin0, Vin0, RotL0, RotE0, RetE0, DiE0, |
|
1170 RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, |
|
1171 ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0) |
|
1172 Etax0 = K0 * Eta0 |
|
1173 Etapx0 = IoutRp0 / IoutTp0 |
|
1174 Etamx0 = IoutRm0 / IoutTm0 |
|
1175 ''' |
|
1176 if(PrintToOutputFile): |
|
1177 f = open('output_ver7.dat', 'w') |
|
1178 old_target = sys.stdout |
|
1179 sys.stdout = f |
|
1180 |
|
1181 print("something") |
|
1182 |
|
1183 if(PrintToOutputFile): |
|
1184 sys.stdout.flush() |
|
1185 f.close |
|
1186 sys.stdout = old_target |
|
1187 ''' |
|
1188 if (Error_Calc): |
|
1189 # --- CALC again assumed truth with LDRCal0 and with assumed true parameters in input file to reset all 0-values |
|
1190 LDRtrue = LDRtrue2 |
|
1191 IoutTp0, IoutTm0, IoutRp0, IoutRm0, It0, Ir0, dIoutTp0, dIoutTm0, dIoutRp0, dIoutRm0, dIt0, dIr0, \ |
|
1192 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \ |
|
1193 Calc(TCalT0, TCalR0, NCalT, NCalR, Qin0, Vin0, RotL0, RotE0, RetE0, DiE0, |
|
1194 RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, |
|
1195 ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0) |
|
1196 Etax0 = K0 * Eta0 |
|
1197 Etapx0 = IoutRp0 / IoutTp0 |
|
1198 Etamx0 = IoutRm0 / IoutTm0 |
|
1199 |
|
1200 # --- Start Error calculation with variable parameters ------------------------------------------------------------------ |
|
1201 # error nNCal: one-sigma in steps to left and right for calibration signals |
|
1202 # error nNI: one-sigma in steps to left and right for 0° signals |
|
1203 |
|
1204 iN = -1 |
|
1205 N = ((nTCalT * 2 + 1) * (nTCalR * 2 + 1) * |
|
1206 (nNCal * 2 + 1) ** 4 * (nNI * 2 + 1) ** 2 * |
|
1207 (nQin * 2 + 1) * (nVin * 2 + 1) * (nRotL * 2 + 1) * |
|
1208 (nRotE * 2 + 1) * (nRetE * 2 + 1) * (nDiE * 2 + 1) * |
|
1209 (nRotO * 2 + 1) * (nRetO * 2 + 1) * (nDiO * 2 + 1) * |
|
1210 (nRotC * 2 + 1) * (nRetC * 2 + 1) * (nDiC * 2 + 1) * |
|
1211 (nTP * 2 + 1) * (nTS * 2 + 1) * (nRP * 2 + 1) * (nRS * 2 + 1) * (nERaT * 2 + 1) * (nERaR * 2 + 1) * |
|
1212 (nRotaT * 2 + 1) * (nRotaR * 2 + 1) * (nRetT * 2 + 1) * (nRetR * 2 + 1) * (nLDRCal * 2 + 1)) |
|
1213 print("number of system variations N = ", N, " ", end="") |
|
1214 |
|
1215 if N > 1e6: |
|
1216 if user_yes_no_query('Warning: processing ' + str( |
|
1217 N) + ' samples will take very long. Do you want to proceed?') == 0: sys.exit() |
|
1218 if N > 5e6: |
|
1219 if user_yes_no_query('Warning: the memory required for ' + str(N) + ' samples might be ' + '{0:5.1f}'.format( |
|
1220 N / 4e6) + ' GB. Do you anyway want to proceed?') == 0: sys.exit() |
|
1221 |
|
1222 # if user_yes_no_query('Warning: processing' + str(N) + ' samples will take very long. Do you want to proceed?') == 0: sys.exit() |
|
1223 |
|
1224 # --- Arrays for plotting ------ |
|
1225 LDRmin = np.zeros(5) |
|
1226 LDRmax = np.zeros(5) |
|
1227 LDRstd = np.zeros(5) |
|
1228 LDRmean = np.zeros(5) |
|
1229 LDRmedian = np.zeros(5) |
|
1230 LDRskew = np.zeros(5) |
|
1231 LDRkurt = np.zeros(5) |
|
1232 LDRsimmin = np.zeros(5) |
|
1233 LDRsimmax = np.zeros(5) |
|
1234 LDRsimmean = np.zeros(5) |
|
1235 |
|
1236 F11min = np.zeros(5) |
|
1237 F11max = np.zeros(5) |
|
1238 Etaxmin = np.zeros(5) |
|
1239 Etaxmax = np.zeros(5) |
|
1240 |
|
1241 aQin = np.zeros(N) |
|
1242 aVin = np.zeros(N) |
|
1243 aERaT = np.zeros(N) |
|
1244 aERaR = np.zeros(N) |
|
1245 aRotaT = np.zeros(N) |
|
1246 aRotaR = np.zeros(N) |
|
1247 aRetT = np.zeros(N) |
|
1248 aRetR = np.zeros(N) |
|
1249 aTP = np.zeros(N) |
|
1250 aTS = np.zeros(N) |
|
1251 aRP = np.zeros(N) |
|
1252 aRS = np.zeros(N) |
|
1253 aDiE = np.zeros(N) |
|
1254 aDiO = np.zeros(N) |
|
1255 aDiC = np.zeros(N) |
|
1256 aRotC = np.zeros(N) |
|
1257 aRetC = np.zeros(N) |
|
1258 aRotL = np.zeros(N) |
|
1259 aRetE = np.zeros(N) |
|
1260 aRotE = np.zeros(N) |
|
1261 aRetO = np.zeros(N) |
|
1262 aRotO = np.zeros(N) |
|
1263 aLDRCal = np.zeros(N) |
|
1264 aNCalTp = np.zeros(N) |
|
1265 aNCalTm = np.zeros(N) |
|
1266 aNCalRp = np.zeros(N) |
|
1267 aNCalRm = np.zeros(N) |
|
1268 aNIt = np.zeros(N) |
|
1269 aNIr = np.zeros(N) |
|
1270 aTCalT = np.zeros(N) |
|
1271 aTCalR = np.zeros(N) |
|
1272 |
|
1273 # each np.zeros((LDRrange, N)) array has the same N-dependency |
|
1274 aLDRcorr = np.zeros((5, N)) |
|
1275 aLDRsim = np.zeros((5, N)) |
|
1276 aF11corr = np.zeros((5, N)) |
|
1277 aPLDR = np.zeros((5, N)) |
|
1278 aEtax = np.zeros((5, N)) |
|
1279 aEtapx = np.zeros((5, N)) |
|
1280 aEtamx = np.zeros((5, N)) |
|
1281 |
|
1282 # np.zeros((GHKs, N)) |
|
1283 aGHK = np.zeros((5, N)) |
|
1284 |
|
1285 atime = clock() |
|
1286 dtime = clock() |
|
1287 |
|
1288 # --- Calc Error signals |
|
1289 # ---- Do the calculations of bra-ket vectors |
|
1290 h = -1. if TypeC == 2 else 1 |
|
1291 |
|
1292 for iLDRCal in range(-nLDRCal, nLDRCal + 1): |
|
1293 # from input file: LDRCal for calibration measurements |
|
1294 LDRCal = LDRCal0 |
|
1295 if nLDRCal > 0: |
|
1296 LDRCal = LDRCal0 + iLDRCal * dLDRCal / nLDRCal |
|
1297 # provides the intensities of the calibration measurements at various LDRCal for signal noise errors |
|
1298 # IoutTp, IoutTm, IoutRp, IoutRm, dIoutTp, dIoutTm, dIoutRp, dIoutRm |
|
1299 |
|
1300 aCal = (1. - LDRCal) / (1. + LDRCal) |
|
1301 for iQin, iVin, iRotL, iRotE, iRetE, iDiE \ |
|
1302 in [(iQin, iVin, iRotL, iRotE, iRetE, iDiE) |
|
1303 for iQin in range(-nQin, nQin + 1) |
|
1304 for iVin in range(-nVin, nVin + 1) |
|
1305 for iRotL in range(-nRotL, nRotL + 1) |
|
1306 for iRotE in range(-nRotE, nRotE + 1) |
|
1307 for iRetE in range(-nRetE, nRetE + 1) |
|
1308 for iDiE in range(-nDiE, nDiE + 1)]: |
|
1309 |
|
1310 if nQin > 0: Qin = Qin0 + iQin * dQin / nQin |
|
1311 if nVin > 0: Vin = Vin0 + iVin * dVin / nVin |
|
1312 if nRotL > 0: RotL = RotL0 + iRotL * dRotL / nRotL |
|
1313 if nRotE > 0: RotE = RotE0 + iRotE * dRotE / nRotE |
|
1314 if nRetE > 0: RetE = RetE0 + iRetE * dRetE / nRetE |
|
1315 if nDiE > 0: DiE = DiE0 + iDiE * dDiE / nDiE |
|
1316 |
|
1317 if ((Qin ** 2 + Vin ** 2) ** 0.5) > 1.0: |
|
1318 print("Error: degree of polarisation of laser > 1. Check Qin and Vin! ") |
|
1319 sys.exit() |
|
1320 # angles of emitter and laser and calibrator and receiver optics |
|
1321 # RotL = alpha, RotE = beta, RotO = gamma, RotC = epsilon |
|
1322 S2a = np.sin(2 * np.deg2rad(RotL)) |
|
1323 C2a = np.cos(2 * np.deg2rad(RotL)) |
|
1324 S2b = np.sin(2 * np.deg2rad(RotE)) |
|
1325 C2b = np.cos(2 * np.deg2rad(RotE)) |
|
1326 S2ab = np.sin(np.deg2rad(2 * RotL - 2 * RotE)) |
|
1327 C2ab = np.cos(np.deg2rad(2 * RotL - 2 * RotE)) |
|
1328 |
|
1329 # Laser with Degree of linear polarization DOLP |
|
1330 IinL = 1. |
|
1331 QinL = Qin |
|
1332 UinL = 0. |
|
1333 VinL = Vin |
|
1334 # VinL = (1. - DOLP ** 2) ** 0.5 |
|
1335 |
|
1336 # Stokes Input Vector rotation Eq. E.4 |
|
1337 A = C2a * QinL - S2a * UinL |
|
1338 B = S2a * QinL + C2a * UinL |
|
1339 # Stokes Input Vector rotation Eq. E.9 |
|
1340 C = C2ab * QinL - S2ab * UinL |
|
1341 D = S2ab * QinL + C2ab * UinL |
|
1342 |
|
1343 # emitter optics |
|
1344 CosE = np.cos(np.deg2rad(RetE)) |
|
1345 SinE = np.sin(np.deg2rad(RetE)) |
|
1346 ZiE = (1. - DiE ** 2) ** 0.5 |
|
1347 WiE = (1. - ZiE * CosE) |
|
1348 |
|
1349 # Stokes Input Vector after emitter optics equivalent to Eq. E.9 with already rotated input vector from Eq. E.4 |
|
1350 # b = beta |
|
1351 IinE = (IinL + DiE * C) |
|
1352 QinE = (C2b * DiE * IinL + A + S2b * (WiE * D - ZiE * SinE * VinL)) |
|
1353 UinE = (S2b * DiE * IinL + B - C2b * (WiE * D - ZiE * SinE * VinL)) |
|
1354 VinE = (-ZiE * SinE * D + ZiE * CosE * VinL) |
|
1355 |
|
1356 # ------------------------- |
|
1357 # F11 assuemd to be = 1 => measured: F11m = IinP / IinE with atrue |
|
1358 # F11sim = (IinE + DiO*atrue*(C2g*QinE - S2g*UinE))/IinE |
|
1359 # ------------------------- |
|
1360 |
|
1361 for iRotO, iRetO, iDiO, iRotC, iRetC, iDiC, iTP, iTS, iRP, iRS, iERaT, iRotaT, iRetT, iERaR, iRotaR, iRetR \ |
|
1362 in [ |
|
1363 (iRotO, iRetO, iDiO, iRotC, iRetC, iDiC, iTP, iTS, iRP, iRS, iERaT, iRotaT, iRetT, iERaR, iRotaR, iRetR) |
|
1364 for iRotO in range(-nRotO, nRotO + 1) |
|
1365 for iRetO in range(-nRetO, nRetO + 1) |
|
1366 for iDiO in range(-nDiO, nDiO + 1) |
|
1367 for iRotC in range(-nRotC, nRotC + 1) |
|
1368 for iRetC in range(-nRetC, nRetC + 1) |
|
1369 for iDiC in range(-nDiC, nDiC + 1) |
|
1370 for iTP in range(-nTP, nTP + 1) |
|
1371 for iTS in range(-nTS, nTS + 1) |
|
1372 for iRP in range(-nRP, nRP + 1) |
|
1373 for iRS in range(-nRS, nRS + 1) |
|
1374 for iERaT in range(-nERaT, nERaT + 1) |
|
1375 for iRotaT in range(-nRotaT, nRotaT + 1) |
|
1376 for iRetT in range(-nRetT, nRetT + 1) |
|
1377 for iERaR in range(-nERaR, nERaR + 1) |
|
1378 for iRotaR in range(-nRotaR, nRotaR + 1) |
|
1379 for iRetR in range(-nRetR, nRetR + 1)]: |
|
1380 |
|
1381 if nRotO > 0: RotO = RotO0 + iRotO * dRotO / nRotO |
|
1382 if nRetO > 0: RetO = RetO0 + iRetO * dRetO / nRetO |
|
1383 if nDiO > 0: DiO = DiO0 + iDiO * dDiO / nDiO |
|
1384 if nRotC > 0: RotC = RotC0 + iRotC * dRotC / nRotC |
|
1385 if nRetC > 0: RetC = RetC0 + iRetC * dRetC / nRetC |
|
1386 if nDiC > 0: DiC = DiC0 + iDiC * dDiC / nDiC |
|
1387 if nTP > 0: TP = TP0 + iTP * dTP / nTP |
|
1388 if nTS > 0: TS = TS0 + iTS * dTS / nTS |
|
1389 if nRP > 0: RP = RP0 + iRP * dRP / nRP |
|
1390 if nRS > 0: RS = RS0 + iRS * dRS / nRS |
|
1391 if nERaT > 0: ERaT = ERaT0 + iERaT * dERaT / nERaT |
|
1392 if nRotaT > 0: RotaT = RotaT0 + iRotaT * dRotaT / nRotaT |
|
1393 if nRetT > 0: RetT = RetT0 + iRetT * dRetT / nRetT |
|
1394 if nERaR > 0: ERaR = ERaR0 + iERaR * dERaR / nERaR |
|
1395 if nRotaR > 0: RotaR = RotaR0 + iRotaR * dRotaR / nRotaR |
|
1396 if nRetR > 0: RetR = RetR0 + iRetR * dRetR / nRetR |
|
1397 |
|
1398 # print("{0:5.2f}, {1:5.2f}, {2:5.2f}, {3:10d}".format(RotL, RotE, RotO, iN)) |
|
1399 |
|
1400 # receiver optics |
|
1401 CosO = np.cos(np.deg2rad(RetO)) |
|
1402 SinO = np.sin(np.deg2rad(RetO)) |
|
1403 ZiO = (1. - DiO ** 2) ** 0.5 |
|
1404 WiO = (1. - ZiO * CosO) |
|
1405 S2g = np.sin(np.deg2rad(2 * RotO)) |
|
1406 C2g = np.cos(np.deg2rad(2 * RotO)) |
|
1407 # calibrator |
|
1408 CosC = np.cos(np.deg2rad(RetC)) |
|
1409 SinC = np.sin(np.deg2rad(RetC)) |
|
1410 ZiC = (1. - DiC ** 2) ** 0.5 |
|
1411 WiC = (1. - ZiC * CosC) |
|
1412 |
|
1413 # analyser |
|
1414 # For POLLY_XTs |
|
1415 if (RS_RP_depend_on_TS_TP): |
|
1416 RS = 1.0 - TS |
|
1417 RP = 1.0 - TP |
|
1418 TiT = 0.5 * (TP + TS) |
|
1419 DiT = (TP - TS) / (TP + TS) |
|
1420 ZiT = (1. - DiT ** 2.) ** 0.5 |
|
1421 TiR = 0.5 * (RP + RS) |
|
1422 DiR = (RP - RS) / (RP + RS) |
|
1423 ZiR = (1. - DiR ** 2.) ** 0.5 |
|
1424 CosT = np.cos(np.deg2rad(RetT)) |
|
1425 SinT = np.sin(np.deg2rad(RetT)) |
|
1426 CosR = np.cos(np.deg2rad(RetR)) |
|
1427 SinR = np.sin(np.deg2rad(RetR)) |
|
1428 |
|
1429 # cleaning pol-filter |
|
1430 DaT = (1.0 - ERaT) / (1.0 + ERaT) |
|
1431 DaR = (1.0 - ERaR) / (1.0 + ERaR) |
|
1432 TaT = 0.5 * (1.0 + ERaT) |
|
1433 TaR = 0.5 * (1.0 + ERaR) |
|
1434 |
|
1435 S2aT = np.sin(np.deg2rad(h * 2.0 * RotaT)) |
|
1436 C2aT = np.cos(np.deg2rad(2.0 * RotaT)) |
|
1437 S2aR = np.sin(np.deg2rad(h * 2.0 * RotaR)) |
|
1438 C2aR = np.cos(np.deg2rad(2.0 * RotaR)) |
|
1439 |
|
1440 # Analyzer As before the PBS Eq. D.5; combined PBS and cleaning pol-filter |
|
1441 ATPT = (1 + C2aT * DaT * DiT) # unpolarized transmission correction |
|
1442 TTa = TiT * TaT * ATPT # unpolarized transmission |
|
1443 ATP1 = 1.0 |
|
1444 ATP2 = Y * (DiT + C2aT * DaT) / ATPT |
|
1445 ATP3 = Y * S2aT * DaT * ZiT * CosT / ATPT |
|
1446 ATP4 = S2aT * DaT * ZiT * SinT / ATPT |
|
1447 ATP = np.array([ATP1, ATP2, ATP3, ATP4]) |
|
1448 DTa = ATP2 * Y |
|
1449 |
|
1450 ARPT = (1 + C2aR * DaR * DiR) # unpolarized transmission correction |
|
1451 TRa = TiR * TaR * ARPT # unpolarized transmission |
|
1452 ARP1 = 1 |
|
1453 ARP2 = Y * (DiR + C2aR * DaR) / ARPT |
|
1454 ARP3 = Y * S2aR * DaR * ZiR * CosR / ARPT |
|
1455 ARP4 = S2aR * DaR * ZiR * SinR / ARPT |
|
1456 ARP = np.array([ARP1, ARP2, ARP3, ARP4]) |
|
1457 DRa = ARP2 * Y |
|
1458 |
|
1459 # ---- Calculate signals and correction parameters for diffeent locations and calibrators |
|
1460 if LocC == 4: # Calibrator before the PBS |
|
1461 # print("Calibrator location not implemented yet") |
|
1462 |
|
1463 # S2ge = np.sin(np.deg2rad(2*RotO + h*2*RotC)) |
|
1464 # C2ge = np.cos(np.deg2rad(2*RotO + h*2*RotC)) |
|
1465 S2e = np.sin(np.deg2rad(h * 2 * RotC)) |
|
1466 C2e = np.cos(np.deg2rad(2 * RotC)) |
|
1467 # rotated AinP by epsilon Eq. C.3 |
|
1468 ATP2e = C2e * ATP2 + S2e * ATP3 |
|
1469 ATP3e = C2e * ATP3 - S2e * ATP2 |
|
1470 ARP2e = C2e * ARP2 + S2e * ARP3 |
|
1471 ARP3e = C2e * ARP3 - S2e * ARP2 |
|
1472 ATPe = np.array([ATP1, ATP2e, ATP3e, ATP4]) |
|
1473 ARPe = np.array([ARP1, ARP2e, ARP3e, ARP4]) |
|
1474 # Stokes Input Vector before the polarising beam splitter Eq. E.31 |
|
1475 A = C2g * QinE - S2g * UinE |
|
1476 B = S2g * QinE + C2g * UinE |
|
1477 # C = (WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE) |
|
1478 Co = ZiO * SinO * VinE |
|
1479 Ca = (WiO * B - 2 * ZiO * SinO * VinE) |
|
1480 # C = Co + aCal*Ca |
|
1481 # IinP = (IinE + DiO*aCal*A) |
|
1482 # QinP = (C2g*DiO*IinE + aCal*QinE - S2g*C) |
|
1483 # UinP = (S2g*DiO*IinE - aCal*UinE + C2g*C) |
|
1484 # VinP = (ZiO*SinO*aCal*B + ZiO*CosO*(1-2*aCal)*VinE) |
|
1485 IinPo = IinE |
|
1486 QinPo = (C2g * DiO * IinE - S2g * Co) |
|
1487 UinPo = (S2g * DiO * IinE + C2g * Co) |
|
1488 VinPo = ZiO * CosO * VinE |
|
1489 |
|
1490 IinPa = DiO * A |
|
1491 QinPa = QinE - S2g * Ca |
|
1492 UinPa = -UinE + C2g * Ca |
|
1493 VinPa = ZiO * (SinO * B - 2 * CosO * VinE) |
|
1494 |
|
1495 IinP = IinPo + aCal * IinPa |
|
1496 QinP = QinPo + aCal * QinPa |
|
1497 UinP = UinPo + aCal * UinPa |
|
1498 VinP = VinPo + aCal * VinPa |
|
1499 # Stokes Input Vector before the polarising beam splitter rotated by epsilon Eq. C.3 |
|
1500 # QinPe = C2e*QinP + S2e*UinP |
|
1501 # UinPe = C2e*UinP - S2e*QinP |
|
1502 QinPoe = C2e * QinPo + S2e * UinPo |
|
1503 UinPoe = C2e * UinPo - S2e * QinPo |
|
1504 QinPae = C2e * QinPa + S2e * UinPa |
|
1505 UinPae = C2e * UinPa - S2e * QinPa |
|
1506 QinPe = C2e * QinP + S2e * UinP |
|
1507 UinPe = C2e * UinP - S2e * QinP |
|
1508 |
|
1509 # Calibration signals and Calibration correction K from measurements with LDRCal / aCal |
|
1510 if (TypeC == 2) or (TypeC == 1): # rotator calibration Eq. C.4 |
|
1511 # parameters for calibration with aCal |
|
1512 AT = ATP1 * IinP + h * ATP4 * VinP |
|
1513 BT = ATP3e * QinP - h * ATP2e * UinP |
|
1514 AR = ARP1 * IinP + h * ARP4 * VinP |
|
1515 BR = ARP3e * QinP - h * ARP2e * UinP |
|
1516 # Correction parameters for normal measurements; they are independent of LDR |
|
1517 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out |
|
1518 IS1 = np.array([IinPo, QinPo, UinPo, VinPo]) |
|
1519 IS2 = np.array([IinPa, QinPa, UinPa, VinPa]) |
|
1520 GT = np.dot(ATP, IS1) |
|
1521 GR = np.dot(ARP, IS1) |
|
1522 HT = np.dot(ATP, IS2) |
|
1523 HR = np.dot(ARP, IS2) |
|
1524 else: |
|
1525 IS1 = np.array([IinPo, QinPo, UinPo, VinPo]) |
|
1526 IS2 = np.array([IinPa, QinPa, UinPa, VinPa]) |
|
1527 GT = np.dot(ATPe, IS1) |
|
1528 GR = np.dot(ARPe, IS1) |
|
1529 HT = np.dot(ATPe, IS2) |
|
1530 HR = np.dot(ARPe, IS2) |
|
1531 elif (TypeC == 3) or (TypeC == 4): # linear polariser calibration Eq. C.5 |
|
1532 # parameters for calibration with aCal |
|
1533 AT = ATP1 * IinP + ATP3e * UinPe + ZiC * CosC * (ATP2e * QinPe + ATP4 * VinP) |
|
1534 BT = DiC * (ATP1 * UinPe + ATP3e * IinP) - ZiC * SinC * (ATP2e * VinP - ATP4 * QinPe) |
|
1535 AR = ARP1 * IinP + ARP3e * UinPe + ZiC * CosC * (ARP2e * QinPe + ARP4 * VinP) |
|
1536 BR = DiC * (ARP1 * UinPe + ARP3e * IinP) - ZiC * SinC * (ARP2e * VinP - ARP4 * QinPe) |
|
1537 # Correction parameters for normal measurements; they are independent of LDR |
|
1538 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out |
|
1539 IS1 = np.array([IinPo, QinPo, UinPo, VinPo]) |
|
1540 IS2 = np.array([IinPa, QinPa, UinPa, VinPa]) |
|
1541 GT = np.dot(ATP, IS1) |
|
1542 GR = np.dot(ARP, IS1) |
|
1543 HT = np.dot(ATP, IS2) |
|
1544 HR = np.dot(ARP, IS2) |
|
1545 else: |
|
1546 IS1e = np.array( |
|
1547 [IinPo + DiC * QinPoe, DiC * IinPo + QinPoe, ZiC * (CosC * UinPoe + SinC * VinPo), |
|
1548 -ZiC * (SinC * UinPoe - CosC * VinPo)]) |
|
1549 IS2e = np.array( |
|
1550 [IinPa + DiC * QinPae, DiC * IinPa + QinPae, ZiC * (CosC * UinPae + SinC * VinPa), |
|
1551 -ZiC * (SinC * UinPae - CosC * VinPa)]) |
|
1552 GT = np.dot(ATPe, IS1e) |
|
1553 GR = np.dot(ARPe, IS1e) |
|
1554 HT = np.dot(ATPe, IS2e) |
|
1555 HR = np.dot(ARPe, IS2e) |
|
1556 elif (TypeC == 6): # diattenuator calibration +-22.5° rotated_diattenuator_X22x5deg.odt |
|
1557 # parameters for calibration with aCal |
|
1558 AT = ATP1 * IinP + sqr05 * DiC * (ATP1 * QinPe + ATP2e * IinP) + (1 - 0.5 * WiC) * ( |
|
1559 ATP2e * QinPe + ATP3e * UinPe) + ZiC * ( |
|
1560 sqr05 * SinC * (ATP3e * VinP - ATP4 * UinPe) + ATP4 * CosC * VinP) |
|
1561 BT = sqr05 * DiC * (ATP1 * UinPe + ATP3e * IinP) + 0.5 * WiC * ( |
|
1562 ATP2e * UinPe + ATP3e * QinPe) - sqr05 * ZiC * SinC * (ATP2e * VinP - ATP4 * QinPe) |
|
1563 AR = ARP1 * IinP + sqr05 * DiC * (ARP1 * QinPe + ARP2e * IinP) + (1 - 0.5 * WiC) * ( |
|
1564 ARP2e * QinPe + ARP3e * UinPe) + ZiC * ( |
|
1565 sqr05 * SinC * (ARP3e * VinP - ARP4 * UinPe) + ARP4 * CosC * VinP) |
|
1566 BR = sqr05 * DiC * (ARP1 * UinPe + ARP3e * IinP) + 0.5 * WiC * ( |
|
1567 ARP2e * UinPe + ARP3e * QinPe) - sqr05 * ZiC * SinC * (ARP2e * VinP - ARP4 * QinPe) |
|
1568 # Correction parameters for normal measurements; they are independent of LDR |
|
1569 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out |
|
1570 IS1 = np.array([IinPo, QinPo, UinPo, VinPo]) |
|
1571 IS2 = np.array([IinPa, QinPa, UinPa, VinPa]) |
|
1572 GT = np.dot(ATP, IS1) |
|
1573 GR = np.dot(ARP, IS1) |
|
1574 HT = np.dot(ATP, IS2) |
|
1575 HR = np.dot(ARP, IS2) |
|
1576 else: |
|
1577 IS1e = np.array( |
|
1578 [IinPo + DiC * QinPoe, DiC * IinPo + QinPoe, ZiC * (CosC * UinPoe + SinC * VinPo), |
|
1579 -ZiC * (SinC * UinPoe - CosC * VinPo)]) |
|
1580 IS2e = np.array( |
|
1581 [IinPa + DiC * QinPae, DiC * IinPa + QinPae, ZiC * (CosC * UinPae + SinC * VinPa), |
|
1582 -ZiC * (SinC * UinPae - CosC * VinPa)]) |
|
1583 GT = np.dot(ATPe, IS1e) |
|
1584 GR = np.dot(ARPe, IS1e) |
|
1585 HT = np.dot(ATPe, IS2e) |
|
1586 HR = np.dot(ARPe, IS2e) |
|
1587 else: |
|
1588 print("Calibrator not implemented yet") |
|
1589 sys.exit() |
|
1590 |
|
1591 elif LocC == 3: # C before receiver optics Eq.57 |
|
1592 |
|
1593 # S2ge = np.sin(np.deg2rad(2*RotO - 2*RotC)) |
|
1594 # C2ge = np.cos(np.deg2rad(2*RotO - 2*RotC)) |
|
1595 S2e = np.sin(np.deg2rad(2 * RotC)) |
|
1596 C2e = np.cos(np.deg2rad(2 * RotC)) |
|
1597 |
|
1598 # AS with C before the receiver optics (see document rotated_diattenuator_X22x5deg.odt) |
|
1599 AF1 = np.array([1, C2g * DiO, S2g * DiO, 0]) |
|
1600 AF2 = np.array([C2g * DiO, 1 - S2g ** 2 * WiO, S2g * C2g * WiO, -S2g * ZiO * SinO]) |
|
1601 AF3 = np.array([S2g * DiO, S2g * C2g * WiO, 1 - C2g ** 2 * WiO, C2g * ZiO * SinO]) |
|
1602 AF4 = np.array([0, S2g * SinO, -C2g * SinO, CosO]) |
|
1603 |
|
1604 ATF = (ATP1 * AF1 + ATP2 * AF2 + ATP3 * AF3 + ATP4 * AF4) |
|
1605 ARF = (ARP1 * AF1 + ARP2 * AF2 + ARP3 * AF3 + ARP4 * AF4) |
|
1606 ATF1 = ATF[0] |
|
1607 ATF2 = ATF[1] |
|
1608 ATF3 = ATF[2] |
|
1609 ATF4 = ATF[3] |
|
1610 ARF1 = ARF[0] |
|
1611 ARF2 = ARF[1] |
|
1612 ARF3 = ARF[2] |
|
1613 ARF4 = ARF[3] |
|
1614 |
|
1615 # rotated AinF by epsilon |
|
1616 ATF2e = C2e * ATF[1] + S2e * ATF[2] |
|
1617 ATF3e = C2e * ATF[2] - S2e * ATF[1] |
|
1618 ARF2e = C2e * ARF[1] + S2e * ARF[2] |
|
1619 ARF3e = C2e * ARF[2] - S2e * ARF[1] |
|
1620 |
|
1621 ATFe = np.array([ATF1, ATF2e, ATF3e, ATF4]) |
|
1622 ARFe = np.array([ARF1, ARF2e, ARF3e, ARF4]) |
|
1623 |
|
1624 QinEe = C2e * QinE + S2e * UinE |
|
1625 UinEe = C2e * UinE - S2e * QinE |
|
1626 |
|
1627 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) |
|
1628 IinF = IinE |
|
1629 QinF = aCal * QinE |
|
1630 UinF = -aCal * UinE |
|
1631 VinF = (1. - 2. * aCal) * VinE |
|
1632 |
|
1633 IinFo = IinE |
|
1634 QinFo = 0. |
|
1635 UinFo = 0. |
|
1636 VinFo = VinE |
|
1637 |
|
1638 IinFa = 0. |
|
1639 QinFa = QinE |
|
1640 UinFa = -UinE |
|
1641 VinFa = -2. * VinE |
|
1642 |
|
1643 # Stokes Input Vector before receiver optics rotated by epsilon Eq. C.3 |
|
1644 QinFe = C2e * QinF + S2e * UinF |
|
1645 UinFe = C2e * UinF - S2e * QinF |
|
1646 QinFoe = C2e * QinFo + S2e * UinFo |
|
1647 UinFoe = C2e * UinFo - S2e * QinFo |
|
1648 QinFae = C2e * QinFa + S2e * UinFa |
|
1649 UinFae = C2e * UinFa - S2e * QinFa |
|
1650 |
|
1651 # Calibration signals and Calibration correction K from measurements with LDRCal / aCal |
|
1652 if (TypeC == 2) or (TypeC == 1): # rotator calibration Eq. C.4 |
|
1653 AT = ATF1 * IinF + ATF4 * h * VinF |
|
1654 BT = ATF3e * QinF - ATF2e * h * UinF |
|
1655 AR = ARF1 * IinF + ARF4 * h * VinF |
|
1656 BR = ARF3e * QinF - ARF2e * h * UinF |
|
1657 |
|
1658 # Correction parameters for normal measurements; they are independent of LDR |
|
1659 if (not RotationErrorEpsilonForNormalMeasurements): |
|
1660 GT = ATF1 * IinE + ATF4 * VinE |
|
1661 GR = ARF1 * IinE + ARF4 * VinE |
|
1662 HT = ATF2 * QinE - ATF3 * UinE - ATF4 * 2 * VinE |
|
1663 HR = ARF2 * QinE - ARF3 * UinE - ARF4 * 2 * VinE |
|
1664 else: |
|
1665 GT = ATF1 * IinE + ATF4 * h * VinE |
|
1666 GR = ARF1 * IinE + ARF4 * h * VinE |
|
1667 HT = ATF2e * QinE - ATF3e * h * UinE - ATF4 * h * 2 * VinE |
|
1668 HR = ARF2e * QinE - ARF3e * h * UinE - ARF4 * h * 2 * VinE |
|
1669 |
|
1670 elif (TypeC == 3) or (TypeC == 4): # linear polariser calibration Eq. C.5 |
|
1671 # p = +45°, m = -45° |
|
1672 IF1e = np.array([IinF, ZiC * CosC * QinFe, UinFe, ZiC * CosC * VinF]) |
|
1673 IF2e = np.array([DiC * UinFe, -ZiC * SinC * VinF, DiC * IinF, ZiC * SinC * QinFe]) |
|
1674 |
|
1675 AT = np.dot(ATFe, IF1e) |
|
1676 AR = np.dot(ARFe, IF1e) |
|
1677 BT = np.dot(ATFe, IF2e) |
|
1678 BR = np.dot(ARFe, IF2e) |
|
1679 |
|
1680 # Correction parameters for normal measurements; they are independent of LDR --- the same as for TypeC = 6 |
|
1681 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out |
|
1682 IS1 = np.array([IinE, 0, 0, VinE]) |
|
1683 IS2 = np.array([0, QinE, -UinE, -2 * VinE]) |
|
1684 |
|
1685 GT = np.dot(ATF, IS1) |
|
1686 GR = np.dot(ARF, IS1) |
|
1687 HT = np.dot(ATF, IS2) |
|
1688 HR = np.dot(ARF, IS2) |
|
1689 else: |
|
1690 IS1e = np.array( |
|
1691 [IinFo + DiC * QinFoe, DiC * IinFo + QinFoe, ZiC * (CosC * UinFoe + SinC * VinFo), |
|
1692 -ZiC * (SinC * UinFoe - CosC * VinFo)]) |
|
1693 IS2e = np.array( |
|
1694 [IinFa + DiC * QinFae, DiC * IinFa + QinFae, ZiC * (CosC * UinFae + SinC * VinFa), |
|
1695 -ZiC * (SinC * UinFae - CosC * VinFa)]) |
|
1696 GT = np.dot(ATFe, IS1e) |
|
1697 GR = np.dot(ARFe, IS1e) |
|
1698 HT = np.dot(ATFe, IS2e) |
|
1699 HR = np.dot(ARFe, IS2e) |
|
1700 |
|
1701 elif (TypeC == 6): # diattenuator calibration +-22.5° rotated_diattenuator_X22x5deg.odt |
|
1702 # p = +22.5°, m = -22.5° |
|
1703 IF1e = np.array([IinF + sqr05 * DiC * QinFe, sqr05 * DiC * IinF + (1 - 0.5 * WiC) * QinFe, |
|
1704 (1 - 0.5 * WiC) * UinFe + sqr05 * ZiC * SinC * VinF, |
|
1705 -sqr05 * ZiC * SinC * UinFe + ZiC * CosC * VinF]) |
|
1706 IF2e = np.array([sqr05 * DiC * UinFe, 0.5 * WiC * UinFe - sqr05 * ZiC * SinC * VinF, |
|
1707 sqr05 * DiC * IinF + 0.5 * WiC * QinFe, sqr05 * ZiC * SinC * QinFe]) |
|
1708 |
|
1709 AT = np.dot(ATFe, IF1e) |
|
1710 AR = np.dot(ARFe, IF1e) |
|
1711 BT = np.dot(ATFe, IF2e) |
|
1712 BR = np.dot(ARFe, IF2e) |
|
1713 |
|
1714 # Correction parameters for normal measurements; they are independent of LDR |
|
1715 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out |
|
1716 # IS1 = np.array([IinE,0,0,VinE]) |
|
1717 # IS2 = np.array([0,QinE,-UinE,-2*VinE]) |
|
1718 IS1 = np.array([IinFo, 0, 0, VinFo]) |
|
1719 IS2 = np.array([0, QinFa, UinFa, VinFa]) |
|
1720 GT = np.dot(ATF, IS1) |
|
1721 GR = np.dot(ARF, IS1) |
|
1722 HT = np.dot(ATF, IS2) |
|
1723 HR = np.dot(ARF, IS2) |
|
1724 else: |
|
1725 # IS1e = np.array([IinE,DiC*IinE,ZiC*SinC*VinE,ZiC*CosC*VinE]) |
|
1726 # IS2e = np.array([DiC*QinEe,QinEe,-ZiC*(CosC*UinEe+2*SinC*VinE),ZiC*(SinC*UinEe-2*CosC*VinE)]) |
|
1727 IS1e = np.array( |
|
1728 [IinFo + DiC * QinFoe, DiC * IinFo + QinFoe, ZiC * (CosC * UinFoe + SinC * VinFo), |
|
1729 -ZiC * (SinC * UinFoe - CosC * VinFo)]) |
|
1730 IS2e = np.array( |
|
1731 [IinFa + DiC * QinFae, DiC * IinFa + QinFae, ZiC * (CosC * UinFae + SinC * VinFa), |
|
1732 -ZiC * (SinC * UinFae - CosC * VinFa)]) |
|
1733 GT = np.dot(ATFe, IS1e) |
|
1734 GR = np.dot(ARFe, IS1e) |
|
1735 HT = np.dot(ATFe, IS2e) |
|
1736 HR = np.dot(ARFe, IS2e) |
|
1737 |
|
1738 |
|
1739 else: |
|
1740 print('Calibrator not implemented yet') |
|
1741 sys.exit() |
|
1742 |
|
1743 elif LocC == 2: # C behind emitter optics Eq.57 |
|
1744 # print("Calibrator location not implemented yet") |
|
1745 S2e = np.sin(np.deg2rad(2 * RotC)) |
|
1746 C2e = np.cos(np.deg2rad(2 * RotC)) |
|
1747 |
|
1748 # AS with C before the receiver optics (see document rotated_diattenuator_X22x5deg.odt) |
|
1749 AF1 = np.array([1, C2g * DiO, S2g * DiO, 0]) |
|
1750 AF2 = np.array([C2g * DiO, 1 - S2g ** 2 * WiO, S2g * C2g * WiO, -S2g * ZiO * SinO]) |
|
1751 AF3 = np.array([S2g * DiO, S2g * C2g * WiO, 1 - C2g ** 2 * WiO, C2g * ZiO * SinO]) |
|
1752 AF4 = np.array([0, S2g * SinO, -C2g * SinO, CosO]) |
|
1753 |
|
1754 ATF = (ATP1 * AF1 + ATP2 * AF2 + ATP3 * AF3 + ATP4 * AF4) |
|
1755 ARF = (ARP1 * AF1 + ARP2 * AF2 + ARP3 * AF3 + ARP4 * AF4) |
|
1756 ATF1 = ATF[0] |
|
1757 ATF2 = ATF[1] |
|
1758 ATF3 = ATF[2] |
|
1759 ATF4 = ATF[3] |
|
1760 ARF1 = ARF[0] |
|
1761 ARF2 = ARF[1] |
|
1762 ARF3 = ARF[2] |
|
1763 ARF4 = ARF[3] |
|
1764 |
|
1765 # AS with C behind the emitter -------------------------------------------- |
|
1766 # terms without aCal |
|
1767 ATE1o, ARE1o = ATF1, ARF1 |
|
1768 ATE2o, ARE2o = 0., 0. |
|
1769 ATE3o, ARE3o = 0., 0. |
|
1770 ATE4o, ARE4o = ATF4, ARF4 |
|
1771 # terms with aCal |
|
1772 ATE1a, ARE1a = 0., 0. |
|
1773 ATE2a, ARE2a = ATF2, ARF2 |
|
1774 ATE3a, ARE3a = -ATF3, -ARF3 |
|
1775 ATE4a, ARE4a = -2 * ATF4, -2 * ARF4 |
|
1776 # rotated AinEa by epsilon |
|
1777 ATE2ae = C2e * ATF2 + S2e * ATF3 |
|
1778 ATE3ae = -S2e * ATF2 - C2e * ATF3 |
|
1779 ARE2ae = C2e * ARF2 + S2e * ARF3 |
|
1780 ARE3ae = -S2e * ARF2 - C2e * ARF3 |
|
1781 |
|
1782 ATE1 = ATE1o |
|
1783 ATE2e = aCal * ATE2ae |
|
1784 ATE3e = aCal * ATE3ae |
|
1785 ATE4 = (1 - 2 * aCal) * ATF4 |
|
1786 ARE1 = ARE1o |
|
1787 ARE2e = aCal * ARE2ae |
|
1788 ARE3e = aCal * ARE3ae |
|
1789 ARE4 = (1. - 2. * aCal) * ARF4 |
|
1790 |
|
1791 # rotated IinE |
|
1792 QinEe = C2e * QinE + S2e * UinE |
|
1793 UinEe = C2e * UinE - S2e * QinE |
|
1794 |
|
1795 # --- Calibration signals and Calibration correction K from measurements with LDRCal / aCal |
|
1796 if (TypeC == 2) or (TypeC == 1): # +++++++++ rotator calibration Eq. C.4 |
|
1797 AT = ATE1o * IinE + (ATE4o + aCal * ATE4a) * h * VinE |
|
1798 BT = aCal * (ATE3ae * QinEe - ATE2ae * h * UinEe) |
|
1799 AR = ARE1o * IinE + (ARE4o + aCal * ARE4a) * h * VinE |
|
1800 BR = aCal * (ARE3ae * QinEe - ARE2ae * h * UinEe) |
|
1801 |
|
1802 # Correction parameters for normal measurements; they are independent of LDR |
|
1803 if (not RotationErrorEpsilonForNormalMeasurements): |
|
1804 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) |
|
1805 GT = ATE1o * IinE + ATE4o * h * VinE |
|
1806 GR = ARE1o * IinE + ARE4o * h * VinE |
|
1807 HT = ATE2a * QinE + ATE3a * h * UinEe + ATE4a * h * VinE |
|
1808 HR = ARE2a * QinE + ARE3a * h * UinEe + ARE4a * h * VinE |
|
1809 else: |
|
1810 GT = ATE1o * IinE + ATE4o * h * VinE |
|
1811 GR = ARE1o * IinE + ARE4o * h * VinE |
|
1812 HT = ATE2ae * QinE + ATE3ae * h * UinEe + ATE4a * h * VinE |
|
1813 HR = ARE2ae * QinE + ARE3ae * h * UinEe + ARE4a * h * VinE |
|
1814 |
|
1815 elif (TypeC == 3) or (TypeC == 4): # +++++++++ linear polariser calibration Eq. C.5 |
|
1816 # p = +45°, m = -45° |
|
1817 AT = ATE1 * IinE + ZiC * CosC * (ATE2e * QinEe + ATE4 * VinE) + ATE3e * UinEe |
|
1818 BT = DiC * (ATE1 * UinEe + ATE3e * IinE) + ZiC * SinC * (ATE4 * QinEe - ATE2e * VinE) |
|
1819 AR = ARE1 * IinE + ZiC * CosC * (ARE2e * QinEe + ARE4 * VinE) + ARE3e * UinEe |
|
1820 BR = DiC * (ARE1 * UinEe + ARE3e * IinE) + ZiC * SinC * (ARE4 * QinEe - ARE2e * VinE) |
|
1821 |
|
1822 # Correction parameters for normal measurements; they are independent of LDR |
|
1823 if (not RotationErrorEpsilonForNormalMeasurements): |
|
1824 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) |
|
1825 GT = ATE1o * IinE + ATE4o * VinE |
|
1826 GR = ARE1o * IinE + ARE4o * VinE |
|
1827 HT = ATE2a * QinE + ATE3a * UinE + ATE4a * VinE |
|
1828 HR = ARE2a * QinE + ARE3a * UinE + ARE4a * VinE |
|
1829 else: |
|
1830 D = IinE + DiC * QinEe |
|
1831 A = DiC * IinE + QinEe |
|
1832 B = ZiC * (CosC * UinEe + SinC * VinE) |
|
1833 C = -ZiC * (SinC * UinEe - CosC * VinE) |
|
1834 GT = ATE1o * D + ATE4o * C |
|
1835 GR = ARE1o * D + ARE4o * C |
|
1836 HT = ATE2a * A + ATE3a * B + ATE4a * C |
|
1837 HR = ARE2a * A + ARE3a * B + ARE4a * C |
|
1838 |
|
1839 elif (TypeC == 6): # real HWP calibration +-22.5° rotated_diattenuator_X22x5deg.odt |
|
1840 # p = +22.5°, m = -22.5° |
|
1841 IE1e = np.array([IinE + sqr05 * DiC * QinEe, sqr05 * DiC * IinE + (1 - 0.5 * WiC) * QinEe, |
|
1842 (1. - 0.5 * WiC) * UinEe + sqr05 * ZiC * SinC * VinE, |
|
1843 -sqr05 * ZiC * SinC * UinEe + ZiC * CosC * VinE]) |
|
1844 IE2e = np.array([sqr05 * DiC * UinEe, 0.5 * WiC * UinEe - sqr05 * ZiC * SinC * VinE, |
|
1845 sqr05 * DiC * IinE + 0.5 * WiC * QinEe, sqr05 * ZiC * SinC * QinEe]) |
|
1846 ATEe = np.array([ATE1, ATE2e, ATE3e, ATE4]) |
|
1847 AREe = np.array([ARE1, ARE2e, ARE3e, ARE4]) |
|
1848 AT = np.dot(ATEe, IE1e) |
|
1849 AR = np.dot(AREe, IE1e) |
|
1850 BT = np.dot(ATEe, IE2e) |
|
1851 BR = np.dot(AREe, IE2e) |
|
1852 |
|
1853 # Correction parameters for normal measurements; they are independent of LDR |
|
1854 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out |
|
1855 GT = ATE1o * IinE + ATE4o * VinE |
|
1856 GR = ARE1o * IinE + ARE4o * VinE |
|
1857 HT = ATE2a * QinE + ATE3a * UinE + ATE4a * VinE |
|
1858 HR = ARE2a * QinE + ARE3a * UinE + ARE4a * VinE |
|
1859 else: |
|
1860 D = IinE + DiC * QinEe |
|
1861 A = DiC * IinE + QinEe |
|
1862 B = ZiC * (CosC * UinEe + SinC * VinE) |
|
1863 C = -ZiC * (SinC * UinEe - CosC * VinE) |
|
1864 GT = ATE1o * D + ATE4o * C |
|
1865 GR = ARE1o * D + ARE4o * C |
|
1866 HT = ATE2a * A + ATE3a * B + ATE4a * C |
|
1867 HR = ARE2a * A + ARE3a * B + ARE4a * C |
|
1868 else: |
|
1869 print('Calibrator not implemented yet') |
|
1870 sys.exit() |
|
1871 |
|
1872 for iTCalT, iTCalR, iNCalTp, iNCalTm, iNCalRp, iNCalRm, iNIt, iNIr \ |
|
1873 in [ |
|
1874 (iTCalT, iTCalR, iNCalTp, iNCalTm, iNCalRp, iNCalRm, iNIt, iNIr) |
|
1875 for iTCalT in range(-nTCalT, nTCalT + 1) # Etax |
|
1876 for iTCalR in range(-nTCalR, nTCalR + 1) # Etax |
|
1877 for iNCalTp in range(-nNCal, nNCal + 1) # noise error of calibration signals => Etax |
|
1878 for iNCalTm in range(-nNCal, nNCal + 1) # noise error of calibration signals => Etax |
|
1879 for iNCalRp in range(-nNCal, nNCal + 1) # noise error of calibration signals => Etax |
|
1880 for iNCalRm in range(-nNCal, nNCal + 1) # noise error of calibration signals => Etax |
|
1881 for iNIt in range(-nNI, nNI + 1) |
|
1882 for iNIr in range(-nNI, nNI + 1)]: |
|
1883 |
|
1884 # Calibration signals with aCal => Determination of the correction K of the real calibration factor |
|
1885 IoutTp = TTa * TiC * TiO * TiE * (AT + BT) |
|
1886 IoutTm = TTa * TiC * TiO * TiE * (AT - BT) |
|
1887 IoutRp = TRa * TiC * TiO * TiE * (AR + BR) |
|
1888 IoutRm = TRa * TiC * TiO * TiE * (AR - BR) |
|
1889 |
|
1890 if nTCalT > 0: TCalT = TCalT0 + iTCalT * dTCalT / nTCalT |
|
1891 if nTCalR > 0: TCalR = TCalR0 + iTCalR * dTCalR / nTCalR |
|
1892 # signal noise errors |
|
1893 # ----- random error calculation ---------- |
|
1894 # noise must be calculated from/with the actually measured signals; influence of TCalT, TCalR errors on noise are not considered ? |
|
1895 # actually measured signal counts are in input file and don't change |
|
1896 # relative standard deviation of calibration signals with LDRcal; assumed to be statisitcally independent |
|
1897 # error nNCal: one-sigma in steps to left and right for calibration signals |
|
1898 if nNCal > 0: |
|
1899 if (CalcFrom0deg): |
|
1900 dIoutTp = (NCalT * IoutTp) ** -0.5 |
|
1901 dIoutTm = (NCalT * IoutTm) ** -0.5 |
|
1902 dIoutRp = (NCalR * IoutRp) ** -0.5 |
|
1903 dIoutRm = (NCalR * IoutRm) ** -0.5 |
|
1904 else: |
|
1905 dIoutTp = dIoutTp0 * (IoutTp / IoutTp0) |
|
1906 dIoutTm = dIoutTm0 * (IoutTm / IoutTm0) |
|
1907 dIoutRp = dIoutRp0 * (IoutRp / IoutRp0) |
|
1908 dIoutRm = dIoutRm0 * (IoutRm / IoutRm0) |
|
1909 # print(iTCalT, iTCalR, iNCalTp, iNCalTm, iNCalRp, iNCalRm, iNIt, iNIr, IoutTp, dIoutTp) |
|
1910 IoutTp = IoutTp * (1. + iNCalTp * dIoutTp / nNCal) |
|
1911 IoutTm = IoutTm * (1. + iNCalTm * dIoutTm / nNCal) |
|
1912 IoutRp = IoutRp * (1. + iNCalRp * dIoutRp / nNCal) |
|
1913 IoutRm = IoutRm * (1. + iNCalRm * dIoutRm / nNCal) |
|
1914 |
|
1915 IoutTp = IoutTp * TCalT / TCalT0 |
|
1916 IoutTm = IoutTm * TCalT / TCalT0 |
|
1917 IoutRp = IoutRp * TCalR / TCalR0 |
|
1918 IoutRm = IoutRm * TCalR / TCalR0 |
|
1919 # --- Results and Corrections; electronic etaR and etaT are assumed to be 1 for true and assumed true systems |
|
1920 # calibration factor |
|
1921 Eta = (TRa / TTa) # = TRa / TTa; Eta = Eta*/K Eq. 84; corrected according to the papers supplement Eqs. (S.10.10.1) ff |
|
1922 # possibly real calibration factor |
|
1923 Etapx = IoutRp / IoutTp |
|
1924 Etamx = IoutRm / IoutTm |
|
1925 Etax = (Etapx * Etamx) ** 0.5 |
|
1926 K = Etax / Eta |
|
1927 # 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)) |
|
1928 # print("{0:6.3f},{1:6.3f},{2:6.3f},{3:6.3f}".format(DiC, ZiC, Kp, Km)) |
|
1929 |
|
1930 # For comparison with Volkers Libreoffice Müller Matrix spreadsheet |
|
1931 # Eta_test_p = (IoutRp/IoutTp) |
|
1932 # Eta_test_m = (IoutRm/IoutTm) |
|
1933 # Eta_test = (Eta_test_p*Eta_test_m)**0.5 |
|
1934 ''' |
|
1935 for iIt, iIr \ |
|
1936 in [(iIt, iIr) |
|
1937 for iIt in range(-nNI, nNI + 1) |
|
1938 for iIr in range(-nNI, nNI + 1)]: |
|
1939 ''' |
|
1940 |
|
1941 iN = iN + 1 |
|
1942 if (iN == 10001): |
|
1943 ctime = clock() |
|
1944 print(" estimated time ", "{0:4.2f}".format(N / 10000 * (ctime - atime)), "sec ") # , end="") |
|
1945 print("\r elapsed time ", "{0:5.0f}".format((ctime - atime)), "sec ", end="\r") |
|
1946 ctime = clock() |
|
1947 if ((ctime - dtime) > 10): |
|
1948 print("\r elapsed time ", "{0:5.0f}".format((ctime - atime)), "sec ", end="\r") |
|
1949 dtime = ctime |
|
1950 |
|
1951 # *** loop for different real LDRs ********************************************************************** |
|
1952 iLDR = -1 |
|
1953 for LDRTrue in LDRrange: |
|
1954 iLDR = iLDR + 1 |
|
1955 atrue = (1. - LDRTrue) / (1. + LDRTrue) |
|
1956 # ----- Forward simulated signals and LDRsim with atrue; from input file; not considering TiC. |
|
1957 It = TTa * TiO * TiE * (GT + atrue * HT) # TaT*TiT*TiC*TiO*IinL*(GT+atrue*HT) |
|
1958 Ir = TRa * TiO * TiE * (GR + atrue * HR) # TaR*TiR*TiC*TiO*IinL*(GR+atrue*HR) |
|
1959 # # signal noise errors; standard deviation of signals; assumed to be statisitcally independent |
|
1960 # because the signals depend on LDRtrue, the errors dIt and dIr must be calculated for each LDRtrue |
|
1961 if (CalcFrom0deg): |
|
1962 ''' |
|
1963 dIt = ((NCalT * It / IoutTp * NILfac / TCalT) ** -0.5) |
|
1964 dIr = ((NCalR * Ir / IoutRp * NILfac / TCalR) ** -0.5) |
|
1965 ''' |
|
1966 dIt = ((It * NI * eFacT) ** -0.5) |
|
1967 dIr = ((Ir * NI * eFacR) ** -0.5) |
|
1968 else: |
|
1969 dIt = ((It * NI * eFacT) ** -0.5) |
|
1970 dIr = ((Ir * NI * eFacR) ** -0.5) |
|
1971 ''' |
|
1972 # does this work? Why not as above? |
|
1973 dIt = ((NCalT * 2. * NILfac / TCalT ) ** -0.5) |
|
1974 dIr = ((NCalR * 2. * NILfac / TCalR) ** -0.5) |
|
1975 ''' |
|
1976 # error nNI: one-sigma in steps to left and right for 0° signals |
|
1977 if nNI > 0: |
|
1978 It = It * (1. + iNIt * dIt / nNI) |
|
1979 Ir = Ir * (1. + iNIr * dIr / nNI) |
|
1980 |
|
1981 # LDRsim = 1/Eta*Ir/It # simulated LDR* with Y from input file |
|
1982 LDRsim = Ir / It # simulated uncorrected LDR with Y from input file |
|
1983 |
|
1984 # ----- Backward correction |
|
1985 # Corrected LDRCorr with assumed true G0,H0,K0,Eta0 from forward simulated (real) LDRsim(atrue) |
|
1986 LDRCorr = (LDRsim / (Etax / K0) * (GT0 + HT0) - (GR0 + HR0)) / ((GR0 - HR0) - LDRsim / (Etax / K0) * (GT0 - HT0)) |
|
1987 |
|
1988 # The following is a test whether the equations for calibration Etax and normal signal (GHK, LDRsim) are consistent |
|
1989 # LDRCorr = (LDRsim / Eta * (GT + HT) - (GR + HR)) / ((GR - HR) - LDRsim / Eta * (GT - HT)) |
|
1990 # Without any correction |
|
1991 LDRunCorr = LDRsim / Etax |
|
1992 # 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))) |
|
1993 |
|
1994 |
|
1995 ''' |
|
1996 # -- F11corr from It and Ir and calibration EtaX |
|
1997 Text1 = "!!! EXPERIMENTAL !!! F11corr from It and Ir with calibration EtaX: x-axis: F11corr(LDRtrue) / F11corr(LDRtrue = 0.004) - 1" |
|
1998 F11corr = 1 / (TiO * TiE) * ( |
|
1999 (HR0 * Etax / K0 * It / TTa - HT0 * Ir / TRa) / (HR0 * GT0 - HT0 * GR0)) # IL = 1 Eq.(64); Etax/K0 = Eta0. |
|
2000 ''' |
|
2001 # Corrected F11corr with assumed true G0,H0,K0 from forward simulated (real) It and Ir (atrue) |
|
2002 Text1 = "!!! EXPERIMENTAL !!! F11corr from real It and Ir with real calibration EtaX: x-axis: F11corr(LDRtrue) / aF11sim0(LDRtrue) - 1" |
|
2003 F11corr = 1 / (TiO * TiE) * ( |
|
2004 (HR0 * Etax / K0 * It / TTa - HT0 * Ir / TRa) / (HR0 * GT0 - HT0 * GR0)) # IL = 1 Eq.(64); Etax/K0 = Eta0. |
|
2005 |
|
2006 # Text1 = "F11corr from It and Ir without corrections but with calibration EtaX: x-axis: F11corr(LDRtrue) devided by F11corr(LDRtrue = 0.004)" |
|
2007 # F11corr = 0.5/(TiO*TiE)*(Etax*It/TTa+Ir/TRa) # IL = 1 Eq.(64) |
|
2008 |
|
2009 # -- It from It only with atrue without corrections - for BERTHA (and PollyXTs) |
|
2010 # Text1 = " x-axis: IT(LDRtrue) / IT(LDRtrue = 0.004) - 1" |
|
2011 # F11corr = It/(TaT*TiT*TiO*TiE) #/(TaT*TiT*TiO*TiE*(GT0+atrue*HT0)) |
|
2012 # ! see below line 1673ff |
|
2013 |
|
2014 aF11corr[iLDR, iN] = F11corr |
|
2015 aLDRcorr[iLDR, iN] = LDRCorr # LDRCorr # LDRsim # for test only |
|
2016 aLDRsim[iLDR, iN] = LDRsim # LDRCorr # LDRsim # for test only |
|
2017 # aPLDR[iLDR, iN] = CalcPLDR(LDRCorr, BSR[iLDR], LDRm0) |
|
2018 aEtax[iLDR, iN] = Etax |
|
2019 aEtapx[iLDR, iN] = Etapx |
|
2020 aEtamx[iLDR, iN] = Etamx |
|
2021 |
|
2022 aGHK[0, iN] = GR |
|
2023 aGHK[1, iN] = GT |
|
2024 aGHK[2, iN] = HR |
|
2025 aGHK[3, iN] = HT |
|
2026 aGHK[4, iN] = K |
|
2027 |
|
2028 aLDRCal[iN] = iLDRCal |
|
2029 aQin[iN] = iQin |
|
2030 aVin[iN] = iVin |
|
2031 aERaT[iN] = iERaT |
|
2032 aERaR[iN] = iERaR |
|
2033 aRotaT[iN] = iRotaT |
|
2034 aRotaR[iN] = iRotaR |
|
2035 aRetT[iN] = iRetT |
|
2036 aRetR[iN] = iRetR |
|
2037 |
|
2038 aRotL[iN] = iRotL |
|
2039 aRotE[iN] = iRotE |
|
2040 aRetE[iN] = iRetE |
|
2041 aRotO[iN] = iRotO |
|
2042 aRetO[iN] = iRetO |
|
2043 aRotC[iN] = iRotC |
|
2044 aRetC[iN] = iRetC |
|
2045 aDiO[iN] = iDiO |
|
2046 aDiE[iN] = iDiE |
|
2047 aDiC[iN] = iDiC |
|
2048 aTP[iN] = iTP |
|
2049 aTS[iN] = iTS |
|
2050 aRP[iN] = iRP |
|
2051 aRS[iN] = iRS |
|
2052 aTCalT[iN] = iTCalT |
|
2053 aTCalR[iN] = iTCalR |
|
2054 |
|
2055 aNCalTp[iN] = iNCalTp # IoutTp, IoutTm, IoutRp, IoutRm => Etax |
|
2056 aNCalTm[iN] = iNCalTm # IoutTp, IoutTm, IoutRp, IoutRm => Etax |
|
2057 aNCalRp[iN] = iNCalRp # IoutTp, IoutTm, IoutRp, IoutRm => Etax |
|
2058 aNCalRm[iN] = iNCalRm # IoutTp, IoutTm, IoutRp, IoutRm => Etax |
|
2059 aNIt[iN] = iNIt # It, Tr |
|
2060 aNIr[iN] = iNIr # It, Tr |
|
2061 |
|
2062 # --- END loop |
|
2063 btime = clock() |
|
2064 # print("\r done in ", "{0:5.0f}".format(btime - atime), "sec. => producing plots now .... some more seconds ..."), # , end="\r"); |
|
2065 print(" done in ", "{0:5.0f}".format(btime - atime), "sec. => producing plots now .... some more seconds ...") |
|
2066 # --- Plot ----------------------------------------------------------------- |
|
2067 print("Errors from GHK correction uncertainties:") |
|
2068 if (sns_loaded): |
|
2069 sns.set_style("whitegrid") |
|
2070 sns.set_palette("bright6", 6) |
|
2071 # for older seaborn versions use: |
|
2072 # sns.set_palette("bright", 6) |
|
2073 |
|
2074 ''' |
|
2075 fig2 = plt.figure() |
|
2076 plt.plot(aLDRcorr[2,:],'b.') |
|
2077 plt.plot(aLDRcorr[3,:],'r.') |
|
2078 plt.plot(aLDRcorr[4,:],'g.') |
|
2079 #plt.plot(aLDRcorr[6,:],'c.') |
|
2080 plt.show |
|
2081 ''' |
|
2082 |
|
2083 # Plot LDR |
|
2084 def PlotSubHist(aVar, aX, X0, daX, iaX, naX): |
|
2085 # aVar is the name of the parameter and aX is the subset of aLDRcorr which is coloured in the plot |
|
2086 # example: PlotSubHist("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP) |
|
2087 fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) |
|
2088 iLDR = -1 |
|
2089 for LDRTrue in LDRrange: |
|
2090 aXmean = np.zeros(2 * naX + 1) |
|
2091 iLDR = iLDR + 1 |
|
2092 LDRmin[iLDR] = np.amin(aLDRcorr[iLDR, :]) |
|
2093 LDRmax[iLDR] = np.amax(aLDRcorr[iLDR, :]) |
|
2094 if (LDRmax[iLDR] > 10): LDRmax[iLDR] = 10 |
|
2095 if (LDRmin[iLDR] < -10): LDRmin[iLDR] = -10 |
|
2096 Rmin = LDRmin[iLDR] * 0.995 # np.min(aLDRcorr[iLDR,:]) * 0.995 |
|
2097 Rmax = LDRmax[iLDR] * 1.005 # np.max(aLDRcorr[iLDR,:]) * 1.005 |
|
2098 |
|
2099 # Determine mean distance of all aXmean from each other for each iLDR |
|
2100 meanDist = 0.0 |
|
2101 for iaX in range(-naX, naX + 1): |
|
2102 # mean LDRCorr value for certain error (iaX) of parameter aVar |
|
2103 aXmean[iaX + naX] = np.mean(aLDRcorr[iLDR, aX == iaX]) |
|
2104 # relative to absolute spread of LDRCorrs |
|
2105 meanDist = (np.max(aXmean) - np.min(aXmean)) / (LDRmax[iLDR] - LDRmin[iLDR]) * 100 |
|
2106 |
|
2107 plt.subplot(1, 5, iLDR + 1) |
|
2108 (n, bins, patches) = plt.hist(aLDRcorr[iLDR, :], |
|
2109 bins=100, log=False, |
|
2110 range=[Rmin, Rmax], |
|
2111 alpha=0.5, density=False, color='0.5', histtype='stepfilled') |
|
2112 |
|
2113 for iaX in range(-naX, naX + 1): |
|
2114 # mean LDRCorr value for certain error (iaX) of parameter aVar |
|
2115 plt.hist(aLDRcorr[iLDR, aX == iaX], |
|
2116 range=[Rmin, Rmax], |
|
2117 bins=100, log=False, alpha=0.3, density=False, histtype='stepfilled', |
|
2118 label=str(round(X0 + iaX * daX / naX, 5))) |
|
2119 |
|
2120 if (iLDR == 2): |
|
2121 leg = plt.legend() |
|
2122 leg.get_frame().set_alpha(0.1) |
|
2123 |
|
2124 plt.tick_params(axis='both', labelsize=10) |
|
2125 plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2) |
|
2126 plt.gca().set_title("{0:3.0f}%".format(meanDist)) |
|
2127 plt.gca().set_xlabel('LDRtrue', color="red") |
|
2128 |
|
2129 # plt.ylabel('frequency', fontsize=10) |
|
2130 # plt.xlabel('LDRCorr', fontsize=10) |
|
2131 # fig.tight_layout() |
|
2132 fig.suptitle(LID + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar + ' error contribution', fontsize=14, y=1.10) |
|
2133 # plt.show() |
|
2134 # fig.savefig(LID + '_' + aVar + '.png', dpi=150, bbox_inches='tight', pad_inches=0) |
|
2135 # plt.close |
|
2136 return |
|
2137 |
|
2138 def PlotLDRsim(aVar, aX, X0, daX, iaX, naX): |
|
2139 # aVar is the name of the parameter and aX is the subset of aLDRsim which is coloured in the plot |
|
2140 # example: PlotSubHist("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP) |
|
2141 fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) |
|
2142 iLDR = -1 |
|
2143 for LDRTrue in LDRrange: |
|
2144 aXmean = np.zeros(2 * naX + 1) |
|
2145 iLDR = iLDR + 1 |
|
2146 LDRsimmin[iLDR] = np.amin(aLDRsim[iLDR, :]) |
|
2147 LDRsimmax[iLDR] = np.amax(aLDRsim[iLDR, :]) |
|
2148 # print("LDRsimmin[iLDR], LDRsimmax[iLDR] = ", LDRsimmin[iLDR], LDRsimmax[iLDR]) |
|
2149 # if (LDRsimmax[iLDR] > 10): LDRsimmax[iLDR] = 10 |
|
2150 # if (LDRsimmin[iLDR] < -10): LDRsimmin[iLDR] = -10 |
|
2151 Rmin = LDRsimmin[iLDR] * 0.995 # np.min(aLDRsim[iLDR,:]) * 0.995 |
|
2152 Rmax = LDRsimmax[iLDR] * 1.005 # np.max(aLDRsim[iLDR,:]) * 1.005 |
|
2153 # print("Rmin, Rmax = ", Rmin, Rmax) |
|
2154 |
|
2155 # Determine mean distance of all aXmean from each other for each iLDR |
|
2156 meanDist = 0.0 |
|
2157 for iaX in range(-naX, naX + 1): |
|
2158 # mean LDRCorr value for certain error (iaX) of parameter aVar |
|
2159 aXmean[iaX + naX] = np.mean(aLDRsim[iLDR, aX == iaX]) |
|
2160 # relative to absolute spread of LDRCorrs |
|
2161 meanDist = (np.max(aXmean) - np.min(aXmean)) / (LDRsimmax[iLDR] - LDRsimmin[iLDR]) * 100 |
|
2162 |
|
2163 plt.subplot(1, 5, iLDR + 1) |
|
2164 (n, bins, patches) = plt.hist(aLDRsim[iLDR, :], |
|
2165 bins=100, log=False, |
|
2166 range=[Rmin, Rmax], |
|
2167 alpha=0.5, density=False, color='0.5', histtype='stepfilled') |
|
2168 |
|
2169 for iaX in range(-naX, naX + 1): |
|
2170 # mean LDRCorr value for certain error (iaX) of parameter aVar |
|
2171 plt.hist(aLDRsim[iLDR, aX == iaX], |
|
2172 range=[Rmin, Rmax], |
|
2173 bins=100, log=False, alpha=0.3, density=False, histtype='stepfilled', |
|
2174 label=str(round(X0 + iaX * daX / naX, 5))) |
|
2175 |
|
2176 if (iLDR == 2): |
|
2177 leg = plt.legend() |
|
2178 leg.get_frame().set_alpha(0.1) |
|
2179 |
|
2180 plt.tick_params(axis='both', labelsize=10) |
|
2181 plt.plot([LDRsim0[iLDR], LDRsim0[iLDR]], [0, np.max(n)], 'r-', lw=2) |
|
2182 plt.gca().set_title("{0:3.0f}%".format(meanDist)) |
|
2183 plt.gca().set_xlabel('LDRsim0', color="red") |
|
2184 |
|
2185 fig.suptitle('LDRsim - ' +LID + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar + ' error contribution', fontsize=14, y=1.10) |
|
2186 return |
|
2187 |
|
2188 |
|
2189 # Plot Etax |
|
2190 def PlotEtax(aVar, aX, X0, daX, iaX, naX): |
|
2191 # aVar is the name of the parameter and aX is the subset of aLDRcorr which is coloured in the plot |
|
2192 # example: PlotSubHist("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP) |
|
2193 fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) |
|
2194 iLDR = -1 |
|
2195 for LDRTrue in LDRrange: |
|
2196 aXmean = np.zeros(2 * naX + 1) |
|
2197 iLDR = iLDR + 1 |
|
2198 Etaxmin = np.amin(aEtax[iLDR, :]) |
|
2199 Etaxmax = np.amax(aEtax[iLDR, :]) |
|
2200 Rmin = Etaxmin * 0.995 # np.min(aLDRcorr[iLDR,:]) * 0.995 |
|
2201 Rmax = Etaxmax * 1.005 # np.max(aLDRcorr[iLDR,:]) * 1.005 |
|
2202 |
|
2203 # Determine mean distance of all aXmean from each other for each iLDR |
|
2204 meanDist = 0.0 |
|
2205 for iaX in range(-naX, naX + 1): |
|
2206 # mean Etax value for certain error (iaX) of parameter aVar |
|
2207 aXmean[iaX + naX] = np.mean(aEtax[iLDR, aX == iaX]) |
|
2208 # relative to absolute spread of Etax |
|
2209 meanDist = (np.max(aXmean) - np.min(aXmean)) / (Etaxmax - Etaxmin) * 100 |
|
2210 |
|
2211 plt.subplot(1, 5, iLDR + 1) |
|
2212 (n, bins, patches) = plt.hist(aEtax[iLDR, :], |
|
2213 bins=50, log=False, |
|
2214 range=[Rmin, Rmax], |
|
2215 alpha=0.5, density=False, color='0.5', histtype='stepfilled') |
|
2216 for iaX in range(-naX, naX + 1): |
|
2217 plt.hist(aEtax[iLDR, aX == iaX], |
|
2218 range=[Rmin, Rmax], |
|
2219 bins=50, log=False, alpha=0.3, density=False, histtype='stepfilled', |
|
2220 label=str(round(X0 + iaX * daX / naX, 5))) |
|
2221 if (iLDR == 2): |
|
2222 leg = plt.legend() |
|
2223 leg.get_frame().set_alpha(0.1) |
|
2224 plt.tick_params(axis='both', labelsize=10) |
|
2225 plt.plot([Etax0, Etax0], [0, np.max(n)], 'r-', lw=2) |
|
2226 plt.gca().set_title("{0:3.0f}%".format(meanDist)) |
|
2227 plt.gca().set_xlabel('Etax0', color="red") |
|
2228 fig.suptitle('Etax - ' + LID + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar + ' error contribution', fontsize=14, y=1.10) |
|
2229 return |
|
2230 |
|
2231 def PlotEtapx(aVar, aX, X0, daX, iaX, naX): |
|
2232 # aVar is the name of the parameter and aX is the subset of aLDRcorr which is coloured in the plot |
|
2233 # example: PlotSubHist("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP) |
|
2234 fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) |
|
2235 iLDR = -1 |
|
2236 for LDRTrue in LDRrange: |
|
2237 aXmean = np.zeros(2 * naX + 1) |
|
2238 iLDR = iLDR + 1 |
|
2239 Etapxmin = np.amin(aEtapx[iLDR, :]) |
|
2240 Etapxmax = np.amax(aEtapx[iLDR, :]) |
|
2241 Rmin = Etapxmin * 0.995 # np.min(aLDRcorr[iLDR,:]) * 0.995 |
|
2242 Rmax = Etapxmax * 1.005 # np.max(aLDRcorr[iLDR,:]) * 1.005 |
|
2243 |
|
2244 # Determine mean distance of all aXmean from each other for each iLDR |
|
2245 meanDist = 0.0 |
|
2246 for iaX in range(-naX, naX + 1): |
|
2247 # mean Etapx value for certain error (iaX) of parameter aVar |
|
2248 aXmean[iaX + naX] = np.mean(aEtapx[iLDR, aX == iaX]) |
|
2249 # relative to absolute spread of Etapx |
|
2250 meanDist = (np.max(aXmean) - np.min(aXmean)) / (Etapxmax - Etapxmin) * 100 |
|
2251 |
|
2252 plt.subplot(1, 5, iLDR + 1) |
|
2253 (n, bins, patches) = plt.hist(aEtapx[iLDR, :], |
|
2254 bins=50, log=False, |
|
2255 range=[Rmin, Rmax], |
|
2256 alpha=0.5, density=False, color='0.5', histtype='stepfilled') |
|
2257 for iaX in range(-naX, naX + 1): |
|
2258 plt.hist(aEtapx[iLDR, aX == iaX], |
|
2259 range=[Rmin, Rmax], |
|
2260 bins=50, log=False, alpha=0.3, density=False, histtype='stepfilled', |
|
2261 label=str(round(X0 + iaX * daX / naX, 5))) |
|
2262 if (iLDR == 2): |
|
2263 leg = plt.legend() |
|
2264 leg.get_frame().set_alpha(0.1) |
|
2265 plt.tick_params(axis='both', labelsize=10) |
|
2266 plt.plot([Etapx0, Etapx0], [0, np.max(n)], 'r-', lw=2) |
|
2267 plt.gca().set_title("{0:3.0f}%".format(meanDist)) |
|
2268 plt.gca().set_xlabel('Etapx0', color="red") |
|
2269 fig.suptitle('Etapx - ' + LID + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar + ' error contribution', fontsize=14, y=1.10) |
|
2270 return |
|
2271 |
|
2272 def PlotEtamx(aVar, aX, X0, daX, iaX, naX): |
|
2273 # aVar is the name of the parameter and aX is the subset of aLDRcorr which is coloured in the plot |
|
2274 # example: PlotSubHist("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP) |
|
2275 fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) |
|
2276 iLDR = -1 |
|
2277 for LDRTrue in LDRrange: |
|
2278 aXmean = np.zeros(2 * naX + 1) |
|
2279 iLDR = iLDR + 1 |
|
2280 Etamxmin = np.amin(aEtamx[iLDR, :]) |
|
2281 Etamxmax = np.amax(aEtamx[iLDR, :]) |
|
2282 Rmin = Etamxmin * 0.995 # np.min(aLDRcorr[iLDR,:]) * 0.995 |
|
2283 Rmax = Etamxmax * 1.005 # np.max(aLDRcorr[iLDR,:]) * 1.005 |
|
2284 |
|
2285 # Determine mean distance of all aXmean from each other for each iLDR |
|
2286 meanDist = 0.0 |
|
2287 for iaX in range(-naX, naX + 1): |
|
2288 # mean Etamx value for certain error (iaX) of parameter aVar |
|
2289 aXmean[iaX + naX] = np.mean(aEtamx[iLDR, aX == iaX]) |
|
2290 # relative to absolute spread of Etamx |
|
2291 meanDist = (np.max(aXmean) - np.min(aXmean)) / (Etamxmax - Etamxmin) * 100 |
|
2292 |
|
2293 plt.subplot(1, 5, iLDR + 1) |
|
2294 (n, bins, patches) = plt.hist(aEtamx[iLDR, :], |
|
2295 bins=50, log=False, |
|
2296 range=[Rmin, Rmax], |
|
2297 alpha=0.5, density=False, color='0.5', histtype='stepfilled') |
|
2298 for iaX in range(-naX, naX + 1): |
|
2299 plt.hist(aEtamx[iLDR, aX == iaX], |
|
2300 range=[Rmin, Rmax], |
|
2301 bins=50, log=False, alpha=0.3, density=False, histtype='stepfilled', |
|
2302 label=str(round(X0 + iaX * daX / naX, 5))) |
|
2303 if (iLDR == 2): |
|
2304 leg = plt.legend() |
|
2305 leg.get_frame().set_alpha(0.1) |
|
2306 plt.tick_params(axis='both', labelsize=10) |
|
2307 plt.plot([Etamx0, Etamx0], [0, np.max(n)], 'r-', lw=2) |
|
2308 plt.gca().set_title("{0:3.0f}%".format(meanDist)) |
|
2309 plt.gca().set_xlabel('Etamx0', color="red") |
|
2310 fig.suptitle('Etamx - ' + LID + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar + ' error contribution', fontsize=14, y=1.10) |
|
2311 return |
|
2312 |
|
2313 # calc contribution of the error of aVar = aX to aY for each LDRtrue |
|
2314 def Contribution(aVar, aX, X0, daX, iaX, naX, aY, Ysum, widthSum): |
|
2315 # aVar is the name of the parameter and aX is the subset of aY which is coloured in the plot |
|
2316 # example: Contribution("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP, aLDRcorr, DOLPcontr) |
|
2317 iLDR = -1 |
|
2318 # Ysum, widthSum = np.zeros(5) |
|
2319 meanDist = np.zeros(5) # iLDR |
|
2320 widthDist = np.zeros(5) # iLDR |
|
2321 for LDRTrue in LDRrange: |
|
2322 aXmean = np.zeros(2 * naX + 1) |
|
2323 aXwidth = np.zeros(2 * naX + 1) |
|
2324 iLDR = iLDR + 1 |
|
2325 # total width of distribution |
|
2326 aYmin = np.amin(aY[iLDR, :]) |
|
2327 aYmax = np.amax(aY[iLDR, :]) |
|
2328 aYwidth = aYmax - aYmin |
|
2329 # Determine mean distance of all aXmean from each other for each iLDR |
|
2330 for iaX in range(-naX, naX + 1): |
|
2331 # mean LDRCorr value for all errors iaX of parameter aVar |
|
2332 aXmean[iaX + naX] = np.mean(aY[iLDR, aX == iaX]) |
|
2333 aXwidth[iaX + naX] = np.max(aY[iLDR, aX == iaX]) - np.min(aY[iLDR, aX == iaX]) |
|
2334 # relative to absolute spread of LDRCorrs |
|
2335 meanDist[iLDR] = (np.max(aXmean) - np.min(aXmean)) / aYwidth * 1000 |
|
2336 # meanDist[iLDR] = (aYwidth - aXwidth[naX]) / aYwidth * 1000 |
|
2337 widthDist[iLDR] = (np.max(aXwidth) - aXwidth[naX]) / aYwidth * 1000 |
|
2338 |
|
2339 print("{:12}{:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f}"\ |
|
2340 .format(aVar,meanDist[0],meanDist[1],meanDist[2],meanDist[3],meanDist[4],widthDist[0],widthDist[1],widthDist[2],widthDist[3],widthDist[4])) |
|
2341 Ysum = Ysum + meanDist |
|
2342 widthSum = widthSum + widthDist |
|
2343 return(Ysum, widthSum) |
|
2344 |
|
2345 # print(.format(LDRrangeA[iLDR],)) |
|
2346 |
|
2347 # error contributions to a certain output aY; loop over all variables |
|
2348 def Contribution_aY(aYvar, aY): |
|
2349 Ysum = np.zeros(5) |
|
2350 widthSum = np.zeros(5) |
|
2351 # meanDist = np.zeros(5) # iLDR |
|
2352 LDRrangeA = np.array(LDRrange) |
|
2353 print() |
|
2354 print(aYvar + ": contribution to the total error (per mill)") |
|
2355 print(" of individual parameter errors of combined parameter errors") |
|
2356 print(" at LDRtrue {:5.3f} {:5.3f} {:5.3f} {:5.3f} {:5.3f} {:5.3f} {:5.3f} {:5.3f} {:5.3f} {:5.3f}"\ |
|
2357 .format(LDRrangeA[0],LDRrangeA[1],LDRrangeA[2],LDRrangeA[3],LDRrangeA[4],LDRrangeA[0],LDRrangeA[1],LDRrangeA[2],LDRrangeA[3],LDRrangeA[4])) |
|
2358 print() |
|
2359 if (nQin > 0): Ysum, widthSum = Contribution("Qin", aQin, Qin0, dQin, iQin, nQin, aY, Ysum, widthSum) |
|
2360 if (nVin > 0): Ysum, widthSum = Contribution("Vin", aVin, Vin0, dVin, iVin, nVin, aY, Ysum, widthSum) |
|
2361 if (nRotL > 0): Ysum, widthSum = Contribution("RotL", aRotL, RotL0, dRotL, iRotL, nRotL, aY, Ysum, widthSum) |
|
2362 if (nRetE > 0): Ysum, widthSum = Contribution("RetE", aRetE, RetE0, dRetE, iRetE, nRetE, aY, Ysum, widthSum) |
|
2363 if (nRotE > 0): Ysum, widthSum = Contribution("RotE", aRotE, RotE0, dRotE, iRotE, nRotE, aY, Ysum, widthSum) |
|
2364 if (nDiE > 0): Ysum, widthSum = Contribution("DiE", aDiE, DiE0, dDiE, iDiE, nDiE, aY, Ysum, widthSum) |
|
2365 if (nRetO > 0): Ysum, widthSum = Contribution("RetO", aRetO, RetO0, dRetO, iRetO, nRetO, aY, Ysum, widthSum) |
|
2366 if (nRotO > 0): Ysum, widthSum = Contribution("RotO", aRotO, RotO0, dRotO, iRotO, nRotO, aY, Ysum, widthSum) |
|
2367 if (nDiO > 0): Ysum, widthSum = Contribution("DiO", aDiO, DiO0, dDiO, iDiO, nDiO, aY, Ysum, widthSum) |
|
2368 if (nDiC > 0): Ysum, widthSum = Contribution("DiC", aDiC, DiC0, dDiC, iDiC, nDiC, aY, Ysum, widthSum) |
|
2369 if (nRotC > 0): Ysum, widthSum = Contribution("RotC", aRotC, RotC0, dRotC, iRotC, nRotC, aY, Ysum, widthSum) |
|
2370 if (nRetC > 0): Ysum, widthSum = Contribution("RetC", aRetC, RetC0, dRetC, iRetC, nRetC, aY, Ysum, widthSum) |
|
2371 if (nTP > 0): Ysum, widthSum = Contribution("TP", aTP, TP0, dTP, iTP, nTP, aY, Ysum, widthSum) |
|
2372 if (nTS > 0): Ysum, widthSum = Contribution("TS", aTS, TS0, dTS, iTS, nTS, aY, Ysum, widthSum) |
|
2373 if (nRP > 0): Ysum, widthSum = Contribution("RP", aRP, RP0, dRP, iRP, nRP, aY, Ysum, widthSum) |
|
2374 if (nRS > 0): Ysum, widthSum = Contribution("RS", aRS, RS0, dRS, iRS, nRS, aY, Ysum, widthSum) |
|
2375 if (nRetT > 0): Ysum, widthSum = Contribution("RetT", aRetT, RetT0, dRetT, iRetT, nRetT, aY, Ysum, widthSum) |
|
2376 if (nRetR > 0): Ysum, widthSum = Contribution("RetR", aRetR, RetR0, dRetR, iRetR, nRetR, aY, Ysum, widthSum) |
|
2377 if (nERaT > 0): Ysum, widthSum = Contribution("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT, aY, Ysum, widthSum) |
|
2378 if (nERaR > 0): Ysum, widthSum = Contribution("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR, aY, Ysum, widthSum) |
|
2379 if (nRotaT > 0): Ysum, widthSum = Contribution("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT, aY, Ysum, widthSum) |
|
2380 if (nRotaR > 0): Ysum, widthSum = Contribution("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR, aY, Ysum, widthSum) |
|
2381 if (nLDRCal > 0): Ysum, widthSum = Contribution("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal, aY, Ysum, widthSum) |
|
2382 if (nTCalT > 0): Ysum, widthSum = Contribution("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT, aY, Ysum, widthSum) |
|
2383 if (nTCalR > 0): Ysum, widthSum = Contribution("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR, aY, Ysum, widthSum) |
|
2384 if (nNCal > 0): Ysum, widthSum = Contribution("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal, aY, Ysum, widthSum) |
|
2385 if (nNCal > 0): Ysum, widthSum = Contribution("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal, aY, Ysum, widthSum) |
|
2386 if (nNCal > 0): Ysum, widthSum = Contribution("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal, aY, Ysum, widthSum) |
|
2387 if (nNCal > 0): Ysum, widthSum = Contribution("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal, aY, Ysum, widthSum) |
|
2388 if (nNI > 0): Ysum, widthSum = Contribution("SigNoiseIt", aNIt, 0, 1, iNIt, nNI, aY, Ysum, widthSum) |
|
2389 if (nNI > 0): Ysum, widthSum = Contribution("SigNoiseIr", aNIr, 0, 1, iNIr, nNI, aY, Ysum, widthSum) |
|
2390 print("{:12}{:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f}"\ |
|
2391 .format("Sum ",Ysum[0],Ysum[1],Ysum[2],Ysum[3],Ysum[4],widthSum[0],widthSum[1],widthSum[2],widthSum[3],widthSum[4])) |
|
2392 |
|
2393 |
|
2394 # Plot LDR histograms |
|
2395 if (nQin > 0): PlotSubHist("Qin", aQin, Qin0, dQin, iQin, nQin) |
|
2396 if (nVin > 0): PlotSubHist("Vin", aVin, Vin0, dVin, iVin, nVin) |
|
2397 if (nRotL > 0): PlotSubHist("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) |
|
2398 if (nRetE > 0): PlotSubHist("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) |
|
2399 if (nRotE > 0): PlotSubHist("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) |
|
2400 if (nDiE > 0): PlotSubHist("DiE", aDiE, DiE0, dDiE, iDiE, nDiE) |
|
2401 if (nRetO > 0): PlotSubHist("RetO", aRetO, RetO0, dRetO, iRetO, nRetO) |
|
2402 if (nRotO > 0): PlotSubHist("RotO", aRotO, RotO0, dRotO, iRotO, nRotO) |
|
2403 if (nDiO > 0): PlotSubHist("DiO", aDiO, DiO0, dDiO, iDiO, nDiO) |
|
2404 if (nDiC > 0): PlotSubHist("DiC", aDiC, DiC0, dDiC, iDiC, nDiC) |
|
2405 if (nRotC > 0): PlotSubHist("RotC", aRotC, RotC0, dRotC, iRotC, nRotC) |
|
2406 if (nRetC > 0): PlotSubHist("RetC", aRetC, RetC0, dRetC, iRetC, nRetC) |
|
2407 if (nTP > 0): PlotSubHist("TP", aTP, TP0, dTP, iTP, nTP) |
|
2408 if (nTS > 0): PlotSubHist("TS", aTS, TS0, dTS, iTS, nTS) |
|
2409 if (nRP > 0): PlotSubHist("RP", aRP, RP0, dRP, iRP, nRP) |
|
2410 if (nRS > 0): PlotSubHist("RS", aRS, RS0, dRS, iRS, nRS) |
|
2411 if (nRetT > 0): PlotSubHist("RetT", aRetT, RetT0, dRetT, iRetT, nRetT) |
|
2412 if (nRetR > 0): PlotSubHist("RetR", aRetR, RetR0, dRetR, iRetR, nRetR) |
|
2413 if (nERaT > 0): PlotSubHist("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) |
|
2414 if (nERaR > 0): PlotSubHist("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) |
|
2415 if (nRotaT > 0): PlotSubHist("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) |
|
2416 if (nRotaR > 0): PlotSubHist("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) |
|
2417 if (nLDRCal > 0): PlotSubHist("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) |
|
2418 if (nTCalT > 0): PlotSubHist("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT) |
|
2419 if (nTCalR > 0): PlotSubHist("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR) |
|
2420 if (nNCal > 0): PlotSubHist("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal) |
|
2421 if (nNCal > 0): PlotSubHist("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal) |
|
2422 if (nNCal > 0): PlotSubHist("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal) |
|
2423 if (nNCal > 0): PlotSubHist("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal) |
|
2424 if (nNI > 0): PlotSubHist("SigNoiseIt", aNIt, 0, 1, iNIt, nNI) |
|
2425 if (nNI > 0): PlotSubHist("SigNoiseIr", aNIr, 0, 1, iNIr, nNI) |
|
2426 plt.show() |
|
2427 plt.close |
|
2428 |
|
2429 |
|
2430 |
|
2431 # --- Plot LDRmin, LDRmax |
|
2432 iLDR = -1 |
|
2433 for LDRTrue in LDRrange: |
|
2434 iLDR = iLDR + 1 |
|
2435 LDRmin[iLDR] = np.amin(aLDRcorr[iLDR, :]) |
|
2436 LDRmax[iLDR] = np.amax(aLDRcorr[iLDR, :]) |
|
2437 LDRstd[iLDR] = np.std(aLDRcorr[iLDR, :]) |
|
2438 LDRmean[iLDR] = np.mean(aLDRcorr[iLDR, :]) |
|
2439 LDRmedian[iLDR] = np.median(aLDRcorr[iLDR, :]) |
|
2440 LDRskew[iLDR] = skew(aLDRcorr[iLDR, :],bias=False) |
|
2441 LDRkurt[iLDR] = kurtosis(aLDRcorr[iLDR, :],fisher=True,bias=False) |
|
2442 |
|
2443 fig2 = plt.figure() |
|
2444 LDRrangeA = np.array(LDRrange) |
|
2445 if((np.amax(LDRmax - LDRrangeA)-np.amin(LDRmin - LDRrangeA)) < 0.001): |
|
2446 plt.ylim(-0.001,0.001) |
|
2447 plt.plot(LDRrangeA, LDRmax - LDRrangeA, linewidth=2.0, color='b') |
|
2448 plt.plot(LDRrangeA, LDRmin - LDRrangeA, linewidth=2.0, color='g') |
|
2449 |
|
2450 plt.xlabel('LDRtrue', fontsize=18) |
|
2451 plt.ylabel('LDRTrue-LDRmin, LDRTrue-LDRmax', fontsize=14) |
|
2452 plt.title(LID + ' ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]), fontsize=18) |
|
2453 # plt.ylimit(-0.07, 0.07) |
|
2454 plt.show() |
|
2455 plt.close |
|
2456 |
|
2457 # --- Save LDRmin, LDRmax to file |
|
2458 # http://stackoverflow.com/questions/4675728/redirect-stdout-to-a-file-in-python |
|
2459 with open('output_files\\' + OutputFile, 'a') as f: |
|
2460 # with open('output_files\\' + LID + '-' + InputFile[0:-3] + '-LDR_min_max.dat', 'w') as f: |
|
2461 with redirect_stdout(f): |
|
2462 print("Lidar ID: " + LID) |
|
2463 print() |
|
2464 print("minimum and maximum values of the distributions of possibly measured LDR for different LDRtrue") |
|
2465 print("LDRtrue , LDRmin, LDRmax") |
|
2466 for i in range(len(LDRrangeA)): |
|
2467 print("{0:7.4f},{1:7.4f},{2:7.4f}".format(LDRrangeA[i], LDRmin[i], LDRmax[i])) |
|
2468 print() |
|
2469 # Print LDR statistics |
|
2470 print("LDRtrue , mean , median, max-mean, min-mean, std, excess_kurtosis, skewness") |
|
2471 iLDR = -1 |
|
2472 LDRrangeA = np.array(LDRrange) |
|
2473 for LDRTrue in LDRrange: |
|
2474 iLDR = iLDR + 1 |
|
2475 print("{0:8.5f},{1:8.5f},{2:8.5f}, {3:8.5f},{4:8.5f},{5:8.5f}, {6:8.5f},{7:8.5f}"\ |
|
2476 .format(LDRrangeA[iLDR], LDRmean[iLDR], LDRmedian[iLDR], LDRmax[iLDR]-LDRrangeA[iLDR], \ |
|
2477 LDRmin[iLDR]-LDRrangeA[iLDR], LDRstd[iLDR], LDRkurt[iLDR], LDRskew[iLDR])) |
|
2478 print() |
|
2479 # Calculate and print statistics for calibration factors |
|
2480 print("minimum and maximum values of the distributions of signal ratios and calibration factors for different LDRtrue") |
|
2481 iLDR = -1 |
|
2482 LDRrangeA = np.array(LDRrange) |
|
2483 print("LDRtrue , LDRsim, (max-min)/2, relerr") |
|
2484 for LDRTrue in LDRrange: |
|
2485 iLDR = iLDR + 1 |
|
2486 LDRsimmin[iLDR] = np.amin(aLDRsim[iLDR, :]) |
|
2487 LDRsimmax[iLDR] = np.amax(aLDRsim[iLDR, :]) |
|
2488 # LDRsimstd = np.std(aLDRsim[iLDR, :]) |
|
2489 LDRsimmean[iLDR] = np.mean(aLDRsim[iLDR, :]) |
|
2490 # LDRsimmedian = np.median(aLDRsim[iLDR, :]) |
|
2491 print("{0:8.5f}, {1:8.5f}, {2:8.5f}, {3:8.5f}".format(LDRrangeA[iLDR],LDRsimmean[iLDR],(LDRsimmax[iLDR]-LDRsimmin[iLDR])/2,(LDRsimmax[iLDR]-LDRsimmin[iLDR])/2/LDRsimmean[iLDR])) |
|
2492 iLDR = -1 |
|
2493 print("LDRtrue , Etax , (max-min)/2, relerr") |
|
2494 for LDRTrue in LDRrange: |
|
2495 iLDR = iLDR + 1 |
|
2496 Etaxmin = np.amin(aEtax[iLDR, :]) |
|
2497 Etaxmax = np.amax(aEtax[iLDR, :]) |
|
2498 # Etaxstd = np.std(aEtax[iLDR, :]) |
|
2499 Etaxmean = np.mean(aEtax[iLDR, :]) |
|
2500 # Etaxmedian = np.median(aEtax[iLDR, :]) |
|
2501 print("{0:8.5f}, {1:8.5f}, {2:8.5f}, {3:8.5f}".format(LDRrangeA[iLDR], Etaxmean, (Etaxmax-Etaxmin)/2, (Etaxmax-Etaxmin)/2/Etaxmean)) |
|
2502 iLDR = -1 |
|
2503 print("LDRtrue , Etapx , (max-min)/2, relerr") |
|
2504 for LDRTrue in LDRrange: |
|
2505 iLDR = iLDR + 1 |
|
2506 Etapxmin = np.amin(aEtapx[iLDR, :]) |
|
2507 Etapxmax = np.amax(aEtapx[iLDR, :]) |
|
2508 # Etapxstd = np.std(aEtapx[iLDR, :]) |
|
2509 Etapxmean = np.mean(aEtapx[iLDR, :]) |
|
2510 # Etapxmedian = np.median(aEtapx[iLDR, :]) |
|
2511 print("{0:8.5f}, {1:8.5f}, {2:8.5f}, {3:8.5f}".format(LDRrangeA[iLDR], Etapxmean, (Etapxmax-Etapxmin)/2, (Etapxmax-Etapxmin)/2/Etapxmean)) |
|
2512 iLDR = -1 |
|
2513 print("LDRtrue , Etamx , (max-min)/2, relerr") |
|
2514 for LDRTrue in LDRrange: |
|
2515 iLDR = iLDR + 1 |
|
2516 Etamxmin = np.amin(aEtamx[iLDR, :]) |
|
2517 Etamxmax = np.amax(aEtamx[iLDR, :]) |
|
2518 # Etamxstd = np.std(aEtamx[iLDR, :]) |
|
2519 Etamxmean = np.mean(aEtamx[iLDR, :]) |
|
2520 # Etamxmedian = np.median(aEtamx[iLDR, :]) |
|
2521 print("{0:8.5f}, {1:8.5f}, {2:8.5f}, {3:8.5f}".format(LDRrangeA[iLDR], Etamxmean, (Etamxmax-Etamxmin)/2, (Etamxmax-Etamxmin)/2/Etamxmean)) |
|
2522 |
|
2523 # Print LDR statistics |
|
2524 print("LDRtrue , mean , median, max-mean, min-mean, std, excess_kurtosis, skewness") |
|
2525 iLDR = -1 |
|
2526 LDRrangeA = np.array(LDRrange) |
|
2527 for LDRTrue in LDRrange: |
|
2528 iLDR = iLDR + 1 |
|
2529 print("{0:8.5f},{1:8.5f},{2:8.5f}, {3:8.5f},{4:8.5f},{5:8.5f}, {6:8.5f},{7:8.5f}".format(LDRrangeA[iLDR], LDRmean[iLDR], LDRmedian[iLDR], LDRmax[iLDR]-LDRrangeA[iLDR], LDRmin[iLDR]-LDRrangeA[iLDR], LDRstd[iLDR],LDRkurt[iLDR],LDRskew[iLDR])) |
|
2530 |
|
2531 |
|
2532 with open('output_files\\' + OutputFile, 'a') as f: |
|
2533 # with open('output_files\\' + LID + '-' + InputFile[0:-3] + '-LDR_min_max.dat', 'a') as f: |
|
2534 with redirect_stdout(f): |
|
2535 Contribution_aY("LDRCorr", aLDRcorr) |
|
2536 Contribution_aY("LDRsim", aLDRsim) |
|
2537 Contribution_aY("EtaX, D90", aEtax) |
|
2538 Contribution_aY("Etapx, +45°", aEtapx) |
|
2539 Contribution_aY("Etamx -45°", aEtamx) |
|
2540 |
|
2541 |
|
2542 # Plot other histograms |
|
2543 if (bPlotEtax): |
|
2544 |
|
2545 if (nQin > 0): PlotLDRsim("Qin", aQin, Qin0, dQin, iQin, nQin) |
|
2546 if (nVin > 0): PlotLDRsim("Vin", aVin, Vin0, dVin, iVin, nVin) |
|
2547 if (nRotL > 0): PlotLDRsim("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) |
|
2548 if (nRetE > 0): PlotLDRsim("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) |
|
2549 if (nRotE > 0): PlotLDRsim("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) |
|
2550 if (nDiE > 0): PlotLDRsim("DiE", aDiE, DiE0, dDiE, iDiE, nDiE) |
|
2551 if (nRetO > 0): PlotLDRsim("RetO", aRetO, RetO0, dRetO, iRetO, nRetO) |
|
2552 if (nRotO > 0): PlotLDRsim("RotO", aRotO, RotO0, dRotO, iRotO, nRotO) |
|
2553 if (nDiO > 0): PlotLDRsim("DiO", aDiO, DiO0, dDiO, iDiO, nDiO) |
|
2554 if (nDiC > 0): PlotLDRsim("DiC", aDiC, DiC0, dDiC, iDiC, nDiC) |
|
2555 if (nRotC > 0): PlotLDRsim("RotC", aRotC, RotC0, dRotC, iRotC, nRotC) |
|
2556 if (nRetC > 0): PlotLDRsim("RetC", aRetC, RetC0, dRetC, iRetC, nRetC) |
|
2557 if (nTP > 0): PlotLDRsim("TP", aTP, TP0, dTP, iTP, nTP) |
|
2558 if (nTS > 0): PlotLDRsim("TS", aTS, TS0, dTS, iTS, nTS) |
|
2559 if (nRP > 0): PlotLDRsim("RP", aRP, RP0, dRP, iRP, nRP) |
|
2560 if (nRS > 0): PlotLDRsim("RS", aRS, RS0, dRS, iRS, nRS) |
|
2561 if (nRetT > 0): PlotLDRsim("RetT", aRetT, RetT0, dRetT, iRetT, nRetT) |
|
2562 if (nRetR > 0): PlotLDRsim("RetR", aRetR, RetR0, dRetR, iRetR, nRetR) |
|
2563 if (nERaT > 0): PlotLDRsim("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) |
|
2564 if (nERaR > 0): PlotLDRsim("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) |
|
2565 if (nRotaT > 0): PlotLDRsim("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) |
|
2566 if (nRotaR > 0): PlotLDRsim("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) |
|
2567 if (nLDRCal > 0): PlotLDRsim("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) |
|
2568 if (nTCalT > 0): PlotLDRsim("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT) |
|
2569 if (nTCalR > 0): PlotLDRsim("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR) |
|
2570 if (nNCal > 0): PlotLDRsim("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal) |
|
2571 if (nNCal > 0): PlotLDRsim("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal) |
|
2572 if (nNCal > 0): PlotLDRsim("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal) |
|
2573 if (nNCal > 0): PlotLDRsim("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal) |
|
2574 if (nNI > 0): PlotLDRsim("SigNoiseIt", aNIt, 0, 1, iNIt, nNI) |
|
2575 if (nNI > 0): PlotLDRsim("SigNoiseIr", aNIr, 0, 1, iNIr, nNI) |
|
2576 plt.show() |
|
2577 plt.close |
|
2578 print("---------------------------------------...producing more plots...------------------------------------------------------------------") |
|
2579 |
|
2580 if (nQin > 0): PlotEtax("Qin", aQin, Qin0, dQin, iQin, nQin) |
|
2581 if (nVin > 0): PlotEtax("Vin", aVin, Vin0, dVin, iVin, nVin) |
|
2582 if (nRotL > 0): PlotEtax("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) |
|
2583 if (nRetE > 0): PlotEtax("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) |
|
2584 if (nRotE > 0): PlotEtax("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) |
|
2585 if (nDiE > 0): PlotEtax("DiE", aDiE, DiE0, dDiE, iDiE, nDiE) |
|
2586 if (nRetO > 0): PlotEtax("RetO", aRetO, RetO0, dRetO, iRetO, nRetO) |
|
2587 if (nRotO > 0): PlotEtax("RotO", aRotO, RotO0, dRotO, iRotO, nRotO) |
|
2588 if (nDiO > 0): PlotEtax("DiO", aDiO, DiO0, dDiO, iDiO, nDiO) |
|
2589 if (nDiC > 0): PlotEtax("DiC", aDiC, DiC0, dDiC, iDiC, nDiC) |
|
2590 if (nRotC > 0): PlotEtax("RotC", aRotC, RotC0, dRotC, iRotC, nRotC) |
|
2591 if (nRetC > 0): PlotEtax("RetC", aRetC, RetC0, dRetC, iRetC, nRetC) |
|
2592 if (nTP > 0): PlotEtax("TP", aTP, TP0, dTP, iTP, nTP) |
|
2593 if (nTS > 0): PlotEtax("TS", aTS, TS0, dTS, iTS, nTS) |
|
2594 if (nRP > 0): PlotEtax("RP", aRP, RP0, dRP, iRP, nRP) |
|
2595 if (nRS > 0): PlotEtax("RS", aRS, RS0, dRS, iRS, nRS) |
|
2596 if (nRetT > 0): PlotEtax("RetT", aRetT, RetT0, dRetT, iRetT, nRetT) |
|
2597 if (nRetR > 0): PlotEtax("RetR", aRetR, RetR0, dRetR, iRetR, nRetR) |
|
2598 if (nERaT > 0): PlotEtax("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) |
|
2599 if (nERaR > 0): PlotEtax("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) |
|
2600 if (nRotaT > 0): PlotEtax("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) |
|
2601 if (nRotaR > 0): PlotEtax("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) |
|
2602 if (nLDRCal > 0): PlotEtax("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) |
|
2603 if (nTCalT > 0): PlotEtax("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT) |
|
2604 if (nTCalR > 0): PlotEtax("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR) |
|
2605 if (nNCal > 0): PlotEtax("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal) |
|
2606 if (nNCal > 0): PlotEtax("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal) |
|
2607 if (nNCal > 0): PlotEtax("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal) |
|
2608 if (nNCal > 0): PlotEtax("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal) |
|
2609 if (nNI > 0): PlotEtax("SigNoiseIt", aNIt, 0, 1, iNIt, nNI) |
|
2610 if (nNI > 0): PlotEtax("SigNoiseIr", aNIr, 0, 1, iNIr, nNI) |
|
2611 plt.show() |
|
2612 plt.close |
|
2613 print("---------------------------------------...producing more plots...------------------------------------------------------------------") |
|
2614 |
|
2615 if (nQin > 0): PlotEtapx("Qin", aQin, Qin0, dQin, iQin, nQin) |
|
2616 if (nVin > 0): PlotEtapx("Vin", aVin, Vin0, dVin, iVin, nVin) |
|
2617 if (nRotL > 0): PlotEtapx("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) |
|
2618 if (nRetE > 0): PlotEtapx("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) |
|
2619 if (nRotE > 0): PlotEtapx("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) |
|
2620 if (nDiE > 0): PlotEtapx("DiE", aDiE, DiE0, dDiE, iDiE, nDiE) |
|
2621 if (nRetO > 0): PlotEtapx("RetO", aRetO, RetO0, dRetO, iRetO, nRetO) |
|
2622 if (nRotO > 0): PlotEtapx("RotO", aRotO, RotO0, dRotO, iRotO, nRotO) |
|
2623 if (nDiO > 0): PlotEtapx("DiO", aDiO, DiO0, dDiO, iDiO, nDiO) |
|
2624 if (nDiC > 0): PlotEtapx("DiC", aDiC, DiC0, dDiC, iDiC, nDiC) |
|
2625 if (nRotC > 0): PlotEtapx("RotC", aRotC, RotC0, dRotC, iRotC, nRotC) |
|
2626 if (nRetC > 0): PlotEtapx("RetC", aRetC, RetC0, dRetC, iRetC, nRetC) |
|
2627 if (nTP > 0): PlotEtapx("TP", aTP, TP0, dTP, iTP, nTP) |
|
2628 if (nTS > 0): PlotEtapx("TS", aTS, TS0, dTS, iTS, nTS) |
|
2629 if (nRP > 0): PlotEtapx("RP", aRP, RP0, dRP, iRP, nRP) |
|
2630 if (nRS > 0): PlotEtapx("RS", aRS, RS0, dRS, iRS, nRS) |
|
2631 if (nRetT > 0): PlotEtapx("RetT", aRetT, RetT0, dRetT, iRetT, nRetT) |
|
2632 if (nRetR > 0): PlotEtapx("RetR", aRetR, RetR0, dRetR, iRetR, nRetR) |
|
2633 if (nERaT > 0): PlotEtapx("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) |
|
2634 if (nERaR > 0): PlotEtapx("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) |
|
2635 if (nRotaT > 0): PlotEtapx("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) |
|
2636 if (nRotaR > 0): PlotEtapx("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) |
|
2637 if (nLDRCal > 0): PlotEtapx("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) |
|
2638 if (nTCalT > 0): PlotEtapx("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT) |
|
2639 if (nTCalR > 0): PlotEtapx("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR) |
|
2640 if (nNCal > 0): PlotEtapx("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal) |
|
2641 if (nNCal > 0): PlotEtapx("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal) |
|
2642 if (nNCal > 0): PlotEtapx("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal) |
|
2643 if (nNCal > 0): PlotEtapx("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal) |
|
2644 if (nNI > 0): PlotEtapx("SigNoiseIt", aNIt, 0, 1, iNIt, nNI) |
|
2645 if (nNI > 0): PlotEtapx("SigNoiseIr", aNIr, 0, 1, iNIr, nNI) |
|
2646 plt.show() |
|
2647 plt.close |
|
2648 print("---------------------------------------...producing more plots...------------------------------------------------------------------") |
|
2649 |
|
2650 if (nQin > 0): PlotEtamx("Qin", aQin, Qin0, dQin, iQin, nQin) |
|
2651 if (nVin > 0): PlotEtamx("Vin", aVin, Vin0, dVin, iVin, nVin) |
|
2652 if (nRotL > 0): PlotEtamx("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) |
|
2653 if (nRetE > 0): PlotEtamx("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) |
|
2654 if (nRotE > 0): PlotEtamx("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) |
|
2655 if (nDiE > 0): PlotEtamx("DiE", aDiE, DiE0, dDiE, iDiE, nDiE) |
|
2656 if (nRetO > 0): PlotEtamx("RetO", aRetO, RetO0, dRetO, iRetO, nRetO) |
|
2657 if (nRotO > 0): PlotEtamx("RotO", aRotO, RotO0, dRotO, iRotO, nRotO) |
|
2658 if (nDiO > 0): PlotEtamx("DiO", aDiO, DiO0, dDiO, iDiO, nDiO) |
|
2659 if (nDiC > 0): PlotEtamx("DiC", aDiC, DiC0, dDiC, iDiC, nDiC) |
|
2660 if (nRotC > 0): PlotEtamx("RotC", aRotC, RotC0, dRotC, iRotC, nRotC) |
|
2661 if (nRetC > 0): PlotEtamx("RetC", aRetC, RetC0, dRetC, iRetC, nRetC) |
|
2662 if (nTP > 0): PlotEtamx("TP", aTP, TP0, dTP, iTP, nTP) |
|
2663 if (nTS > 0): PlotEtamx("TS", aTS, TS0, dTS, iTS, nTS) |
|
2664 if (nRP > 0): PlotEtamx("RP", aRP, RP0, dRP, iRP, nRP) |
|
2665 if (nRS > 0): PlotEtamx("RS", aRS, RS0, dRS, iRS, nRS) |
|
2666 if (nRetT > 0): PlotEtamx("RetT", aRetT, RetT0, dRetT, iRetT, nRetT) |
|
2667 if (nRetR > 0): PlotEtamx("RetR", aRetR, RetR0, dRetR, iRetR, nRetR) |
|
2668 if (nERaT > 0): PlotEtamx("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) |
|
2669 if (nERaR > 0): PlotEtamx("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) |
|
2670 if (nRotaT > 0): PlotEtamx("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) |
|
2671 if (nRotaR > 0): PlotEtamx("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) |
|
2672 if (nLDRCal > 0): PlotEtamx("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) |
|
2673 if (nTCalT > 0): PlotEtamx("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT) |
|
2674 if (nTCalR > 0): PlotEtamx("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR) |
|
2675 if (nNCal > 0): PlotEtamx("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal) |
|
2676 if (nNCal > 0): PlotEtamx("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal) |
|
2677 if (nNCal > 0): PlotEtamx("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal) |
|
2678 if (nNCal > 0): PlotEtamx("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal) |
|
2679 if (nNI > 0): PlotEtamx("SigNoiseIt", aNIt, 0, 1, iNIt, nNI) |
|
2680 if (nNI > 0): PlotEtamx("SigNoiseIr", aNIr, 0, 1, iNIr, nNI) |
|
2681 plt.show() |
|
2682 plt.close |
|
2683 |
|
2684 # Print Etax statistics |
|
2685 Etaxmin = np.amin(aEtax[1, :]) |
|
2686 Etaxmax = np.amax(aEtax[1, :]) |
|
2687 Etaxstd = np.std(aEtax[1, :]) |
|
2688 Etaxmean = np.mean(aEtax[1, :]) |
|
2689 Etaxmedian = np.median(aEtax[1, :]) |
|
2690 print("Etax , max-mean, min-mean, median, mean ± std, eta") |
|
2691 print("{0:8.5f} ±({1:8.5f},{2:8.5f}),{3:8.5f},{4:8.5f}±{5:8.5f},{6:8.5f}".format(Etax0, Etaxmax-Etax0, Etaxmin-Etax0, Etaxmedian, Etaxmean, Etaxstd, Etax0 / K0)) |
|
2692 print() |
|
2693 |
|
2694 # Calculate and print statistics for calibration factors |
|
2695 iLDR = -1 |
|
2696 LDRrangeA = np.array(LDRrange) |
|
2697 print("LDR...., LDRsim, (max-min)/2, relerr") |
|
2698 for LDRTrue in LDRrange: |
|
2699 iLDR = iLDR + 1 |
|
2700 LDRsimmin[iLDR] = np.amin(aLDRsim[iLDR, :]) |
|
2701 LDRsimmax[iLDR] = np.amax(aLDRsim[iLDR, :]) |
|
2702 # LDRsimstd = np.std(aLDRsim[iLDR, :]) |
|
2703 LDRsimmean[iLDR] = np.mean(aLDRsim[iLDR, :]) |
|
2704 # LDRsimmedian = np.median(aLDRsim[iLDR, :]) |
|
2705 print("{0:8.5f}, {1:8.5f}, {2:8.5f}, {3:8.5f}".format(LDRrangeA[iLDR], LDRsimmean[iLDR], (LDRsimmax[iLDR]-LDRsimmin[iLDR])/2, (LDRsimmax[iLDR]-LDRsimmin[iLDR])/2/LDRsimmean[iLDR])) |
|
2706 iLDR = -1 |
|
2707 print("LDR...., Etax , (max-min)/2, relerr") |
|
2708 for LDRTrue in LDRrange: |
|
2709 iLDR = iLDR + 1 |
|
2710 Etaxmin = np.amin(aEtax[iLDR, :]) |
|
2711 Etaxmax = np.amax(aEtax[iLDR, :]) |
|
2712 # Etaxstd = np.std(aEtax[iLDR, :]) |
|
2713 Etaxmean = np.mean(aEtax[iLDR, :]) |
|
2714 # Etaxmedian = np.median(aEtax[iLDR, :]) |
|
2715 print("{0:8.5f}, {1:8.5f}, {2:8.5f}, {3:8.5f}".format(LDRrangeA[iLDR], Etaxmean, (Etaxmax-Etaxmin)/2, (Etaxmax-Etaxmin)/2/Etaxmean)) |
|
2716 iLDR = -1 |
|
2717 print("LDR...., Etapx , (max-min)/2, relerr") |
|
2718 for LDRTrue in LDRrange: |
|
2719 iLDR = iLDR + 1 |
|
2720 Etapxmin = np.amin(aEtapx[iLDR, :]) |
|
2721 Etapxmax = np.amax(aEtapx[iLDR, :]) |
|
2722 # Etapxstd = np.std(aEtapx[iLDR, :]) |
|
2723 Etapxmean = np.mean(aEtapx[iLDR, :]) |
|
2724 # Etapxmedian = np.median(aEtapx[iLDR, :]) |
|
2725 print("{0:8.5f}, {1:8.5f}, {2:8.5f}, {3:8.5f}".format(LDRrangeA[iLDR], Etapxmean, (Etapxmax-Etapxmin)/2, (Etapxmax-Etapxmin)/2/Etapxmean)) |
|
2726 iLDR = -1 |
|
2727 print("LDR...., Etamx , (max-min)/2, relerr") |
|
2728 for LDRTrue in LDRrange: |
|
2729 iLDR = iLDR + 1 |
|
2730 Etamxmin = np.amin(aEtamx[iLDR, :]) |
|
2731 Etamxmax = np.amax(aEtamx[iLDR, :]) |
|
2732 # Etamxstd = np.std(aEtamx[iLDR, :]) |
|
2733 Etamxmean = np.mean(aEtamx[iLDR, :]) |
|
2734 # Etamxmedian = np.median(aEtamx[iLDR, :]) |
|
2735 print("{0:8.5f}, {1:8.5f}, {2:8.5f}, {3:8.5f}".format(LDRrangeA[iLDR], Etamxmean, (Etamxmax-Etamxmin)/2, (Etamxmax-Etamxmin)/2/Etamxmean)) |
|
2736 |
|
2737 f.close() |
|
2738 |
|
2739 |
|
2740 ''' |
|
2741 # --- Plot F11 histograms |
|
2742 print() |
|
2743 print(" ############################################################################## ") |
|
2744 print(Text1) |
|
2745 print() |
|
2746 |
|
2747 iLDR = 5 |
|
2748 for LDRTrue in LDRrange: |
|
2749 iLDR = iLDR - 1 |
|
2750 #aF11corr[iLDR,:] = aF11corr[iLDR,:] / aF11corr[0,:] - 1.0 |
|
2751 aF11corr[iLDR,:] = aF11corr[iLDR,:] / aF11sim0[iLDR] - 1.0 |
|
2752 # Plot F11 |
|
2753 def PlotSubHistF11(aVar, aX, X0, daX, iaX, naX): |
|
2754 fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) |
|
2755 iLDR = -1 |
|
2756 for LDRTrue in LDRrange: |
|
2757 iLDR = iLDR + 1 |
|
2758 |
|
2759 #F11min[iLDR] = np.min(aF11corr[iLDR,:]) |
|
2760 #F11max[iLDR] = np.max(aF11corr[iLDR,:]) |
|
2761 #Rmin = F11min[iLDR] * 0.995 # np.min(aLDRcorr[iLDR,:]) * 0.995 |
|
2762 #Rmax = F11max[iLDR] * 1.005 # np.max(aLDRcorr[iLDR,:]) * 1.005 |
|
2763 |
|
2764 #Rmin = 0.8 |
|
2765 #Rmax = 1.2 |
|
2766 |
|
2767 #plt.subplot(5,2,iLDR+1) |
|
2768 plt.subplot(1,5,iLDR+1) |
|
2769 (n, bins, patches) = plt.hist(aF11corr[iLDR,:], |
|
2770 bins=100, log=False, |
|
2771 alpha=0.5, density=False, color = '0.5', histtype='stepfilled') |
|
2772 |
|
2773 for iaX in range(-naX,naX+1): |
|
2774 plt.hist(aF11corr[iLDR,aX == iaX], |
|
2775 bins=100, log=False, alpha=0.3, density=False, histtype='stepfilled', label = str(round(X0 + iaX*daX/naX,5))) |
|
2776 |
|
2777 if (iLDR == 2): plt.legend() |
|
2778 |
|
2779 plt.tick_params(axis='both', labelsize=9) |
|
2780 #plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2) |
|
2781 |
|
2782 #plt.title(LID + ' ' + aVar, fontsize=18) |
|
2783 #plt.ylabel('frequency', fontsize=10) |
|
2784 #plt.xlabel('LDRCorr', fontsize=10) |
|
2785 #fig.tight_layout() |
|
2786 fig.suptitle(LID + ' ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar, fontsize=14, y=1.05) |
|
2787 #plt.show() |
|
2788 #fig.savefig(LID + '_' + aVar + '.png', dpi=150, bbox_inches='tight', pad_inches=0) |
|
2789 #plt.close |
|
2790 return |
|
2791 |
|
2792 if (nQin > 0): PlotSubHistF11("Qin", aQin, Qin0, dQin, iQin, nQin) |
|
2793 if (nVin > 0): PlotSubHistF11("Vin", aVin, Vin0, dVin, iVin, nVin) |
|
2794 if (nRotL > 0): PlotSubHistF11("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) |
|
2795 if (nRetE > 0): PlotSubHistF11("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) |
|
2796 if (nRotE > 0): PlotSubHistF11("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) |
|
2797 if (nDiE > 0): PlotSubHistF11("DiE", aDiE, DiE0, dDiE, iDiE, nDiE) |
|
2798 if (nRetO > 0): PlotSubHistF11("RetO", aRetO, RetO0, dRetO, iRetO, nRetO) |
|
2799 if (nRotO > 0): PlotSubHistF11("RotO", aRotO, RotO0, dRotO, iRotO, nRotO) |
|
2800 if (nDiO > 0): PlotSubHistF11("DiO", aDiO, DiO0, dDiO, iDiO, nDiO) |
|
2801 if (nDiC > 0): PlotSubHistF11("DiC", aDiC, DiC0, dDiC, iDiC, nDiC) |
|
2802 if (nRotC > 0): PlotSubHistF11("RotC", aRotC, RotC0, dRotC, iRotC, nRotC) |
|
2803 if (nRetC > 0): PlotSubHistF11("RetC", aRetC, RetC0, dRetC, iRetC, nRetC) |
|
2804 if (nTP > 0): PlotSubHistF11("TP", aTP, TP0, dTP, iTP, nTP) |
|
2805 if (nTS > 0): PlotSubHistF11("TS", aTS, TS0, dTS, iTS, nTS) |
|
2806 if (nRP > 0): PlotSubHistF11("RP", aRP, RP0, dRP, iRP, nRP) |
|
2807 if (nRS > 0): PlotSubHistF11("RS", aRS, RS0, dRS, iRS, nRS) |
|
2808 if (nRetT > 0): PlotSubHistF11("RetT", aRetT, RetT0, dRetT, iRetT, nRetT) |
|
2809 if (nRetR > 0): PlotSubHistF11("RetR", aRetR, RetR0, dRetR, iRetR, nRetR) |
|
2810 if (nERaT > 0): PlotSubHistF11("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) |
|
2811 if (nERaR > 0): PlotSubHistF11("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) |
|
2812 if (nRotaT > 0): PlotSubHistF11("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) |
|
2813 if (nRotaR > 0): PlotSubHistF11("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) |
|
2814 if (nLDRCal > 0): PlotSubHistF11("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) |
|
2815 if (nTCalT > 0): PlotSubHistF11("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT) |
|
2816 if (nTCalR > 0): PlotSubHistF11("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR) |
|
2817 if (nNCal > 0): PlotSubHistF11("CalNoise", aNCal, 0, 1/nNCal, iNCal, nNCal) |
|
2818 if (nNI > 0): PlotSubHistF11("SigNoise", aNI, 0, 1/nNI, iNI, nNI) |
|
2819 |
|
2820 |
|
2821 plt.show() |
|
2822 plt.close |
|
2823 |
|
2824 ''' |
|
2825 ''' |
|
2826 # only histogram |
|
2827 #print("******************* " + aVar + " *******************") |
|
2828 fig, ax = plt.subplots(nrows=5, ncols=2, sharex=True, sharey=True, figsize=(10, 10)) |
|
2829 iLDR = -1 |
|
2830 for LDRTrue in LDRrange: |
|
2831 iLDR = iLDR + 1 |
|
2832 LDRmin[iLDR] = np.min(aLDRcorr[iLDR,:]) |
|
2833 LDRmax[iLDR] = np.max(aLDRcorr[iLDR,:]) |
|
2834 Rmin = np.min(aLDRcorr[iLDR,:]) * 0.999 |
|
2835 Rmax = np.max(aLDRcorr[iLDR,:]) * 1.001 |
|
2836 plt.subplot(5,2,iLDR+1) |
|
2837 (n, bins, patches) = plt.hist(aLDRcorr[iLDR,:], |
|
2838 range=[Rmin, Rmax], |
|
2839 bins=200, log=False, alpha=0.2, density=False, color = '0.5', histtype='stepfilled') |
|
2840 plt.tick_params(axis='both', labelsize=9) |
|
2841 plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2) |
|
2842 plt.show() |
|
2843 plt.close |
|
2844 # --- End of Plot F11 histograms |
|
2845 ''' |
|
2846 |
|
2847 |
|
2848 ''' |
|
2849 # --- Plot K over LDRCal |
|
2850 fig3 = plt.figure() |
|
2851 plt.plot(LDRCal0+aLDRCal*dLDRCal/nLDRCal,aGHK[4,:], linewidth=2.0, color='b') |
|
2852 |
|
2853 plt.xlabel('LDRCal', fontsize=18) |
|
2854 plt.ylabel('K', fontsize=14) |
|
2855 plt.title(LID, fontsize=18) |
|
2856 plt.show() |
|
2857 plt.close |
|
2858 ''' |
|
2859 |
|
2860 # Additional plot routines ======> |
|
2861 ''' |
|
2862 #****************************************************************************** |
|
2863 # 1. Plot LDRCorrected - LDR(measured Icross/Iparallel) |
|
2864 LDRa = np.arange(1.,100.)*0.005 |
|
2865 LDRCorra = np.arange(1.,100.) |
|
2866 if Y == - 1.: LDRa = 1./LDRa |
|
2867 LDRCorra = (1./Eta*LDRa*(GT+HT)-(GR+HR))/((GR-HR)-1./Eta*LDRa*(GT-HT)) |
|
2868 if Y == - 1.: LDRa = 1./LDRa |
|
2869 # |
|
2870 #fig = plt.figure() |
|
2871 plt.plot(LDRa,LDRCorra-LDRa) |
|
2872 plt.plot([0.,0.5],[0.,0.5]) |
|
2873 plt.suptitle('LDRCorrected - LDR(measured Icross/Iparallel)', fontsize=16) |
|
2874 plt.xlabel('LDR', fontsize=18) |
|
2875 plt.ylabel('LDRCorr - LDR', fontsize=16) |
|
2876 #plt.savefig('test.png') |
|
2877 # |
|
2878 ''' |
|
2879 ''' |
|
2880 #****************************************************************************** |
|
2881 # 2. Plot LDRsim (simulated measurements without corrections = Icross/Iparallel) over LDRtrue |
|
2882 LDRa = np.arange(1.,100.)*0.005 |
|
2883 LDRsima = np.arange(1.,100.) |
|
2884 |
|
2885 atruea = (1.-LDRa)/(1+LDRa) |
|
2886 Ita = TiT*TiO*IinL*(GT+atruea*HT) |
|
2887 Ira = TiR*TiO*IinL*(GR+atruea*HR) |
|
2888 LDRsima = Ira/Ita # simulated uncorrected LDR with Y from input file |
|
2889 if Y == -1.: LDRsima = 1./LDRsima |
|
2890 # |
|
2891 #fig = plt.figure() |
|
2892 plt.plot(LDRa,LDRsima) |
|
2893 plt.plot([0.,0.5],[0.,0.5]) |
|
2894 plt.suptitle('LDRsim (simulated measurements without corrections = Icross/Iparallel) over LDRtrue', fontsize=10) |
|
2895 plt.xlabel('LDRtrue', fontsize=18) |
|
2896 plt.ylabel('LDRsim', fontsize=16) |
|
2897 #plt.savefig('test.png') |
|
2898 # |
|
2899 ''' |
|