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