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