Lidar polarisation GHK-parameter and error calculation: lidar_correction

lidar_correction_ghk.py@79fa4a41420f (annotated)

lidar_correction_ghk.py

Thu, 24 Jan 2019 21:24:25 +0100

author: Volker Freudenthaler <volker.freudenthaler@lmu.de>
date: Thu, 24 Jan 2019 21:24:25 +0100
changeset 28: 79fa4a41420f
parent 26: 28b5510492ba
permissions: -rw-r--r--

Bug fixes and further options. Read Improvements_of_lidar_correction_ghk_ver.0.9.8_190124.pdf

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

Mercurial > public > atmospheric_lidar_ghk / annotate

lidar_correction_ghk.py@79fa4a41420f (annotated)

lidar_correction_ghk.py