Volker@43: # -*- coding: utf-8 -*- Volker@43: """ Volker@43: Copyright 2016, 2019 Volker Freudenthaler Volker@43: Volker@43: Licensed under the EUPL, Version 1.1 only (the "Licence"). Volker@43: Volker@43: You may not use this work except in compliance with the Licence. Volker@43: A copy of the licence is distributed with the code. Alternatively, you may obtain Volker@43: a copy of the Licence at: Volker@43: Volker@43: https://joinup.ec.europa.eu/community/eupl/og_page/eupl Volker@43: Volker@43: Unless required by applicable law or agreed to in writing, software distributed Volker@43: under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR CONDITIONS Volker@43: OF ANY KIND, either express or implied. See the Licence for the specific language governing Volker@43: permissions and limitations under the Licence. Volker@43: Volker@43: Equation reference: http://www.atmos-meas-tech-discuss.net/amt-2015-338/amt-2015-338.pdf Volker@43: With equations code from Appendix C Volker@43: Python 3.7, seaborn 0.9.0 Volker@43: Volker@43: Code description: Volker@43: Volker@43: From measured lidar signals we cannot directly determine the desired backscatter coefficient (F11) and the linear depolarization ratio (LDR) Volker@43: because of the cross talk between the channles and systematic errors of a lidar system. Volker@43: 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@43: with which the measured signals can be corrected. Volker@43: 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@43: 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@43: 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@43: 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@43: GHK parameters (GT0, GR0, HT0, HR0, and K0) to derive finally the distributions of "possibly real" linear depolarization ratios (LDRCorr), Volker@43: which are plotted for five different input linear depolarization ratios (LDRtrue). The red bars in the plots represent the input values of LDRtrue. Volker@43: 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@43: 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@43: "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@43: 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@43: 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@43: Volker@43: 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@43: Volker@43: Ver. 0.9.7: includes the random error (signal noise) of the calibration and standard measurements Volker@43: Changes: Volker@43: Line 1687 Eta = (TaR * TiR) / (TaT * TiT) Volker@43: Line 1691 K = Etax / Eta # K of the real system; but correction in Line 1721 with K0 / Etax Volker@43: should work with nTCalT = nTCalR = 0 Volker@43: Ver. 0.9.7b: Volker@43: ToDo: include error due to TCalT und TCalR => determination of NCalT and NCalR etc. in error calculation line 1741ff Volker@43: combined error loops iNI and INCal for signals Volker@43: Ver. 0.9.7c: individual error loops for each of the six signals Volker@43: Ver. 0.9.7c2: different calculation of the signal noise errors Volker@43: Ver. 0.9.7c3: n.a.different calculation of the signal noise errors Volker@43: 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@43: Ver. 0.9.8: Volker@43: - correct calculation of Eta for cleaned anaylsers considering the combined transmission Eta = (TaT* TiT)(1 + cos2RotaT * DaT * DiT) and (TaR * TiR)(1 + cos2RotaR * DaR * DiR) according to the papers supplement Eqs. (S.10.10.1) ff Volker@43: - calculation of the PLDR from LDR and BSR, BSR, and LDRm Volker@43: - 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@43: 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@43: - the strong Tp dependence of the errors Volker@43: - the factor 2 in the GH parameters Volker@43: - see c:\technik\Optik\Polarizers\DepCal\ApplOpt\GH-parameters-190114.odt Volker@43: Ver. 0.9.8c: includes error of Etax Volker@43: Ver. 0.9.8d: Eta0, K0 etc in error loop replaced by Eta0y, K0y etc. Changes in signal noise calculations Volker@43: Ver. 0.9.8e: ambiguous laser spec. DOLP (no discrimination between left and right circular polarisation) replaced by Stokes parameters Qin, Uin Volker@43: Ver. 0.9.8e2: Added plot of LDRsim, Etax, Etapx, Etamx; LDRCorr and aLDRcorr consistently named Volker@43: 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@43: Ver. 0.9.8e4: text changed for y=+-1 (see line 274 ff and line 1044 ff Volker@43: Ver. 0.9.8e5: changed: LDRunCorr = LDRsim / Etax Volker@43: Volker@43: ======================================================== Volker@43: simulation: LDRsim = Ir / It with variable parameters (possible truths) Volker@43: G,H,Eta,Etax,K Volker@43: It = TaT * TiT * ATP1 * TiO * TiE * (GT + atrue * HT) Volker@43: LDRsim = Ir / It Volker@43: consistency test: is forward simulation and correction consistent? Volker@43: LDRCorr = (LDRsim / Eta * (GT + HT) - (GR + HR)) / ((GR - HR) - LDRsim / Eta * (GT - HT)) => atrue? Volker@43: assumed true: G0,H0,Eta0,Etax0,K0 => actual retrievals of LDRCorr Volker@43: => correct possible truths with assumed true G0,H0,Eta0 Volker@43: measure: It, Ir, EtaX Volker@43: LDRunCorr = LDRsim / Etax Volker@43: correct it with G0,H0,K0: Volker@43: LDRCorr = (LDRsim / (Etax / K0) * (GT0 + HT0) - (GR0 + HR0)) / ((GR0 - HR0) - LDRsim0 / (Etax / K0) * (GT0 - HT0)) Volker@43: """ Volker@43: # 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@43: from __future__ import print_function Volker@43: # !/usr/bin/env python3 Volker@43: Volker@43: import os Volker@43: import sys Volker@43: Volker@43: from scipy.stats import kurtosis Volker@43: from scipy.stats import skew Volker@43: # use: kurtosis(data, fisher=True,bias=False) => 0; skew(data,bias=False) => 0 Volker@43: # Comment: the seaborn library makes nicer plots, but the code works also without it. Volker@43: import numpy as np Volker@43: import matplotlib.pyplot as plt Volker@43: Volker@43: try: Volker@43: import seaborn as sns Volker@43: Volker@43: sns_loaded = True Volker@43: except ImportError: Volker@43: sns_loaded = False Volker@43: Volker@43: # from time import clock # python 2 Volker@43: from timeit import default_timer as clock Volker@43: Volker@43: # from matplotlib.backends.backend_pdf import PdfPages Volker@43: # pdffile = '{}.pdf'.format('path') Volker@43: # pp = PdfPages(pdffile) Volker@43: ## pp.savefig can be called multiple times to save to multiple pages Volker@43: # pp.savefig() Volker@43: # pp.close() Volker@43: Volker@43: from contextlib import contextmanager Volker@43: Volker@43: @contextmanager Volker@43: def redirect_stdout(new_target): Volker@43: old_target, sys.stdout = sys.stdout, new_target # replace sys.stdout Volker@43: try: Volker@43: yield new_target # run some code with the replaced stdout Volker@43: finally: Volker@43: sys.stdout.flush() Volker@43: sys.stdout = old_target # restore to the previous value Volker@43: Volker@43: ''' Volker@43: real_raw_input = vars(__builtins__).get('raw_input',input) Volker@43: ''' Volker@43: try: Volker@43: import __builtin__ Volker@43: Volker@43: input = getattr(__builtin__, 'raw_input') Volker@43: except (ImportError, AttributeError): Volker@43: pass Volker@43: Volker@43: from distutils.util import strtobool Volker@43: Volker@43: Volker@43: def user_yes_no_query(question): Volker@43: sys.stdout.write('%s [y/n]\n' % question) Volker@43: while True: Volker@43: try: Volker@43: return strtobool(input().lower()) Volker@43: except ValueError: Volker@43: sys.stdout.write('Please respond with \'y\' or \'n\'.\n') Volker@43: Volker@43: Volker@43: # if user_yes_no_query('want to exit?') == 1: sys.exit() Volker@43: Volker@43: abspath = os.path.abspath(__file__) Volker@43: dname = os.path.dirname(abspath) Volker@43: fname = os.path.basename(abspath) Volker@43: os.chdir(dname) Volker@43: Volker@43: # PrintToOutputFile = True Volker@43: Volker@43: sqr05 = 0.5 ** 0.5 Volker@43: Volker@43: # ---- Initial definition of variables; the actual values will be read in with exec(open('./optic_input.py').read()) below Volker@43: # Do you want to calculate the errors? If not, just the GHK-parameters are determined. Volker@43: Error_Calc = True Volker@43: LID = "internal" Volker@43: EID = "internal" Volker@43: # --- IL Laser IL and +-Uncertainty Volker@43: Qin, dQin, nQin = 1., 0.0, 0 # second Stokes vector parameter; default 1 => linear polarization Volker@43: Vin, dVin, nVin = 0., 0.0, 0 # fourth Stokes vector parameter Volker@43: RotL, dRotL, nRotL = 0.0, 0.0, 1 # alpha; rotation of laser polarization in degrees; default 0 Volker@43: # IL = 1e5 #photons in the laser beam, including detection efficiency of the telescope, atmodspheric and r^2 attenuation Volker@43: # --- ME Emitter and +-Uncertainty Volker@43: DiE, dDiE, nDiE = 0., 0.00, 1 # Diattenuation Volker@43: TiE = 1. # Unpolarized transmittance Volker@43: RetE, dRetE, nRetE = 0., 180.0, 0 # Retardance in degrees Volker@43: RotE, dRotE, nRotE = 0., 0.0, 0 # beta: Rotation of optical element in degrees Volker@43: # --- MO Receiver Optics including telescope Volker@43: DiO, dDiO, nDiO = -0.055, 0.003, 1 Volker@43: TiO = 0.9 Volker@43: RetO, dRetO, nRetO = 0., 180.0, 2 Volker@43: RotO, dRotO, nRotO = 0., 0.1, 1 # gamma Volker@43: # --- PBS MT transmitting path defined with (TS,TP); and +-Uncertainty Volker@43: TP, dTP, nTP = 0.98, 0.02, 1 Volker@43: TS, dTS, nTS = 0.001, 0.001, 1 Volker@43: TiT = 0.5 * (TP + TS) Volker@43: DiT = (TP - TS) / (TP + TS) Volker@43: # PolFilter Volker@43: RetT, dRetT, nRetT = 0., 180., 0 Volker@43: ERaT, dERaT, nERaT = 0.001, 0.001, 1 Volker@43: RotaT, dRotaT, nRotaT = 0., 3., 1 Volker@43: DaT = (1 - ERaT) / (1 + ERaT) Volker@43: TaT = 0.5 * (1 + ERaT) Volker@43: # --- PBS MR reflecting path defined with (RS,RP); and +-Uncertainty Volker@43: RS_RP_depend_on_TS_TP = False Volker@43: if (RS_RP_depend_on_TS_TP): Volker@43: RP, dRP, nRP = 1 - TP, 0.0, 0 Volker@43: RS, dRS, nRS = 1 - TS, 0.0, 0 Volker@43: else: Volker@43: RP, dRP, nRP = 0.05, 0.01, 1 Volker@43: RS, dRS, nRS = 0.98, 0.01, 1 Volker@43: TiR = 0.5 * (RP + RS) Volker@43: DiR = (RP - RS) / (RP + RS) Volker@43: # PolFilter Volker@43: RetR, dRetR, nRetR = 0., 180., 0 Volker@43: ERaR, dERaR, nERaR = 0.001, 0.001, 1 Volker@43: RotaR, dRotaR, nRotaR = 90., 3., 1 Volker@43: DaR = (1 - ERaR) / (1 + ERaR) Volker@43: TaR = 0.5 * (1 + ERaR) Volker@43: Volker@43: # +++ Orientation of the PBS with respect to the reference plane (see Polarisation-orientation.png and Polarisation-orientation-2.png in /system_settings) Volker@43: # Y = +1: PBS incidence plane is parallel to reference plane and polarisation in reference plane is finally transmitted. Volker@43: # Y = -1: PBS incidence plane is perpendicular to reference plane and polarisation in reference plane is finally reflected. Volker@43: Y = 1. Volker@43: Volker@43: # Calibrator = type defined by matrix values Volker@43: LocC = 4 # location of calibrator: behind laser = 1; behind emitter = 2; before receiver = 3; before PBS = 4 Volker@43: Volker@43: # --- Additional attenuation (transmission of the ND-filter) during the calibration Volker@43: TCalT, dTCalT, nTCalT = 1, 0., 0 # transmitting path; error calc not working yet Volker@43: TCalR, dTCalR, nTCalR = 1, 0., 0 # reflecting path; error calc not working yet Volker@43: Volker@43: # *** signal noise error calculation Volker@43: # --- number of photon counts in the signal summed up in the calibration range during the calibration measurements Volker@43: NCalT = 1e6 # default 1e6, assumed the same in +45° and -45° signals Volker@43: NCalR = 1e6 # default 1e6, assumed the same in +45° and -45° signals Volker@43: NILfac = 1.0 # duration of standard (0°) measurement relative to calibration measurements Volker@43: nNCal = 0 # error nNCal: one-sigma in steps to left and right for calibration signals Volker@43: nNI = 0 # error nNI: one-sigma in steps to left and right for 0° signals Volker@43: NI = 50000 #number of photon counts in the parallel 0°-signal Volker@43: eFacT = 1.0 # rel. amplification of transmitted channel, approximate values are sufficient; def. = 1 Volker@43: eFacR = 10.0 Volker@43: IoutTp0, IoutTp, dIoutTp0 = 0.5, 0.5, 0.0 Volker@43: IoutTm0, IoutTm, dIoutTm0 = 0.5, 0.5, 0.0 Volker@43: IoutRp0, IoutRp, dIoutRp0 = 0.5, 0.5, 0.0 Volker@43: IoutRm0, IoutRm, dIoutRm0 = 0.5, 0.5, 0.0 Volker@43: It0, It, dIt0 = 1 , 1, 0 Volker@43: Ir0, Ir, dTr0 = 1 , 1, 0 Volker@43: CalcFrom0deg = True Volker@43: Volker@43: TypeC = 3 # linear polarizer calibrator Volker@43: # example with extinction ratio 0.001 Volker@43: DiC, dDiC, nDiC = 1.0, 0., 0 # ideal 1.0 Volker@43: TiC = 0.5 # ideal 0.5 Volker@43: RetC, dRetC, nRetC = 0.0, 0.0, 0 Volker@43: RotC, dRotC, nRotC = 0.0, 0.1, 0 # constant calibrator offset epsilon Volker@43: RotationErrorEpsilonForNormalMeasurements = False # is in general False for TypeC == 3 calibrator Volker@43: Volker@43: # Rotation error without calibrator: if False, then epsilon = 0 for normal measurements Volker@43: RotationErrorEpsilonForNormalMeasurements = True Volker@43: # BSR backscatter ratio Volker@43: # BSR, dBSR, nBSR = 10, 0.05, 1 Volker@43: BSR = np.zeros(5) Volker@43: BSR = [1.1, 2, 5, 10., 50.] Volker@43: # theoretical molecular LDR LDRm Volker@43: LDRm, dLDRm, nLDRm = 0.004, 0.001, 1 Volker@43: # LDRCal assumed atmospheric linear depolarization ratio during the calibration measurements (first guess) Volker@43: LDRCal0, dLDRCal, nLDRCal = 0.25, 0.04, 1 Volker@43: LDRCal = LDRCal0 Volker@43: # measured LDRm will be corrected with calculated parameters Volker@43: LDRmeas = 0.015 Volker@43: # LDRtrue for simulation of measurement => LDRsim Volker@43: LDRtrue = 0.004 Volker@43: LDRtrue2 = 0.004 Volker@43: LDRunCorr = 1. Volker@43: # Initialize other values to 0 Volker@43: ER, nER, dER = 0.001, 0, 0.001 Volker@43: K = 0. Volker@43: Km = 0. Volker@43: Kp = 0. Volker@43: LDRCorr = 0. Volker@43: Eta = 0. Volker@43: Ir = 0. Volker@43: It = 0. Volker@43: h = 1. Volker@43: Volker@43: Loc = ['', 'behind laser', 'behind emitter', 'before receiver', 'before PBS'] Volker@43: Type = ['', 'mechanical rotator', 'hwp rotator', 'linear polarizer', 'qwp rotator', 'circular polarizer', Volker@43: 'real HWP +-22.5°'] Volker@43: Volker@43: bPlotEtax = False Volker@43: Volker@43: # end of initial definition of variables Volker@43: # ******************************************************************************************************************************* Volker@43: # --- Read actual lidar system parameters from optic_input.py (must be in the programs sub-directory 'system_settings') Volker@43: # ******************************************************************************************************************************* Volker@43: Volker@43: # InputFile = 'optic_input_0.9.8e4-PollyXT_Lacros.py' Volker@43: InputFile = 'optic_input_example_lidar_ver0.9.8e.py' Volker@43: Volker@43: # ******************************************************************************************************************************* Volker@43: Volker@43: ''' Volker@43: print("From ", dname) Volker@43: print("Running ", fname) Volker@43: print("Reading input file ", InputFile, " for") Volker@43: ''' Volker@43: input_path = os.path.join('.', 'system_settings', InputFile) Volker@43: # this works with Python 2 and 3! Volker@43: exec(open(input_path).read(), globals()) Volker@43: # end of read actual system parameters Volker@43: Volker@43: Volker@43: # --- Manual Parameter Change --- Volker@43: # (use for quick parameter changes without changing the input file ) Volker@43: # DiO = 0. Volker@43: # LDRtrue = 0.45 Volker@43: # LDRtrue2 = 0.004 Volker@43: # Y = -1 Volker@43: # LocC = 4 #location of calibrator: 1 = behind laser; 2 = behind emitter; 3 = before receiver; 4 = before PBS Volker@43: # #TypeC = 6 Don't change the TypeC here Volker@43: # RotationErrorEpsilonForNormalMeasurements = True Volker@43: # LDRCal = 0.25 Volker@43: # # --- Errors Volker@43: Qin0, dQin, nQin = Qin, dQin, nQin Volker@43: Vin0, dVin, nVin = Vin, dVin, nVin Volker@43: RotL0, dRotL, nRotL = RotL, dRotL, nRotL Volker@43: Volker@43: DiE0, dDiE, nDiE = DiE, dDiE, nDiE Volker@43: RetE0, dRetE, nRetE = RetE, dRetE, nRetE Volker@43: RotE0, dRotE, nRotE = RotE, dRotE, nRotE Volker@43: Volker@43: DiO0, dDiO, nDiO = DiO, dDiO, nDiO Volker@43: RetO0, dRetO, nRetO = RetO, dRetO, nRetO Volker@43: RotO0, dRotO, nRotO = RotO, dRotO, nRotO Volker@43: Volker@43: DiC0, dDiC, nDiC = DiC, dDiC, nDiC Volker@43: RetC0, dRetC, nRetC = RetC, dRetC, nRetC Volker@43: RotC0, dRotC, nRotC = RotC, dRotC, nRotC Volker@43: Volker@43: TP0, dTP, nTP = TP, dTP, nTP Volker@43: TS0, dTS, nTS = TS, dTS, nTS Volker@43: RetT0, dRetT, nRetT = RetT, dRetT, nRetT Volker@43: Volker@43: ERaT0, dERaT, nERaT = ERaT, dERaT, nERaT Volker@43: RotaT0, dRotaT, nRotaT = RotaT, dRotaT, nRotaT Volker@43: Volker@43: RP0, dRP, nRP = RP, dRP, nRP Volker@43: RS0, dRS, nRS = RS, dRS, nRS Volker@43: RetR0, dRetR, nRetR = RetR, dRetR, nRetR Volker@43: Volker@43: ERaR0, dERaR, nERaR = ERaR, dERaR, nERaR Volker@43: RotaR0, dRotaR, nRotaR = RotaR, dRotaR, nRotaR Volker@43: Volker@43: LDRCal0, dLDRCal, nLDRCal = LDRCal, dLDRCal, nLDRCal Volker@43: Volker@43: # BSR0, dBSR, nBSR = BSR, dBSR, nBSR Volker@43: LDRm0, dLDRm, nLDRm = LDRm, dLDRm, nLDRm Volker@43: # ---------- End of manual parameter change Volker@43: Volker@43: RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC = RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0 Volker@43: TP, TS, RP, RS, ERaT, RotaT, RetT, ERaR, RotaR, RetR = TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0 Volker@43: LDRCal = LDRCal0 Volker@43: DTa0, TTa0, DRa0, TRa0, LDRsimx, LDRCorr = 0., 0., 0., 0., 0., 0. Volker@43: TCalT0, TCalR0 = TCalT, TCalR Volker@43: Volker@43: TiT = 0.5 * (TP + TS) Volker@43: DiT = (TP - TS) / (TP + TS) Volker@43: ZiT = (1. - DiT ** 2) ** 0.5 Volker@43: TiR = 0.5 * (RP + RS) Volker@43: DiR = (RP - RS) / (RP + RS) Volker@43: ZiR = (1. - DiR ** 2) ** 0.5 Volker@43: Volker@43: C2aT = np.cos(np.deg2rad(2. * RotaT)) Volker@43: C2aR = np.cos(np.deg2rad(2. * RotaR)) Volker@43: ATPT = float(1. + C2aT * DaT * DiT) Volker@43: ARPT = float(1. + C2aR * DaR * DiR) Volker@43: TTa = TiT * TaT * ATPT # unpolarized transmission Volker@43: TRa = TiR * TaR * ARPT # unpolarized transmission Volker@43: Eta0 = TRa / TTa Volker@43: Volker@43: # --- alternative texts for output Volker@43: dY = ['perpendicular', '', 'parallel'] Volker@43: dY2 = ['reflected', '', 'transmitted'] Volker@43: if ((abs(RotL) < 45 and Y == 1) or (abs(RotL) >= 45 and Y == -1)): Volker@44: dY3 = "the parallel laser polarisation is detected in the transmitted channel." Volker@43: else: Volker@44: dY3 = "the parallel laser polarisation is detected in the reflected channel." Volker@43: Volker@43: # --- check input errors Volker@43: if ((Qin ** 2 + Vin ** 2) ** 0.5) > 1: Volker@43: print("Error: degree of polarisation of laser > 1. Check Qin and Vin! ") Volker@43: sys.exit() Volker@43: Volker@43: # --- this subroutine is for the calculation of the PLDR from LDR, BSR, and LDRm ------------------- Volker@43: def CalcPLDR(LDR, BSR, LDRm): Volker@43: PLDR = (BSR * (1. + LDRm) * LDR - LDRm * (1. + LDR)) / (BSR * (1. + LDRm) - (1. + LDR)) Volker@43: return (PLDR) Volker@43: # --- this subroutine is for the calculation with certain fixed parameters ------------------------ Volker@43: def Calc(TCalT, TCalR, NCalT, NCalR, Qin, Vin, RotL, RotE, RetE, DiE, RotO, RetO, DiO, Volker@43: RotC, RetC, DiC, TP, TS, RP, RS, Volker@43: ERaT, RotaT, RetT, ERaR, RotaR, RetR, LDRCal): Volker@43: # ---- Do the calculations of bra-ket vectors Volker@43: h = -1. if TypeC == 2 else 1 Volker@43: # from input file: assumed LDRCal for calibration measurements Volker@43: aCal = (1. - LDRCal) / (1. + LDRCal) Volker@43: atrue = (1. - LDRtrue) / (1. + LDRtrue) Volker@43: Volker@43: # angles of emitter and laser and calibrator and receiver optics Volker@43: # RotL = alpha, RotE = beta, RotO = gamma, RotC = epsilon Volker@43: S2a = np.sin(2 * np.deg2rad(RotL)) Volker@43: C2a = np.cos(2 * np.deg2rad(RotL)) Volker@43: S2b = np.sin(2 * np.deg2rad(RotE)) Volker@43: C2b = np.cos(2 * np.deg2rad(RotE)) Volker@43: S2ab = np.sin(np.deg2rad(2 * RotL - 2 * RotE)) Volker@43: C2ab = np.cos(np.deg2rad(2 * RotL - 2 * RotE)) Volker@43: S2g = np.sin(np.deg2rad(2 * RotO)) Volker@43: C2g = np.cos(np.deg2rad(2 * RotO)) Volker@43: Volker@43: # Laser with Degree of linear polarization DOLP Volker@43: IinL = 1. Volker@43: QinL = Qin Volker@43: UinL = 0. Volker@43: VinL = Vin Volker@43: # VinL = (1. - DOLP ** 2) ** 0.5 Volker@43: Volker@43: # Stokes Input Vector rotation Eq. E.4 Volker@43: A = C2a * QinL - S2a * UinL Volker@43: B = S2a * QinL + C2a * UinL Volker@43: # Stokes Input Vector rotation Eq. E.9 Volker@43: C = C2ab * QinL - S2ab * UinL Volker@43: D = S2ab * QinL + C2ab * UinL Volker@43: Volker@43: # emitter optics Volker@43: CosE = np.cos(np.deg2rad(RetE)) Volker@43: SinE = np.sin(np.deg2rad(RetE)) Volker@43: ZiE = (1. - DiE ** 2) ** 0.5 Volker@43: WiE = (1. - ZiE * CosE) Volker@43: Volker@43: # Stokes Input Vector after emitter optics equivalent to Eq. E.9 with already rotated input vector from Eq. E.4 Volker@43: # b = beta Volker@43: IinE = (IinL + DiE * C) Volker@43: QinE = (C2b * DiE * IinL + A + S2b * (WiE * D - ZiE * SinE * VinL)) Volker@43: UinE = (S2b * DiE * IinL + B - C2b * (WiE * D - ZiE * SinE * VinL)) Volker@43: VinE = (-ZiE * SinE * D + ZiE * CosE * VinL) Volker@43: Volker@43: # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) Volker@43: IinF = IinE Volker@43: QinF = aCal * QinE Volker@43: UinF = -aCal * UinE Volker@43: VinF = (1. - 2. * aCal) * VinE Volker@43: Volker@43: # receiver optics Volker@43: CosO = np.cos(np.deg2rad(RetO)) Volker@43: SinO = np.sin(np.deg2rad(RetO)) Volker@43: ZiO = (1. - DiO ** 2) ** 0.5 Volker@43: WiO = (1. - ZiO * CosO) Volker@43: Volker@43: # calibrator Volker@43: CosC = np.cos(np.deg2rad(RetC)) Volker@43: SinC = np.sin(np.deg2rad(RetC)) Volker@43: ZiC = (1. - DiC ** 2) ** 0.5 Volker@43: WiC = (1. - ZiC * CosC) Volker@43: Volker@43: # Stokes Input Vector before the polarising beam splitter Eq. E.31 Volker@43: A = C2g * QinE - S2g * UinE Volker@43: B = S2g * QinE + C2g * UinE Volker@43: Volker@43: IinP = (IinE + DiO * aCal * A) Volker@43: QinP = (C2g * DiO * IinE + aCal * QinE - S2g * (WiO * aCal * B + ZiO * SinO * (1. - 2. * aCal) * VinE)) Volker@43: UinP = (S2g * DiO * IinE - aCal * UinE + C2g * (WiO * aCal * B + ZiO * SinO * (1. - 2. * aCal) * VinE)) Volker@43: VinP = (ZiO * SinO * aCal * B + ZiO * CosO * (1. - 2. * aCal) * VinE) Volker@43: Volker@43: # ------------------------- Volker@43: # F11 assuemd to be = 1 => measured: F11m = IinP / IinE with atrue Volker@43: # F11sim = TiO*(IinE + DiO*atrue*A)/IinE Volker@43: # ------------------------- Volker@43: Volker@43: # analyser Volker@43: if (RS_RP_depend_on_TS_TP): Volker@43: RS = 1. - TS Volker@43: RP = 1. - TP Volker@43: Volker@43: TiT = 0.5 * (TP + TS) Volker@43: DiT = (TP - TS) / (TP + TS) Volker@43: ZiT = (1. - DiT ** 2) ** 0.5 Volker@43: TiR = 0.5 * (RP + RS) Volker@43: DiR = (RP - RS) / (RP + RS) Volker@43: ZiR = (1. - DiR ** 2) ** 0.5 Volker@43: CosT = np.cos(np.deg2rad(RetT)) Volker@43: SinT = np.sin(np.deg2rad(RetT)) Volker@43: CosR = np.cos(np.deg2rad(RetR)) Volker@43: SinR = np.sin(np.deg2rad(RetR)) Volker@43: Volker@43: DaT = (1. - ERaT) / (1. + ERaT) Volker@43: DaR = (1. - ERaR) / (1. + ERaR) Volker@43: TaT = 0.5 * (1. + ERaT) Volker@43: TaR = 0.5 * (1. + ERaR) Volker@43: Volker@43: S2aT = np.sin(np.deg2rad(h * 2 * RotaT)) Volker@43: C2aT = np.cos(np.deg2rad(2 * RotaT)) Volker@43: S2aR = np.sin(np.deg2rad(h * 2 * RotaR)) Volker@43: C2aR = np.cos(np.deg2rad(2 * RotaR)) Volker@43: Volker@43: # Analyzer As before the PBS Eq. D.5; combined PBS and cleaning pol-filter Volker@43: ATPT = (1. + C2aT * DaT * DiT) # unpolarized transmission correction Volker@43: TTa = TiT * TaT * ATPT # unpolarized transmission Volker@43: ATP1 = 1. Volker@43: ATP2 = Y * (DiT + C2aT * DaT) / ATPT Volker@43: ATP3 = Y * S2aT * DaT * ZiT * CosT / ATPT Volker@43: ATP4 = S2aT * DaT * ZiT * SinT / ATPT Volker@43: ATP = np.array([ATP1, ATP2, ATP3, ATP4]) Volker@43: DTa = ATP2 * Y Volker@43: Volker@43: ARPT = (1 + C2aR * DaR * DiR) # unpolarized transmission correction Volker@43: TRa = TiR * TaR * ARPT # unpolarized transmission Volker@43: ARP1 = 1 Volker@43: ARP2 = Y * (DiR + C2aR * DaR) / ARPT Volker@43: ARP3 = Y * S2aR * DaR * ZiR * CosR / ARPT Volker@43: ARP4 = S2aR * DaR * ZiR * SinR / ARPT Volker@43: ARP = np.array([ARP1, ARP2, ARP3, ARP4]) Volker@43: DRa = ARP2 * Y Volker@43: Volker@43: Volker@43: # ---- Calculate signals and correction parameters for diffeent locations and calibrators Volker@43: if LocC == 4: # Calibrator before the PBS Volker@43: # print("Calibrator location not implemented yet") Volker@43: Volker@43: # S2ge = np.sin(np.deg2rad(2*RotO + h*2*RotC)) Volker@43: # C2ge = np.cos(np.deg2rad(2*RotO + h*2*RotC)) Volker@43: S2e = np.sin(np.deg2rad(h * 2 * RotC)) Volker@43: C2e = np.cos(np.deg2rad(2 * RotC)) Volker@43: # rotated AinP by epsilon Eq. C.3 Volker@43: ATP2e = C2e * ATP2 + S2e * ATP3 Volker@43: ATP3e = C2e * ATP3 - S2e * ATP2 Volker@43: ARP2e = C2e * ARP2 + S2e * ARP3 Volker@43: ARP3e = C2e * ARP3 - S2e * ARP2 Volker@43: ATPe = np.array([ATP1, ATP2e, ATP3e, ATP4]) Volker@43: ARPe = np.array([ARP1, ARP2e, ARP3e, ARP4]) Volker@43: # Stokes Input Vector before the polarising beam splitter Eq. E.31 Volker@43: A = C2g * QinE - S2g * UinE Volker@43: B = S2g * QinE + C2g * UinE Volker@43: # C = (WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE) Volker@43: Co = ZiO * SinO * VinE Volker@43: Ca = (WiO * B - 2 * ZiO * SinO * VinE) Volker@43: # C = Co + aCal*Ca Volker@43: # IinP = (IinE + DiO*aCal*A) Volker@43: # QinP = (C2g*DiO*IinE + aCal*QinE - S2g*C) Volker@43: # UinP = (S2g*DiO*IinE - aCal*UinE + C2g*C) Volker@43: # VinP = (ZiO*SinO*aCal*B + ZiO*CosO*(1-2*aCal)*VinE) Volker@43: IinPo = IinE Volker@43: QinPo = (C2g * DiO * IinE - S2g * Co) Volker@43: UinPo = (S2g * DiO * IinE + C2g * Co) Volker@43: VinPo = ZiO * CosO * VinE Volker@43: Volker@43: IinPa = DiO * A Volker@43: QinPa = QinE - S2g * Ca Volker@43: UinPa = -UinE + C2g * Ca Volker@43: VinPa = ZiO * (SinO * B - 2 * CosO * VinE) Volker@43: Volker@43: IinP = IinPo + aCal * IinPa Volker@43: QinP = QinPo + aCal * QinPa Volker@43: UinP = UinPo + aCal * UinPa Volker@43: VinP = VinPo + aCal * VinPa Volker@43: # Stokes Input Vector before the polarising beam splitter rotated by epsilon Eq. C.3 Volker@43: # QinPe = C2e*QinP + S2e*UinP Volker@43: # UinPe = C2e*UinP - S2e*QinP Volker@43: QinPoe = C2e * QinPo + S2e * UinPo Volker@43: UinPoe = C2e * UinPo - S2e * QinPo Volker@43: QinPae = C2e * QinPa + S2e * UinPa Volker@43: UinPae = C2e * UinPa - S2e * QinPa Volker@43: QinPe = C2e * QinP + S2e * UinP Volker@43: UinPe = C2e * UinP - S2e * QinP Volker@43: Volker@43: # Calibration signals and Calibration correction K from measurements with LDRCal / aCal Volker@43: if (TypeC == 2) or (TypeC == 1): # rotator calibration Eq. C.4 Volker@43: # parameters for calibration with aCal Volker@43: AT = ATP1 * IinP + h * ATP4 * VinP Volker@43: BT = ATP3e * QinP - h * ATP2e * UinP Volker@43: AR = ARP1 * IinP + h * ARP4 * VinP Volker@43: BR = ARP3e * QinP - h * ARP2e * UinP Volker@43: # Correction parameters for normal measurements; they are independent of LDR Volker@43: if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out Volker@43: IS1 = np.array([IinPo, QinPo, UinPo, VinPo]) Volker@43: IS2 = np.array([IinPa, QinPa, UinPa, VinPa]) Volker@43: GT = np.dot(ATP, IS1) Volker@43: GR = np.dot(ARP, IS1) Volker@43: HT = np.dot(ATP, IS2) Volker@43: HR = np.dot(ARP, IS2) Volker@43: else: Volker@43: IS1 = np.array([IinPo, QinPo, UinPo, VinPo]) Volker@43: IS2 = np.array([IinPa, QinPa, UinPa, VinPa]) Volker@43: GT = np.dot(ATPe, IS1) Volker@43: GR = np.dot(ARPe, IS1) Volker@43: HT = np.dot(ATPe, IS2) Volker@43: HR = np.dot(ARPe, IS2) Volker@43: elif (TypeC == 3) or (TypeC == 4): # linear polariser calibration Eq. C.5 Volker@43: # parameters for calibration with aCal Volker@43: AT = ATP1 * IinP + ATP3e * UinPe + ZiC * CosC * (ATP2e * QinPe + ATP4 * VinP) Volker@43: BT = DiC * (ATP1 * UinPe + ATP3e * IinP) - ZiC * SinC * (ATP2e * VinP - ATP4 * QinPe) Volker@43: AR = ARP1 * IinP + ARP3e * UinPe + ZiC * CosC * (ARP2e * QinPe + ARP4 * VinP) Volker@43: BR = DiC * (ARP1 * UinPe + ARP3e * IinP) - ZiC * SinC * (ARP2e * VinP - ARP4 * QinPe) Volker@43: # Correction parameters for normal measurements; they are independent of LDR Volker@43: if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out Volker@43: IS1 = np.array([IinPo, QinPo, UinPo, VinPo]) Volker@43: IS2 = np.array([IinPa, QinPa, UinPa, VinPa]) Volker@43: GT = np.dot(ATP, IS1) Volker@43: GR = np.dot(ARP, IS1) Volker@43: HT = np.dot(ATP, IS2) Volker@43: HR = np.dot(ARP, IS2) Volker@43: else: Volker@43: IS1e = np.array([IinPo + DiC * QinPoe, DiC * IinPo + QinPoe, ZiC * (CosC * UinPoe + SinC * VinPo), Volker@43: -ZiC * (SinC * UinPoe - CosC * VinPo)]) Volker@43: IS2e = np.array([IinPa + DiC * QinPae, DiC * IinPa + QinPae, ZiC * (CosC * UinPae + SinC * VinPa), Volker@43: -ZiC * (SinC * UinPae - CosC * VinPa)]) Volker@43: GT = np.dot(ATPe, IS1e) Volker@43: GR = np.dot(ARPe, IS1e) Volker@43: HT = np.dot(ATPe, IS2e) Volker@43: HR = np.dot(ARPe, IS2e) Volker@43: elif (TypeC == 6): # diattenuator calibration +-22.5° rotated_diattenuator_X22x5deg.odt Volker@43: # parameters for calibration with aCal Volker@43: AT = ATP1 * IinP + sqr05 * DiC * (ATP1 * QinPe + ATP2e * IinP) + (1. - 0.5 * WiC) * ( Volker@43: ATP2e * QinPe + ATP3e * UinPe) + ZiC * (sqr05 * SinC * (ATP3e * VinP - ATP4 * UinPe) + ATP4 * CosC * VinP) Volker@43: BT = sqr05 * DiC * (ATP1 * UinPe + ATP3e * IinP) + 0.5 * WiC * ( Volker@43: ATP2e * UinPe + ATP3e * QinPe) - sqr05 * ZiC * SinC * (ATP2e * VinP - ATP4 * QinPe) Volker@43: AR = ARP1 * IinP + sqr05 * DiC * (ARP1 * QinPe + ARP2e * IinP) + (1. - 0.5 * WiC) * ( Volker@43: ARP2e * QinPe + ARP3e * UinPe) + ZiC * (sqr05 * SinC * (ARP3e * VinP - ARP4 * UinPe) + ARP4 * CosC * VinP) Volker@43: BR = sqr05 * DiC * (ARP1 * UinPe + ARP3e * IinP) + 0.5 * WiC * ( Volker@43: ARP2e * UinPe + ARP3e * QinPe) - sqr05 * ZiC * SinC * (ARP2e * VinP - ARP4 * QinPe) Volker@43: # Correction parameters for normal measurements; they are independent of LDR Volker@43: if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out Volker@43: IS1 = np.array([IinPo, QinPo, UinPo, VinPo]) Volker@43: IS2 = np.array([IinPa, QinPa, UinPa, VinPa]) Volker@43: GT = np.dot(ATP, IS1) Volker@43: GR = np.dot(ARP, IS1) Volker@43: HT = np.dot(ATP, IS2) Volker@43: HR = np.dot(ARP, IS2) Volker@43: else: Volker@43: IS1e = np.array([IinPo + DiC * QinPoe, DiC * IinPo + QinPoe, ZiC * (CosC * UinPoe + SinC * VinPo), Volker@43: -ZiC * (SinC * UinPoe - CosC * VinPo)]) Volker@43: IS2e = np.array([IinPa + DiC * QinPae, DiC * IinPa + QinPae, ZiC * (CosC * UinPae + SinC * VinPa), Volker@43: -ZiC * (SinC * UinPae - CosC * VinPa)]) Volker@43: GT = np.dot(ATPe, IS1e) Volker@43: GR = np.dot(ARPe, IS1e) Volker@43: HT = np.dot(ATPe, IS2e) Volker@43: HR = np.dot(ARPe, IS2e) Volker@43: else: Volker@43: print("Calibrator not implemented yet") Volker@43: sys.exit() Volker@43: Volker@43: elif LocC == 3: # C before receiver optics Eq.57 Volker@43: Volker@43: # S2ge = np.sin(np.deg2rad(2*RotO - 2*RotC)) Volker@43: # C2ge = np.cos(np.deg2rad(2*RotO - 2*RotC)) Volker@43: S2e = np.sin(np.deg2rad(2. * RotC)) Volker@43: C2e = np.cos(np.deg2rad(2. * RotC)) Volker@43: Volker@43: # As with C before the receiver optics (rotated_diattenuator_X22x5deg.odt) Volker@43: AF1 = np.array([1., C2g * DiO, S2g * DiO, 0.]) Volker@43: AF2 = np.array([C2g * DiO, 1. - S2g ** 2 * WiO, S2g * C2g * WiO, -S2g * ZiO * SinO]) Volker@43: AF3 = np.array([S2g * DiO, S2g * C2g * WiO, 1. - C2g ** 2 * WiO, C2g * ZiO * SinO]) Volker@43: AF4 = np.array([0., S2g * SinO, -C2g * SinO, CosO]) Volker@43: Volker@43: ATF = (ATP1 * AF1 + ATP2 * AF2 + ATP3 * AF3 + ATP4 * AF4) Volker@43: ARF = (ARP1 * AF1 + ARP2 * AF2 + ARP3 * AF3 + ARP4 * AF4) Volker@43: ATF2 = ATF[1] Volker@43: ATF3 = ATF[2] Volker@43: ARF2 = ARF[1] Volker@43: ARF3 = ARF[2] Volker@43: Volker@43: # rotated AinF by epsilon Volker@43: ATF1 = ATF[0] Volker@43: ATF4 = ATF[3] Volker@43: ATF2e = C2e * ATF[1] + S2e * ATF[2] Volker@43: ATF3e = C2e * ATF[2] - S2e * ATF[1] Volker@43: ARF1 = ARF[0] Volker@43: ARF4 = ARF[3] Volker@43: ARF2e = C2e * ARF[1] + S2e * ARF[2] Volker@43: ARF3e = C2e * ARF[2] - S2e * ARF[1] Volker@43: Volker@43: ATFe = np.array([ATF1, ATF2e, ATF3e, ATF4]) Volker@43: ARFe = np.array([ARF1, ARF2e, ARF3e, ARF4]) Volker@43: Volker@43: QinEe = C2e * QinE + S2e * UinE Volker@43: UinEe = C2e * UinE - S2e * QinE Volker@43: Volker@43: # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) Volker@43: IinF = IinE Volker@43: QinF = aCal * QinE Volker@43: UinF = -aCal * UinE Volker@43: VinF = (1. - 2. * aCal) * VinE Volker@43: Volker@43: IinFo = IinE Volker@43: QinFo = 0. Volker@43: UinFo = 0. Volker@43: VinFo = VinE Volker@43: Volker@43: IinFa = 0. Volker@43: QinFa = QinE Volker@43: UinFa = -UinE Volker@43: VinFa = -2. * VinE Volker@43: Volker@43: # Stokes Input Vector before receiver optics rotated by epsilon Eq. C.3 Volker@43: QinFe = C2e * QinF + S2e * UinF Volker@43: UinFe = C2e * UinF - S2e * QinF Volker@43: QinFoe = C2e * QinFo + S2e * UinFo Volker@43: UinFoe = C2e * UinFo - S2e * QinFo Volker@43: QinFae = C2e * QinFa + S2e * UinFa Volker@43: UinFae = C2e * UinFa - S2e * QinFa Volker@43: Volker@43: # Calibration signals and Calibration correction K from measurements with LDRCal / aCal Volker@43: if (TypeC == 2) or (TypeC == 1): # rotator calibration Eq. C.4 Volker@43: # parameters for calibration with aCal Volker@43: AT = ATF1 * IinF + ATF4 * h * VinF Volker@43: BT = ATF3e * QinF - ATF2e * h * UinF Volker@43: AR = ARF1 * IinF + ARF4 * h * VinF Volker@43: BR = ARF3e * QinF - ARF2e * h * UinF Volker@43: # Correction parameters for normal measurements; they are independent of LDR Volker@43: if (not RotationErrorEpsilonForNormalMeasurements): Volker@43: GT = ATF1 * IinE + ATF4 * VinE Volker@43: GR = ARF1 * IinE + ARF4 * VinE Volker@43: HT = ATF2 * QinE - ATF3 * UinE - ATF4 * 2 * VinE Volker@43: HR = ARF2 * QinE - ARF3 * UinE - ARF4 * 2 * VinE Volker@43: else: Volker@43: GT = ATF1 * IinE + ATF4 * h * VinE Volker@43: GR = ARF1 * IinE + ARF4 * h * VinE Volker@43: HT = ATF2e * QinE - ATF3e * h * UinE - ATF4 * h * 2 * VinE Volker@43: HR = ARF2e * QinE - ARF3e * h * UinE - ARF4 * h * 2 * VinE Volker@43: elif (TypeC == 3) or (TypeC == 4): # linear polariser calibration Eq. C.5 Volker@43: # p = +45°, m = -45° Volker@43: IF1e = np.array([IinF, ZiC * CosC * QinFe, UinFe, ZiC * CosC * VinF]) Volker@43: IF2e = np.array([DiC * UinFe, -ZiC * SinC * VinF, DiC * IinF, ZiC * SinC * QinFe]) Volker@43: AT = np.dot(ATFe, IF1e) Volker@43: AR = np.dot(ARFe, IF1e) Volker@43: BT = np.dot(ATFe, IF2e) Volker@43: BR = np.dot(ARFe, IF2e) Volker@43: Volker@43: # Correction parameters for normal measurements; they are independent of LDR --- the same as for TypeC = 6 Volker@43: if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out Volker@43: IS1 = np.array([IinE, 0., 0., VinE]) Volker@43: IS2 = np.array([0., QinE, -UinE, -2. * VinE]) Volker@43: GT = np.dot(ATF, IS1) Volker@43: GR = np.dot(ARF, IS1) Volker@43: HT = np.dot(ATF, IS2) Volker@43: HR = np.dot(ARF, IS2) Volker@43: else: Volker@43: IS1e = np.array([IinFo + DiC * QinFoe, DiC * IinFo + QinFoe, ZiC * (CosC * UinFoe + SinC * VinFo), Volker@43: -ZiC * (SinC * UinFoe - CosC * VinFo)]) Volker@43: IS2e = np.array([IinFa + DiC * QinFae, DiC * IinFa + QinFae, ZiC * (CosC * UinFae + SinC * VinFa), Volker@43: -ZiC * (SinC * UinFae - CosC * VinFa)]) Volker@43: GT = np.dot(ATFe, IS1e) Volker@43: GR = np.dot(ARFe, IS1e) Volker@43: HT = np.dot(ATFe, IS2e) Volker@43: HR = np.dot(ARFe, IS2e) Volker@43: Volker@43: elif (TypeC == 6): # diattenuator calibration +-22.5° rotated_diattenuator_X22x5deg.odt Volker@43: # parameters for calibration with aCal Volker@43: IF1e = np.array([IinF + sqr05 * DiC * QinFe, sqr05 * DiC * IinF + (1. - 0.5 * WiC) * QinFe, Volker@43: (1. - 0.5 * WiC) * UinFe + sqr05 * ZiC * SinC * VinF, Volker@43: -sqr05 * ZiC * SinC * UinFe + ZiC * CosC * VinF]) Volker@43: IF2e = np.array([sqr05 * DiC * UinFe, 0.5 * WiC * UinFe - sqr05 * ZiC * SinC * VinF, Volker@43: sqr05 * DiC * IinF + 0.5 * WiC * QinFe, sqr05 * ZiC * SinC * QinFe]) Volker@43: AT = np.dot(ATFe, IF1e) Volker@43: AR = np.dot(ARFe, IF1e) Volker@43: BT = np.dot(ATFe, IF2e) Volker@43: BR = np.dot(ARFe, IF2e) Volker@43: Volker@43: # Correction parameters for normal measurements; they are independent of LDR Volker@43: if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out Volker@43: # IS1 = np.array([IinE,0,0,VinE]) Volker@43: # IS2 = np.array([0,QinE,-UinE,-2*VinE]) Volker@43: IS1 = np.array([IinFo, 0., 0., VinFo]) Volker@43: IS2 = np.array([0., QinFa, UinFa, VinFa]) Volker@43: GT = np.dot(ATF, IS1) Volker@43: GR = np.dot(ARF, IS1) Volker@43: HT = np.dot(ATF, IS2) Volker@43: HR = np.dot(ARF, IS2) Volker@43: else: Volker@43: IS1e = np.array([IinFo + DiC * QinFoe, DiC * IinFo + QinFoe, ZiC * (CosC * UinFoe + SinC * VinFo), Volker@43: -ZiC * (SinC * UinFoe - CosC * VinFo)]) Volker@43: IS2e = np.array([IinFa + DiC * QinFae, DiC * IinFa + QinFae, ZiC * (CosC * UinFae + SinC * VinFa), Volker@43: -ZiC * (SinC * UinFae - CosC * VinFa)]) Volker@43: # IS1e = np.array([IinFo,0,0,VinFo]) Volker@43: # IS2e = np.array([0,QinFae,UinFae,VinFa]) Volker@43: GT = np.dot(ATFe, IS1e) Volker@43: GR = np.dot(ARFe, IS1e) Volker@43: HT = np.dot(ATFe, IS2e) Volker@43: HR = np.dot(ARFe, IS2e) Volker@43: Volker@43: else: Volker@43: print('Calibrator not implemented yet') Volker@43: sys.exit() Volker@43: Volker@43: elif LocC == 2: # C behind emitter optics Eq.57 ------------------------------------------------------- Volker@43: # print("Calibrator location not implemented yet") Volker@43: S2e = np.sin(np.deg2rad(2. * RotC)) Volker@43: C2e = np.cos(np.deg2rad(2. * RotC)) Volker@43: Volker@43: # AS with C before the receiver optics (see document rotated_diattenuator_X22x5deg.odt) Volker@43: AF1 = np.array([1, C2g * DiO, S2g * DiO, 0.]) Volker@43: AF2 = np.array([C2g * DiO, 1. - S2g ** 2 * WiO, S2g * C2g * WiO, -S2g * ZiO * SinO]) Volker@43: AF3 = np.array([S2g * DiO, S2g * C2g * WiO, 1. - C2g ** 2 * WiO, C2g * ZiO * SinO]) Volker@43: AF4 = np.array([0., S2g * SinO, -C2g * SinO, CosO]) Volker@43: Volker@43: ATF = (ATP1 * AF1 + ATP2 * AF2 + ATP3 * AF3 + ATP4 * AF4) Volker@43: ARF = (ARP1 * AF1 + ARP2 * AF2 + ARP3 * AF3 + ARP4 * AF4) Volker@43: ATF1 = ATF[0] Volker@43: ATF2 = ATF[1] Volker@43: ATF3 = ATF[2] Volker@43: ATF4 = ATF[3] Volker@43: ARF1 = ARF[0] Volker@43: ARF2 = ARF[1] Volker@43: ARF3 = ARF[2] Volker@43: ARF4 = ARF[3] Volker@43: Volker@43: # AS with C behind the emitter Volker@43: # terms without aCal Volker@43: ATE1o, ARE1o = ATF1, ARF1 Volker@43: ATE2o, ARE2o = 0., 0. Volker@43: ATE3o, ARE3o = 0., 0. Volker@43: ATE4o, ARE4o = ATF4, ARF4 Volker@43: # terms with aCal Volker@43: ATE1a, ARE1a = 0., 0. Volker@43: ATE2a, ARE2a = ATF2, ARF2 Volker@43: ATE3a, ARE3a = -ATF3, -ARF3 Volker@43: ATE4a, ARE4a = -2. * ATF4, -2. * ARF4 Volker@43: # rotated AinEa by epsilon Volker@43: ATE2ae = C2e * ATF2 + S2e * ATF3 Volker@43: ATE3ae = -S2e * ATF2 - C2e * ATF3 Volker@43: ARE2ae = C2e * ARF2 + S2e * ARF3 Volker@43: ARE3ae = -S2e * ARF2 - C2e * ARF3 Volker@43: Volker@43: ATE1 = ATE1o Volker@43: ATE2e = aCal * ATE2ae Volker@43: ATE3e = aCal * ATE3ae Volker@43: ATE4 = (1 - 2 * aCal) * ATF4 Volker@43: ARE1 = ARE1o Volker@43: ARE2e = aCal * ARE2ae Volker@43: ARE3e = aCal * ARE3ae Volker@43: ARE4 = (1 - 2 * aCal) * ARF4 Volker@43: Volker@43: # rotated IinE Volker@43: QinEe = C2e * QinE + S2e * UinE Volker@43: UinEe = C2e * UinE - S2e * QinE Volker@43: Volker@43: # Calibration signals and Calibration correction K from measurements with LDRCal / aCal Volker@43: if (TypeC == 2) or (TypeC == 1): # +++++++++ rotator calibration Eq. C.4 Volker@43: AT = ATE1o * IinE + (ATE4o + aCal * ATE4a) * h * VinE Volker@43: BT = aCal * (ATE3ae * QinEe - ATE2ae * h * UinEe) Volker@43: AR = ARE1o * IinE + (ARE4o + aCal * ARE4a) * h * VinE Volker@43: BR = aCal * (ARE3ae * QinEe - ARE2ae * h * UinEe) Volker@43: Volker@43: # Correction parameters for normal measurements; they are independent of LDR Volker@43: if (not RotationErrorEpsilonForNormalMeasurements): Volker@43: # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) Volker@43: GT = ATE1o * IinE + ATE4o * h * VinE Volker@43: GR = ARE1o * IinE + ARE4o * h * VinE Volker@43: HT = ATE2a * QinE + ATE3a * h * UinEe + ATE4a * h * VinE Volker@43: HR = ARE2a * QinE + ARE3a * h * UinEe + ARE4a * h * VinE Volker@43: else: Volker@43: GT = ATE1o * IinE + ATE4o * h * VinE Volker@43: GR = ARE1o * IinE + ARE4o * h * VinE Volker@43: HT = ATE2ae * QinE + ATE3ae * h * UinEe + ATE4a * h * VinE Volker@43: HR = ARE2ae * QinE + ARE3ae * h * UinEe + ARE4a * h * VinE Volker@43: Volker@43: elif (TypeC == 3) or (TypeC == 4): # +++++++++ linear polariser calibration Eq. C.5 Volker@43: # p = +45°, m = -45° Volker@43: AT = ATE1 * IinE + ZiC * CosC * (ATE2e * QinEe + ATE4 * VinE) + ATE3e * UinEe Volker@43: BT = DiC * (ATE1 * UinEe + ATE3e * IinE) + ZiC * SinC * (ATE4 * QinEe - ATE2e * VinE) Volker@43: AR = ARE1 * IinE + ZiC * CosC * (ARE2e * QinEe + ARE4 * VinE) + ARE3e * UinEe Volker@43: BR = DiC * (ARE1 * UinEe + ARE3e * IinE) + ZiC * SinC * (ARE4 * QinEe - ARE2e * VinE) Volker@43: Volker@43: # Correction parameters for normal measurements; they are independent of LDR Volker@43: if (not RotationErrorEpsilonForNormalMeasurements): Volker@43: # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) Volker@43: GT = ATE1o * IinE + ATE4o * VinE Volker@43: GR = ARE1o * IinE + ARE4o * VinE Volker@43: HT = ATE2a * QinE + ATE3a * UinE + ATE4a * VinE Volker@43: HR = ARE2a * QinE + ARE3a * UinE + ARE4a * VinE Volker@43: else: Volker@43: D = IinE + DiC * QinEe Volker@43: A = DiC * IinE + QinEe Volker@43: B = ZiC * (CosC * UinEe + SinC * VinE) Volker@43: C = -ZiC * (SinC * UinEe - CosC * VinE) Volker@43: GT = ATE1o * D + ATE4o * C Volker@43: GR = ARE1o * D + ARE4o * C Volker@43: HT = ATE2a * A + ATE3a * B + ATE4a * C Volker@43: HR = ARE2a * A + ARE3a * B + ARE4a * C Volker@43: Volker@43: elif (TypeC == 6): # real HWP calibration +-22.5° rotated_diattenuator_X22x5deg.odt Volker@43: # p = +22.5°, m = -22.5° Volker@43: IE1e = np.array([IinE + sqr05 * DiC * QinEe, sqr05 * DiC * IinE + (1 - 0.5 * WiC) * QinEe, Volker@43: (1 - 0.5 * WiC) * UinEe + sqr05 * ZiC * SinC * VinE, Volker@43: -sqr05 * ZiC * SinC * UinEe + ZiC * CosC * VinE]) Volker@43: IE2e = np.array([sqr05 * DiC * UinEe, 0.5 * WiC * UinEe - sqr05 * ZiC * SinC * VinE, Volker@43: sqr05 * DiC * IinE + 0.5 * WiC * QinEe, sqr05 * ZiC * SinC * QinEe]) Volker@43: ATEe = np.array([ATE1, ATE2e, ATE3e, ATE4]) Volker@43: AREe = np.array([ARE1, ARE2e, ARE3e, ARE4]) Volker@43: AT = np.dot(ATEe, IE1e) Volker@43: AR = np.dot(AREe, IE1e) Volker@43: BT = np.dot(ATEe, IE2e) Volker@43: BR = np.dot(AREe, IE2e) Volker@43: Volker@43: # Correction parameters for normal measurements; they are independent of LDR Volker@43: if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out Volker@43: GT = ATE1o * IinE + ATE4o * VinE Volker@43: GR = ARE1o * IinE + ARE4o * VinE Volker@43: HT = ATE2a * QinE + ATE3a * UinE + ATE4a * VinE Volker@43: HR = ARE2a * QinE + ARE3a * UinE + ARE4a * VinE Volker@43: else: Volker@43: D = IinE + DiC * QinEe Volker@43: A = DiC * IinE + QinEe Volker@43: B = ZiC * (CosC * UinEe + SinC * VinE) Volker@43: C = -ZiC * (SinC * UinEe - CosC * VinE) Volker@43: GT = ATE1o * D + ATE4o * C Volker@43: GR = ARE1o * D + ARE4o * C Volker@43: HT = ATE2a * A + ATE3a * B + ATE4a * C Volker@43: HR = ARE2a * A + ARE3a * B + ARE4a * C Volker@43: Volker@43: else: Volker@43: print('Calibrator not implemented yet') Volker@43: sys.exit() Volker@43: Volker@43: else: Volker@43: print("Calibrator location not implemented yet") Volker@43: sys.exit() Volker@43: Volker@43: # Determination of the correction K of the calibration factor. Volker@43: IoutTp = TTa * TiC * TiO * TiE * (AT + BT) Volker@43: IoutTm = TTa * TiC * TiO * TiE * (AT - BT) Volker@43: IoutRp = TRa * TiC * TiO * TiE * (AR + BR) Volker@43: IoutRm = TRa * TiC * TiO * TiE * (AR - BR) Volker@43: # --- Results and Corrections; electronic etaR and etaT are assumed to be 1 Volker@43: Etapx = IoutRp / IoutTp Volker@43: Etamx = IoutRm / IoutTm Volker@43: Etax = (Etapx * Etamx) ** 0.5 Volker@43: Volker@43: 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@43: K = Etax / Eta Volker@43: Volker@43: # For comparison with Volkers Libreoffice Müller Matrix spreadsheet Volker@43: # Eta_test_p = (IoutRp/IoutTp) Volker@43: # Eta_test_m = (IoutRm/IoutTm) Volker@43: # Eta_test = (Eta_test_p*Eta_test_m)**0.5 Volker@43: Volker@43: # ----- random error calculation ---------- Volker@43: # noise must be calculated with the photon counts of measured signals; Volker@43: # relative standard deviation of calibration signals with LDRcal; assumed to be statisitcally independent Volker@43: # normalised noise errors Volker@43: if (CalcFrom0deg): Volker@43: dIoutTp = (NCalT * IoutTp) ** -0.5 Volker@43: dIoutTm = (NCalT * IoutTm) ** -0.5 Volker@43: dIoutRp = (NCalR * IoutRp) ** -0.5 Volker@43: dIoutRm = (NCalR * IoutRm) ** -0.5 Volker@43: else: Volker@43: dIoutTp = (NCalT ** -0.5) Volker@43: dIoutTm = (NCalT ** -0.5) Volker@43: dIoutRp = (NCalR ** -0.5) Volker@43: dIoutRm = (NCalR ** -0.5) Volker@43: # Forward simulated 0°-signals with LDRCal with atrue; from input file Volker@43: Volker@43: It = TTa * TiO * TiE * (GT + atrue * HT) Volker@43: Ir = TRa * TiO * TiE * (GR + atrue * HR) Volker@43: # relative standard deviation of standard signals with LDRmeas; assumed to be statisitcally independent Volker@43: if (CalcFrom0deg): # this works! Volker@43: dIt = ((It * NI * eFacT) ** -0.5) Volker@43: dIr = ((Ir * NI * eFacR) ** -0.5) Volker@43: ''' Volker@43: dIt = ((NCalT * It / IoutTp * NILfac / TCalT) ** -0.5) Volker@43: dIr = ((NCalR * Ir / IoutRp * NILfac / TCalR) ** -0.5) Volker@43: ''' Volker@43: else: # does this work? Why not as above? Volker@43: dIt = ((NCalT * 2 * NILfac / TCalT ) ** -0.5) Volker@43: dIr = ((NCalR * 2 * NILfac / TCalR) ** -0.5) Volker@43: Volker@43: # ----- Forward simulated LDRsim = 1/Eta*Ir/It # simulated LDR* with Y from input file Volker@43: LDRsim = Ir / It # simulated uncorrected LDR with Y from input file Volker@43: # Corrected LDRsimCorr from forward simulated LDRsim (atrue) Volker@43: # LDRsimCorr = (1./Eta*LDRsim*(GT+HT)-(GR+HR))/((GR-HR)-1./Eta*LDRsim*(GT-HT)) Volker@43: ''' Volker@43: if ((Y == -1.) and (abs(RotL0) < 45)) or ((Y == +1.) and (abs(RotL0) > 45)): Volker@43: LDRsimx = 1. / LDRsim / Etax Volker@43: else: Volker@43: LDRsimx = LDRsim / Etax Volker@43: ''' Volker@43: LDRsimx = LDRsim Volker@43: Volker@43: # The following is correct without doubt Volker@43: # LDRCorr = (LDRsim/(Etax/K)*(GT+HT)-(GR+HR))/((GR-HR)-LDRsim/(Etax/K)*(GT-HT)) Volker@43: Volker@43: # The following is a test whether the equations for calibration Etax and normal signal (GHK, LDRsim) are consistent Volker@43: LDRCorr = (LDRsim / (Etax / K) * (GT + HT) - (GR + HR)) / ((GR - HR) - LDRsim / (Etax / K) * (GT - HT)) Volker@43: # here we could also use Eta instead of Etax / K => how to test whether Etax is correct? => comparison with MüllerMatrix simulation! Volker@43: # Without any correction: only measured It, Ir, EtaX are used Volker@43: LDRunCorr = LDRsim / Etax Volker@43: # 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@43: Volker@43: #LDRCorr = LDRsimx # for test only Volker@43: Volker@43: 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@43: Volker@43: return (IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr, Volker@43: GT, HT, GR, HR, K, Eta, LDRsimx, LDRCorr, DTa, DRa, TTa, TRa, F11sim, LDRunCorr) Volker@43: Volker@43: Volker@43: Volker@43: # ******************************************************************************************************************************* Volker@43: Volker@43: # --- CALC with assumed true parameters from the input file Volker@43: LDRtrue = LDRtrue2 Volker@43: IoutTp0, IoutTm0, IoutRp0, IoutRm0, It0, Ir0, dIoutTp0, dIoutTm0, dIoutRp0, dIoutRm0, dIt0, dIr0, \ Volker@43: GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \ Volker@43: Calc(TCalT, TCalR, NCalT, NCalR, Qin0, Vin0, RotL0, RotE0, RetE0, DiE0, Volker@43: RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, Volker@43: ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0) Volker@43: Etax0 = K0 * Eta0 Volker@43: Etapx0 = IoutRp0 / IoutTp0 Volker@43: Etamx0 = IoutRm0 / IoutTm0 Volker@43: # --- Print parameters to console and output file Volker@43: OutputFile = 'output_' + InputFile[0:-3] + '_' + fname[0:-3] +'.dat' Volker@43: with open('output_files\\' + OutputFile, 'w') as f: Volker@43: with redirect_stdout(f): Volker@43: print("From ", dname) Volker@43: print("Running ", fname) Volker@43: print("Reading input file ", InputFile) # , " for Lidar system :", EID, ", ", LID) Volker@43: print("for Lidar system: ", EID, ", ", LID) Volker@43: # --- Print iput information********************************* Volker@43: print(" --- Input parameters: value ±error / ±steps ----------------------") Volker@43: print("{0:7}{1:17} {2:6.4f}±{3:7.4f}/{4:2d}".format("Laser: ", "Qin =", Qin0, dQin, nQin)) Volker@43: print("{0:7}{1:17} {2:6.4f}±{3:7.4f}/{4:2d}".format("", "Vin =", Vin0, dVin, nVin)) Volker@43: print("{0:7}{1:17} {2:6.4f}±{3:7.4f}/{4:2d}".format("", "Rotation alpha = ", RotL0, dRotL, nRotL)) Volker@43: print("{0:7}{1:15} {2:8.4f} {3:17}".format("", "=> DOP", ((Qin ** 2 + Vin ** 2) ** 0.5), " (degree of polarisation)")) Volker@43: Volker@43: print("Optic: Diatt., Tunpol, Retard., Rotation (deg)") Volker@43: 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@43: "Emitter ", DiE0, dDiE, TiE, RetE0, dRetE, RotE0, dRotE, nDiE, nRetE, nRotE)) Volker@43: 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@43: "Receiver ", DiO0, dDiO, TiO, RetO0, dRetO, RotO0, dRotO, nDiO, nRetO, nRotO)) Volker@43: 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@43: "Calibrator ", DiC0, dDiC, TiC, RetC0, dRetC, RotC0, dRotC, nDiC, nRetC, nRotC)) Volker@43: print("{0:12}".format(" Pol.-filter ------ ")) Volker@43: print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format( Volker@43: "ERT, RotT :", ERaT0, dERaT, nERaT, RotaT0, dRotaT, nRotaT)) Volker@43: print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format( Volker@43: "ERR, RotR :", ERaR0, dERaR, nERaR, RotaR0, dRotaR, nRotaR)) Volker@43: print("{0:12}".format(" PBS ------ ")) Volker@43: print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format( Volker@43: "TP,TS :", TP0, dTP, nTP, TS0, dTS, nTS)) Volker@43: print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format( Volker@43: "RP,RS :", RP0, dRP, nRP, RS0, dRS, nRS)) Volker@43: print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}, {5:1.0f}".format( Volker@43: "DT,TT,DR,TR,Y :", DiT, TiT, DiR, TiR, Y)) Volker@43: print("{0:12}".format(" Combined PBS + Pol.-filter ------ ")) Volker@43: print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}".format( Volker@43: "DT,TT,DR,TR :", DTa0, TTa0, DRa0, TRa0)) Volker@43: print("{0:26}: {1:6.3f}± {2:5.3f}/{3:2d}".format( Volker@43: "LDRCal during calibration in calibration range", LDRCal0, dLDRCal, nLDRCal)) Volker@43: print("{0:12}".format(" --- Additional ND filter attenuation (transmission) during the calibration ---")) Volker@43: print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format( Volker@43: "TCalT,TCalR :", TCalT0, dTCalT, nTCalT, TCalR0, dTCalR, nTCalR)) Volker@43: print() Volker@44: print(Type[TypeC],"calibrator is located", Loc[LocC]) Volker@44: print("Rotation error epsilon is considered also for normal measurements = ", RotationErrorEpsilonForNormalMeasurements) Volker@44: print("The PBS incidence plane is ", dY[int(Y + 1)], "to the reference plane" ) Volker@44: print("The laser polarisation in the reference plane is finally", dY2[int(Y + 1)], "=>", dY3) Volker@43: print("RS_RP_depend_on_TS_TP = ", RS_RP_depend_on_TS_TP) Volker@43: # end of print actual system parameters Volker@43: # ****************************************************************************** Volker@43: Volker@43: Volker@43: print() Volker@43: Volker@43: K0List = np.zeros(7) Volker@43: LDRsimxList = np.zeros(7) Volker@43: LDRCalList = 0.0, 0.004, 0.02, 0.1, 0.2, 0.3, 0.45 Volker@43: # 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@43: # Still with assumed true parameters in input file Volker@43: Volker@43: ''' Volker@43: facIt = NCalT / TCalT0 * NILfac Volker@43: facIr = NCalR / TCalR0 * NILfac Volker@43: ''' Volker@43: facIt = NI * eFacT Volker@43: facIr = NI * eFacR Volker@43: if (bPlotEtax): Volker@43: # check error signals Volker@43: # dIs are relative stdevs Volker@43: print("LDRCal, IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp,dIoutTm,dIoutRp,dIoutRm,dIt, dIr") Volker@43: Volker@43: for i, LDRCal in enumerate(LDRCalList): Volker@43: IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr, \ Volker@43: GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \ Volker@43: Calc(TCalT0, TCalR0, NCalT, NCalR, Qin0, Vin0, RotL0, RotE0, RetE0, DiE0, Volker@43: RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, Volker@43: ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal) Volker@43: K0List[i] = K0 Volker@43: LDRsimxList[i] = LDRsimx Volker@43: Volker@43: if (bPlotEtax): Volker@43: # check error signals Volker@43: 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@43: #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@43: # end check error signals Volker@43: print('===========================================================================================================') Volker@43: print("{0:8},{1:8},{2:8},{3:8},{4:9},{5:8},{6:9},{7:9},{8:9},{9:9},{10:9}".format( Volker@43: " GR", " GT", " HR", " HT", " K(0.000)", " K(0.004)", " K(0.02)", " K(0.1)", " K(0.2)", " K(0.3)", " K(0.45)")) Volker@43: print("{0:8.5f},{1:8.5f},{2:8.5f},{3:8.5f},{4:9.5f},{5:9.5f},{6:9.5f},{7:9.5f},{8:9.5f},{9:9.5f},{10:9.5f}".format( Volker@43: GR0, GT0, HR0, HT0, K0List[0], K0List[1], K0List[2], K0List[3], K0List[4], K0List[5], K0List[6])) Volker@43: print('===========================================================================================================') Volker@43: print() Volker@43: print("Errors from neglecting GHK corrections and/or calibration:") Volker@43: print("{0:>10},{1:>10},{2:>10},{3:>10},{4:>10},{5:>10}".format( Volker@43: "LDRtrue", "LDRunCorr", "1/LDRunCorr", "LDRsimx", "1/LDRsimx", "LDRCorr")) Volker@43: Volker@43: aF11sim0 = np.zeros(5) Volker@43: LDRrange = np.zeros(5) Volker@43: LDRsim0 = np.zeros(5) Volker@43: LDRrange = [0.004, 0.02, 0.1, 0.3, 0.45] # list Volker@43: LDRrange[0] = LDRtrue2 # value in the input file; default 0.004 Volker@43: Volker@43: # The loop over LDRtrueList is only for checking how much the uncorrected LDRsimx deviates from LDRtrue ... and whether the corrections work. Volker@43: # LDRsimx = LDRsim = Ir / It or 1/LDRsim Volker@43: # Still with assumed true parameters in input file Volker@43: for i, LDRtrue in enumerate(LDRrange): Volker@43: #for LDRtrue in LDRrange: Volker@43: IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr, \ Volker@43: GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \ Volker@43: Calc(TCalT0, TCalR0, NCalT, NCalR, Qin0, Vin0, RotL0, RotE0, RetE0, DiE0, Volker@43: RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, Volker@43: ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0) Volker@43: 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@43: aF11sim0[i] = F11sim0 Volker@43: LDRsim0[i] = Ir / It Volker@43: # the assumed true aF11sim0 results will be used below to calc the deviation from the real signals Volker@43: print("LDRsimx = LDR of the nominal system directly from measured signals without calibration and GHK-corrections") Volker@43: print("LDRunCorr = LDR of the nominal system directly from measured signals with calibration but without GHK-corrections; electronic amplifications = 1 assumed") Volker@43: print("LDRCorr = LDR calibrated and GHK-corrected") Volker@43: print() Volker@43: print("Errors from signal noise:") Volker@43: 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@43: NI, NCalT, NCalR, NILfac, nNCal, nNI, 1.0 / NI ** 0.5)) Volker@43: print() Volker@43: print() Volker@43: '''# das muß wieder weg Volker@43: print("IoutTp, IoutTm, IoutRp, IoutRm, It , Ir , dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr") Volker@43: LDRCal = 0.01 Volker@43: for i, LDRtrue in enumerate(LDRrange): Volker@43: IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr, \ Volker@43: GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \ Volker@43: Calc(TCalT0, TCalR0, NCalT, NCalR, DOLP0, RotL0, RotE0, RetE0, DiE0, Volker@43: RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, Volker@43: ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0) Volker@43: 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@43: IoutTp * NCalT, IoutTm * NCalT, IoutRp * NCalR, IoutRm * NCalR, It * facIt, Ir * facIr, Volker@43: dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr)) Volker@43: aF11sim0[i] = F11sim0 Volker@43: # the assumed true aF11sim0 results will be used below to calc the deviation from the real signals Volker@43: # bis hierher weg Volker@43: ''' Volker@43: Volker@43: file = open('output_files\\' + OutputFile, 'r') Volker@43: print(file.read()) Volker@43: file.close() Volker@43: Volker@43: # --- CALC again assumed truth with LDRCal0 and with assumed true parameters in input file to reset all 0-values Volker@43: LDRtrue = LDRtrue2 Volker@43: IoutTp0, IoutTm0, IoutRp0, IoutRm0, It0, Ir0, dIoutTp0, dIoutTm0, dIoutRp0, dIoutRm0, dIt0, dIr0, \ Volker@43: GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \ Volker@43: Calc(TCalT0, TCalR0, NCalT, NCalR, Qin0, Vin0, RotL0, RotE0, RetE0, DiE0, Volker@43: RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, Volker@43: ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0) Volker@43: Etax0 = K0 * Eta0 Volker@43: Etapx0 = IoutRp0 / IoutTp0 Volker@43: Etamx0 = IoutRm0 / IoutTm0 Volker@43: ''' Volker@43: if(PrintToOutputFile): Volker@43: f = open('output_ver7.dat', 'w') Volker@43: old_target = sys.stdout Volker@43: sys.stdout = f Volker@43: Volker@43: print("something") Volker@43: Volker@43: if(PrintToOutputFile): Volker@43: sys.stdout.flush() Volker@43: f.close Volker@43: sys.stdout = old_target Volker@43: ''' Volker@43: if (Error_Calc): Volker@43: # --- CALC again assumed truth with LDRCal0 and with assumed true parameters in input file to reset all 0-values Volker@43: LDRtrue = LDRtrue2 Volker@43: IoutTp0, IoutTm0, IoutRp0, IoutRm0, It0, Ir0, dIoutTp0, dIoutTm0, dIoutRp0, dIoutRm0, dIt0, dIr0, \ Volker@43: GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \ Volker@43: Calc(TCalT0, TCalR0, NCalT, NCalR, Qin0, Vin0, RotL0, RotE0, RetE0, DiE0, Volker@43: RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, Volker@43: ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0) Volker@43: Etax0 = K0 * Eta0 Volker@43: Etapx0 = IoutRp0 / IoutTp0 Volker@43: Etamx0 = IoutRm0 / IoutTm0 Volker@43: Volker@43: # --- Start Error calculation with variable parameters ------------------------------------------------------------------ Volker@43: # error nNCal: one-sigma in steps to left and right for calibration signals Volker@43: # error nNI: one-sigma in steps to left and right for 0° signals Volker@43: Volker@43: iN = -1 Volker@43: N = ((nTCalT * 2 + 1) * (nTCalR * 2 + 1) * Volker@43: (nNCal * 2 + 1) ** 4 * (nNI * 2 + 1) ** 2 * Volker@43: (nQin * 2 + 1) * (nVin * 2 + 1) * (nRotL * 2 + 1) * Volker@43: (nRotE * 2 + 1) * (nRetE * 2 + 1) * (nDiE * 2 + 1) * Volker@43: (nRotO * 2 + 1) * (nRetO * 2 + 1) * (nDiO * 2 + 1) * Volker@43: (nRotC * 2 + 1) * (nRetC * 2 + 1) * (nDiC * 2 + 1) * Volker@43: (nTP * 2 + 1) * (nTS * 2 + 1) * (nRP * 2 + 1) * (nRS * 2 + 1) * (nERaT * 2 + 1) * (nERaR * 2 + 1) * Volker@43: (nRotaT * 2 + 1) * (nRotaR * 2 + 1) * (nRetT * 2 + 1) * (nRetR * 2 + 1) * (nLDRCal * 2 + 1)) Volker@43: print("number of system variations N = ", N, " ", end="") Volker@43: Volker@43: if N > 1e6: Volker@43: if user_yes_no_query('Warning: processing ' + str( Volker@43: N) + ' samples will take very long. Do you want to proceed?') == 0: sys.exit() Volker@43: if N > 5e6: Volker@43: if user_yes_no_query('Warning: the memory required for ' + str(N) + ' samples might be ' + '{0:5.1f}'.format( Volker@43: N / 4e6) + ' GB. Do you anyway want to proceed?') == 0: sys.exit() Volker@43: Volker@43: # if user_yes_no_query('Warning: processing' + str(N) + ' samples will take very long. Do you want to proceed?') == 0: sys.exit() Volker@43: Volker@43: # --- Arrays for plotting ------ Volker@43: LDRmin = np.zeros(5) Volker@43: LDRmax = np.zeros(5) Volker@43: LDRstd = np.zeros(5) Volker@43: LDRmean = np.zeros(5) Volker@43: LDRmedian = np.zeros(5) Volker@43: LDRskew = np.zeros(5) Volker@43: LDRkurt = np.zeros(5) Volker@43: LDRsimmin = np.zeros(5) Volker@43: LDRsimmax = np.zeros(5) Volker@43: LDRsimmean = np.zeros(5) Volker@43: Volker@43: F11min = np.zeros(5) Volker@43: F11max = np.zeros(5) Volker@43: Etaxmin = np.zeros(5) Volker@43: Etaxmax = np.zeros(5) Volker@43: Volker@43: aQin = np.zeros(N) Volker@43: aVin = np.zeros(N) Volker@43: aERaT = np.zeros(N) Volker@43: aERaR = np.zeros(N) Volker@43: aRotaT = np.zeros(N) Volker@43: aRotaR = np.zeros(N) Volker@43: aRetT = np.zeros(N) Volker@43: aRetR = np.zeros(N) Volker@43: aTP = np.zeros(N) Volker@43: aTS = np.zeros(N) Volker@43: aRP = np.zeros(N) Volker@43: aRS = np.zeros(N) Volker@43: aDiE = np.zeros(N) Volker@43: aDiO = np.zeros(N) Volker@43: aDiC = np.zeros(N) Volker@43: aRotC = np.zeros(N) Volker@43: aRetC = np.zeros(N) Volker@43: aRotL = np.zeros(N) Volker@43: aRetE = np.zeros(N) Volker@43: aRotE = np.zeros(N) Volker@43: aRetO = np.zeros(N) Volker@43: aRotO = np.zeros(N) Volker@43: aLDRCal = np.zeros(N) Volker@43: aNCalTp = np.zeros(N) Volker@43: aNCalTm = np.zeros(N) Volker@43: aNCalRp = np.zeros(N) Volker@43: aNCalRm = np.zeros(N) Volker@43: aNIt = np.zeros(N) Volker@43: aNIr = np.zeros(N) Volker@43: aTCalT = np.zeros(N) Volker@43: aTCalR = np.zeros(N) Volker@43: Volker@43: # each np.zeros((LDRrange, N)) array has the same N-dependency Volker@43: aLDRcorr = np.zeros((5, N)) Volker@43: aLDRsim = np.zeros((5, N)) Volker@43: aF11corr = np.zeros((5, N)) Volker@43: aPLDR = np.zeros((5, N)) Volker@43: aEtax = np.zeros((5, N)) Volker@43: aEtapx = np.zeros((5, N)) Volker@43: aEtamx = np.zeros((5, N)) Volker@43: Volker@43: # np.zeros((GHKs, N)) Volker@43: aGHK = np.zeros((5, N)) Volker@43: Volker@43: atime = clock() Volker@43: dtime = clock() Volker@43: Volker@43: # --- Calc Error signals Volker@43: # ---- Do the calculations of bra-ket vectors Volker@43: h = -1. if TypeC == 2 else 1 Volker@43: Volker@43: for iLDRCal in range(-nLDRCal, nLDRCal + 1): Volker@43: # from input file: LDRCal for calibration measurements Volker@43: LDRCal = LDRCal0 Volker@43: if nLDRCal > 0: Volker@43: LDRCal = LDRCal0 + iLDRCal * dLDRCal / nLDRCal Volker@43: # provides the intensities of the calibration measurements at various LDRCal for signal noise errors Volker@43: # IoutTp, IoutTm, IoutRp, IoutRm, dIoutTp, dIoutTm, dIoutRp, dIoutRm Volker@43: Volker@43: aCal = (1. - LDRCal) / (1. + LDRCal) Volker@43: for iQin, iVin, iRotL, iRotE, iRetE, iDiE \ Volker@43: in [(iQin, iVin, iRotL, iRotE, iRetE, iDiE) Volker@43: for iQin in range(-nQin, nQin + 1) Volker@43: for iVin in range(-nVin, nVin + 1) Volker@43: for iRotL in range(-nRotL, nRotL + 1) Volker@43: for iRotE in range(-nRotE, nRotE + 1) Volker@43: for iRetE in range(-nRetE, nRetE + 1) Volker@43: for iDiE in range(-nDiE, nDiE + 1)]: Volker@43: Volker@43: if nQin > 0: Qin = Qin0 + iQin * dQin / nQin Volker@43: if nVin > 0: Vin = Vin0 + iVin * dVin / nVin Volker@43: if nRotL > 0: RotL = RotL0 + iRotL * dRotL / nRotL Volker@43: if nRotE > 0: RotE = RotE0 + iRotE * dRotE / nRotE Volker@43: if nRetE > 0: RetE = RetE0 + iRetE * dRetE / nRetE Volker@43: if nDiE > 0: DiE = DiE0 + iDiE * dDiE / nDiE Volker@43: Volker@43: if ((Qin ** 2 + Vin ** 2) ** 0.5) > 1.0: Volker@43: print("Error: degree of polarisation of laser > 1. Check Qin and Vin! ") Volker@43: sys.exit() Volker@43: # angles of emitter and laser and calibrator and receiver optics Volker@43: # RotL = alpha, RotE = beta, RotO = gamma, RotC = epsilon Volker@43: S2a = np.sin(2 * np.deg2rad(RotL)) Volker@43: C2a = np.cos(2 * np.deg2rad(RotL)) Volker@43: S2b = np.sin(2 * np.deg2rad(RotE)) Volker@43: C2b = np.cos(2 * np.deg2rad(RotE)) Volker@43: S2ab = np.sin(np.deg2rad(2 * RotL - 2 * RotE)) Volker@43: C2ab = np.cos(np.deg2rad(2 * RotL - 2 * RotE)) Volker@43: Volker@43: # Laser with Degree of linear polarization DOLP Volker@43: IinL = 1. Volker@43: QinL = Qin Volker@43: UinL = 0. Volker@43: VinL = Vin Volker@43: # VinL = (1. - DOLP ** 2) ** 0.5 Volker@43: Volker@43: # Stokes Input Vector rotation Eq. E.4 Volker@43: A = C2a * QinL - S2a * UinL Volker@43: B = S2a * QinL + C2a * UinL Volker@43: # Stokes Input Vector rotation Eq. E.9 Volker@43: C = C2ab * QinL - S2ab * UinL Volker@43: D = S2ab * QinL + C2ab * UinL Volker@43: Volker@43: # emitter optics Volker@43: CosE = np.cos(np.deg2rad(RetE)) Volker@43: SinE = np.sin(np.deg2rad(RetE)) Volker@43: ZiE = (1. - DiE ** 2) ** 0.5 Volker@43: WiE = (1. - ZiE * CosE) Volker@43: Volker@43: # Stokes Input Vector after emitter optics equivalent to Eq. E.9 with already rotated input vector from Eq. E.4 Volker@43: # b = beta Volker@43: IinE = (IinL + DiE * C) Volker@43: QinE = (C2b * DiE * IinL + A + S2b * (WiE * D - ZiE * SinE * VinL)) Volker@43: UinE = (S2b * DiE * IinL + B - C2b * (WiE * D - ZiE * SinE * VinL)) Volker@43: VinE = (-ZiE * SinE * D + ZiE * CosE * VinL) Volker@43: Volker@43: # ------------------------- Volker@43: # F11 assuemd to be = 1 => measured: F11m = IinP / IinE with atrue Volker@43: # F11sim = (IinE + DiO*atrue*(C2g*QinE - S2g*UinE))/IinE Volker@43: # ------------------------- Volker@43: Volker@43: for iRotO, iRetO, iDiO, iRotC, iRetC, iDiC, iTP, iTS, iRP, iRS, iERaT, iRotaT, iRetT, iERaR, iRotaR, iRetR \ Volker@43: in [ Volker@43: (iRotO, iRetO, iDiO, iRotC, iRetC, iDiC, iTP, iTS, iRP, iRS, iERaT, iRotaT, iRetT, iERaR, iRotaR, iRetR) Volker@43: for iRotO in range(-nRotO, nRotO + 1) Volker@43: for iRetO in range(-nRetO, nRetO + 1) Volker@43: for iDiO in range(-nDiO, nDiO + 1) Volker@43: for iRotC in range(-nRotC, nRotC + 1) Volker@43: for iRetC in range(-nRetC, nRetC + 1) Volker@43: for iDiC in range(-nDiC, nDiC + 1) Volker@43: for iTP in range(-nTP, nTP + 1) Volker@43: for iTS in range(-nTS, nTS + 1) Volker@43: for iRP in range(-nRP, nRP + 1) Volker@43: for iRS in range(-nRS, nRS + 1) Volker@43: for iERaT in range(-nERaT, nERaT + 1) Volker@43: for iRotaT in range(-nRotaT, nRotaT + 1) Volker@43: for iRetT in range(-nRetT, nRetT + 1) Volker@43: for iERaR in range(-nERaR, nERaR + 1) Volker@43: for iRotaR in range(-nRotaR, nRotaR + 1) Volker@43: for iRetR in range(-nRetR, nRetR + 1)]: Volker@43: Volker@43: if nRotO > 0: RotO = RotO0 + iRotO * dRotO / nRotO Volker@43: if nRetO > 0: RetO = RetO0 + iRetO * dRetO / nRetO Volker@43: if nDiO > 0: DiO = DiO0 + iDiO * dDiO / nDiO Volker@43: if nRotC > 0: RotC = RotC0 + iRotC * dRotC / nRotC Volker@43: if nRetC > 0: RetC = RetC0 + iRetC * dRetC / nRetC Volker@43: if nDiC > 0: DiC = DiC0 + iDiC * dDiC / nDiC Volker@43: if nTP > 0: TP = TP0 + iTP * dTP / nTP Volker@43: if nTS > 0: TS = TS0 + iTS * dTS / nTS Volker@43: if nRP > 0: RP = RP0 + iRP * dRP / nRP Volker@43: if nRS > 0: RS = RS0 + iRS * dRS / nRS Volker@43: if nERaT > 0: ERaT = ERaT0 + iERaT * dERaT / nERaT Volker@43: if nRotaT > 0: RotaT = RotaT0 + iRotaT * dRotaT / nRotaT Volker@43: if nRetT > 0: RetT = RetT0 + iRetT * dRetT / nRetT Volker@43: if nERaR > 0: ERaR = ERaR0 + iERaR * dERaR / nERaR Volker@43: if nRotaR > 0: RotaR = RotaR0 + iRotaR * dRotaR / nRotaR Volker@43: if nRetR > 0: RetR = RetR0 + iRetR * dRetR / nRetR Volker@43: Volker@43: # print("{0:5.2f}, {1:5.2f}, {2:5.2f}, {3:10d}".format(RotL, RotE, RotO, iN)) Volker@43: Volker@43: # receiver optics Volker@43: CosO = np.cos(np.deg2rad(RetO)) Volker@43: SinO = np.sin(np.deg2rad(RetO)) Volker@43: ZiO = (1. - DiO ** 2) ** 0.5 Volker@43: WiO = (1. - ZiO * CosO) Volker@43: S2g = np.sin(np.deg2rad(2 * RotO)) Volker@43: C2g = np.cos(np.deg2rad(2 * RotO)) Volker@43: # calibrator Volker@43: CosC = np.cos(np.deg2rad(RetC)) Volker@43: SinC = np.sin(np.deg2rad(RetC)) Volker@43: ZiC = (1. - DiC ** 2) ** 0.5 Volker@43: WiC = (1. - ZiC * CosC) Volker@43: Volker@43: # analyser Volker@43: # For POLLY_XTs Volker@43: if (RS_RP_depend_on_TS_TP): Volker@43: RS = 1.0 - TS Volker@43: RP = 1.0 - TP Volker@43: TiT = 0.5 * (TP + TS) Volker@43: DiT = (TP - TS) / (TP + TS) Volker@43: ZiT = (1. - DiT ** 2.) ** 0.5 Volker@43: TiR = 0.5 * (RP + RS) Volker@43: DiR = (RP - RS) / (RP + RS) Volker@43: ZiR = (1. - DiR ** 2.) ** 0.5 Volker@43: CosT = np.cos(np.deg2rad(RetT)) Volker@43: SinT = np.sin(np.deg2rad(RetT)) Volker@43: CosR = np.cos(np.deg2rad(RetR)) Volker@43: SinR = np.sin(np.deg2rad(RetR)) Volker@43: Volker@43: # cleaning pol-filter Volker@43: DaT = (1.0 - ERaT) / (1.0 + ERaT) Volker@43: DaR = (1.0 - ERaR) / (1.0 + ERaR) Volker@43: TaT = 0.5 * (1.0 + ERaT) Volker@43: TaR = 0.5 * (1.0 + ERaR) Volker@43: Volker@43: S2aT = np.sin(np.deg2rad(h * 2.0 * RotaT)) Volker@43: C2aT = np.cos(np.deg2rad(2.0 * RotaT)) Volker@43: S2aR = np.sin(np.deg2rad(h * 2.0 * RotaR)) Volker@43: C2aR = np.cos(np.deg2rad(2.0 * RotaR)) Volker@43: Volker@43: # Analyzer As before the PBS Eq. D.5; combined PBS and cleaning pol-filter Volker@43: ATPT = (1 + C2aT * DaT * DiT) # unpolarized transmission correction Volker@43: TTa = TiT * TaT * ATPT # unpolarized transmission Volker@43: ATP1 = 1.0 Volker@43: ATP2 = Y * (DiT + C2aT * DaT) / ATPT Volker@43: ATP3 = Y * S2aT * DaT * ZiT * CosT / ATPT Volker@43: ATP4 = S2aT * DaT * ZiT * SinT / ATPT Volker@43: ATP = np.array([ATP1, ATP2, ATP3, ATP4]) Volker@43: DTa = ATP2 * Y Volker@43: Volker@43: ARPT = (1 + C2aR * DaR * DiR) # unpolarized transmission correction Volker@43: TRa = TiR * TaR * ARPT # unpolarized transmission Volker@43: ARP1 = 1 Volker@43: ARP2 = Y * (DiR + C2aR * DaR) / ARPT Volker@43: ARP3 = Y * S2aR * DaR * ZiR * CosR / ARPT Volker@43: ARP4 = S2aR * DaR * ZiR * SinR / ARPT Volker@43: ARP = np.array([ARP1, ARP2, ARP3, ARP4]) Volker@43: DRa = ARP2 * Y Volker@43: Volker@43: # ---- Calculate signals and correction parameters for diffeent locations and calibrators Volker@43: if LocC == 4: # Calibrator before the PBS Volker@43: # print("Calibrator location not implemented yet") Volker@43: Volker@43: # S2ge = np.sin(np.deg2rad(2*RotO + h*2*RotC)) Volker@43: # C2ge = np.cos(np.deg2rad(2*RotO + h*2*RotC)) Volker@43: S2e = np.sin(np.deg2rad(h * 2 * RotC)) Volker@43: C2e = np.cos(np.deg2rad(2 * RotC)) Volker@43: # rotated AinP by epsilon Eq. C.3 Volker@43: ATP2e = C2e * ATP2 + S2e * ATP3 Volker@43: ATP3e = C2e * ATP3 - S2e * ATP2 Volker@43: ARP2e = C2e * ARP2 + S2e * ARP3 Volker@43: ARP3e = C2e * ARP3 - S2e * ARP2 Volker@43: ATPe = np.array([ATP1, ATP2e, ATP3e, ATP4]) Volker@43: ARPe = np.array([ARP1, ARP2e, ARP3e, ARP4]) Volker@43: # Stokes Input Vector before the polarising beam splitter Eq. E.31 Volker@43: A = C2g * QinE - S2g * UinE Volker@43: B = S2g * QinE + C2g * UinE Volker@43: # C = (WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE) Volker@43: Co = ZiO * SinO * VinE Volker@43: Ca = (WiO * B - 2 * ZiO * SinO * VinE) Volker@43: # C = Co + aCal*Ca Volker@43: # IinP = (IinE + DiO*aCal*A) Volker@43: # QinP = (C2g*DiO*IinE + aCal*QinE - S2g*C) Volker@43: # UinP = (S2g*DiO*IinE - aCal*UinE + C2g*C) Volker@43: # VinP = (ZiO*SinO*aCal*B + ZiO*CosO*(1-2*aCal)*VinE) Volker@43: IinPo = IinE Volker@43: QinPo = (C2g * DiO * IinE - S2g * Co) Volker@43: UinPo = (S2g * DiO * IinE + C2g * Co) Volker@43: VinPo = ZiO * CosO * VinE Volker@43: Volker@43: IinPa = DiO * A Volker@43: QinPa = QinE - S2g * Ca Volker@43: UinPa = -UinE + C2g * Ca Volker@43: VinPa = ZiO * (SinO * B - 2 * CosO * VinE) Volker@43: Volker@43: IinP = IinPo + aCal * IinPa Volker@43: QinP = QinPo + aCal * QinPa Volker@43: UinP = UinPo + aCal * UinPa Volker@43: VinP = VinPo + aCal * VinPa Volker@43: # Stokes Input Vector before the polarising beam splitter rotated by epsilon Eq. C.3 Volker@43: # QinPe = C2e*QinP + S2e*UinP Volker@43: # UinPe = C2e*UinP - S2e*QinP Volker@43: QinPoe = C2e * QinPo + S2e * UinPo Volker@43: UinPoe = C2e * UinPo - S2e * QinPo Volker@43: QinPae = C2e * QinPa + S2e * UinPa Volker@43: UinPae = C2e * UinPa - S2e * QinPa Volker@43: QinPe = C2e * QinP + S2e * UinP Volker@43: UinPe = C2e * UinP - S2e * QinP Volker@43: Volker@43: # Calibration signals and Calibration correction K from measurements with LDRCal / aCal Volker@43: if (TypeC == 2) or (TypeC == 1): # rotator calibration Eq. C.4 Volker@43: # parameters for calibration with aCal Volker@43: AT = ATP1 * IinP + h * ATP4 * VinP Volker@43: BT = ATP3e * QinP - h * ATP2e * UinP Volker@43: AR = ARP1 * IinP + h * ARP4 * VinP Volker@43: BR = ARP3e * QinP - h * ARP2e * UinP Volker@43: # Correction parameters for normal measurements; they are independent of LDR Volker@43: if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out Volker@43: IS1 = np.array([IinPo, QinPo, UinPo, VinPo]) Volker@43: IS2 = np.array([IinPa, QinPa, UinPa, VinPa]) Volker@43: GT = np.dot(ATP, IS1) Volker@43: GR = np.dot(ARP, IS1) Volker@43: HT = np.dot(ATP, IS2) Volker@43: HR = np.dot(ARP, IS2) Volker@43: else: Volker@43: IS1 = np.array([IinPo, QinPo, UinPo, VinPo]) Volker@43: IS2 = np.array([IinPa, QinPa, UinPa, VinPa]) Volker@43: GT = np.dot(ATPe, IS1) Volker@43: GR = np.dot(ARPe, IS1) Volker@43: HT = np.dot(ATPe, IS2) Volker@43: HR = np.dot(ARPe, IS2) Volker@43: elif (TypeC == 3) or (TypeC == 4): # linear polariser calibration Eq. C.5 Volker@43: # parameters for calibration with aCal Volker@43: AT = ATP1 * IinP + ATP3e * UinPe + ZiC * CosC * (ATP2e * QinPe + ATP4 * VinP) Volker@43: BT = DiC * (ATP1 * UinPe + ATP3e * IinP) - ZiC * SinC * (ATP2e * VinP - ATP4 * QinPe) Volker@43: AR = ARP1 * IinP + ARP3e * UinPe + ZiC * CosC * (ARP2e * QinPe + ARP4 * VinP) Volker@43: BR = DiC * (ARP1 * UinPe + ARP3e * IinP) - ZiC * SinC * (ARP2e * VinP - ARP4 * QinPe) Volker@43: # Correction parameters for normal measurements; they are independent of LDR Volker@43: if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out Volker@43: IS1 = np.array([IinPo, QinPo, UinPo, VinPo]) Volker@43: IS2 = np.array([IinPa, QinPa, UinPa, VinPa]) Volker@43: GT = np.dot(ATP, IS1) Volker@43: GR = np.dot(ARP, IS1) Volker@43: HT = np.dot(ATP, IS2) Volker@43: HR = np.dot(ARP, IS2) Volker@43: else: Volker@43: IS1e = np.array( Volker@43: [IinPo + DiC * QinPoe, DiC * IinPo + QinPoe, ZiC * (CosC * UinPoe + SinC * VinPo), Volker@43: -ZiC * (SinC * UinPoe - CosC * VinPo)]) Volker@43: IS2e = np.array( Volker@43: [IinPa + DiC * QinPae, DiC * IinPa + QinPae, ZiC * (CosC * UinPae + SinC * VinPa), Volker@43: -ZiC * (SinC * UinPae - CosC * VinPa)]) Volker@43: GT = np.dot(ATPe, IS1e) Volker@43: GR = np.dot(ARPe, IS1e) Volker@43: HT = np.dot(ATPe, IS2e) Volker@43: HR = np.dot(ARPe, IS2e) Volker@43: elif (TypeC == 6): # diattenuator calibration +-22.5° rotated_diattenuator_X22x5deg.odt Volker@43: # parameters for calibration with aCal Volker@43: AT = ATP1 * IinP + sqr05 * DiC * (ATP1 * QinPe + ATP2e * IinP) + (1 - 0.5 * WiC) * ( Volker@43: ATP2e * QinPe + ATP3e * UinPe) + ZiC * ( Volker@43: sqr05 * SinC * (ATP3e * VinP - ATP4 * UinPe) + ATP4 * CosC * VinP) Volker@43: BT = sqr05 * DiC * (ATP1 * UinPe + ATP3e * IinP) + 0.5 * WiC * ( Volker@43: ATP2e * UinPe + ATP3e * QinPe) - sqr05 * ZiC * SinC * (ATP2e * VinP - ATP4 * QinPe) Volker@43: AR = ARP1 * IinP + sqr05 * DiC * (ARP1 * QinPe + ARP2e * IinP) + (1 - 0.5 * WiC) * ( Volker@43: ARP2e * QinPe + ARP3e * UinPe) + ZiC * ( Volker@43: sqr05 * SinC * (ARP3e * VinP - ARP4 * UinPe) + ARP4 * CosC * VinP) Volker@43: BR = sqr05 * DiC * (ARP1 * UinPe + ARP3e * IinP) + 0.5 * WiC * ( Volker@43: ARP2e * UinPe + ARP3e * QinPe) - sqr05 * ZiC * SinC * (ARP2e * VinP - ARP4 * QinPe) Volker@43: # Correction parameters for normal measurements; they are independent of LDR Volker@43: if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out Volker@43: IS1 = np.array([IinPo, QinPo, UinPo, VinPo]) Volker@43: IS2 = np.array([IinPa, QinPa, UinPa, VinPa]) Volker@43: GT = np.dot(ATP, IS1) Volker@43: GR = np.dot(ARP, IS1) Volker@43: HT = np.dot(ATP, IS2) Volker@43: HR = np.dot(ARP, IS2) Volker@43: else: Volker@43: IS1e = np.array( Volker@43: [IinPo + DiC * QinPoe, DiC * IinPo + QinPoe, ZiC * (CosC * UinPoe + SinC * VinPo), Volker@43: -ZiC * (SinC * UinPoe - CosC * VinPo)]) Volker@43: IS2e = np.array( Volker@43: [IinPa + DiC * QinPae, DiC * IinPa + QinPae, ZiC * (CosC * UinPae + SinC * VinPa), Volker@43: -ZiC * (SinC * UinPae - CosC * VinPa)]) Volker@43: GT = np.dot(ATPe, IS1e) Volker@43: GR = np.dot(ARPe, IS1e) Volker@43: HT = np.dot(ATPe, IS2e) Volker@43: HR = np.dot(ARPe, IS2e) Volker@43: else: Volker@43: print("Calibrator not implemented yet") Volker@43: sys.exit() Volker@43: Volker@43: elif LocC == 3: # C before receiver optics Eq.57 Volker@43: Volker@43: # S2ge = np.sin(np.deg2rad(2*RotO - 2*RotC)) Volker@43: # C2ge = np.cos(np.deg2rad(2*RotO - 2*RotC)) Volker@43: S2e = np.sin(np.deg2rad(2 * RotC)) Volker@43: C2e = np.cos(np.deg2rad(2 * RotC)) Volker@43: Volker@43: # AS with C before the receiver optics (see document rotated_diattenuator_X22x5deg.odt) Volker@43: AF1 = np.array([1, C2g * DiO, S2g * DiO, 0]) Volker@43: AF2 = np.array([C2g * DiO, 1 - S2g ** 2 * WiO, S2g * C2g * WiO, -S2g * ZiO * SinO]) Volker@43: AF3 = np.array([S2g * DiO, S2g * C2g * WiO, 1 - C2g ** 2 * WiO, C2g * ZiO * SinO]) Volker@43: AF4 = np.array([0, S2g * SinO, -C2g * SinO, CosO]) Volker@43: Volker@43: ATF = (ATP1 * AF1 + ATP2 * AF2 + ATP3 * AF3 + ATP4 * AF4) Volker@43: ARF = (ARP1 * AF1 + ARP2 * AF2 + ARP3 * AF3 + ARP4 * AF4) Volker@43: ATF1 = ATF[0] Volker@43: ATF2 = ATF[1] Volker@43: ATF3 = ATF[2] Volker@43: ATF4 = ATF[3] Volker@43: ARF1 = ARF[0] Volker@43: ARF2 = ARF[1] Volker@43: ARF3 = ARF[2] Volker@43: ARF4 = ARF[3] Volker@43: Volker@43: # rotated AinF by epsilon Volker@43: ATF2e = C2e * ATF[1] + S2e * ATF[2] Volker@43: ATF3e = C2e * ATF[2] - S2e * ATF[1] Volker@43: ARF2e = C2e * ARF[1] + S2e * ARF[2] Volker@43: ARF3e = C2e * ARF[2] - S2e * ARF[1] Volker@43: Volker@43: ATFe = np.array([ATF1, ATF2e, ATF3e, ATF4]) Volker@43: ARFe = np.array([ARF1, ARF2e, ARF3e, ARF4]) Volker@43: Volker@43: QinEe = C2e * QinE + S2e * UinE Volker@43: UinEe = C2e * UinE - S2e * QinE Volker@43: Volker@43: # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) Volker@43: IinF = IinE Volker@43: QinF = aCal * QinE Volker@43: UinF = -aCal * UinE Volker@43: VinF = (1. - 2. * aCal) * VinE Volker@43: Volker@43: IinFo = IinE Volker@43: QinFo = 0. Volker@43: UinFo = 0. Volker@43: VinFo = VinE Volker@43: Volker@43: IinFa = 0. Volker@43: QinFa = QinE Volker@43: UinFa = -UinE Volker@43: VinFa = -2. * VinE Volker@43: Volker@43: # Stokes Input Vector before receiver optics rotated by epsilon Eq. C.3 Volker@43: QinFe = C2e * QinF + S2e * UinF Volker@43: UinFe = C2e * UinF - S2e * QinF Volker@43: QinFoe = C2e * QinFo + S2e * UinFo Volker@43: UinFoe = C2e * UinFo - S2e * QinFo Volker@43: QinFae = C2e * QinFa + S2e * UinFa Volker@43: UinFae = C2e * UinFa - S2e * QinFa Volker@43: Volker@43: # Calibration signals and Calibration correction K from measurements with LDRCal / aCal Volker@43: if (TypeC == 2) or (TypeC == 1): # rotator calibration Eq. C.4 Volker@43: AT = ATF1 * IinF + ATF4 * h * VinF Volker@43: BT = ATF3e * QinF - ATF2e * h * UinF Volker@43: AR = ARF1 * IinF + ARF4 * h * VinF Volker@43: BR = ARF3e * QinF - ARF2e * h * UinF Volker@43: Volker@43: # Correction parameters for normal measurements; they are independent of LDR Volker@43: if (not RotationErrorEpsilonForNormalMeasurements): Volker@43: GT = ATF1 * IinE + ATF4 * VinE Volker@43: GR = ARF1 * IinE + ARF4 * VinE Volker@43: HT = ATF2 * QinE - ATF3 * UinE - ATF4 * 2 * VinE Volker@43: HR = ARF2 * QinE - ARF3 * UinE - ARF4 * 2 * VinE Volker@43: else: Volker@43: GT = ATF1 * IinE + ATF4 * h * VinE Volker@43: GR = ARF1 * IinE + ARF4 * h * VinE Volker@43: HT = ATF2e * QinE - ATF3e * h * UinE - ATF4 * h * 2 * VinE Volker@43: HR = ARF2e * QinE - ARF3e * h * UinE - ARF4 * h * 2 * VinE Volker@43: Volker@43: elif (TypeC == 3) or (TypeC == 4): # linear polariser calibration Eq. C.5 Volker@43: # p = +45°, m = -45° Volker@43: IF1e = np.array([IinF, ZiC * CosC * QinFe, UinFe, ZiC * CosC * VinF]) Volker@43: IF2e = np.array([DiC * UinFe, -ZiC * SinC * VinF, DiC * IinF, ZiC * SinC * QinFe]) Volker@43: Volker@43: AT = np.dot(ATFe, IF1e) Volker@43: AR = np.dot(ARFe, IF1e) Volker@43: BT = np.dot(ATFe, IF2e) Volker@43: BR = np.dot(ARFe, IF2e) Volker@43: Volker@43: # Correction parameters for normal measurements; they are independent of LDR --- the same as for TypeC = 6 Volker@43: if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out Volker@43: IS1 = np.array([IinE, 0, 0, VinE]) Volker@43: IS2 = np.array([0, QinE, -UinE, -2 * VinE]) Volker@43: Volker@43: GT = np.dot(ATF, IS1) Volker@43: GR = np.dot(ARF, IS1) Volker@43: HT = np.dot(ATF, IS2) Volker@43: HR = np.dot(ARF, IS2) Volker@43: else: Volker@43: IS1e = np.array( Volker@43: [IinFo + DiC * QinFoe, DiC * IinFo + QinFoe, ZiC * (CosC * UinFoe + SinC * VinFo), Volker@43: -ZiC * (SinC * UinFoe - CosC * VinFo)]) Volker@43: IS2e = np.array( Volker@43: [IinFa + DiC * QinFae, DiC * IinFa + QinFae, ZiC * (CosC * UinFae + SinC * VinFa), Volker@43: -ZiC * (SinC * UinFae - CosC * VinFa)]) Volker@43: GT = np.dot(ATFe, IS1e) Volker@43: GR = np.dot(ARFe, IS1e) Volker@43: HT = np.dot(ATFe, IS2e) Volker@43: HR = np.dot(ARFe, IS2e) Volker@43: Volker@43: elif (TypeC == 6): # diattenuator calibration +-22.5° rotated_diattenuator_X22x5deg.odt Volker@43: # p = +22.5°, m = -22.5° Volker@43: IF1e = np.array([IinF + sqr05 * DiC * QinFe, sqr05 * DiC * IinF + (1 - 0.5 * WiC) * QinFe, Volker@43: (1 - 0.5 * WiC) * UinFe + sqr05 * ZiC * SinC * VinF, Volker@43: -sqr05 * ZiC * SinC * UinFe + ZiC * CosC * VinF]) Volker@43: IF2e = np.array([sqr05 * DiC * UinFe, 0.5 * WiC * UinFe - sqr05 * ZiC * SinC * VinF, Volker@43: sqr05 * DiC * IinF + 0.5 * WiC * QinFe, sqr05 * ZiC * SinC * QinFe]) Volker@43: Volker@43: AT = np.dot(ATFe, IF1e) Volker@43: AR = np.dot(ARFe, IF1e) Volker@43: BT = np.dot(ATFe, IF2e) Volker@43: BR = np.dot(ARFe, IF2e) Volker@43: Volker@43: # Correction parameters for normal measurements; they are independent of LDR Volker@43: if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out Volker@43: # IS1 = np.array([IinE,0,0,VinE]) Volker@43: # IS2 = np.array([0,QinE,-UinE,-2*VinE]) Volker@43: IS1 = np.array([IinFo, 0, 0, VinFo]) Volker@43: IS2 = np.array([0, QinFa, UinFa, VinFa]) Volker@43: GT = np.dot(ATF, IS1) Volker@43: GR = np.dot(ARF, IS1) Volker@43: HT = np.dot(ATF, IS2) Volker@43: HR = np.dot(ARF, IS2) Volker@43: else: Volker@43: # IS1e = np.array([IinE,DiC*IinE,ZiC*SinC*VinE,ZiC*CosC*VinE]) Volker@43: # IS2e = np.array([DiC*QinEe,QinEe,-ZiC*(CosC*UinEe+2*SinC*VinE),ZiC*(SinC*UinEe-2*CosC*VinE)]) Volker@43: IS1e = np.array( Volker@43: [IinFo + DiC * QinFoe, DiC * IinFo + QinFoe, ZiC * (CosC * UinFoe + SinC * VinFo), Volker@43: -ZiC * (SinC * UinFoe - CosC * VinFo)]) Volker@43: IS2e = np.array( Volker@43: [IinFa + DiC * QinFae, DiC * IinFa + QinFae, ZiC * (CosC * UinFae + SinC * VinFa), Volker@43: -ZiC * (SinC * UinFae - CosC * VinFa)]) Volker@43: GT = np.dot(ATFe, IS1e) Volker@43: GR = np.dot(ARFe, IS1e) Volker@43: HT = np.dot(ATFe, IS2e) Volker@43: HR = np.dot(ARFe, IS2e) Volker@43: Volker@43: Volker@43: else: Volker@43: print('Calibrator not implemented yet') Volker@43: sys.exit() Volker@43: Volker@43: elif LocC == 2: # C behind emitter optics Eq.57 Volker@43: # print("Calibrator location not implemented yet") Volker@43: S2e = np.sin(np.deg2rad(2 * RotC)) Volker@43: C2e = np.cos(np.deg2rad(2 * RotC)) Volker@43: Volker@43: # AS with C before the receiver optics (see document rotated_diattenuator_X22x5deg.odt) Volker@43: AF1 = np.array([1, C2g * DiO, S2g * DiO, 0]) Volker@43: AF2 = np.array([C2g * DiO, 1 - S2g ** 2 * WiO, S2g * C2g * WiO, -S2g * ZiO * SinO]) Volker@43: AF3 = np.array([S2g * DiO, S2g * C2g * WiO, 1 - C2g ** 2 * WiO, C2g * ZiO * SinO]) Volker@43: AF4 = np.array([0, S2g * SinO, -C2g * SinO, CosO]) Volker@43: Volker@43: ATF = (ATP1 * AF1 + ATP2 * AF2 + ATP3 * AF3 + ATP4 * AF4) Volker@43: ARF = (ARP1 * AF1 + ARP2 * AF2 + ARP3 * AF3 + ARP4 * AF4) Volker@43: ATF1 = ATF[0] Volker@43: ATF2 = ATF[1] Volker@43: ATF3 = ATF[2] Volker@43: ATF4 = ATF[3] Volker@43: ARF1 = ARF[0] Volker@43: ARF2 = ARF[1] Volker@43: ARF3 = ARF[2] Volker@43: ARF4 = ARF[3] Volker@43: Volker@43: # AS with C behind the emitter -------------------------------------------- Volker@43: # terms without aCal Volker@43: ATE1o, ARE1o = ATF1, ARF1 Volker@43: ATE2o, ARE2o = 0., 0. Volker@43: ATE3o, ARE3o = 0., 0. Volker@43: ATE4o, ARE4o = ATF4, ARF4 Volker@43: # terms with aCal Volker@43: ATE1a, ARE1a = 0., 0. Volker@43: ATE2a, ARE2a = ATF2, ARF2 Volker@43: ATE3a, ARE3a = -ATF3, -ARF3 Volker@43: ATE4a, ARE4a = -2 * ATF4, -2 * ARF4 Volker@43: # rotated AinEa by epsilon Volker@43: ATE2ae = C2e * ATF2 + S2e * ATF3 Volker@43: ATE3ae = -S2e * ATF2 - C2e * ATF3 Volker@43: ARE2ae = C2e * ARF2 + S2e * ARF3 Volker@43: ARE3ae = -S2e * ARF2 - C2e * ARF3 Volker@43: Volker@43: ATE1 = ATE1o Volker@43: ATE2e = aCal * ATE2ae Volker@43: ATE3e = aCal * ATE3ae Volker@43: ATE4 = (1 - 2 * aCal) * ATF4 Volker@43: ARE1 = ARE1o Volker@43: ARE2e = aCal * ARE2ae Volker@43: ARE3e = aCal * ARE3ae Volker@43: ARE4 = (1. - 2. * aCal) * ARF4 Volker@43: Volker@43: # rotated IinE Volker@43: QinEe = C2e * QinE + S2e * UinE Volker@43: UinEe = C2e * UinE - S2e * QinE Volker@43: Volker@43: # --- Calibration signals and Calibration correction K from measurements with LDRCal / aCal Volker@43: if (TypeC == 2) or (TypeC == 1): # +++++++++ rotator calibration Eq. C.4 Volker@43: AT = ATE1o * IinE + (ATE4o + aCal * ATE4a) * h * VinE Volker@43: BT = aCal * (ATE3ae * QinEe - ATE2ae * h * UinEe) Volker@43: AR = ARE1o * IinE + (ARE4o + aCal * ARE4a) * h * VinE Volker@43: BR = aCal * (ARE3ae * QinEe - ARE2ae * h * UinEe) Volker@43: Volker@43: # Correction parameters for normal measurements; they are independent of LDR Volker@43: if (not RotationErrorEpsilonForNormalMeasurements): Volker@43: # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) Volker@43: GT = ATE1o * IinE + ATE4o * h * VinE Volker@43: GR = ARE1o * IinE + ARE4o * h * VinE Volker@43: HT = ATE2a * QinE + ATE3a * h * UinEe + ATE4a * h * VinE Volker@43: HR = ARE2a * QinE + ARE3a * h * UinEe + ARE4a * h * VinE Volker@43: else: Volker@43: GT = ATE1o * IinE + ATE4o * h * VinE Volker@43: GR = ARE1o * IinE + ARE4o * h * VinE Volker@43: HT = ATE2ae * QinE + ATE3ae * h * UinEe + ATE4a * h * VinE Volker@43: HR = ARE2ae * QinE + ARE3ae * h * UinEe + ARE4a * h * VinE Volker@43: Volker@43: elif (TypeC == 3) or (TypeC == 4): # +++++++++ linear polariser calibration Eq. C.5 Volker@43: # p = +45°, m = -45° Volker@43: AT = ATE1 * IinE + ZiC * CosC * (ATE2e * QinEe + ATE4 * VinE) + ATE3e * UinEe Volker@43: BT = DiC * (ATE1 * UinEe + ATE3e * IinE) + ZiC * SinC * (ATE4 * QinEe - ATE2e * VinE) Volker@43: AR = ARE1 * IinE + ZiC * CosC * (ARE2e * QinEe + ARE4 * VinE) + ARE3e * UinEe Volker@43: BR = DiC * (ARE1 * UinEe + ARE3e * IinE) + ZiC * SinC * (ARE4 * QinEe - ARE2e * VinE) Volker@43: Volker@43: # Correction parameters for normal measurements; they are independent of LDR Volker@43: if (not RotationErrorEpsilonForNormalMeasurements): Volker@43: # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) Volker@43: GT = ATE1o * IinE + ATE4o * VinE Volker@43: GR = ARE1o * IinE + ARE4o * VinE Volker@43: HT = ATE2a * QinE + ATE3a * UinE + ATE4a * VinE Volker@43: HR = ARE2a * QinE + ARE3a * UinE + ARE4a * VinE Volker@43: else: Volker@43: D = IinE + DiC * QinEe Volker@43: A = DiC * IinE + QinEe Volker@43: B = ZiC * (CosC * UinEe + SinC * VinE) Volker@43: C = -ZiC * (SinC * UinEe - CosC * VinE) Volker@43: GT = ATE1o * D + ATE4o * C Volker@43: GR = ARE1o * D + ARE4o * C Volker@43: HT = ATE2a * A + ATE3a * B + ATE4a * C Volker@43: HR = ARE2a * A + ARE3a * B + ARE4a * C Volker@43: Volker@43: elif (TypeC == 6): # real HWP calibration +-22.5° rotated_diattenuator_X22x5deg.odt Volker@43: # p = +22.5°, m = -22.5° Volker@43: IE1e = np.array([IinE + sqr05 * DiC * QinEe, sqr05 * DiC * IinE + (1 - 0.5 * WiC) * QinEe, Volker@43: (1. - 0.5 * WiC) * UinEe + sqr05 * ZiC * SinC * VinE, Volker@43: -sqr05 * ZiC * SinC * UinEe + ZiC * CosC * VinE]) Volker@43: IE2e = np.array([sqr05 * DiC * UinEe, 0.5 * WiC * UinEe - sqr05 * ZiC * SinC * VinE, Volker@43: sqr05 * DiC * IinE + 0.5 * WiC * QinEe, sqr05 * ZiC * SinC * QinEe]) Volker@43: ATEe = np.array([ATE1, ATE2e, ATE3e, ATE4]) Volker@43: AREe = np.array([ARE1, ARE2e, ARE3e, ARE4]) Volker@43: AT = np.dot(ATEe, IE1e) Volker@43: AR = np.dot(AREe, IE1e) Volker@43: BT = np.dot(ATEe, IE2e) Volker@43: BR = np.dot(AREe, IE2e) Volker@43: Volker@43: # Correction parameters for normal measurements; they are independent of LDR Volker@43: if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out Volker@43: GT = ATE1o * IinE + ATE4o * VinE Volker@43: GR = ARE1o * IinE + ARE4o * VinE Volker@43: HT = ATE2a * QinE + ATE3a * UinE + ATE4a * VinE Volker@43: HR = ARE2a * QinE + ARE3a * UinE + ARE4a * VinE Volker@43: else: Volker@43: D = IinE + DiC * QinEe Volker@43: A = DiC * IinE + QinEe Volker@43: B = ZiC * (CosC * UinEe + SinC * VinE) Volker@43: C = -ZiC * (SinC * UinEe - CosC * VinE) Volker@43: GT = ATE1o * D + ATE4o * C Volker@43: GR = ARE1o * D + ARE4o * C Volker@43: HT = ATE2a * A + ATE3a * B + ATE4a * C Volker@43: HR = ARE2a * A + ARE3a * B + ARE4a * C Volker@43: else: Volker@43: print('Calibrator not implemented yet') Volker@43: sys.exit() Volker@43: Volker@43: for iTCalT, iTCalR, iNCalTp, iNCalTm, iNCalRp, iNCalRm, iNIt, iNIr \ Volker@43: in [ Volker@43: (iTCalT, iTCalR, iNCalTp, iNCalTm, iNCalRp, iNCalRm, iNIt, iNIr) Volker@43: for iTCalT in range(-nTCalT, nTCalT + 1) # Etax Volker@43: for iTCalR in range(-nTCalR, nTCalR + 1) # Etax Volker@43: for iNCalTp in range(-nNCal, nNCal + 1) # noise error of calibration signals => Etax Volker@43: for iNCalTm in range(-nNCal, nNCal + 1) # noise error of calibration signals => Etax Volker@43: for iNCalRp in range(-nNCal, nNCal + 1) # noise error of calibration signals => Etax Volker@43: for iNCalRm in range(-nNCal, nNCal + 1) # noise error of calibration signals => Etax Volker@43: for iNIt in range(-nNI, nNI + 1) Volker@43: for iNIr in range(-nNI, nNI + 1)]: Volker@43: Volker@43: # Calibration signals with aCal => Determination of the correction K of the real calibration factor Volker@43: IoutTp = TTa * TiC * TiO * TiE * (AT + BT) Volker@43: IoutTm = TTa * TiC * TiO * TiE * (AT - BT) Volker@43: IoutRp = TRa * TiC * TiO * TiE * (AR + BR) Volker@43: IoutRm = TRa * TiC * TiO * TiE * (AR - BR) Volker@43: Volker@43: if nTCalT > 0: TCalT = TCalT0 + iTCalT * dTCalT / nTCalT Volker@43: if nTCalR > 0: TCalR = TCalR0 + iTCalR * dTCalR / nTCalR Volker@43: # signal noise errors Volker@43: # ----- random error calculation ---------- Volker@43: # noise must be calculated from/with the actually measured signals; influence of TCalT, TCalR errors on noise are not considered ? Volker@43: # actually measured signal counts are in input file and don't change Volker@43: # relative standard deviation of calibration signals with LDRcal; assumed to be statisitcally independent Volker@43: # error nNCal: one-sigma in steps to left and right for calibration signals Volker@43: if nNCal > 0: Volker@43: if (CalcFrom0deg): Volker@43: dIoutTp = (NCalT * IoutTp) ** -0.5 Volker@43: dIoutTm = (NCalT * IoutTm) ** -0.5 Volker@43: dIoutRp = (NCalR * IoutRp) ** -0.5 Volker@43: dIoutRm = (NCalR * IoutRm) ** -0.5 Volker@43: else: Volker@43: dIoutTp = dIoutTp0 * (IoutTp / IoutTp0) Volker@43: dIoutTm = dIoutTm0 * (IoutTm / IoutTm0) Volker@43: dIoutRp = dIoutRp0 * (IoutRp / IoutRp0) Volker@43: dIoutRm = dIoutRm0 * (IoutRm / IoutRm0) Volker@43: # print(iTCalT, iTCalR, iNCalTp, iNCalTm, iNCalRp, iNCalRm, iNIt, iNIr, IoutTp, dIoutTp) Volker@43: IoutTp = IoutTp * (1. + iNCalTp * dIoutTp / nNCal) Volker@43: IoutTm = IoutTm * (1. + iNCalTm * dIoutTm / nNCal) Volker@43: IoutRp = IoutRp * (1. + iNCalRp * dIoutRp / nNCal) Volker@43: IoutRm = IoutRm * (1. + iNCalRm * dIoutRm / nNCal) Volker@43: Volker@43: IoutTp = IoutTp * TCalT / TCalT0 Volker@43: IoutTm = IoutTm * TCalT / TCalT0 Volker@43: IoutRp = IoutRp * TCalR / TCalR0 Volker@43: IoutRm = IoutRm * TCalR / TCalR0 Volker@43: # --- Results and Corrections; electronic etaR and etaT are assumed to be 1 for true and assumed true systems Volker@43: # calibration factor Volker@43: Eta = (TRa / TTa) # = TRa / TTa; Eta = Eta*/K Eq. 84; corrected according to the papers supplement Eqs. (S.10.10.1) ff Volker@43: # possibly real calibration factor Volker@43: Etapx = IoutRp / IoutTp Volker@43: Etamx = IoutRm / IoutTm Volker@43: Etax = (Etapx * Etamx) ** 0.5 Volker@43: K = Etax / Eta Volker@43: # 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@43: # print("{0:6.3f},{1:6.3f},{2:6.3f},{3:6.3f}".format(DiC, ZiC, Kp, Km)) Volker@43: Volker@43: # For comparison with Volkers Libreoffice Müller Matrix spreadsheet Volker@43: # Eta_test_p = (IoutRp/IoutTp) Volker@43: # Eta_test_m = (IoutRm/IoutTm) Volker@43: # Eta_test = (Eta_test_p*Eta_test_m)**0.5 Volker@43: ''' Volker@43: for iIt, iIr \ Volker@43: in [(iIt, iIr) Volker@43: for iIt in range(-nNI, nNI + 1) Volker@43: for iIr in range(-nNI, nNI + 1)]: Volker@43: ''' Volker@43: Volker@43: iN = iN + 1 Volker@43: if (iN == 10001): Volker@43: ctime = clock() Volker@43: print(" estimated time ", "{0:4.2f}".format(N / 10000 * (ctime - atime)), "sec ") # , end="") Volker@43: print("\r elapsed time ", "{0:5.0f}".format((ctime - atime)), "sec ", end="\r") Volker@43: ctime = clock() Volker@43: if ((ctime - dtime) > 10): Volker@43: print("\r elapsed time ", "{0:5.0f}".format((ctime - atime)), "sec ", end="\r") Volker@43: dtime = ctime Volker@43: Volker@43: # *** loop for different real LDRs ********************************************************************** Volker@43: iLDR = -1 Volker@43: for LDRTrue in LDRrange: Volker@43: iLDR = iLDR + 1 Volker@43: atrue = (1. - LDRTrue) / (1. + LDRTrue) Volker@43: # ----- Forward simulated signals and LDRsim with atrue; from input file; not considering TiC. Volker@43: It = TTa * TiO * TiE * (GT + atrue * HT) # TaT*TiT*TiC*TiO*IinL*(GT+atrue*HT) Volker@43: Ir = TRa * TiO * TiE * (GR + atrue * HR) # TaR*TiR*TiC*TiO*IinL*(GR+atrue*HR) Volker@43: # # signal noise errors; standard deviation of signals; assumed to be statisitcally independent Volker@43: # because the signals depend on LDRtrue, the errors dIt and dIr must be calculated for each LDRtrue Volker@43: if (CalcFrom0deg): Volker@43: ''' Volker@43: dIt = ((NCalT * It / IoutTp * NILfac / TCalT) ** -0.5) Volker@43: dIr = ((NCalR * Ir / IoutRp * NILfac / TCalR) ** -0.5) Volker@43: ''' Volker@43: dIt = ((It * NI * eFacT) ** -0.5) Volker@43: dIr = ((Ir * NI * eFacR) ** -0.5) Volker@43: else: Volker@43: dIt = ((It * NI * eFacT) ** -0.5) Volker@43: dIr = ((Ir * NI * eFacR) ** -0.5) Volker@43: ''' Volker@43: # does this work? Why not as above? Volker@43: dIt = ((NCalT * 2. * NILfac / TCalT ) ** -0.5) Volker@43: dIr = ((NCalR * 2. * NILfac / TCalR) ** -0.5) Volker@43: ''' Volker@43: # error nNI: one-sigma in steps to left and right for 0° signals Volker@43: if nNI > 0: Volker@43: It = It * (1. + iNIt * dIt / nNI) Volker@43: Ir = Ir * (1. + iNIr * dIr / nNI) Volker@43: Volker@43: # LDRsim = 1/Eta*Ir/It # simulated LDR* with Y from input file Volker@43: LDRsim = Ir / It # simulated uncorrected LDR with Y from input file Volker@43: Volker@43: # ----- Backward correction Volker@43: # Corrected LDRCorr with assumed true G0,H0,K0,Eta0 from forward simulated (real) LDRsim(atrue) Volker@43: LDRCorr = (LDRsim / (Etax / K0) * (GT0 + HT0) - (GR0 + HR0)) / ((GR0 - HR0) - LDRsim / (Etax / K0) * (GT0 - HT0)) Volker@43: Volker@43: # The following is a test whether the equations for calibration Etax and normal signal (GHK, LDRsim) are consistent Volker@43: # LDRCorr = (LDRsim / Eta * (GT + HT) - (GR + HR)) / ((GR - HR) - LDRsim / Eta * (GT - HT)) Volker@43: # Without any correction Volker@43: LDRunCorr = LDRsim / Etax Volker@43: # 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@43: Volker@43: Volker@43: ''' Volker@43: # -- F11corr from It and Ir and calibration EtaX Volker@43: Text1 = "!!! EXPERIMENTAL !!! F11corr from It and Ir with calibration EtaX: x-axis: F11corr(LDRtrue) / F11corr(LDRtrue = 0.004) - 1" Volker@43: F11corr = 1 / (TiO * TiE) * ( Volker@43: (HR0 * Etax / K0 * It / TTa - HT0 * Ir / TRa) / (HR0 * GT0 - HT0 * GR0)) # IL = 1 Eq.(64); Etax/K0 = Eta0. Volker@43: ''' Volker@43: # Corrected F11corr with assumed true G0,H0,K0 from forward simulated (real) It and Ir (atrue) Volker@43: Text1 = "!!! EXPERIMENTAL !!! F11corr from real It and Ir with real calibration EtaX: x-axis: F11corr(LDRtrue) / aF11sim0(LDRtrue) - 1" Volker@43: F11corr = 1 / (TiO * TiE) * ( Volker@43: (HR0 * Etax / K0 * It / TTa - HT0 * Ir / TRa) / (HR0 * GT0 - HT0 * GR0)) # IL = 1 Eq.(64); Etax/K0 = Eta0. Volker@43: Volker@43: # Text1 = "F11corr from It and Ir without corrections but with calibration EtaX: x-axis: F11corr(LDRtrue) devided by F11corr(LDRtrue = 0.004)" Volker@43: # F11corr = 0.5/(TiO*TiE)*(Etax*It/TTa+Ir/TRa) # IL = 1 Eq.(64) Volker@43: Volker@43: # -- It from It only with atrue without corrections - for BERTHA (and PollyXTs) Volker@43: # Text1 = " x-axis: IT(LDRtrue) / IT(LDRtrue = 0.004) - 1" Volker@43: # F11corr = It/(TaT*TiT*TiO*TiE) #/(TaT*TiT*TiO*TiE*(GT0+atrue*HT0)) Volker@43: # ! see below line 1673ff Volker@43: Volker@43: aF11corr[iLDR, iN] = F11corr Volker@43: aLDRcorr[iLDR, iN] = LDRCorr # LDRCorr # LDRsim # for test only Volker@43: aLDRsim[iLDR, iN] = LDRsim # LDRCorr # LDRsim # for test only Volker@43: # aPLDR[iLDR, iN] = CalcPLDR(LDRCorr, BSR[iLDR], LDRm0) Volker@43: aEtax[iLDR, iN] = Etax Volker@43: aEtapx[iLDR, iN] = Etapx Volker@43: aEtamx[iLDR, iN] = Etamx Volker@43: Volker@43: aGHK[0, iN] = GR Volker@43: aGHK[1, iN] = GT Volker@43: aGHK[2, iN] = HR Volker@43: aGHK[3, iN] = HT Volker@43: aGHK[4, iN] = K Volker@43: Volker@43: aLDRCal[iN] = iLDRCal Volker@43: aQin[iN] = iQin Volker@43: aVin[iN] = iVin Volker@43: aERaT[iN] = iERaT Volker@43: aERaR[iN] = iERaR Volker@43: aRotaT[iN] = iRotaT Volker@43: aRotaR[iN] = iRotaR Volker@43: aRetT[iN] = iRetT Volker@43: aRetR[iN] = iRetR Volker@43: Volker@43: aRotL[iN] = iRotL Volker@43: aRotE[iN] = iRotE Volker@43: aRetE[iN] = iRetE Volker@43: aRotO[iN] = iRotO Volker@43: aRetO[iN] = iRetO Volker@43: aRotC[iN] = iRotC Volker@43: aRetC[iN] = iRetC Volker@43: aDiO[iN] = iDiO Volker@43: aDiE[iN] = iDiE Volker@43: aDiC[iN] = iDiC Volker@43: aTP[iN] = iTP Volker@43: aTS[iN] = iTS Volker@43: aRP[iN] = iRP Volker@43: aRS[iN] = iRS Volker@43: aTCalT[iN] = iTCalT Volker@43: aTCalR[iN] = iTCalR Volker@43: Volker@43: aNCalTp[iN] = iNCalTp # IoutTp, IoutTm, IoutRp, IoutRm => Etax Volker@43: aNCalTm[iN] = iNCalTm # IoutTp, IoutTm, IoutRp, IoutRm => Etax Volker@43: aNCalRp[iN] = iNCalRp # IoutTp, IoutTm, IoutRp, IoutRm => Etax Volker@43: aNCalRm[iN] = iNCalRm # IoutTp, IoutTm, IoutRp, IoutRm => Etax Volker@43: aNIt[iN] = iNIt # It, Tr Volker@43: aNIr[iN] = iNIr # It, Tr Volker@43: Volker@43: # --- END loop Volker@43: btime = clock() Volker@43: # print("\r done in ", "{0:5.0f}".format(btime - atime), "sec. => producing plots now .... some more seconds ..."), # , end="\r"); Volker@43: print(" done in ", "{0:5.0f}".format(btime - atime), "sec. => producing plots now .... some more seconds ...") Volker@43: # --- Plot ----------------------------------------------------------------- Volker@43: print("Errors from GHK correction uncertainties:") Volker@43: if (sns_loaded): Volker@43: sns.set_style("whitegrid") Volker@43: sns.set_palette("bright6", 6) Volker@43: # for older seaborn versions use: Volker@43: # sns.set_palette("bright", 6) Volker@43: Volker@43: ''' Volker@43: fig2 = plt.figure() Volker@43: plt.plot(aLDRcorr[2,:],'b.') Volker@43: plt.plot(aLDRcorr[3,:],'r.') Volker@43: plt.plot(aLDRcorr[4,:],'g.') Volker@43: #plt.plot(aLDRcorr[6,:],'c.') Volker@43: plt.show Volker@43: ''' Volker@43: Volker@43: # Plot LDR Volker@43: def PlotSubHist(aVar, aX, X0, daX, iaX, naX): Volker@43: # aVar is the name of the parameter and aX is the subset of aLDRcorr which is coloured in the plot Volker@43: # example: PlotSubHist("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP) Volker@43: fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) Volker@43: iLDR = -1 Volker@43: for LDRTrue in LDRrange: Volker@43: aXmean = np.zeros(2 * naX + 1) Volker@43: iLDR = iLDR + 1 Volker@43: LDRmin[iLDR] = np.amin(aLDRcorr[iLDR, :]) Volker@43: LDRmax[iLDR] = np.amax(aLDRcorr[iLDR, :]) Volker@43: if (LDRmax[iLDR] > 10): LDRmax[iLDR] = 10 Volker@43: if (LDRmin[iLDR] < -10): LDRmin[iLDR] = -10 Volker@43: Rmin = LDRmin[iLDR] * 0.995 # np.min(aLDRcorr[iLDR,:]) * 0.995 Volker@43: Rmax = LDRmax[iLDR] * 1.005 # np.max(aLDRcorr[iLDR,:]) * 1.005 Volker@43: Volker@43: # Determine mean distance of all aXmean from each other for each iLDR Volker@43: meanDist = 0.0 Volker@43: for iaX in range(-naX, naX + 1): Volker@43: # mean LDRCorr value for certain error (iaX) of parameter aVar Volker@43: aXmean[iaX + naX] = np.mean(aLDRcorr[iLDR, aX == iaX]) Volker@43: # relative to absolute spread of LDRCorrs Volker@43: meanDist = (np.max(aXmean) - np.min(aXmean)) / (LDRmax[iLDR] - LDRmin[iLDR]) * 100 Volker@43: Volker@43: plt.subplot(1, 5, iLDR + 1) Volker@43: (n, bins, patches) = plt.hist(aLDRcorr[iLDR, :], Volker@43: bins=100, log=False, Volker@43: range=[Rmin, Rmax], Volker@43: alpha=0.5, density=False, color='0.5', histtype='stepfilled') Volker@43: Volker@43: for iaX in range(-naX, naX + 1): Volker@43: # mean LDRCorr value for certain error (iaX) of parameter aVar Volker@43: plt.hist(aLDRcorr[iLDR, aX == iaX], Volker@43: range=[Rmin, Rmax], Volker@43: bins=100, log=False, alpha=0.3, density=False, histtype='stepfilled', Volker@43: label=str(round(X0 + iaX * daX / naX, 5))) Volker@43: Volker@43: if (iLDR == 2): Volker@43: leg = plt.legend() Volker@43: leg.get_frame().set_alpha(0.1) Volker@43: Volker@43: plt.tick_params(axis='both', labelsize=10) Volker@43: plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2) Volker@43: plt.gca().set_title("{0:3.0f}%".format(meanDist)) Volker@43: plt.gca().set_xlabel('LDRtrue', color="red") Volker@43: Volker@43: # plt.ylabel('frequency', fontsize=10) Volker@43: # plt.xlabel('LDRCorr', fontsize=10) Volker@43: # fig.tight_layout() Volker@43: fig.suptitle(LID + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar + ' error contribution', fontsize=14, y=1.10) Volker@43: # plt.show() Volker@43: # fig.savefig(LID + '_' + aVar + '.png', dpi=150, bbox_inches='tight', pad_inches=0) Volker@43: # plt.close Volker@43: return Volker@43: Volker@43: def PlotLDRsim(aVar, aX, X0, daX, iaX, naX): Volker@43: # aVar is the name of the parameter and aX is the subset of aLDRsim which is coloured in the plot Volker@43: # example: PlotSubHist("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP) Volker@43: fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) Volker@43: iLDR = -1 Volker@43: for LDRTrue in LDRrange: Volker@43: aXmean = np.zeros(2 * naX + 1) Volker@43: iLDR = iLDR + 1 Volker@43: LDRsimmin[iLDR] = np.amin(aLDRsim[iLDR, :]) Volker@43: LDRsimmax[iLDR] = np.amax(aLDRsim[iLDR, :]) Volker@43: # print("LDRsimmin[iLDR], LDRsimmax[iLDR] = ", LDRsimmin[iLDR], LDRsimmax[iLDR]) Volker@43: # if (LDRsimmax[iLDR] > 10): LDRsimmax[iLDR] = 10 Volker@43: # if (LDRsimmin[iLDR] < -10): LDRsimmin[iLDR] = -10 Volker@43: Rmin = LDRsimmin[iLDR] * 0.995 # np.min(aLDRsim[iLDR,:]) * 0.995 Volker@43: Rmax = LDRsimmax[iLDR] * 1.005 # np.max(aLDRsim[iLDR,:]) * 1.005 Volker@43: # print("Rmin, Rmax = ", Rmin, Rmax) Volker@43: Volker@43: # Determine mean distance of all aXmean from each other for each iLDR Volker@43: meanDist = 0.0 Volker@43: for iaX in range(-naX, naX + 1): Volker@43: # mean LDRCorr value for certain error (iaX) of parameter aVar Volker@43: aXmean[iaX + naX] = np.mean(aLDRsim[iLDR, aX == iaX]) Volker@43: # relative to absolute spread of LDRCorrs Volker@43: meanDist = (np.max(aXmean) - np.min(aXmean)) / (LDRsimmax[iLDR] - LDRsimmin[iLDR]) * 100 Volker@43: Volker@43: plt.subplot(1, 5, iLDR + 1) Volker@43: (n, bins, patches) = plt.hist(aLDRsim[iLDR, :], Volker@43: bins=100, log=False, Volker@43: range=[Rmin, Rmax], Volker@43: alpha=0.5, density=False, color='0.5', histtype='stepfilled') Volker@43: Volker@43: for iaX in range(-naX, naX + 1): Volker@43: # mean LDRCorr value for certain error (iaX) of parameter aVar Volker@43: plt.hist(aLDRsim[iLDR, aX == iaX], Volker@43: range=[Rmin, Rmax], Volker@43: bins=100, log=False, alpha=0.3, density=False, histtype='stepfilled', Volker@43: label=str(round(X0 + iaX * daX / naX, 5))) Volker@43: Volker@43: if (iLDR == 2): Volker@43: leg = plt.legend() Volker@43: leg.get_frame().set_alpha(0.1) Volker@43: Volker@43: plt.tick_params(axis='both', labelsize=10) Volker@43: plt.plot([LDRsim0[iLDR], LDRsim0[iLDR]], [0, np.max(n)], 'r-', lw=2) Volker@43: plt.gca().set_title("{0:3.0f}%".format(meanDist)) Volker@43: plt.gca().set_xlabel('LDRsim0', color="red") Volker@43: Volker@43: fig.suptitle('LDRsim - ' +LID + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar + ' error contribution', fontsize=14, y=1.10) Volker@43: return Volker@43: Volker@43: Volker@43: # Plot Etax Volker@43: def PlotEtax(aVar, aX, X0, daX, iaX, naX): Volker@43: # aVar is the name of the parameter and aX is the subset of aLDRcorr which is coloured in the plot Volker@43: # example: PlotSubHist("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP) Volker@43: fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) Volker@43: iLDR = -1 Volker@43: for LDRTrue in LDRrange: Volker@43: aXmean = np.zeros(2 * naX + 1) Volker@43: iLDR = iLDR + 1 Volker@43: Etaxmin = np.amin(aEtax[iLDR, :]) Volker@43: Etaxmax = np.amax(aEtax[iLDR, :]) Volker@43: Rmin = Etaxmin * 0.995 # np.min(aLDRcorr[iLDR,:]) * 0.995 Volker@43: Rmax = Etaxmax * 1.005 # np.max(aLDRcorr[iLDR,:]) * 1.005 Volker@43: Volker@43: # Determine mean distance of all aXmean from each other for each iLDR Volker@43: meanDist = 0.0 Volker@43: for iaX in range(-naX, naX + 1): Volker@43: # mean Etax value for certain error (iaX) of parameter aVar Volker@43: aXmean[iaX + naX] = np.mean(aEtax[iLDR, aX == iaX]) Volker@43: # relative to absolute spread of Etax Volker@43: meanDist = (np.max(aXmean) - np.min(aXmean)) / (Etaxmax - Etaxmin) * 100 Volker@43: Volker@43: plt.subplot(1, 5, iLDR + 1) Volker@43: (n, bins, patches) = plt.hist(aEtax[iLDR, :], Volker@43: bins=50, log=False, Volker@43: range=[Rmin, Rmax], Volker@43: alpha=0.5, density=False, color='0.5', histtype='stepfilled') Volker@43: for iaX in range(-naX, naX + 1): Volker@43: plt.hist(aEtax[iLDR, aX == iaX], Volker@43: range=[Rmin, Rmax], Volker@43: bins=50, log=False, alpha=0.3, density=False, histtype='stepfilled', Volker@43: label=str(round(X0 + iaX * daX / naX, 5))) Volker@43: if (iLDR == 2): Volker@43: leg = plt.legend() Volker@43: leg.get_frame().set_alpha(0.1) Volker@43: plt.tick_params(axis='both', labelsize=10) Volker@43: plt.plot([Etax0, Etax0], [0, np.max(n)], 'r-', lw=2) Volker@43: plt.gca().set_title("{0:3.0f}%".format(meanDist)) Volker@43: plt.gca().set_xlabel('Etax0', color="red") Volker@43: fig.suptitle('Etax - ' + LID + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar + ' error contribution', fontsize=14, y=1.10) Volker@43: return Volker@43: Volker@43: def PlotEtapx(aVar, aX, X0, daX, iaX, naX): Volker@43: # aVar is the name of the parameter and aX is the subset of aLDRcorr which is coloured in the plot Volker@43: # example: PlotSubHist("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP) Volker@43: fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) Volker@43: iLDR = -1 Volker@43: for LDRTrue in LDRrange: Volker@43: aXmean = np.zeros(2 * naX + 1) Volker@43: iLDR = iLDR + 1 Volker@43: Etapxmin = np.amin(aEtapx[iLDR, :]) Volker@43: Etapxmax = np.amax(aEtapx[iLDR, :]) Volker@43: Rmin = Etapxmin * 0.995 # np.min(aLDRcorr[iLDR,:]) * 0.995 Volker@43: Rmax = Etapxmax * 1.005 # np.max(aLDRcorr[iLDR,:]) * 1.005 Volker@43: Volker@43: # Determine mean distance of all aXmean from each other for each iLDR Volker@43: meanDist = 0.0 Volker@43: for iaX in range(-naX, naX + 1): Volker@43: # mean Etapx value for certain error (iaX) of parameter aVar Volker@43: aXmean[iaX + naX] = np.mean(aEtapx[iLDR, aX == iaX]) Volker@43: # relative to absolute spread of Etapx Volker@43: meanDist = (np.max(aXmean) - np.min(aXmean)) / (Etapxmax - Etapxmin) * 100 Volker@43: Volker@43: plt.subplot(1, 5, iLDR + 1) Volker@43: (n, bins, patches) = plt.hist(aEtapx[iLDR, :], Volker@43: bins=50, log=False, Volker@43: range=[Rmin, Rmax], Volker@43: alpha=0.5, density=False, color='0.5', histtype='stepfilled') Volker@43: for iaX in range(-naX, naX + 1): Volker@43: plt.hist(aEtapx[iLDR, aX == iaX], Volker@43: range=[Rmin, Rmax], Volker@43: bins=50, log=False, alpha=0.3, density=False, histtype='stepfilled', Volker@43: label=str(round(X0 + iaX * daX / naX, 5))) Volker@43: if (iLDR == 2): Volker@43: leg = plt.legend() Volker@43: leg.get_frame().set_alpha(0.1) Volker@43: plt.tick_params(axis='both', labelsize=10) Volker@43: plt.plot([Etapx0, Etapx0], [0, np.max(n)], 'r-', lw=2) Volker@43: plt.gca().set_title("{0:3.0f}%".format(meanDist)) Volker@43: plt.gca().set_xlabel('Etapx0', color="red") Volker@43: fig.suptitle('Etapx - ' + LID + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar + ' error contribution', fontsize=14, y=1.10) Volker@43: return Volker@43: Volker@43: def PlotEtamx(aVar, aX, X0, daX, iaX, naX): Volker@43: # aVar is the name of the parameter and aX is the subset of aLDRcorr which is coloured in the plot Volker@43: # example: PlotSubHist("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP) Volker@43: fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) Volker@43: iLDR = -1 Volker@43: for LDRTrue in LDRrange: Volker@43: aXmean = np.zeros(2 * naX + 1) Volker@43: iLDR = iLDR + 1 Volker@43: Etamxmin = np.amin(aEtamx[iLDR, :]) Volker@43: Etamxmax = np.amax(aEtamx[iLDR, :]) Volker@43: Rmin = Etamxmin * 0.995 # np.min(aLDRcorr[iLDR,:]) * 0.995 Volker@43: Rmax = Etamxmax * 1.005 # np.max(aLDRcorr[iLDR,:]) * 1.005 Volker@43: Volker@43: # Determine mean distance of all aXmean from each other for each iLDR Volker@43: meanDist = 0.0 Volker@43: for iaX in range(-naX, naX + 1): Volker@43: # mean Etamx value for certain error (iaX) of parameter aVar Volker@43: aXmean[iaX + naX] = np.mean(aEtamx[iLDR, aX == iaX]) Volker@43: # relative to absolute spread of Etamx Volker@43: meanDist = (np.max(aXmean) - np.min(aXmean)) / (Etamxmax - Etamxmin) * 100 Volker@43: Volker@43: plt.subplot(1, 5, iLDR + 1) Volker@43: (n, bins, patches) = plt.hist(aEtamx[iLDR, :], Volker@43: bins=50, log=False, Volker@43: range=[Rmin, Rmax], Volker@43: alpha=0.5, density=False, color='0.5', histtype='stepfilled') Volker@43: for iaX in range(-naX, naX + 1): Volker@43: plt.hist(aEtamx[iLDR, aX == iaX], Volker@43: range=[Rmin, Rmax], Volker@43: bins=50, log=False, alpha=0.3, density=False, histtype='stepfilled', Volker@43: label=str(round(X0 + iaX * daX / naX, 5))) Volker@43: if (iLDR == 2): Volker@43: leg = plt.legend() Volker@43: leg.get_frame().set_alpha(0.1) Volker@43: plt.tick_params(axis='both', labelsize=10) Volker@43: plt.plot([Etamx0, Etamx0], [0, np.max(n)], 'r-', lw=2) Volker@43: plt.gca().set_title("{0:3.0f}%".format(meanDist)) Volker@43: plt.gca().set_xlabel('Etamx0', color="red") Volker@43: fig.suptitle('Etamx - ' + LID + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar + ' error contribution', fontsize=14, y=1.10) Volker@43: return Volker@43: Volker@43: # calc contribution of the error of aVar = aX to aY for each LDRtrue Volker@43: def Contribution(aVar, aX, X0, daX, iaX, naX, aY, Ysum, widthSum): Volker@43: # aVar is the name of the parameter and aX is the subset of aY which is coloured in the plot Volker@43: # example: Contribution("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP, aLDRcorr, DOLPcontr) Volker@43: iLDR = -1 Volker@43: # Ysum, widthSum = np.zeros(5) Volker@43: meanDist = np.zeros(5) # iLDR Volker@43: widthDist = np.zeros(5) # iLDR Volker@43: for LDRTrue in LDRrange: Volker@43: aXmean = np.zeros(2 * naX + 1) Volker@43: aXwidth = np.zeros(2 * naX + 1) Volker@43: iLDR = iLDR + 1 Volker@43: # total width of distribution Volker@43: aYmin = np.amin(aY[iLDR, :]) Volker@43: aYmax = np.amax(aY[iLDR, :]) Volker@43: aYwidth = aYmax - aYmin Volker@43: # Determine mean distance of all aXmean from each other for each iLDR Volker@43: for iaX in range(-naX, naX + 1): Volker@43: # mean LDRCorr value for all errors iaX of parameter aVar Volker@43: aXmean[iaX + naX] = np.mean(aY[iLDR, aX == iaX]) Volker@43: aXwidth[iaX + naX] = np.max(aY[iLDR, aX == iaX]) - np.min(aY[iLDR, aX == iaX]) Volker@43: # relative to absolute spread of LDRCorrs Volker@43: meanDist[iLDR] = (np.max(aXmean) - np.min(aXmean)) / aYwidth * 1000 Volker@43: # meanDist[iLDR] = (aYwidth - aXwidth[naX]) / aYwidth * 1000 Volker@43: widthDist[iLDR] = (np.max(aXwidth) - aXwidth[naX]) / aYwidth * 1000 Volker@43: Volker@43: print("{:12}{:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f}"\ Volker@43: .format(aVar,meanDist[0],meanDist[1],meanDist[2],meanDist[3],meanDist[4],widthDist[0],widthDist[1],widthDist[2],widthDist[3],widthDist[4])) Volker@43: Ysum = Ysum + meanDist Volker@43: widthSum = widthSum + widthDist Volker@43: return(Ysum, widthSum) Volker@43: Volker@43: # print(.format(LDRrangeA[iLDR],)) Volker@43: Volker@43: # error contributions to a certain output aY; loop over all variables Volker@43: def Contribution_aY(aYvar, aY): Volker@43: Ysum = np.zeros(5) Volker@43: widthSum = np.zeros(5) Volker@43: # meanDist = np.zeros(5) # iLDR Volker@43: LDRrangeA = np.array(LDRrange) Volker@43: print() Volker@43: print(aYvar + ": contribution to the total error (per mill)") Volker@43: print(" of individual parameter errors of combined parameter errors") Volker@43: print(" at LDRtrue {:5.3f} {:5.3f} {:5.3f} {:5.3f} {:5.3f} {:5.3f} {:5.3f} {:5.3f} {:5.3f} {:5.3f}"\ Volker@43: .format(LDRrangeA[0],LDRrangeA[1],LDRrangeA[2],LDRrangeA[3],LDRrangeA[4],LDRrangeA[0],LDRrangeA[1],LDRrangeA[2],LDRrangeA[3],LDRrangeA[4])) Volker@43: print() Volker@43: if (nQin > 0): Ysum, widthSum = Contribution("Qin", aQin, Qin0, dQin, iQin, nQin, aY, Ysum, widthSum) Volker@43: if (nVin > 0): Ysum, widthSum = Contribution("Vin", aVin, Vin0, dVin, iVin, nVin, aY, Ysum, widthSum) Volker@43: if (nRotL > 0): Ysum, widthSum = Contribution("RotL", aRotL, RotL0, dRotL, iRotL, nRotL, aY, Ysum, widthSum) Volker@43: if (nRetE > 0): Ysum, widthSum = Contribution("RetE", aRetE, RetE0, dRetE, iRetE, nRetE, aY, Ysum, widthSum) Volker@43: if (nRotE > 0): Ysum, widthSum = Contribution("RotE", aRotE, RotE0, dRotE, iRotE, nRotE, aY, Ysum, widthSum) Volker@43: if (nDiE > 0): Ysum, widthSum = Contribution("DiE", aDiE, DiE0, dDiE, iDiE, nDiE, aY, Ysum, widthSum) Volker@43: if (nRetO > 0): Ysum, widthSum = Contribution("RetO", aRetO, RetO0, dRetO, iRetO, nRetO, aY, Ysum, widthSum) Volker@43: if (nRotO > 0): Ysum, widthSum = Contribution("RotO", aRotO, RotO0, dRotO, iRotO, nRotO, aY, Ysum, widthSum) Volker@43: if (nDiO > 0): Ysum, widthSum = Contribution("DiO", aDiO, DiO0, dDiO, iDiO, nDiO, aY, Ysum, widthSum) Volker@43: if (nDiC > 0): Ysum, widthSum = Contribution("DiC", aDiC, DiC0, dDiC, iDiC, nDiC, aY, Ysum, widthSum) Volker@43: if (nRotC > 0): Ysum, widthSum = Contribution("RotC", aRotC, RotC0, dRotC, iRotC, nRotC, aY, Ysum, widthSum) Volker@43: if (nRetC > 0): Ysum, widthSum = Contribution("RetC", aRetC, RetC0, dRetC, iRetC, nRetC, aY, Ysum, widthSum) Volker@43: if (nTP > 0): Ysum, widthSum = Contribution("TP", aTP, TP0, dTP, iTP, nTP, aY, Ysum, widthSum) Volker@43: if (nTS > 0): Ysum, widthSum = Contribution("TS", aTS, TS0, dTS, iTS, nTS, aY, Ysum, widthSum) Volker@43: if (nRP > 0): Ysum, widthSum = Contribution("RP", aRP, RP0, dRP, iRP, nRP, aY, Ysum, widthSum) Volker@43: if (nRS > 0): Ysum, widthSum = Contribution("RS", aRS, RS0, dRS, iRS, nRS, aY, Ysum, widthSum) Volker@43: if (nRetT > 0): Ysum, widthSum = Contribution("RetT", aRetT, RetT0, dRetT, iRetT, nRetT, aY, Ysum, widthSum) Volker@43: if (nRetR > 0): Ysum, widthSum = Contribution("RetR", aRetR, RetR0, dRetR, iRetR, nRetR, aY, Ysum, widthSum) Volker@43: if (nERaT > 0): Ysum, widthSum = Contribution("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT, aY, Ysum, widthSum) Volker@43: if (nERaR > 0): Ysum, widthSum = Contribution("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR, aY, Ysum, widthSum) Volker@43: if (nRotaT > 0): Ysum, widthSum = Contribution("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT, aY, Ysum, widthSum) Volker@43: if (nRotaR > 0): Ysum, widthSum = Contribution("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR, aY, Ysum, widthSum) Volker@43: if (nLDRCal > 0): Ysum, widthSum = Contribution("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal, aY, Ysum, widthSum) Volker@43: if (nTCalT > 0): Ysum, widthSum = Contribution("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT, aY, Ysum, widthSum) Volker@43: if (nTCalR > 0): Ysum, widthSum = Contribution("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR, aY, Ysum, widthSum) Volker@43: if (nNCal > 0): Ysum, widthSum = Contribution("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal, aY, Ysum, widthSum) Volker@43: if (nNCal > 0): Ysum, widthSum = Contribution("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal, aY, Ysum, widthSum) Volker@43: if (nNCal > 0): Ysum, widthSum = Contribution("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal, aY, Ysum, widthSum) Volker@43: if (nNCal > 0): Ysum, widthSum = Contribution("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal, aY, Ysum, widthSum) Volker@43: if (nNI > 0): Ysum, widthSum = Contribution("SigNoiseIt", aNIt, 0, 1, iNIt, nNI, aY, Ysum, widthSum) Volker@43: if (nNI > 0): Ysum, widthSum = Contribution("SigNoiseIr", aNIr, 0, 1, iNIr, nNI, aY, Ysum, widthSum) Volker@43: print("{:12}{:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f}"\ Volker@43: .format("Sum ",Ysum[0],Ysum[1],Ysum[2],Ysum[3],Ysum[4],widthSum[0],widthSum[1],widthSum[2],widthSum[3],widthSum[4])) Volker@43: Volker@43: Volker@43: # Plot LDR histograms Volker@43: if (nQin > 0): PlotSubHist("Qin", aQin, Qin0, dQin, iQin, nQin) Volker@43: if (nVin > 0): PlotSubHist("Vin", aVin, Vin0, dVin, iVin, nVin) Volker@43: if (nRotL > 0): PlotSubHist("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) Volker@43: if (nRetE > 0): PlotSubHist("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) Volker@43: if (nRotE > 0): PlotSubHist("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) Volker@43: if (nDiE > 0): PlotSubHist("DiE", aDiE, DiE0, dDiE, iDiE, nDiE) Volker@43: if (nRetO > 0): PlotSubHist("RetO", aRetO, RetO0, dRetO, iRetO, nRetO) Volker@43: if (nRotO > 0): PlotSubHist("RotO", aRotO, RotO0, dRotO, iRotO, nRotO) Volker@43: if (nDiO > 0): PlotSubHist("DiO", aDiO, DiO0, dDiO, iDiO, nDiO) Volker@43: if (nDiC > 0): PlotSubHist("DiC", aDiC, DiC0, dDiC, iDiC, nDiC) Volker@43: if (nRotC > 0): PlotSubHist("RotC", aRotC, RotC0, dRotC, iRotC, nRotC) Volker@43: if (nRetC > 0): PlotSubHist("RetC", aRetC, RetC0, dRetC, iRetC, nRetC) Volker@43: if (nTP > 0): PlotSubHist("TP", aTP, TP0, dTP, iTP, nTP) Volker@43: if (nTS > 0): PlotSubHist("TS", aTS, TS0, dTS, iTS, nTS) Volker@43: if (nRP > 0): PlotSubHist("RP", aRP, RP0, dRP, iRP, nRP) Volker@43: if (nRS > 0): PlotSubHist("RS", aRS, RS0, dRS, iRS, nRS) Volker@43: if (nRetT > 0): PlotSubHist("RetT", aRetT, RetT0, dRetT, iRetT, nRetT) Volker@43: if (nRetR > 0): PlotSubHist("RetR", aRetR, RetR0, dRetR, iRetR, nRetR) Volker@43: if (nERaT > 0): PlotSubHist("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) Volker@43: if (nERaR > 0): PlotSubHist("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) Volker@43: if (nRotaT > 0): PlotSubHist("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) Volker@43: if (nRotaR > 0): PlotSubHist("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) Volker@43: if (nLDRCal > 0): PlotSubHist("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) Volker@43: if (nTCalT > 0): PlotSubHist("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT) Volker@43: if (nTCalR > 0): PlotSubHist("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR) Volker@43: if (nNCal > 0): PlotSubHist("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal) Volker@43: if (nNCal > 0): PlotSubHist("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal) Volker@43: if (nNCal > 0): PlotSubHist("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal) Volker@43: if (nNCal > 0): PlotSubHist("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal) Volker@43: if (nNI > 0): PlotSubHist("SigNoiseIt", aNIt, 0, 1, iNIt, nNI) Volker@43: if (nNI > 0): PlotSubHist("SigNoiseIr", aNIr, 0, 1, iNIr, nNI) Volker@43: plt.show() Volker@43: plt.close Volker@43: Volker@43: Volker@43: Volker@43: # --- Plot LDRmin, LDRmax Volker@43: iLDR = -1 Volker@43: for LDRTrue in LDRrange: Volker@43: iLDR = iLDR + 1 Volker@43: LDRmin[iLDR] = np.amin(aLDRcorr[iLDR, :]) Volker@43: LDRmax[iLDR] = np.amax(aLDRcorr[iLDR, :]) Volker@43: LDRstd[iLDR] = np.std(aLDRcorr[iLDR, :]) Volker@43: LDRmean[iLDR] = np.mean(aLDRcorr[iLDR, :]) Volker@43: LDRmedian[iLDR] = np.median(aLDRcorr[iLDR, :]) Volker@43: LDRskew[iLDR] = skew(aLDRcorr[iLDR, :],bias=False) Volker@43: LDRkurt[iLDR] = kurtosis(aLDRcorr[iLDR, :],fisher=True,bias=False) Volker@43: Volker@43: fig2 = plt.figure() Volker@43: LDRrangeA = np.array(LDRrange) Volker@43: if((np.amax(LDRmax - LDRrangeA)-np.amin(LDRmin - LDRrangeA)) < 0.001): Volker@43: plt.ylim(-0.001,0.001) Volker@43: plt.plot(LDRrangeA, LDRmax - LDRrangeA, linewidth=2.0, color='b') Volker@43: plt.plot(LDRrangeA, LDRmin - LDRrangeA, linewidth=2.0, color='g') Volker@43: Volker@43: plt.xlabel('LDRtrue', fontsize=18) Volker@43: plt.ylabel('LDRTrue-LDRmin, LDRTrue-LDRmax', fontsize=14) Volker@43: plt.title(LID + ' ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]), fontsize=18) Volker@43: # plt.ylimit(-0.07, 0.07) Volker@43: plt.show() Volker@43: plt.close Volker@43: Volker@43: # --- Save LDRmin, LDRmax to file Volker@43: # http://stackoverflow.com/questions/4675728/redirect-stdout-to-a-file-in-python Volker@43: with open('output_files\\' + OutputFile, 'a') as f: Volker@43: # with open('output_files\\' + LID + '-' + InputFile[0:-3] + '-LDR_min_max.dat', 'w') as f: Volker@43: with redirect_stdout(f): Volker@43: print("Lidar ID: " + LID) Volker@43: print() Volker@43: print("minimum and maximum values of the distributions of possibly measured LDR for different LDRtrue") Volker@43: print("LDRtrue , LDRmin, LDRmax") Volker@43: for i in range(len(LDRrangeA)): Volker@43: print("{0:7.4f},{1:7.4f},{2:7.4f}".format(LDRrangeA[i], LDRmin[i], LDRmax[i])) Volker@43: print() Volker@43: # Print LDR statistics Volker@43: print("LDRtrue , mean , median, max-mean, min-mean, std, excess_kurtosis, skewness") Volker@43: iLDR = -1 Volker@43: LDRrangeA = np.array(LDRrange) Volker@43: for LDRTrue in LDRrange: Volker@43: iLDR = iLDR + 1 Volker@43: 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@43: .format(LDRrangeA[iLDR], LDRmean[iLDR], LDRmedian[iLDR], LDRmax[iLDR]-LDRrangeA[iLDR], \ Volker@43: LDRmin[iLDR]-LDRrangeA[iLDR], LDRstd[iLDR], LDRkurt[iLDR], LDRskew[iLDR])) Volker@43: print() Volker@43: # Calculate and print statistics for calibration factors Volker@43: print("minimum and maximum values of the distributions of signal ratios and calibration factors for different LDRtrue") Volker@43: iLDR = -1 Volker@43: LDRrangeA = np.array(LDRrange) Volker@43: print("LDRtrue , LDRsim, (max-min)/2, relerr") Volker@43: for LDRTrue in LDRrange: Volker@43: iLDR = iLDR + 1 Volker@43: LDRsimmin[iLDR] = np.amin(aLDRsim[iLDR, :]) Volker@43: LDRsimmax[iLDR] = np.amax(aLDRsim[iLDR, :]) Volker@43: # LDRsimstd = np.std(aLDRsim[iLDR, :]) Volker@43: LDRsimmean[iLDR] = np.mean(aLDRsim[iLDR, :]) Volker@43: # LDRsimmedian = np.median(aLDRsim[iLDR, :]) Volker@43: 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@43: iLDR = -1 Volker@43: print("LDRtrue , Etax , (max-min)/2, relerr") Volker@43: for LDRTrue in LDRrange: Volker@43: iLDR = iLDR + 1 Volker@43: Etaxmin = np.amin(aEtax[iLDR, :]) Volker@43: Etaxmax = np.amax(aEtax[iLDR, :]) Volker@43: # Etaxstd = np.std(aEtax[iLDR, :]) Volker@43: Etaxmean = np.mean(aEtax[iLDR, :]) Volker@43: # Etaxmedian = np.median(aEtax[iLDR, :]) Volker@43: 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@43: iLDR = -1 Volker@43: print("LDRtrue , Etapx , (max-min)/2, relerr") Volker@43: for LDRTrue in LDRrange: Volker@43: iLDR = iLDR + 1 Volker@43: Etapxmin = np.amin(aEtapx[iLDR, :]) Volker@43: Etapxmax = np.amax(aEtapx[iLDR, :]) Volker@43: # Etapxstd = np.std(aEtapx[iLDR, :]) Volker@43: Etapxmean = np.mean(aEtapx[iLDR, :]) Volker@43: # Etapxmedian = np.median(aEtapx[iLDR, :]) Volker@43: 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@43: iLDR = -1 Volker@43: print("LDRtrue , Etamx , (max-min)/2, relerr") Volker@43: for LDRTrue in LDRrange: Volker@43: iLDR = iLDR + 1 Volker@43: Etamxmin = np.amin(aEtamx[iLDR, :]) Volker@43: Etamxmax = np.amax(aEtamx[iLDR, :]) Volker@43: # Etamxstd = np.std(aEtamx[iLDR, :]) Volker@43: Etamxmean = np.mean(aEtamx[iLDR, :]) Volker@43: # Etamxmedian = np.median(aEtamx[iLDR, :]) Volker@43: 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@43: Volker@43: # Print LDR statistics Volker@43: print("LDRtrue , mean , median, max-mean, min-mean, std, excess_kurtosis, skewness") Volker@43: iLDR = -1 Volker@43: LDRrangeA = np.array(LDRrange) Volker@43: for LDRTrue in LDRrange: Volker@43: iLDR = iLDR + 1 Volker@43: 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@43: Volker@43: Volker@43: with open('output_files\\' + OutputFile, 'a') as f: Volker@43: # with open('output_files\\' + LID + '-' + InputFile[0:-3] + '-LDR_min_max.dat', 'a') as f: Volker@43: with redirect_stdout(f): Volker@43: Contribution_aY("LDRCorr", aLDRcorr) Volker@43: Contribution_aY("LDRsim", aLDRsim) Volker@43: Contribution_aY("EtaX, D90", aEtax) Volker@43: Contribution_aY("Etapx, +45°", aEtapx) Volker@43: Contribution_aY("Etamx -45°", aEtamx) Volker@43: Volker@43: Volker@43: # Plot other histograms Volker@43: if (bPlotEtax): Volker@43: Volker@43: if (nQin > 0): PlotLDRsim("Qin", aQin, Qin0, dQin, iQin, nQin) Volker@43: if (nVin > 0): PlotLDRsim("Vin", aVin, Vin0, dVin, iVin, nVin) Volker@43: if (nRotL > 0): PlotLDRsim("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) Volker@43: if (nRetE > 0): PlotLDRsim("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) Volker@43: if (nRotE > 0): PlotLDRsim("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) Volker@43: if (nDiE > 0): PlotLDRsim("DiE", aDiE, DiE0, dDiE, iDiE, nDiE) Volker@43: if (nRetO > 0): PlotLDRsim("RetO", aRetO, RetO0, dRetO, iRetO, nRetO) Volker@43: if (nRotO > 0): PlotLDRsim("RotO", aRotO, RotO0, dRotO, iRotO, nRotO) Volker@43: if (nDiO > 0): PlotLDRsim("DiO", aDiO, DiO0, dDiO, iDiO, nDiO) Volker@43: if (nDiC > 0): PlotLDRsim("DiC", aDiC, DiC0, dDiC, iDiC, nDiC) Volker@43: if (nRotC > 0): PlotLDRsim("RotC", aRotC, RotC0, dRotC, iRotC, nRotC) Volker@43: if (nRetC > 0): PlotLDRsim("RetC", aRetC, RetC0, dRetC, iRetC, nRetC) Volker@43: if (nTP > 0): PlotLDRsim("TP", aTP, TP0, dTP, iTP, nTP) Volker@43: if (nTS > 0): PlotLDRsim("TS", aTS, TS0, dTS, iTS, nTS) Volker@43: if (nRP > 0): PlotLDRsim("RP", aRP, RP0, dRP, iRP, nRP) Volker@43: if (nRS > 0): PlotLDRsim("RS", aRS, RS0, dRS, iRS, nRS) Volker@43: if (nRetT > 0): PlotLDRsim("RetT", aRetT, RetT0, dRetT, iRetT, nRetT) Volker@43: if (nRetR > 0): PlotLDRsim("RetR", aRetR, RetR0, dRetR, iRetR, nRetR) Volker@43: if (nERaT > 0): PlotLDRsim("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) Volker@43: if (nERaR > 0): PlotLDRsim("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) Volker@43: if (nRotaT > 0): PlotLDRsim("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) Volker@43: if (nRotaR > 0): PlotLDRsim("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) Volker@43: if (nLDRCal > 0): PlotLDRsim("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) Volker@43: if (nTCalT > 0): PlotLDRsim("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT) Volker@43: if (nTCalR > 0): PlotLDRsim("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR) Volker@43: if (nNCal > 0): PlotLDRsim("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal) Volker@43: if (nNCal > 0): PlotLDRsim("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal) Volker@43: if (nNCal > 0): PlotLDRsim("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal) Volker@43: if (nNCal > 0): PlotLDRsim("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal) Volker@43: if (nNI > 0): PlotLDRsim("SigNoiseIt", aNIt, 0, 1, iNIt, nNI) Volker@43: if (nNI > 0): PlotLDRsim("SigNoiseIr", aNIr, 0, 1, iNIr, nNI) Volker@43: plt.show() Volker@43: plt.close Volker@43: print("---------------------------------------...producing more plots...------------------------------------------------------------------") Volker@43: Volker@43: if (nQin > 0): PlotEtax("Qin", aQin, Qin0, dQin, iQin, nQin) Volker@43: if (nVin > 0): PlotEtax("Vin", aVin, Vin0, dVin, iVin, nVin) Volker@43: if (nRotL > 0): PlotEtax("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) Volker@43: if (nRetE > 0): PlotEtax("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) Volker@43: if (nRotE > 0): PlotEtax("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) Volker@43: if (nDiE > 0): PlotEtax("DiE", aDiE, DiE0, dDiE, iDiE, nDiE) Volker@43: if (nRetO > 0): PlotEtax("RetO", aRetO, RetO0, dRetO, iRetO, nRetO) Volker@43: if (nRotO > 0): PlotEtax("RotO", aRotO, RotO0, dRotO, iRotO, nRotO) Volker@43: if (nDiO > 0): PlotEtax("DiO", aDiO, DiO0, dDiO, iDiO, nDiO) Volker@43: if (nDiC > 0): PlotEtax("DiC", aDiC, DiC0, dDiC, iDiC, nDiC) Volker@43: if (nRotC > 0): PlotEtax("RotC", aRotC, RotC0, dRotC, iRotC, nRotC) Volker@43: if (nRetC > 0): PlotEtax("RetC", aRetC, RetC0, dRetC, iRetC, nRetC) Volker@43: if (nTP > 0): PlotEtax("TP", aTP, TP0, dTP, iTP, nTP) Volker@43: if (nTS > 0): PlotEtax("TS", aTS, TS0, dTS, iTS, nTS) Volker@43: if (nRP > 0): PlotEtax("RP", aRP, RP0, dRP, iRP, nRP) Volker@43: if (nRS > 0): PlotEtax("RS", aRS, RS0, dRS, iRS, nRS) Volker@43: if (nRetT > 0): PlotEtax("RetT", aRetT, RetT0, dRetT, iRetT, nRetT) Volker@43: if (nRetR > 0): PlotEtax("RetR", aRetR, RetR0, dRetR, iRetR, nRetR) Volker@43: if (nERaT > 0): PlotEtax("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) Volker@43: if (nERaR > 0): PlotEtax("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) Volker@43: if (nRotaT > 0): PlotEtax("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) Volker@43: if (nRotaR > 0): PlotEtax("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) Volker@43: if (nLDRCal > 0): PlotEtax("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) Volker@43: if (nTCalT > 0): PlotEtax("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT) Volker@43: if (nTCalR > 0): PlotEtax("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR) Volker@43: if (nNCal > 0): PlotEtax("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal) Volker@43: if (nNCal > 0): PlotEtax("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal) Volker@43: if (nNCal > 0): PlotEtax("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal) Volker@43: if (nNCal > 0): PlotEtax("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal) Volker@43: if (nNI > 0): PlotEtax("SigNoiseIt", aNIt, 0, 1, iNIt, nNI) Volker@43: if (nNI > 0): PlotEtax("SigNoiseIr", aNIr, 0, 1, iNIr, nNI) Volker@43: plt.show() Volker@43: plt.close Volker@43: print("---------------------------------------...producing more plots...------------------------------------------------------------------") Volker@43: Volker@43: if (nQin > 0): PlotEtapx("Qin", aQin, Qin0, dQin, iQin, nQin) Volker@43: if (nVin > 0): PlotEtapx("Vin", aVin, Vin0, dVin, iVin, nVin) Volker@43: if (nRotL > 0): PlotEtapx("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) Volker@43: if (nRetE > 0): PlotEtapx("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) Volker@43: if (nRotE > 0): PlotEtapx("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) Volker@43: if (nDiE > 0): PlotEtapx("DiE", aDiE, DiE0, dDiE, iDiE, nDiE) Volker@43: if (nRetO > 0): PlotEtapx("RetO", aRetO, RetO0, dRetO, iRetO, nRetO) Volker@43: if (nRotO > 0): PlotEtapx("RotO", aRotO, RotO0, dRotO, iRotO, nRotO) Volker@43: if (nDiO > 0): PlotEtapx("DiO", aDiO, DiO0, dDiO, iDiO, nDiO) Volker@43: if (nDiC > 0): PlotEtapx("DiC", aDiC, DiC0, dDiC, iDiC, nDiC) Volker@43: if (nRotC > 0): PlotEtapx("RotC", aRotC, RotC0, dRotC, iRotC, nRotC) Volker@43: if (nRetC > 0): PlotEtapx("RetC", aRetC, RetC0, dRetC, iRetC, nRetC) Volker@43: if (nTP > 0): PlotEtapx("TP", aTP, TP0, dTP, iTP, nTP) Volker@43: if (nTS > 0): PlotEtapx("TS", aTS, TS0, dTS, iTS, nTS) Volker@43: if (nRP > 0): PlotEtapx("RP", aRP, RP0, dRP, iRP, nRP) Volker@43: if (nRS > 0): PlotEtapx("RS", aRS, RS0, dRS, iRS, nRS) Volker@43: if (nRetT > 0): PlotEtapx("RetT", aRetT, RetT0, dRetT, iRetT, nRetT) Volker@43: if (nRetR > 0): PlotEtapx("RetR", aRetR, RetR0, dRetR, iRetR, nRetR) Volker@43: if (nERaT > 0): PlotEtapx("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) Volker@43: if (nERaR > 0): PlotEtapx("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) Volker@43: if (nRotaT > 0): PlotEtapx("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) Volker@43: if (nRotaR > 0): PlotEtapx("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) Volker@43: if (nLDRCal > 0): PlotEtapx("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) Volker@43: if (nTCalT > 0): PlotEtapx("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT) Volker@43: if (nTCalR > 0): PlotEtapx("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR) Volker@43: if (nNCal > 0): PlotEtapx("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal) Volker@43: if (nNCal > 0): PlotEtapx("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal) Volker@43: if (nNCal > 0): PlotEtapx("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal) Volker@43: if (nNCal > 0): PlotEtapx("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal) Volker@43: if (nNI > 0): PlotEtapx("SigNoiseIt", aNIt, 0, 1, iNIt, nNI) Volker@43: if (nNI > 0): PlotEtapx("SigNoiseIr", aNIr, 0, 1, iNIr, nNI) Volker@43: plt.show() Volker@43: plt.close Volker@43: print("---------------------------------------...producing more plots...------------------------------------------------------------------") Volker@43: Volker@43: if (nQin > 0): PlotEtamx("Qin", aQin, Qin0, dQin, iQin, nQin) Volker@43: if (nVin > 0): PlotEtamx("Vin", aVin, Vin0, dVin, iVin, nVin) Volker@43: if (nRotL > 0): PlotEtamx("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) Volker@43: if (nRetE > 0): PlotEtamx("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) Volker@43: if (nRotE > 0): PlotEtamx("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) Volker@43: if (nDiE > 0): PlotEtamx("DiE", aDiE, DiE0, dDiE, iDiE, nDiE) Volker@43: if (nRetO > 0): PlotEtamx("RetO", aRetO, RetO0, dRetO, iRetO, nRetO) Volker@43: if (nRotO > 0): PlotEtamx("RotO", aRotO, RotO0, dRotO, iRotO, nRotO) Volker@43: if (nDiO > 0): PlotEtamx("DiO", aDiO, DiO0, dDiO, iDiO, nDiO) Volker@43: if (nDiC > 0): PlotEtamx("DiC", aDiC, DiC0, dDiC, iDiC, nDiC) Volker@43: if (nRotC > 0): PlotEtamx("RotC", aRotC, RotC0, dRotC, iRotC, nRotC) Volker@43: if (nRetC > 0): PlotEtamx("RetC", aRetC, RetC0, dRetC, iRetC, nRetC) Volker@43: if (nTP > 0): PlotEtamx("TP", aTP, TP0, dTP, iTP, nTP) Volker@43: if (nTS > 0): PlotEtamx("TS", aTS, TS0, dTS, iTS, nTS) Volker@43: if (nRP > 0): PlotEtamx("RP", aRP, RP0, dRP, iRP, nRP) Volker@43: if (nRS > 0): PlotEtamx("RS", aRS, RS0, dRS, iRS, nRS) Volker@43: if (nRetT > 0): PlotEtamx("RetT", aRetT, RetT0, dRetT, iRetT, nRetT) Volker@43: if (nRetR > 0): PlotEtamx("RetR", aRetR, RetR0, dRetR, iRetR, nRetR) Volker@43: if (nERaT > 0): PlotEtamx("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) Volker@43: if (nERaR > 0): PlotEtamx("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) Volker@43: if (nRotaT > 0): PlotEtamx("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) Volker@43: if (nRotaR > 0): PlotEtamx("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) Volker@43: if (nLDRCal > 0): PlotEtamx("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) Volker@43: if (nTCalT > 0): PlotEtamx("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT) Volker@43: if (nTCalR > 0): PlotEtamx("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR) Volker@43: if (nNCal > 0): PlotEtamx("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal) Volker@43: if (nNCal > 0): PlotEtamx("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal) Volker@43: if (nNCal > 0): PlotEtamx("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal) Volker@43: if (nNCal > 0): PlotEtamx("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal) Volker@43: if (nNI > 0): PlotEtamx("SigNoiseIt", aNIt, 0, 1, iNIt, nNI) Volker@43: if (nNI > 0): PlotEtamx("SigNoiseIr", aNIr, 0, 1, iNIr, nNI) Volker@43: plt.show() Volker@43: plt.close Volker@43: Volker@43: # Print Etax statistics Volker@43: Etaxmin = np.amin(aEtax[1, :]) Volker@43: Etaxmax = np.amax(aEtax[1, :]) Volker@43: Etaxstd = np.std(aEtax[1, :]) Volker@43: Etaxmean = np.mean(aEtax[1, :]) Volker@43: Etaxmedian = np.median(aEtax[1, :]) Volker@43: print("Etax , max-mean, min-mean, median, mean ± std, eta") Volker@43: 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@43: print() Volker@43: Volker@43: # Calculate and print statistics for calibration factors Volker@43: iLDR = -1 Volker@43: LDRrangeA = np.array(LDRrange) Volker@43: print("LDR...., LDRsim, (max-min)/2, relerr") Volker@43: for LDRTrue in LDRrange: Volker@43: iLDR = iLDR + 1 Volker@43: LDRsimmin[iLDR] = np.amin(aLDRsim[iLDR, :]) Volker@43: LDRsimmax[iLDR] = np.amax(aLDRsim[iLDR, :]) Volker@43: # LDRsimstd = np.std(aLDRsim[iLDR, :]) Volker@43: LDRsimmean[iLDR] = np.mean(aLDRsim[iLDR, :]) Volker@43: # LDRsimmedian = np.median(aLDRsim[iLDR, :]) Volker@43: 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@43: iLDR = -1 Volker@43: print("LDR...., Etax , (max-min)/2, relerr") Volker@43: for LDRTrue in LDRrange: Volker@43: iLDR = iLDR + 1 Volker@43: Etaxmin = np.amin(aEtax[iLDR, :]) Volker@43: Etaxmax = np.amax(aEtax[iLDR, :]) Volker@43: # Etaxstd = np.std(aEtax[iLDR, :]) Volker@43: Etaxmean = np.mean(aEtax[iLDR, :]) Volker@43: # Etaxmedian = np.median(aEtax[iLDR, :]) Volker@43: 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@43: iLDR = -1 Volker@43: print("LDR...., Etapx , (max-min)/2, relerr") Volker@43: for LDRTrue in LDRrange: Volker@43: iLDR = iLDR + 1 Volker@43: Etapxmin = np.amin(aEtapx[iLDR, :]) Volker@43: Etapxmax = np.amax(aEtapx[iLDR, :]) Volker@43: # Etapxstd = np.std(aEtapx[iLDR, :]) Volker@43: Etapxmean = np.mean(aEtapx[iLDR, :]) Volker@43: # Etapxmedian = np.median(aEtapx[iLDR, :]) Volker@43: 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@43: iLDR = -1 Volker@43: print("LDR...., Etamx , (max-min)/2, relerr") Volker@43: for LDRTrue in LDRrange: Volker@43: iLDR = iLDR + 1 Volker@43: Etamxmin = np.amin(aEtamx[iLDR, :]) Volker@43: Etamxmax = np.amax(aEtamx[iLDR, :]) Volker@43: # Etamxstd = np.std(aEtamx[iLDR, :]) Volker@43: Etamxmean = np.mean(aEtamx[iLDR, :]) Volker@43: # Etamxmedian = np.median(aEtamx[iLDR, :]) Volker@43: 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@43: Volker@43: f.close() Volker@43: Volker@43: Volker@43: ''' Volker@43: # --- Plot F11 histograms Volker@43: print() Volker@43: print(" ############################################################################## ") Volker@43: print(Text1) Volker@43: print() Volker@43: Volker@43: iLDR = 5 Volker@43: for LDRTrue in LDRrange: Volker@43: iLDR = iLDR - 1 Volker@43: #aF11corr[iLDR,:] = aF11corr[iLDR,:] / aF11corr[0,:] - 1.0 Volker@43: aF11corr[iLDR,:] = aF11corr[iLDR,:] / aF11sim0[iLDR] - 1.0 Volker@43: # Plot F11 Volker@43: def PlotSubHistF11(aVar, aX, X0, daX, iaX, naX): Volker@43: fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) Volker@43: iLDR = -1 Volker@43: for LDRTrue in LDRrange: Volker@43: iLDR = iLDR + 1 Volker@43: Volker@43: #F11min[iLDR] = np.min(aF11corr[iLDR,:]) Volker@43: #F11max[iLDR] = np.max(aF11corr[iLDR,:]) Volker@43: #Rmin = F11min[iLDR] * 0.995 # np.min(aLDRcorr[iLDR,:]) * 0.995 Volker@43: #Rmax = F11max[iLDR] * 1.005 # np.max(aLDRcorr[iLDR,:]) * 1.005 Volker@43: Volker@43: #Rmin = 0.8 Volker@43: #Rmax = 1.2 Volker@43: Volker@43: #plt.subplot(5,2,iLDR+1) Volker@43: plt.subplot(1,5,iLDR+1) Volker@43: (n, bins, patches) = plt.hist(aF11corr[iLDR,:], Volker@43: bins=100, log=False, Volker@43: alpha=0.5, density=False, color = '0.5', histtype='stepfilled') Volker@43: Volker@43: for iaX in range(-naX,naX+1): Volker@43: plt.hist(aF11corr[iLDR,aX == iaX], Volker@43: bins=100, log=False, alpha=0.3, density=False, histtype='stepfilled', label = str(round(X0 + iaX*daX/naX,5))) Volker@43: Volker@43: if (iLDR == 2): plt.legend() Volker@43: Volker@43: plt.tick_params(axis='both', labelsize=9) Volker@43: #plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2) Volker@43: Volker@43: #plt.title(LID + ' ' + aVar, fontsize=18) Volker@43: #plt.ylabel('frequency', fontsize=10) Volker@43: #plt.xlabel('LDRCorr', fontsize=10) Volker@43: #fig.tight_layout() Volker@43: fig.suptitle(LID + ' ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar, fontsize=14, y=1.05) Volker@43: #plt.show() Volker@43: #fig.savefig(LID + '_' + aVar + '.png', dpi=150, bbox_inches='tight', pad_inches=0) Volker@43: #plt.close Volker@43: return Volker@43: Volker@43: if (nQin > 0): PlotSubHistF11("Qin", aQin, Qin0, dQin, iQin, nQin) Volker@43: if (nVin > 0): PlotSubHistF11("Vin", aVin, Vin0, dVin, iVin, nVin) Volker@43: if (nRotL > 0): PlotSubHistF11("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) Volker@43: if (nRetE > 0): PlotSubHistF11("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) Volker@43: if (nRotE > 0): PlotSubHistF11("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) Volker@43: if (nDiE > 0): PlotSubHistF11("DiE", aDiE, DiE0, dDiE, iDiE, nDiE) Volker@43: if (nRetO > 0): PlotSubHistF11("RetO", aRetO, RetO0, dRetO, iRetO, nRetO) Volker@43: if (nRotO > 0): PlotSubHistF11("RotO", aRotO, RotO0, dRotO, iRotO, nRotO) Volker@43: if (nDiO > 0): PlotSubHistF11("DiO", aDiO, DiO0, dDiO, iDiO, nDiO) Volker@43: if (nDiC > 0): PlotSubHistF11("DiC", aDiC, DiC0, dDiC, iDiC, nDiC) Volker@43: if (nRotC > 0): PlotSubHistF11("RotC", aRotC, RotC0, dRotC, iRotC, nRotC) Volker@43: if (nRetC > 0): PlotSubHistF11("RetC", aRetC, RetC0, dRetC, iRetC, nRetC) Volker@43: if (nTP > 0): PlotSubHistF11("TP", aTP, TP0, dTP, iTP, nTP) Volker@43: if (nTS > 0): PlotSubHistF11("TS", aTS, TS0, dTS, iTS, nTS) Volker@43: if (nRP > 0): PlotSubHistF11("RP", aRP, RP0, dRP, iRP, nRP) Volker@43: if (nRS > 0): PlotSubHistF11("RS", aRS, RS0, dRS, iRS, nRS) Volker@43: if (nRetT > 0): PlotSubHistF11("RetT", aRetT, RetT0, dRetT, iRetT, nRetT) Volker@43: if (nRetR > 0): PlotSubHistF11("RetR", aRetR, RetR0, dRetR, iRetR, nRetR) Volker@43: if (nERaT > 0): PlotSubHistF11("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) Volker@43: if (nERaR > 0): PlotSubHistF11("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) Volker@43: if (nRotaT > 0): PlotSubHistF11("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) Volker@43: if (nRotaR > 0): PlotSubHistF11("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) Volker@43: if (nLDRCal > 0): PlotSubHistF11("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) Volker@43: if (nTCalT > 0): PlotSubHistF11("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT) Volker@43: if (nTCalR > 0): PlotSubHistF11("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR) Volker@43: if (nNCal > 0): PlotSubHistF11("CalNoise", aNCal, 0, 1/nNCal, iNCal, nNCal) Volker@43: if (nNI > 0): PlotSubHistF11("SigNoise", aNI, 0, 1/nNI, iNI, nNI) Volker@43: Volker@43: Volker@43: plt.show() Volker@43: plt.close Volker@43: Volker@43: ''' Volker@43: ''' Volker@43: # only histogram Volker@43: #print("******************* " + aVar + " *******************") Volker@43: fig, ax = plt.subplots(nrows=5, ncols=2, sharex=True, sharey=True, figsize=(10, 10)) Volker@43: iLDR = -1 Volker@43: for LDRTrue in LDRrange: Volker@43: iLDR = iLDR + 1 Volker@43: LDRmin[iLDR] = np.min(aLDRcorr[iLDR,:]) Volker@43: LDRmax[iLDR] = np.max(aLDRcorr[iLDR,:]) Volker@43: Rmin = np.min(aLDRcorr[iLDR,:]) * 0.999 Volker@43: Rmax = np.max(aLDRcorr[iLDR,:]) * 1.001 Volker@43: plt.subplot(5,2,iLDR+1) Volker@43: (n, bins, patches) = plt.hist(aLDRcorr[iLDR,:], Volker@43: range=[Rmin, Rmax], Volker@43: bins=200, log=False, alpha=0.2, density=False, color = '0.5', histtype='stepfilled') Volker@43: plt.tick_params(axis='both', labelsize=9) Volker@43: plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2) Volker@43: plt.show() Volker@43: plt.close Volker@43: # --- End of Plot F11 histograms Volker@43: ''' Volker@43: Volker@43: Volker@43: ''' Volker@43: # --- Plot K over LDRCal Volker@43: fig3 = plt.figure() Volker@43: plt.plot(LDRCal0+aLDRCal*dLDRCal/nLDRCal,aGHK[4,:], linewidth=2.0, color='b') Volker@43: Volker@43: plt.xlabel('LDRCal', fontsize=18) Volker@43: plt.ylabel('K', fontsize=14) Volker@43: plt.title(LID, fontsize=18) Volker@43: plt.show() Volker@43: plt.close Volker@43: ''' Volker@43: Volker@43: # Additional plot routines ======> Volker@43: ''' Volker@43: #****************************************************************************** Volker@43: # 1. Plot LDRCorrected - LDR(measured Icross/Iparallel) Volker@43: LDRa = np.arange(1.,100.)*0.005 Volker@43: LDRCorra = np.arange(1.,100.) Volker@43: if Y == - 1.: LDRa = 1./LDRa Volker@43: LDRCorra = (1./Eta*LDRa*(GT+HT)-(GR+HR))/((GR-HR)-1./Eta*LDRa*(GT-HT)) Volker@43: if Y == - 1.: LDRa = 1./LDRa Volker@43: # Volker@43: #fig = plt.figure() Volker@43: plt.plot(LDRa,LDRCorra-LDRa) Volker@43: plt.plot([0.,0.5],[0.,0.5]) Volker@43: plt.suptitle('LDRCorrected - LDR(measured Icross/Iparallel)', fontsize=16) Volker@43: plt.xlabel('LDR', fontsize=18) Volker@43: plt.ylabel('LDRCorr - LDR', fontsize=16) Volker@43: #plt.savefig('test.png') Volker@43: # Volker@43: ''' Volker@43: ''' Volker@43: #****************************************************************************** Volker@43: # 2. Plot LDRsim (simulated measurements without corrections = Icross/Iparallel) over LDRtrue Volker@43: LDRa = np.arange(1.,100.)*0.005 Volker@43: LDRsima = np.arange(1.,100.) Volker@43: Volker@43: atruea = (1.-LDRa)/(1+LDRa) Volker@43: Ita = TiT*TiO*IinL*(GT+atruea*HT) Volker@43: Ira = TiR*TiO*IinL*(GR+atruea*HR) Volker@43: LDRsima = Ira/Ita # simulated uncorrected LDR with Y from input file Volker@43: if Y == -1.: LDRsima = 1./LDRsima Volker@43: # Volker@43: #fig = plt.figure() Volker@43: plt.plot(LDRa,LDRsima) Volker@43: plt.plot([0.,0.5],[0.,0.5]) Volker@43: plt.suptitle('LDRsim (simulated measurements without corrections = Icross/Iparallel) over LDRtrue', fontsize=10) Volker@43: plt.xlabel('LDRtrue', fontsize=18) Volker@43: plt.ylabel('LDRsim', fontsize=16) Volker@43: #plt.savefig('test.png') Volker@43: # Volker@43: '''