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