Mercurial > public > atmospheric_lidar_ghk / file revision

# -*- 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.7:  includes the random error (signal noise) of the calibration and standard measurements
Changes:
    Line 1687   Eta = (TaR * TiR) / (TaT * TiT)
    Line 1691   K = Etax / Eta  # K of the real system; but correction in Line 1721 with K0 / Etax
    should work with nTCalT = nTCalR = 0
Ver. 0.9.7b:
    ToDo: include error due to TCalT und TCalR => determination of NCalT and NCalR etc. in error calculation line 1741ff
    combined error loops iNI and INCal for signals
Ver. 0.9.7c: individual error loops for each of the six signals
Ver. 0.9.7c2: different calculation of the signal noise errors
Ver. 0.9.7c3: n.a.different calculation of the signal noise errors
Ver. 0.9.7c4: test to speed up the loops for error calculation by moving them just before the actual calculation: still some code errors
Ver. 0.9.8:
    - 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
    - calculation of the PLDR from LDR and BSR, BSR, and LDRm
    - 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.
Ver. 0.9.8b:  change from  "TTa = TiT * TaT"  to  "TTa = TiT * TaT * ATPT" etc. (compare ver 0.9.8 with 0.9.8b) removes
	- the strong Tp dependence of the errors
	- the factor 2 in the GH parameters
    - see c:\technik\Optik\Polarizers\DepCal\ApplOpt\GH-parameters-190114.odt
