GHK_0.9.8h_Py3.7.py

Wed, 07 Jul 2021 15:42:30 +0200

author
Volker Freudenthaler
date
Wed, 07 Jul 2021 15:42:30 +0200
changeset 49
7c14c9a085ba
permissions
-rw-r--r--

2nd order polynomial fits to
- K(LDRCal) and
- LDR-uncertainty max. and min depending on LDRtrue

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

mercurial