Ver. 0.9.8c:  includes error of Etax
Ver. 0.9.8d:  Eta0, K0 etc in error loop replaced by Eta0y, K0y etc. Changes in signal noise calculations
Ver. 0.9.8e:  ambiguous laser spec. DOLP (no discrimination between left and right circular polarisation) replaced by Stokes parameters Qin, Uin
Ver. 0.9.8e2:  Added plot of LDRsim, Etax, Etapx, Etamx;  LDRCorr and aLDRcorr consistently named
Ver. 0.9.8e3:  Change of OutputFile name; Change of Ir and It noise if (CalcFrom0deg) = False;  (Different calculation of error contributions tested but not implemented)
Ver. 0.9.8e4:  text changed for y=+-1 (see line 274 ff and line 1044 ff
Ver. 0.9.8e5:  changed: LDRunCorr = LDRsim / Etax
Ver. 0.9.8e6:  K(0.05) instead K(0.02)
Ver. 0.9.8f:  Tip from Ioannis Binietoglou for LINUX compatibility: Using os.path.join should work in all operating systems.
                # After line 1007:
                output_path = os.path.join('output_files', OutputFile)
                with open(output_path, 'w') as f:
                # Line 1147
                file = open(output_path, 'r')
------------

 ========================================================
simulation: LDRsim = Ir / It with variable parameters (possible truths)
    G,H,Eta,Etax,K
    It = TaT * TiT * ATP1 * TiO * TiE * (GT + atrue * HT)
    LDRsim = Ir / It
consistency test: is forward simulation and correction consistent?
    LDRCorr = (LDRsim / Eta * (GT + HT) - (GR + HR)) / ((GR - HR) - LDRsim / Eta * (GT - HT)) => atrue?
assumed true: G0,H0,Eta0,Etax0,K0 => actual retrievals of LDRCorr
    => correct possible truths with assumed true G0,H0,Eta0
    measure: It, Ir, EtaX
    LDRunCorr = LDRsim / Etax
    correct it with G0,H0,K0:
    LDRCorr = (LDRsim / (Etax / K0) * (GT0 + HT0) - (GR0 + HR0)) / ((GR0 - HR0) - LDRsim0 / (Etax / K0) * (GT0 - HT0))
"""
# 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

from scipy.stats import kurtosis
from scipy.stats import skew
# use: kurtosis(data, fisher=True,bias=False) => 0; skew(data,bias=False) => 0
# Comment: the seaborn library makes nicer plots, but the code works also without it.
import numpy as np
import matplotlib.pyplot as plt

try:
    import seaborn as sns

    sns_loaded = True
except ImportError:
    sns_loaded = False

# 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.
Error_Calc = True
LID = "internal"
EID = "internal"
# --- IL Laser IL and +-Uncertainty
Qin, dQin, nQin = 1., 0.0,  0	# second Stokes vector parameter; default 1 => linear polarization
Vin, dVin, nVin = 0., 0.0,  0	# fourth Stokes vector parameter
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.0, 0
    RS, dRS, nRS = 1 - TS, 0.0, 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: PBS incidence plane is parallel to the reference plane and polarisation in reference plane is finally transmitted.
#    Y = -1: PBS incidence plane is perpendicular to the reference plane and 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 = 1.0    # 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
NI = 50000 #number of photon counts in the parallel 0°-signal
eFacT = 1.0                     			# rel. amplification of transmitted channel, approximate values are sufficient; def. = 1
eFacR = 10.0
IoutTp0, IoutTp, dIoutTp0 = 0.5, 0.5, 0.0
IoutTm0, IoutTm, dIoutTm0 = 0.5, 0.5, 0.0
IoutRp0, IoutRp, dIoutRp0 = 0.5, 0.5, 0.0
IoutRm0, IoutRm, dIoutRm0 = 0.5, 0.5, 0.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.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.004
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 the laser', 'behind the emitter', 'before the receiver', 'before the PBS']
Type = ['', 'Mechanical rotator', 'HWP rotator', 'Linear polarizer', 'QWP rotator', 'Circular polarizer',
        'real HWP +-22.5°']

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_0.9.8e4-PollyXT_Lacros.py'
InputFile = 'optic_input_example_lidar_ver0.9.8e.py'
InputFile = 'calibrator-test-1-ver0.9.8e.py'
InputFile = 'optic_input_0.9.8e4-BRC-532.py'
InputFile = ''
InputFile = ''
InputFile = ''
InputFile = 'RALI-may2020-X.py'
InputFile = 'mulhacen_run_532xp-bias.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
Qin0, dQin, nQin = Qin, dQin, nQin
Vin0, dVin, nVin = Vin, dVin, nVin
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 = float(1. + C2aT * DaT * DiT)
ARPT = float(1. + C2aR * DaR * DiR)
TTa = TiT * TaT * ATPT  # unpolarized transmission
TRa = TiR * TaR * ARPT  # unpolarized transmission
Eta0 = TRa / TTa

# --- alternative texts for output
dY = ['perpendicular', '', 'parallel']
dY2 = ['reflected', '', 'transmitted']
if ((abs(RotL) < 45 and Y == 1) or (abs(RotL) >= 45 and Y == -1)):
    dY3 = "the parallel laser polarisation is detected in the transmitted channel."
else:
    dY3 = "the parallel laser polarisation is detected in the reflected channel."

# --- check input errors
if ((Qin ** 2 + Vin ** 2) ** 0.5) > 1:
    print("Error: degree of polarisation of laser > 1. Check Qin and Vin! ")
    sys.exit()

# --- 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, Qin, Vin, 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)
    atrue = (1. - LDRtrue) / (1. + LDRtrue)

    # 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 = Qin
    UinL = 0.
    VinL = Vin
    # 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):	# this works!
        dIt = ((It * NI * eFacT) ** -0.5)
        dIr = ((Ir * NI * eFacR) ** -0.5)
        '''
        dIt = ((NCalT * It / IoutTp * NILfac / TCalT) ** -0.5)
        dIr = ((NCalR * Ir / IoutRp * NILfac / TCalR) ** -0.5)
        '''
    else:	# does this work? Why not as above?
        dIt = ((NCalT * 2 * NILfac / TCalT ) ** -0.5)
        dIr = ((NCalR * 2 * NILfac / TCalR) ** -0.5)

        # ----- 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
    # 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
LDRtrue = LDRtrue2
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, Qin0, Vin0, RotL0, RotE0, RetE0, DiE0,
     RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0,
     ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0)
Etax0 = K0 * Eta0
Etapx0 = IoutRp0 / IoutTp0
Etamx0 = IoutRm0 / IoutTm0
# --- Print parameters to console and output file
OutputFile = 'output_' + InputFile[0:-3] + '_' + fname[0:-3] +'.dat'
output_path = os.path.join('output_files', OutputFile)
# with open('output_files\\' + OutputFile, 'w') as f:
with open(output_path, 'w') as f:
    with redirect_stdout(f):
        print("From ", dname)
        print("Running ", fname)
        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:7}{1:17} {2:6.4f}±{3:7.4f}/{4:2d}".format("Laser: ", "Qin =", Qin0, dQin, nQin))
        print("{0:7}{1:17} {2:6.4f}±{3:7.4f}/{4:2d}".format("", "Vin =", Vin0, dVin, nVin))
        print("{0:7}{1:17} {2:6.4f}±{3:7.4f}/{4:2d}".format("", "Rotation alpha = ", RotL0, dRotL, nRotL))
        print("{0:7}{1:15} {2:8.4f} {3:17}".format("", "=> DOP", ((Qin ** 2 + Vin ** 2) ** 0.5), " (degree of polarisation)"))

        print("Optic:        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:9.6f}±{2:9.6f}/{8:2d}, {3:7.4f}, {4:3.0f}±{5:3.0f}/{9:2d}, {6:7.4f}±{7:7.4f}/{10:2d}".format(
            "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(Type[TypeC],"calibrator is located", Loc[LocC])
        print("Rotation error epsilon is considered also for normal measurements = ", RotationErrorEpsilonForNormalMeasurements)
        print("The PBS incidence plane is ", dY[int(Y + 1)], "to the reference plane" )
        print("The laser polarisation in the reference plane is finally", dY2[int(Y + 1)], "=>", dY3)
        print("RS_RP_depend_on_TS_TP = ", RS_RP_depend_on_TS_TP)
        #  end of print actual system parameters
        # ******************************************************************************


        print()

        K0List = np.zeros(7)
        LDRsimxList = np.zeros(7)
        LDRCalList = 0.0, 0.004, 0.05, 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
        '''
        facIt = NI * eFacT
        facIr = NI * eFacR
        if (bPlotEtax):
            # check error signals
            # dIs are relative stdevs
            print("LDRCal, 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, Qin0, Vin0, 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

            if (bPlotEtax):
                # check error signals
                print( "{:0.2f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}".format(LDRCal, IoutTp * NCalT, IoutTm * NCalT, IoutRp * NCalR, IoutRm * NCalR, It * facIt, Ir * facIr, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr))
                #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},{10:9}".format(
            " GR", " GT", " HR", " HT", "  K(0.000)", "  K(0.004)", " K(0.05)", "  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},{10:9.5f}".format(
            GR0, GT0, HR0, HT0, K0List[0], K0List[1], K0List[2], K0List[3], K0List[4], K0List[5], K0List[6]))
        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"))

        aF11sim0 = np.zeros(5)
        LDRrange = np.zeros(5)
        LDRsim0 = np.zeros(5)
        LDRrange = [0.004, 0.02, 0.1, 0.3, 0.45]  # list
        LDRrange[0] = LDRtrue2  # value in the input file; default 0.004

        # 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, Qin0, Vin0, 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
            LDRsim0[i] = Ir / It
            # 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: NI, NCalT, NCalR, NILfac, nNCal, nNI, stdev(NI)/NI = {0:10.0f},{1:10.0f},{2:10.0f},{3:3.0f},{4:2.0f},{5:2.0f},{6:8.5f}".format(
            NI, NCalT, NCalR, NILfac, nNCal, nNI, 1.0 / NI ** 0.5))
        print()
        print()
        '''# 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\\' + OutputFile, 'r')
file = open(output_path, 'r')
print(file.read())
file.close()

# --- CALC again assumed truth with LDRCal0 and with assumed true parameters in input file to reset all 0-values
LDRtrue = LDRtrue2
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, Qin0, Vin0, RotL0, RotE0, RetE0, DiE0,
     RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0,
     ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0)
Etax0 = K0 * Eta0
Etapx0 = IoutRp0 / IoutTp0
Etamx0 = IoutRm0 / IoutTm0
'''
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
    LDRtrue = LDRtrue2
    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, Qin0, Vin0, RotL0, RotE0, RetE0, DiE0,
         RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0,
         ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0)
    Etax0 = K0 * Eta0
    Etapx0 = IoutRp0 / IoutTp0
    Etamx0 = IoutRm0 / IoutTm0

    # --- 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 *
         (nQin * 2 + 1) * (nVin * 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)
    LDRstd = np.zeros(5)
    LDRmean = np.zeros(5)
    LDRmedian = np.zeros(5)
    LDRskew = np.zeros(5)
    LDRkurt = np.zeros(5)
    LDRsimmin = np.zeros(5)
    LDRsimmax = np.zeros(5)
    LDRsimmean = np.zeros(5)

    F11min = np.zeros(5)
    F11max = np.zeros(5)
    Etaxmin = np.zeros(5)
    Etaxmax = np.zeros(5)

    aQin = np.zeros(N)
    aVin = 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))
    aLDRsim = np.zeros((5, N))
    aF11corr = np.zeros((5, N))
    aPLDR = np.zeros((5, N))
    aEtax = np.zeros((5, N))
    aEtapx = np.zeros((5, N))
    aEtamx = 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

    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

        aCal = (1. - LDRCal) / (1. + LDRCal)
        for iQin, iVin, iRotL, iRotE, iRetE, iDiE \
                in [(iQin, iVin, iRotL, iRotE, iRetE, iDiE)
                    for iQin in range(-nQin, nQin + 1)
                    for iVin in range(-nVin, nVin + 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 nQin > 0: Qin = Qin0 + iQin * dQin / nQin
            if nVin > 0: Vin = Vin0 + iVin * dVin / nVin
            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

            if ((Qin ** 2 + Vin ** 2) ** 0.5) > 1.0:
                print("Error: degree of polarisation of laser > 1. Check Qin and Vin! ")
                sys.exit()
            # 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 = Qin
            UinL = 0.
            VinL = Vin
            # 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.0 - TS
                    RP = 1.0 - 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.0 - ERaT) / (1.0 + ERaT)
                DaR = (1.0 - ERaR) / (1.0 + ERaR)
                TaT = 0.5 * (1.0 + ERaT)
                TaR = 0.5 * (1.0 + ERaR)

                S2aT = np.sin(np.deg2rad(h * 2.0 * RotaT))
                C2aT = np.cos(np.deg2rad(2.0 * RotaT))
                S2aR = np.sin(np.deg2rad(h * 2.0 * RotaR))
                C2aR = np.cos(np.deg2rad(2.0 * 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.0
                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 noise are not considered ?
                        # actually measured signal counts 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)
                            '''
                            dIt = ((It * NI * eFacT) ** -0.5)
                            dIr = ((Ir * NI * eFacR) ** -0.5)
                        else:
                            dIt = ((It * NI * eFacT) ** -0.5)
                            dIr = ((Ir * NI * eFacR) ** -0.5)
                            '''
                            # does this work? Why not as above?
                            dIt = ((NCalT * 2. * NILfac / TCalT ) ** -0.5)
                            dIr = ((NCalR * 2. * NILfac / TCalR) ** -0.5)
                            '''
                        # 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
                        # 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
                        aLDRsim[iLDR, iN] = LDRsim # LDRCorr # LDRsim # for test only
                        # aPLDR[iLDR, iN] = CalcPLDR(LDRCorr, BSR[iLDR], LDRm0)
                        aEtax[iLDR, iN] = Etax
                        aEtapx[iLDR, iN] = Etapx
                        aEtamx[iLDR, iN] = Etamx

                        aGHK[0, iN] = GR
                        aGHK[1, iN] = GT
                        aGHK[2, iN] = HR
                        aGHK[3, iN] = HT
                        aGHK[4, iN] = K

                        aLDRCal[iN] = iLDRCal
                        aQin[iN] = iQin
                        aVin[iN] = iVin
                        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 use:
        # 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:
            aXmean = np.zeros(2 * naX + 1)
            iLDR = iLDR + 1
            LDRmin[iLDR] = np.amin(aLDRcorr[iLDR, :])
            LDRmax[iLDR] = np.amax(aLDRcorr[iLDR, :])
            if (LDRmax[iLDR] > 10): LDRmax[iLDR] = 10
            if (LDRmin[iLDR] < -10): LDRmin[iLDR] = -10
            Rmin = LDRmin[iLDR] * 0.995  # np.min(aLDRcorr[iLDR,:])    * 0.995
            Rmax = LDRmax[iLDR] * 1.005  # np.max(aLDRcorr[iLDR,:])    * 1.005

            # Determine mean distance of all aXmean from each other for each iLDR
            meanDist = 0.0
            for iaX in range(-naX, naX + 1):
            # mean LDRCorr value for certain error (iaX) of parameter aVar
                aXmean[iaX + naX] = np.mean(aLDRcorr[iLDR, aX == iaX])
            # relative to absolute spread of LDRCorrs
            meanDist = (np.max(aXmean) - np.min(aXmean)) / (LDRmax[iLDR] - LDRmin[iLDR]) * 100

            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):
                # mean LDRCorr value for certain error (iaX) of parameter aVar
                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=10)
            plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2)
            plt.gca().set_title("{0:3.0f}%".format(meanDist))
            plt.gca().set_xlabel('LDRtrue', color="red")

        # plt.ylabel('frequency', fontsize=10)
        # plt.xlabel('LDRCorr', fontsize=10)
        # fig.tight_layout()
        fig.suptitle(LID + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar + ' error contribution', fontsize=14, y=1.10)
        # plt.show()
        # fig.savefig(LID + '_' + aVar + '.png', dpi=150, bbox_inches='tight', pad_inches=0)
        # plt.close
        return

    def PlotLDRsim(aVar, aX, X0, daX, iaX, naX):
        # aVar is the name of the parameter and aX is the subset of aLDRsim 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:
            aXmean = np.zeros(2 * naX + 1)
            iLDR = iLDR + 1
            LDRsimmin[iLDR] = np.amin(aLDRsim[iLDR, :])
            LDRsimmax[iLDR] = np.amax(aLDRsim[iLDR, :])
            # print("LDRsimmin[iLDR], LDRsimmax[iLDR] = ", LDRsimmin[iLDR], LDRsimmax[iLDR])
            # if (LDRsimmax[iLDR] > 10): LDRsimmax[iLDR] = 10
            # if (LDRsimmin[iLDR] < -10): LDRsimmin[iLDR] = -10
            Rmin = LDRsimmin[iLDR] * 0.995  # np.min(aLDRsim[iLDR,:])    * 0.995
            Rmax = LDRsimmax[iLDR] * 1.005  # np.max(aLDRsim[iLDR,:])    * 1.005
            # print("Rmin, Rmax = ", Rmin, Rmax)

            # Determine mean distance of all aXmean from each other for each iLDR
            meanDist = 0.0
            for iaX in range(-naX, naX + 1):
            # mean LDRCorr value for certain error (iaX) of parameter aVar
                aXmean[iaX + naX] = np.mean(aLDRsim[iLDR, aX == iaX])
            # relative to absolute spread of LDRCorrs
            meanDist = (np.max(aXmean) - np.min(aXmean)) / (LDRsimmax[iLDR] - LDRsimmin[iLDR]) * 100

            plt.subplot(1, 5, iLDR + 1)
            (n, bins, patches) = plt.hist(aLDRsim[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):
                # mean LDRCorr value for certain error (iaX) of parameter aVar
                plt.hist(aLDRsim[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=10)
            plt.plot([LDRsim0[iLDR], LDRsim0[iLDR]], [0, np.max(n)], 'r-', lw=2)
            plt.gca().set_title("{0:3.0f}%".format(meanDist))
            plt.gca().set_xlabel('LDRsim0', color="red")

        fig.suptitle('LDRsim - ' +LID + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar + ' error contribution', fontsize=14, y=1.10)
        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:
            aXmean = np.zeros(2 * naX + 1)
            iLDR = iLDR + 1
            Etaxmin = np.amin(aEtax[iLDR, :])
            Etaxmax = np.amax(aEtax[iLDR, :])
            Rmin = Etaxmin * 0.995  # np.min(aLDRcorr[iLDR,:])    * 0.995
            Rmax = Etaxmax * 1.005  # np.max(aLDRcorr[iLDR,:])    * 1.005

            # Determine mean distance of all aXmean from each other for each iLDR
            meanDist = 0.0
            for iaX in range(-naX, naX + 1):
            # mean Etax value for certain error (iaX) of parameter aVar
                aXmean[iaX + naX] = np.mean(aEtax[iLDR, aX == iaX])
            # relative to absolute spread of Etax
            meanDist = (np.max(aXmean) - np.min(aXmean)) / (Etaxmax - Etaxmin) * 100

            plt.subplot(1, 5, iLDR + 1)
            (n, bins, patches) = plt.hist(aEtax[iLDR, :],
                                          bins=50, 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=50, 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=10)
            plt.plot([Etax0, Etax0], [0, np.max(n)], 'r-', lw=2)
            plt.gca().set_title("{0:3.0f}%".format(meanDist))
            plt.gca().set_xlabel('Etax0', color="red")
        fig.suptitle('Etax - ' + LID + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar + ' error contribution', fontsize=14, y=1.10)
        return

    def PlotEtapx(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:
            aXmean = np.zeros(2 * naX + 1)
            iLDR = iLDR + 1
            Etapxmin = np.amin(aEtapx[iLDR, :])
            Etapxmax = np.amax(aEtapx[iLDR, :])
            Rmin = Etapxmin * 0.995  # np.min(aLDRcorr[iLDR,:])    * 0.995
            Rmax = Etapxmax * 1.005  # np.max(aLDRcorr[iLDR,:])    * 1.005

            # Determine mean distance of all aXmean from each other for each iLDR
            meanDist = 0.0
            for iaX in range(-naX, naX + 1):
            # mean Etapx value for certain error (iaX) of parameter aVar
                aXmean[iaX + naX] = np.mean(aEtapx[iLDR, aX == iaX])
            # relative to absolute spread of Etapx
            meanDist = (np.max(aXmean) - np.min(aXmean)) / (Etapxmax - Etapxmin) * 100

            plt.subplot(1, 5, iLDR + 1)
            (n, bins, patches) = plt.hist(aEtapx[iLDR, :],
                                          bins=50, log=False,
                                          range=[Rmin, Rmax],
                                          alpha=0.5, density=False, color='0.5', histtype='stepfilled')
            for iaX in range(-naX, naX + 1):
                plt.hist(aEtapx[iLDR, aX == iaX],
                         range=[Rmin, Rmax],
                         bins=50, 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=10)
            plt.plot([Etapx0, Etapx0], [0, np.max(n)], 'r-', lw=2)
            plt.gca().set_title("{0:3.0f}%".format(meanDist))
            plt.gca().set_xlabel('Etapx0', color="red")
        fig.suptitle('Etapx - ' + LID + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar + ' error contribution', fontsize=14, y=1.10)
        return

    def PlotEtamx(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:
            aXmean = np.zeros(2 * naX + 1)
            iLDR = iLDR + 1
            Etamxmin = np.amin(aEtamx[iLDR, :])
            Etamxmax = np.amax(aEtamx[iLDR, :])
            Rmin = Etamxmin * 0.995  # np.min(aLDRcorr[iLDR,:])    * 0.995
            Rmax = Etamxmax * 1.005  # np.max(aLDRcorr[iLDR,:])    * 1.005

            # Determine mean distance of all aXmean from each other for each iLDR
            meanDist = 0.0
            for iaX in range(-naX, naX + 1):
            # mean Etamx value for certain error (iaX) of parameter aVar
                aXmean[iaX + naX] = np.mean(aEtamx[iLDR, aX == iaX])
            # relative to absolute spread of Etamx
            meanDist = (np.max(aXmean) - np.min(aXmean)) / (Etamxmax - Etamxmin) * 100

            plt.subplot(1, 5, iLDR + 1)
            (n, bins, patches) = plt.hist(aEtamx[iLDR, :],
                                          bins=50, log=False,
                                          range=[Rmin, Rmax],
                                          alpha=0.5, density=False, color='0.5', histtype='stepfilled')
            for iaX in range(-naX, naX + 1):
                plt.hist(aEtamx[iLDR, aX == iaX],
                         range=[Rmin, Rmax],
                         bins=50, 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=10)
            plt.plot([Etamx0, Etamx0], [0, np.max(n)], 'r-', lw=2)
            plt.gca().set_title("{0:3.0f}%".format(meanDist))
            plt.gca().set_xlabel('Etamx0', color="red")
        fig.suptitle('Etamx - ' + LID + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar + ' error contribution', fontsize=14, y=1.10)
        return

    # calc contribution of the error of aVar = aX  to aY for each LDRtrue
    def Contribution(aVar, aX, X0, daX, iaX, naX, aY, Ysum, widthSum):
        # aVar is the name of the parameter and aX is the subset of aY which is coloured in the plot
        # example: Contribution("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP, aLDRcorr, DOLPcontr)
        iLDR = -1
        # Ysum, widthSum = np.zeros(5)
        meanDist = np.zeros(5) # iLDR
        widthDist = np.zeros(5) # iLDR
        for LDRTrue in LDRrange:
            aXmean = np.zeros(2 * naX + 1)
            aXwidth = np.zeros(2 * naX + 1)
            iLDR = iLDR + 1
            # total width of distribution
            aYmin = np.amin(aY[iLDR, :])
            aYmax = np.amax(aY[iLDR, :])
            aYwidth = aYmax - aYmin
            # Determine mean distance of all aXmean from each other for each iLDR
            for iaX in range(-naX, naX + 1):
            # mean LDRCorr value for all errors iaX of parameter aVar
                aXmean[iaX + naX] = np.mean(aY[iLDR, aX == iaX])
                aXwidth[iaX + naX] = np.max(aY[iLDR, aX == iaX]) - np.min(aY[iLDR, aX == iaX])
            # relative to absolute spread of LDRCorrs
            meanDist[iLDR] = (np.max(aXmean) - np.min(aXmean)) / aYwidth * 1000
            # meanDist[iLDR] = (aYwidth - aXwidth[naX]) / aYwidth * 1000
            widthDist[iLDR] = (np.max(aXwidth) - aXwidth[naX]) / aYwidth * 1000

        print("{:12}{:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f}    {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f}"\
              .format(aVar,meanDist[0],meanDist[1],meanDist[2],meanDist[3],meanDist[4],widthDist[0],widthDist[1],widthDist[2],widthDist[3],widthDist[4]))
        Ysum = Ysum + meanDist
        widthSum = widthSum + widthDist
        return(Ysum, widthSum)

        # print(.format(LDRrangeA[iLDR],))

    # error contributions to a certain output aY; loop over all variables
    def Contribution_aY(aYvar, aY):
        Ysum = np.zeros(5)
        widthSum = np.zeros(5)
        # meanDist = np.zeros(5) # iLDR
        LDRrangeA = np.array(LDRrange)
        print()
        print(aYvar + ": contribution to the total error (per mill)")
        print("          of individual parameter errors        of combined parameter errors")
        print(" at LDRtrue {:5.3f} {:5.3f} {:5.3f} {:5.3f} {:5.3f}    {:5.3f} {:5.3f} {:5.3f} {:5.3f} {:5.3f}"\
              .format(LDRrangeA[0],LDRrangeA[1],LDRrangeA[2],LDRrangeA[3],LDRrangeA[4],LDRrangeA[0],LDRrangeA[1],LDRrangeA[2],LDRrangeA[3],LDRrangeA[4]))
        print()
        if (nQin > 0): Ysum, widthSum = Contribution("Qin", aQin, Qin0, dQin, iQin, nQin, aY, Ysum, widthSum)
        if (nVin > 0): Ysum, widthSum = Contribution("Vin", aVin, Vin0, dVin, iVin, nVin, aY, Ysum, widthSum)
        if (nRotL > 0): Ysum, widthSum = Contribution("RotL", aRotL, RotL0, dRotL, iRotL, nRotL, aY, Ysum, widthSum)
        if (nRetE > 0): Ysum, widthSum = Contribution("RetE", aRetE, RetE0, dRetE, iRetE, nRetE, aY, Ysum, widthSum)
        if (nRotE > 0): Ysum, widthSum = Contribution("RotE", aRotE, RotE0, dRotE, iRotE, nRotE, aY, Ysum, widthSum)
        if (nDiE > 0): Ysum, widthSum = Contribution("DiE", aDiE, DiE0, dDiE, iDiE, nDiE, aY, Ysum, widthSum)
        if (nRetO > 0): Ysum, widthSum = Contribution("RetO", aRetO, RetO0, dRetO, iRetO, nRetO, aY, Ysum, widthSum)
        if (nRotO > 0): Ysum, widthSum = Contribution("RotO", aRotO, RotO0, dRotO, iRotO, nRotO, aY, Ysum, widthSum)
        if (nDiO > 0): Ysum, widthSum = Contribution("DiO", aDiO, DiO0, dDiO, iDiO, nDiO, aY, Ysum, widthSum)
        if (nDiC > 0): Ysum, widthSum = Contribution("DiC", aDiC, DiC0, dDiC, iDiC, nDiC, aY, Ysum, widthSum)
        if (nRotC > 0): Ysum, widthSum = Contribution("RotC", aRotC, RotC0, dRotC, iRotC, nRotC, aY, Ysum, widthSum)
        if (nRetC > 0): Ysum, widthSum = Contribution("RetC", aRetC, RetC0, dRetC, iRetC, nRetC, aY, Ysum, widthSum)
        if (nTP > 0): Ysum, widthSum = Contribution("TP", aTP, TP0, dTP, iTP, nTP, aY, Ysum, widthSum)
        if (nTS > 0): Ysum, widthSum = Contribution("TS", aTS, TS0, dTS, iTS, nTS, aY, Ysum, widthSum)
        if (nRP > 0): Ysum, widthSum = Contribution("RP", aRP, RP0, dRP, iRP, nRP, aY, Ysum, widthSum)
        if (nRS > 0): Ysum, widthSum = Contribution("RS", aRS, RS0, dRS, iRS, nRS, aY, Ysum, widthSum)
        if (nRetT > 0): Ysum, widthSum = Contribution("RetT", aRetT, RetT0, dRetT, iRetT, nRetT, aY, Ysum, widthSum)
        if (nRetR > 0): Ysum, widthSum = Contribution("RetR", aRetR, RetR0, dRetR, iRetR, nRetR, aY, Ysum, widthSum)
        if (nERaT > 0): Ysum, widthSum = Contribution("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT, aY, Ysum, widthSum)
        if (nERaR > 0): Ysum, widthSum = Contribution("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR, aY, Ysum, widthSum)
        if (nRotaT > 0): Ysum, widthSum = Contribution("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT, aY, Ysum, widthSum)
        if (nRotaR > 0): Ysum, widthSum = Contribution("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR, aY, Ysum, widthSum)
        if (nLDRCal > 0): Ysum, widthSum = Contribution("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal, aY, Ysum, widthSum)
        if (nTCalT > 0): Ysum, widthSum = Contribution("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT, aY, Ysum, widthSum)
        if (nTCalR > 0): Ysum, widthSum = Contribution("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR, aY, Ysum, widthSum)
        if (nNCal > 0): Ysum, widthSum = Contribution("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal, aY, Ysum, widthSum)
        if (nNCal > 0): Ysum, widthSum = Contribution("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal, aY, Ysum, widthSum)
        if (nNCal > 0): Ysum, widthSum = Contribution("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal, aY, Ysum, widthSum)
        if (nNCal > 0): Ysum, widthSum = Contribution("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal, aY, Ysum, widthSum)
        if (nNI > 0): Ysum, widthSum = Contribution("SigNoiseIt", aNIt, 0, 1, iNIt, nNI, aY, Ysum, widthSum)
        if (nNI > 0): Ysum, widthSum = Contribution("SigNoiseIr", aNIr, 0, 1, iNIr, nNI, aY, Ysum, widthSum)
        print("{:12}{:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f}    {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f}"\
              .format("Sum ",Ysum[0],Ysum[1],Ysum[2],Ysum[3],Ysum[4],widthSum[0],widthSum[1],widthSum[2],widthSum[3],widthSum[4]))


    # Plot LDR histograms
    if (nQin > 0): PlotSubHist("Qin", aQin, Qin0, dQin, iQin, nQin)
    if (nVin > 0): PlotSubHist("Vin", aVin, Vin0, dVin, iVin, nVin)
    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 LDRmin, LDRmax
    iLDR = -1
    for LDRTrue in LDRrange:
        iLDR = iLDR + 1
        LDRmin[iLDR] = np.amin(aLDRcorr[iLDR, :])
        LDRmax[iLDR] = np.amax(aLDRcorr[iLDR, :])
        LDRstd[iLDR] = np.std(aLDRcorr[iLDR, :])
        LDRmean[iLDR] = np.mean(aLDRcorr[iLDR, :])
        LDRmedian[iLDR] = np.median(aLDRcorr[iLDR, :])
        LDRskew[iLDR] = skew(aLDRcorr[iLDR, :],bias=False)
        LDRkurt[iLDR] = kurtosis(aLDRcorr[iLDR, :],fisher=True,bias=False)

    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\\' + OutputFile, 'a') as f:
    # with open('output_files\\' + LID + '-' + InputFile[0:-3] + '-LDR_min_max.dat', 'w') as f:
        with redirect_stdout(f):
            print("Lidar ID: " + LID)
            print()
            print("minimum and maximum values of the distributions of possibly measured LDR for different LDRtrue")
            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]))
            print()
            # Print LDR statistics
            print("LDRtrue ,  mean  ,  median,    max-mean,  min-mean, std,   excess_kurtosis, skewness")
            iLDR = -1
            LDRrangeA = np.array(LDRrange)
            for LDRTrue in LDRrange:
                iLDR = iLDR + 1
                print("{0:8.5f},{1:8.5f},{2:8.5f},    {3:8.5f},{4:8.5f},{5:8.5f},   {6:8.5f},{7:8.5f}"\
                      .format(LDRrangeA[iLDR], LDRmean[iLDR], LDRmedian[iLDR], LDRmax[iLDR]-LDRrangeA[iLDR], \
                              LDRmin[iLDR]-LDRrangeA[iLDR], LDRstd[iLDR], LDRkurt[iLDR], LDRskew[iLDR]))
            print()
            # Calculate and print statistics for calibration factors
            print("minimum and maximum values of the distributions of signal ratios and calibration factors for different LDRtrue")
            iLDR = -1
            LDRrangeA = np.array(LDRrange)
            print("LDRtrue  , LDRsim, (max-min)/2, relerr")
            for LDRTrue in LDRrange:
                iLDR = iLDR + 1
                LDRsimmin[iLDR] = np.amin(aLDRsim[iLDR, :])
                LDRsimmax[iLDR] = np.amax(aLDRsim[iLDR, :])
                # LDRsimstd = np.std(aLDRsim[iLDR, :])
                LDRsimmean[iLDR] = np.mean(aLDRsim[iLDR, :])
                # LDRsimmedian = np.median(aLDRsim[iLDR, :])
                print("{0:8.5f}, {1:8.5f}, {2:8.5f}, {3:8.5f}".format(LDRrangeA[iLDR],LDRsimmean[iLDR],(LDRsimmax[iLDR]-LDRsimmin[iLDR])/2,(LDRsimmax[iLDR]-LDRsimmin[iLDR])/2/LDRsimmean[iLDR]))
            iLDR = -1
            print("LDRtrue  , Etax   , (max-min)/2, relerr")
            for LDRTrue in LDRrange:
                iLDR = iLDR + 1
                Etaxmin = np.amin(aEtax[iLDR, :])
                Etaxmax = np.amax(aEtax[iLDR, :])
                # Etaxstd = np.std(aEtax[iLDR, :])
                Etaxmean = np.mean(aEtax[iLDR, :])
                # Etaxmedian = np.median(aEtax[iLDR, :])
                print("{0:8.5f}, {1:8.5f}, {2:8.5f}, {3:8.5f}".format(LDRrangeA[iLDR], Etaxmean, (Etaxmax-Etaxmin)/2, (Etaxmax-Etaxmin)/2/Etaxmean))
            iLDR = -1
            print("LDRtrue  , Etapx  , (max-min)/2, relerr")
            for LDRTrue in LDRrange:
                iLDR = iLDR + 1
                Etapxmin = np.amin(aEtapx[iLDR, :])
                Etapxmax = np.amax(aEtapx[iLDR, :])
                # Etapxstd = np.std(aEtapx[iLDR, :])
                Etapxmean = np.mean(aEtapx[iLDR, :])
                # Etapxmedian = np.median(aEtapx[iLDR, :])
                print("{0:8.5f}, {1:8.5f}, {2:8.5f}, {3:8.5f}".format(LDRrangeA[iLDR], Etapxmean, (Etapxmax-Etapxmin)/2, (Etapxmax-Etapxmin)/2/Etapxmean))
            iLDR = -1
            print("LDRtrue  , Etamx  , (max-min)/2, relerr")
            for LDRTrue in LDRrange:
                iLDR = iLDR + 1
                Etamxmin = np.amin(aEtamx[iLDR, :])
                Etamxmax = np.amax(aEtamx[iLDR, :])
                # Etamxstd = np.std(aEtamx[iLDR, :])
                Etamxmean = np.mean(aEtamx[iLDR, :])
                # Etamxmedian = np.median(aEtamx[iLDR, :])
                print("{0:8.5f}, {1:8.5f}, {2:8.5f}, {3:8.5f}".format(LDRrangeA[iLDR], Etamxmean, (Etamxmax-Etamxmin)/2, (Etamxmax-Etamxmin)/2/Etamxmean))

    # Print LDR statistics
    print("LDRtrue ,  mean  ,  median,    max-mean,  min-mean, std,   excess_kurtosis, skewness")
    iLDR = -1
    LDRrangeA = np.array(LDRrange)
    for LDRTrue in LDRrange:
        iLDR = iLDR + 1
        print("{0:8.5f},{1:8.5f},{2:8.5f},    {3:8.5f},{4:8.5f},{5:8.5f},   {6:8.5f},{7:8.5f}".format(LDRrangeA[iLDR], LDRmean[iLDR], LDRmedian[iLDR], LDRmax[iLDR]-LDRrangeA[iLDR], LDRmin[iLDR]-LDRrangeA[iLDR], LDRstd[iLDR],LDRkurt[iLDR],LDRskew[iLDR]))


    with open('output_files\\' + OutputFile, 'a') as f:
    # with open('output_files\\' + LID + '-' + InputFile[0:-3] + '-LDR_min_max.dat', 'a') as f:
        with redirect_stdout(f):
            Contribution_aY("LDRCorr", aLDRcorr)
            Contribution_aY("LDRsim", aLDRsim)
            Contribution_aY("EtaX, D90", aEtax)
            Contribution_aY("Etapx, +45°", aEtapx)
            Contribution_aY("Etamx -45°", aEtamx)


    # Plot other histograms
    if (bPlotEtax):

        if (nQin > 0): PlotLDRsim("Qin", aQin, Qin0, dQin, iQin, nQin)
        if (nVin > 0): PlotLDRsim("Vin", aVin, Vin0, dVin, iVin, nVin)
        if (nRotL > 0): PlotLDRsim("RotL", aRotL, RotL0, dRotL, iRotL, nRotL)
        if (nRetE > 0): PlotLDRsim("RetE", aRetE, RetE0, dRetE, iRetE, nRetE)
        if (nRotE > 0): PlotLDRsim("RotE", aRotE, RotE0, dRotE, iRotE, nRotE)
        if (nDiE > 0): PlotLDRsim("DiE", aDiE, DiE0, dDiE, iDiE, nDiE)
        if (nRetO > 0): PlotLDRsim("RetO", aRetO, RetO0, dRetO, iRetO, nRetO)
        if (nRotO > 0): PlotLDRsim("RotO", aRotO, RotO0, dRotO, iRotO, nRotO)
        if (nDiO > 0): PlotLDRsim("DiO", aDiO, DiO0, dDiO, iDiO, nDiO)
        if (nDiC > 0): PlotLDRsim("DiC", aDiC, DiC0, dDiC, iDiC, nDiC)
        if (nRotC > 0): PlotLDRsim("RotC", aRotC, RotC0, dRotC, iRotC, nRotC)
        if (nRetC > 0): PlotLDRsim("RetC", aRetC, RetC0, dRetC, iRetC, nRetC)
        if (nTP > 0): PlotLDRsim("TP", aTP, TP0, dTP, iTP, nTP)
        if (nTS > 0): PlotLDRsim("TS", aTS, TS0, dTS, iTS, nTS)
        if (nRP > 0): PlotLDRsim("RP", aRP, RP0, dRP, iRP, nRP)
        if (nRS > 0): PlotLDRsim("RS", aRS, RS0, dRS, iRS, nRS)
        if (nRetT > 0): PlotLDRsim("RetT", aRetT, RetT0, dRetT, iRetT, nRetT)
        if (nRetR > 0): PlotLDRsim("RetR", aRetR, RetR0, dRetR, iRetR, nRetR)
        if (nERaT > 0): PlotLDRsim("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT)
        if (nERaR > 0): PlotLDRsim("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR)
        if (nRotaT > 0): PlotLDRsim("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT)
        if (nRotaR > 0): PlotLDRsim("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR)
        if (nLDRCal > 0): PlotLDRsim("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal)
        if (nTCalT > 0): PlotLDRsim("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT)
        if (nTCalR > 0): PlotLDRsim("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR)
        if (nNCal > 0): PlotLDRsim("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal)
        if (nNCal > 0): PlotLDRsim("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal)
        if (nNCal > 0): PlotLDRsim("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal)
        if (nNCal > 0): PlotLDRsim("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal)
        if (nNI > 0): PlotLDRsim("SigNoiseIt", aNIt, 0, 1, iNIt, nNI)
        if (nNI > 0): PlotLDRsim("SigNoiseIr", aNIr, 0, 1, iNIr, nNI)
        plt.show()
        plt.close
        print("---------------------------------------...producing more plots...------------------------------------------------------------------")

        if (nQin > 0): PlotEtax("Qin", aQin, Qin0, dQin, iQin, nQin)
        if (nVin > 0): PlotEtax("Vin", aVin, Vin0, dVin, iVin, nVin)
        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
        print("---------------------------------------...producing more plots...------------------------------------------------------------------")

        if (nQin > 0): PlotEtapx("Qin", aQin, Qin0, dQin, iQin, nQin)
        if (nVin > 0): PlotEtapx("Vin", aVin, Vin0, dVin, iVin, nVin)
        if (nRotL > 0): PlotEtapx("RotL", aRotL, RotL0, dRotL, iRotL, nRotL)
        if (nRetE > 0): PlotEtapx("RetE", aRetE, RetE0, dRetE, iRetE, nRetE)
        if (nRotE > 0): PlotEtapx("RotE", aRotE, RotE0, dRotE, iRotE, nRotE)
        if (nDiE > 0): PlotEtapx("DiE", aDiE, DiE0, dDiE, iDiE, nDiE)
        if (nRetO > 0): PlotEtapx("RetO", aRetO, RetO0, dRetO, iRetO, nRetO)
        if (nRotO > 0): PlotEtapx("RotO", aRotO, RotO0, dRotO, iRotO, nRotO)
        if (nDiO > 0): PlotEtapx("DiO", aDiO, DiO0, dDiO, iDiO, nDiO)
        if (nDiC > 0): PlotEtapx("DiC", aDiC, DiC0, dDiC, iDiC, nDiC)
        if (nRotC > 0): PlotEtapx("RotC", aRotC, RotC0, dRotC, iRotC, nRotC)
        if (nRetC > 0): PlotEtapx("RetC", aRetC, RetC0, dRetC, iRetC, nRetC)
        if (nTP > 0): PlotEtapx("TP", aTP, TP0, dTP, iTP, nTP)
        if (nTS > 0): PlotEtapx("TS", aTS, TS0, dTS, iTS, nTS)
        if (nRP > 0): PlotEtapx("RP", aRP, RP0, dRP, iRP, nRP)
        if (nRS > 0): PlotEtapx("RS", aRS, RS0, dRS, iRS, nRS)
        if (nRetT > 0): PlotEtapx("RetT", aRetT, RetT0, dRetT, iRetT, nRetT)
        if (nRetR > 0): PlotEtapx("RetR", aRetR, RetR0, dRetR, iRetR, nRetR)
        if (nERaT > 0): PlotEtapx("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT)
        if (nERaR > 0): PlotEtapx("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR)
        if (nRotaT > 0): PlotEtapx("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT)
        if (nRotaR > 0): PlotEtapx("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR)
        if (nLDRCal > 0): PlotEtapx("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal)
        if (nTCalT > 0): PlotEtapx("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT)
        if (nTCalR > 0): PlotEtapx("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR)
        if (nNCal > 0): PlotEtapx("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal)
        if (nNCal > 0): PlotEtapx("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal)
        if (nNCal > 0): PlotEtapx("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal)
        if (nNCal > 0): PlotEtapx("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal)
        if (nNI > 0): PlotEtapx("SigNoiseIt", aNIt, 0, 1, iNIt, nNI)
        if (nNI > 0): PlotEtapx("SigNoiseIr", aNIr, 0, 1, iNIr, nNI)
        plt.show()
        plt.close
        print("---------------------------------------...producing more plots...------------------------------------------------------------------")

        if (nQin > 0): PlotEtamx("Qin", aQin, Qin0, dQin, iQin, nQin)
        if (nVin > 0): PlotEtamx("Vin", aVin, Vin0, dVin, iVin, nVin)
        if (nRotL > 0): PlotEtamx("RotL", aRotL, RotL0, dRotL, iRotL, nRotL)
        if (nRetE > 0): PlotEtamx("RetE", aRetE, RetE0, dRetE, iRetE, nRetE)
        if (nRotE > 0): PlotEtamx("RotE", aRotE, RotE0, dRotE, iRotE, nRotE)
        if (nDiE > 0): PlotEtamx("DiE", aDiE, DiE0, dDiE, iDiE, nDiE)
        if (nRetO > 0): PlotEtamx("RetO", aRetO, RetO0, dRetO, iRetO, nRetO)
        if (nRotO > 0): PlotEtamx("RotO", aRotO, RotO0, dRotO, iRotO, nRotO)
        if (nDiO > 0): PlotEtamx("DiO", aDiO, DiO0, dDiO, iDiO, nDiO)
        if (nDiC > 0): PlotEtamx("DiC", aDiC, DiC0, dDiC, iDiC, nDiC)
        if (nRotC > 0): PlotEtamx("RotC", aRotC, RotC0, dRotC, iRotC, nRotC)
        if (nRetC > 0): PlotEtamx("RetC", aRetC, RetC0, dRetC, iRetC, nRetC)
        if (nTP > 0): PlotEtamx("TP", aTP, TP0, dTP, iTP, nTP)
        if (nTS > 0): PlotEtamx("TS", aTS, TS0, dTS, iTS, nTS)
        if (nRP > 0): PlotEtamx("RP", aRP, RP0, dRP, iRP, nRP)
        if (nRS > 0): PlotEtamx("RS", aRS, RS0, dRS, iRS, nRS)
        if (nRetT > 0): PlotEtamx("RetT", aRetT, RetT0, dRetT, iRetT, nRetT)
        if (nRetR > 0): PlotEtamx("RetR", aRetR, RetR0, dRetR, iRetR, nRetR)
        if (nERaT > 0): PlotEtamx("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT)
        if (nERaR > 0): PlotEtamx("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR)
        if (nRotaT > 0): PlotEtamx("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT)
        if (nRotaR > 0): PlotEtamx("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR)
        if (nLDRCal > 0): PlotEtamx("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal)
        if (nTCalT > 0): PlotEtamx("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT)
        if (nTCalR > 0): PlotEtamx("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR)
        if (nNCal > 0): PlotEtamx("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal)
        if (nNCal > 0): PlotEtamx("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal)
        if (nNCal > 0): PlotEtamx("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal)
        if (nNCal > 0): PlotEtamx("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal)
        if (nNI > 0): PlotEtamx("SigNoiseIt", aNIt, 0, 1, iNIt, nNI)
        if (nNI > 0): PlotEtamx("SigNoiseIr", aNIr, 0, 1, iNIr, nNI)
        plt.show()
        plt.close

        # Print Etax statistics
        Etaxmin = np.amin(aEtax[1, :])
        Etaxmax = np.amax(aEtax[1, :])
        Etaxstd = np.std(aEtax[1, :])
        Etaxmean = np.mean(aEtax[1, :])
        Etaxmedian = np.median(aEtax[1, :])
        print("Etax      , max-mean, min-mean, median, mean ± std, eta")
        print("{0:8.5f} ±({1:8.5f},{2:8.5f}),{3:8.5f},{4:8.5f}±{5:8.5f},{6:8.5f}".format(Etax0, Etaxmax-Etax0, Etaxmin-Etax0, Etaxmedian, Etaxmean, Etaxstd, Etax0 / K0))
        print()

        # Calculate and print statistics for calibration factors
        iLDR = -1
        LDRrangeA = np.array(LDRrange)
        print("LDR...., LDRsim, (max-min)/2, relerr")
        for LDRTrue in LDRrange:
            iLDR = iLDR + 1
            LDRsimmin[iLDR] = np.amin(aLDRsim[iLDR, :])
            LDRsimmax[iLDR] = np.amax(aLDRsim[iLDR, :])
            # LDRsimstd = np.std(aLDRsim[iLDR, :])
            LDRsimmean[iLDR] = np.mean(aLDRsim[iLDR, :])
            # LDRsimmedian = np.median(aLDRsim[iLDR, :])
            print("{0:8.5f}, {1:8.5f}, {2:8.5f}, {3:8.5f}".format(LDRrangeA[iLDR], LDRsimmean[iLDR], (LDRsimmax[iLDR]-LDRsimmin[iLDR])/2,  (LDRsimmax[iLDR]-LDRsimmin[iLDR])/2/LDRsimmean[iLDR]))
        iLDR = -1
        print("LDR...., Etax   , (max-min)/2, relerr")
        for LDRTrue in LDRrange:
            iLDR = iLDR + 1
            Etaxmin = np.amin(aEtax[iLDR, :])
            Etaxmax = np.amax(aEtax[iLDR, :])
            # Etaxstd = np.std(aEtax[iLDR, :])
            Etaxmean = np.mean(aEtax[iLDR, :])
            # Etaxmedian = np.median(aEtax[iLDR, :])
            print("{0:8.5f}, {1:8.5f}, {2:8.5f}, {3:8.5f}".format(LDRrangeA[iLDR], Etaxmean, (Etaxmax-Etaxmin)/2, (Etaxmax-Etaxmin)/2/Etaxmean))
        iLDR = -1
        print("LDR...., Etapx  , (max-min)/2, relerr")
        for LDRTrue in LDRrange:
            iLDR = iLDR + 1
            Etapxmin = np.amin(aEtapx[iLDR, :])
            Etapxmax = np.amax(aEtapx[iLDR, :])
            # Etapxstd = np.std(aEtapx[iLDR, :])
            Etapxmean = np.mean(aEtapx[iLDR, :])
            # Etapxmedian = np.median(aEtapx[iLDR, :])
            print("{0:8.5f}, {1:8.5f}, {2:8.5f}, {3:8.5f}".format(LDRrangeA[iLDR], Etapxmean, (Etapxmax-Etapxmin)/2, (Etapxmax-Etapxmin)/2/Etapxmean))
        iLDR = -1
        print("LDR...., Etamx  , (max-min)/2, relerr")
        for LDRTrue in LDRrange:
            iLDR = iLDR + 1
            Etamxmin = np.amin(aEtamx[iLDR, :])
            Etamxmax = np.amax(aEtamx[iLDR, :])
            # Etamxstd = np.std(aEtamx[iLDR, :])
            Etamxmean = np.mean(aEtamx[iLDR, :])
            # Etamxmedian = np.median(aEtamx[iLDR, :])
            print("{0:8.5f}, {1:8.5f}, {2:8.5f}, {3:8.5f}".format(LDRrangeA[iLDR], Etamxmean, (Etamxmax-Etamxmin)/2, (Etamxmax-Etamxmin)/2/Etamxmean))

    f.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 (nQin > 0): PlotSubHistF11("Qin", aQin, Qin0, dQin, iQin, nQin)
    if (nVin > 0): PlotSubHistF11("Vin", aVin, Vin0, dVin, iVin, nVin)
    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 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')
#
'''