Mercurial > public > atmospheric_lidar_ghk / file revision

# -*- coding: utf-8 -*-
"""
Copyright 2016 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.4.2
"""
# !/usr/bin/env python3

# Comment:  The code 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

import os
import sys

import numpy as np

# Comment: the seaborn library makes nicer plots, but the code works also without it.
try:
    import seaborn as sns

    sns_loaded = True
except ImportError:
    sns_loaded = False

import matplotlib.pyplot as plt
from time import clock

# 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

# Do you want to calculate the errors? If not, just the GHK-parameters are determined.
Error_Calc = True
# ---- Initial definition of variables; the actual values will be read in with exec(open('./optic_input.py').read()) below
LID = "internal"
EID = "internal"
# --- IL Laser IL and +-Uncertainty
bL = 1.  # degree of linear polarization; default 1
RotL, dRotL, nRotL = 0.0, 0.0, 1  # alpha; rotation of laser polarization in degrees; default 0
# --- ME Emitter and +-Uncertainty
DiE, dDiE, nDiE = 0., 0.00, 1  # Diattenuation
TiE = 1.  # Unpolarized transmittance
RetE, dRetE, nRetE = 0., 180.0, 0  # Retardance in degrees
RotE, dRotE, nRotE = 0., 0.0, 0  # beta: Rotation of optical element in degrees
# --- MO Receiver Optics including telescope
DiO, dDiO, nDiO = -0.055, 0.003, 1
TiO = 0.9
RetO, dRetO, nRetO = 0., 180.0, 2
RotO, dRotO, nRotO = 0., 0.1, 1  # gamma
# --- PBS MT transmitting path defined with (TS,TP);  and +-Uncertainty
TP, dTP, nTP = 0.98, 0.02, 1
TS, dTS, nTS = 0.001, 0.001, 1
TiT = 0.5 * (TP + TS)
DiT = (TP - TS) / (TP + TS)
# PolFilter
RetT, dRetT, nRetT = 0., 180., 0
ERaT, dERaT, nERaT = 0.001, 0.001, 1
RotaT, dRotaT, nRotaT = 0., 3., 1
DaT = (1 - ERaT) / (1 + ERaT)
TaT = 0.5 * (1 + ERaT)
# --- PBS MR reflecting path defined with (RS,RP);  and +-Uncertainty
RS_RP_depend_on_TS_TP = False
if (RS_RP_depend_on_TS_TP):
    RP, dRP, nRP = 1 - TP, 0.00, 0
    RS, dRS, nRS = 1 - TS, 0.00, 0
else:
    RP, dRP, nRP = 0.05, 0.01, 1
    RS, dRS, nRS = 0.98, 0.01, 1
TiR = 0.5 * (RP + RS)
DiR = (RP - RS) / (RP + RS)
# PolFilter
RetR, dRetR, nRetR = 0., 180., 0
ERaR, dERaR, nERaR = 0.001, 0.001, 1
RotaR, dRotaR, nRotaR = 90., 3., 1
DaR = (1 - ERaR) / (1 + ERaR)
TaR = 0.5 * (1 + ERaR)

# Parellel signal detected in the transmitted channel => Y = 1, or in the reflected channel => Y = -1
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

TypeC = 3  # linear polarizer calibrator
# example with extinction ratio 0.001
DiC, dDiC, nDiC = 1.0, 0., 0  # ideal 1.0
TiC = 0.5  # ideal 0.5
RetC, dRetC, nRetC = 0., 0., 0
RotC, dRotC, nRotC = 0.0, 0.1, 0  # constant calibrator offset epsilon
RotationErrorEpsilonForNormalMeasurements = False  # is in general False for TypeC == 3 calibrator

# Rotation error without calibrator: if False, then epsilon = 0 for normal measurements
RotationErrorEpsilonForNormalMeasurements = True

# LDRCal assumed atmospheric linear depolarization ratio during the calibration measurements (first guess)
LDRCal0, dLDRCal, nLDRCal = 0.25, 0.04, 1
LDRCal = LDRCal0
# measured LDRm will be corrected with calculated parameters
LDRmeas = 0.015
# LDRtrue for simulation of measurement => LDRsim
LDRtrue = 0.5
LDRtrue2 = 0.004

# Initialize other values to 0
ER, nER, dER = 0.001, 0, 0.001
K = 0.
Km = 0.
Kp = 0.
LDRcorr = 0.
Eta = 0.
Ir = 0.
It = 0.
h = 1.

Loc = ['', 'behind laser', 'behind emitter', 'before receiver', 'before PBS']
Type = ['', 'mechanical rotator', 'hwp rotator', 'linear polarizer', 'qwp rotator', 'circular polarizer',
        'real HWP +-22.5°']
dY = ['reflected channel', '', 'transmitted channel']

#  end of initial definition of variables
# *******************************************************************************************************************************

# --- Read actual lidar system parameters from ./optic_input.py  (must be in the same directory)
InputFile = 'optic_input_example_lidar.py'
# InputFile = 'optic_input_ver6e_POLIS_355.py'
# InputFile = 'optic_input_ver6e_POLIS_355_JA.py'
# InputFile = 'optic_input_ver6c_POLIS_532.py'
# InputFile = 'optic_input_ver6e_POLIS_532.py'
# InputFile = 'optic_input_ver8c_POLIS_532.py'
# InputFile = 'optic_input_ver6e_MUSA.py'
# InputFile = 'optic_input_ver6e_MUSA_JA.py'
# InputFile = 'optic_input_ver6e_PollyXTSea.py'
# InputFile = 'optic_input_ver6e_PollyXTSea_JA.py'
# InputFile = 'optic_input_ver6e_PollyXT_RALPH.py'
# InputFile = 'optic_input_ver8c_PollyXT_RALPH.py'
# InputFile = 'optic_input_ver8c_PollyXT_RALPH_2.py'
# InputFile = 'optic_input_ver8c_PollyXT_RALPH_3.py'
# InputFile = 'optic_input_ver8c_PollyXT_RALPH_4.py'
# InputFile = 'optic_input_ver8c_PollyXT_RALPH_5.py'
# InputFile = 'optic_input_ver8c_PollyXT_RALPH_6.py'
# InputFile = 'optic_input_ver8c_PollyXT_RALPH_7.py'
# InputFile = 'optic_input_ver8a_MOHP_DPL_355.py'
# InputFile = 'optic_input_ver9_MOHP_DPL_355.py'
# InputFile = 'optic_input_ver6e_RALI.py'
# InputFile = 'optic_input_ver6e_RALI_JA.py'
# InputFile = 'optic_input_ver6e_RALI_new.py'
# InputFile = 'optic_input_ver6e_RALI_act.py'
# InputFile = 'optic_input_ver6e_MULHACEN.py'
# InputFile = 'optic_input_ver6e_MULHACEN_JA.py'
# InputFile = 'optic_input_ver6e-IPRAL.py'
# InputFile = 'optic_input_ver6e-IPRAL_JA.py'
# InputFile = 'optic_input_ver6e-LB21.py'
# InputFile = 'optic_input_ver6e-LB21_JA.py'
# InputFile = 'optic_input_ver6e_Bertha_b_355.py'
# InputFile = 'optic_input_ver6e_Bertha_b_532.py'
# InputFile = 'optic_input_ver6e_Bertha_b_1064.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
# bL = 0.8
## --- Errors
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
# LDRCal0,dLDRCal,nLDRCal=LDRCal,dLDRCal,0
# ---------- 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

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


# --- this subroutine is for the calculation with certain fixed parameters -----------------------------------------------------
def Calc(RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC, TP, TS, RP, RS, ERaT, RotaT, RetT, ERaR, RotaR, RetR,
         LDRCal):
    # ---- Do the calculations of bra-ket vectors
    h = -1. if TypeC == 2 else 1
    # from input file:  assumed LDRCal for calibration measurements
    aCal = (1. - LDRCal) / (1 + LDRCal)
    # from input file: measured LDRm and true LDRtrue, LDRtrue2  =>
    # ameas = (1.-LDRmeas)/(1+LDRmeas)
    atrue = (1. - LDRtrue) / (1 + LDRtrue)
    # atrue2 = (1.-LDRtrue2)/(1+LDRtrue2)

    # angles of emitter and laser and calibrator and receiver optics
    # RotL = alpha, RotE = beta, RotO = gamma, RotC = epsilon
    S2a = np.sin(2 * np.deg2rad(RotL))
    C2a = np.cos(2 * np.deg2rad(RotL))
    S2b = np.sin(2 * np.deg2rad(RotE))
    C2b = np.cos(2 * np.deg2rad(RotE))
    S2ab = np.sin(np.deg2rad(2 * RotL - 2 * RotE))
    C2ab = np.cos(np.deg2rad(2 * RotL - 2 * RotE))
    S2g = np.sin(np.deg2rad(2 * RotO))
    C2g = np.cos(np.deg2rad(2 * RotO))

    # Laser with Degree of linear polarization DOLP = bL
    IinL = 1.
    QinL = bL
    UinL = 0.
    VinL = (1. - bL ** 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
    # For POLLY_XTs
    if (RS_RP_depend_on_TS_TP):
        RS = 1 - TS
        RP = 1 - TP
    TiT = 0.5 * (TP + TS)
    DiT = (TP - TS) / (TP + TS)
    ZiT = (1. - DiT ** 2) ** 0.5
    TiR = 0.5 * (RP + RS)
    DiR = (RP - RS) / (RP + RS)
    ZiR = (1. - DiR ** 2) ** 0.5
    CosT = np.cos(np.deg2rad(RetT))
    SinT = np.sin(np.deg2rad(RetT))
    CosR = np.cos(np.deg2rad(RetR))
    SinR = np.sin(np.deg2rad(RetR))

    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))

    # Aanalyzer As before the PBS Eq. D.5
    ATP1 = (1 + C2aT * DaT * DiT)
    ATP2 = Y * (DiT + C2aT * DaT)
    ATP3 = Y * S2aT * DaT * ZiT * CosT
    ATP4 = S2aT * DaT * ZiT * SinT
    ATP = np.array([ATP1, ATP2, ATP3, ATP4])

    ARP1 = (1 + C2aR * DaR * DiR)
    ARP2 = Y * (DiR + C2aR * DaR)
    ARP3 = Y * S2aR * DaR * ZiR * CosR
    ARP4 = S2aR * DaR * ZiR * SinR
    ARP = np.array([ARP1, ARP2, ARP3, ARP4])

    DTa = ATP2 * Y / ATP1
    DRa = ARP2 * Y / ARP1

    # ---- 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 paremeters 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 paremeters 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 paremeters 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 paremeters 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 paremeters 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 paremeters 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 paremeters 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 paremeters 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 paremeters 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 = TaT * TiT * TiO * TiE * (AT + BT)
    IoutTm = TaT * TiT * TiO * TiE * (AT - BT)
    IoutRp = TaR * TiR * TiO * TiE * (AR + BR)
    IoutRm = TaR * TiR * 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 = (TaR * TiR) / (TaT * TiT)  # Eta = Eta*/K  Eq. 84
    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

    # ----- Forward simulated signals and LDRsim with atrue; from input file
    It = TaT * TiT * TiO * TiE * (GT + atrue * HT)
    Ir = TaR * TiR * TiO * TiE * (GR + atrue * HR)
    # 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.:
        LDRsimx = 1. / LDRsim
    else:
        LDRsimx = LDRsim

    # The following is correct without doubt
    # LDRCorr = (LDRsim*K/Etax*(GT+HT)-(GR+HR))/((GR-HR)-LDRsim*K/Etax*(GT-HT))

    # 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 * K / Etax * (GT - HT))

    TTa = TiT * TaT  # *ATP1
    TRa = TiR * TaR  # *ARP1

    F11sim = 1 / (TiO * TiE) * (
    (HR * Etax / K * It / TTa - HT * Ir / TRa) / (HR * GT - HT * GR))  # IL = 1, Etat = Etar = 1

    return (GT, HT, GR, HR, K, Eta, LDRsimx, LDRCorr, DTa, DRa, TTa, TRa, F11sim)


# *******************************************************************************************************************************

# --- CALC truth
GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0,
                                                                                       RetO0, DiO0, RotC0, RetC0, DiC0,
                                                                                       TP0, TS0, RP0, RS0, ERaT0,
                                                                                       RotaT0, RetT0, ERaR0, RotaR0,
                                                                                       RetR0, LDRCal0)

# --- Print parameters to console and output file
with open('output_files\output_' + LID + '.dat', '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:8} {1:8} {2:8.5f}; {3:8} {4:7.4f}±{5:7.4f}/{6:2d}".format("Laser: ", "DOLP = ", bL,
                                                                            "        rotation alpha = ", RotL0, dRotL,
                                                                            nRotL))
        print("              Diatt.,             Tunpol,   Retard.,   Rotation (deg)")
        print("{0:12} {1:7.4f}±{2:7.4f}/{8:2d}, {3:7.4f}, {4:3.0f}±{5:3.0f}/{9:2d}, {6:7.4f}±{7:7.4f}/{10:2d}".format(
            "Emitter    ", DiE0, dDiE, TiE, RetE0, dRetE, RotE0, dRotE, nDiE, nRetE, nRotE))
        print("{0:12} {1:7.4f}±{2:7.4f}/{8:2d}, {3:7.4f}, {4:3.0f}±{5:3.0f}/{9:2d}, {6:7.4f}±{7:7.4f}/{10:2d}".format(
            "Receiver   ", DiO0, dDiO, TiO, RetO0, dRetO, RotO0, dRotO, nDiO, nRetO, nRotO))
        print("{0:12} {1:7.4f}±{2:7.4f}/{8:2d}, {3:7.4f}, {4:3.0f}±{5:3.0f}/{9:2d}, {6:7.4f}±{7:7.4f}/{10:2d}".format(
            "Calibrator ", DiC0, dDiC, TiC, RetC0, dRetC, RotC0, dRotC, nDiC, nRetC, nRotC))
        print("{0:12}".format(" --- Pol.-filter ---"))
        print(
            "{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format("ERT, RotT       :", ERaT0, dERaT, nERaT,
                                                                              RotaT0, dRotaT, nRotaT))
        print(
            "{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format("ERR, RotR       :", ERaR0, dERaR, nERaR,
                                                                              RotaR0, dRotaR, nRotaR))
        print("{0:12}".format(" --- PBS ---"))
        print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format("TP,TS           :", TP0, dTP, nTP, TS0,
                                                                                dTS, nTS))
        print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format("RP,RS           :", RP0, dRP, nRP, RS0,
                                                                                dRS, nRS))
        print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}, {5:1.0f}".format("DT,TT,DR,TR,Y   :", DiT, TiT, DiR, TiR, Y))
        print("{0:12}".format(" --- Combined PBS + Pol.-filter ---"))
        print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}".format("DT,TT,DR,TR     :", DTa0, TTa0, DRa0, TRa0))
        print()
        print("Rotation Error Epsilon For Normal Measurements = ", RotationErrorEpsilonForNormalMeasurements)
        print(Type[TypeC], Loc[LocC])
        print("Parallel signal detected in", dY[int(Y + 1)])
        print("RS_RP_depend_on_TS_TP = ", RS_RP_depend_on_TS_TP)
        #  end of print actual system parameters
        # ******************************************************************************

        # print()
        # print(" --- LDRCal during calibration | simulated and corrected LDRs -------------")
        # print("{0:8} |{1:8}->{2:8},{3:9}->{4:9} |{5:8}->{6:8}".format(" LDRCal"," LDRtrue", " LDRsim"," LDRtrue2", " LDRsim2", " LDRmeas", " LDRcorr"))
        # print("{0:8.5f} |{1:8.5f}->{2:8.5f},{3:9.5f}->{4:9.5f} |{5:8.5f}->{6:8.5f}".format(LDRCal, LDRtrue, LDRsim, LDRtrue2, LDRsim2, LDRmeas, LDRCorr))
        # print("{0:8}       |{1:8}->{2:8}->{3:8}".format(" LDRCal"," LDRtrue", " LDRsimx", " LDRcorr"))
        # print("{0:6.3f}±{1:5.3f}/{2:2d}|{3:8.5f}->{4:8.5f}->{5:8.5f}".format(LDRCal0, dLDRCal, nLDRCal, LDRtrue, LDRsimx, LDRCorr))
        # print("{0:8}       |{1:8}->{2:8}->{3:8}".format(" LDRCal"," LDRtrue", " LDRsimx", " LDRcorr"))
        # print(" --- LDRCal during calibration")
        print("{0:26}: {1:6.3f}±{2:5.3f}/{3:2d}".format("LDRCal during calibration in calibration range", LDRCal0,
                                                        dLDRCal, nLDRCal))

        # print("{0:8}={1:8.5f};{2:8}={3:8.5f}".format(" IinP",IinP," F11sim",F11sim))
        print()

        K0List = np.zeros(3)
        LDRsimxList = np.zeros(3)
        LDRCalList = 0.004, 0.2, 0.45
        for i, LDRCal in enumerate(LDRCalList):
            GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(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

        print('========================================================================')
        print("{0:8},{1:8},{2:8},{3:8},{4:9},{5:8},{6:9}".format(" GR", " GT", " HR", " HT", "  K(0.004)", " K(0.2)",
                                                                 "  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}".format(GR0, GT0, HR0, HT0, K0List[0],
                                                                                      K0List[1], K0List[2]))
        print('========================================================================')

        print("{0:9},{1:9},{2:9}".format("  LDRtrue", "  LDRsimx", "  LDRCorr"))
        LDRtrueList = 0.004, 0.02, 0.2, 0.45
        for i, LDRtrue in enumerate(LDRtrueList):
            GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0,
                                                                                                   DiE0, RotO0, RetO0,
                                                                                                   DiO0, RotC0, RetC0,
                                                                                                   DiC0, TP0, TS0, RP0,
                                                                                                   RS0, ERaT0, RotaT0,
                                                                                                   RetT0, ERaR0, RotaR0,
                                                                                                   RetR0, LDRCal0)
            print("{0:9.5f},{1:9.5f},{2:9.5f}".format(LDRtrue, LDRsimx, LDRCorr))

file = open('output_files\output_' + LID + '.dat', 'r')
print(file.read())
file.close()

'''
if(PrintToOutputFile):
    f = open('output_ver7.dat', 'w')
    old_target = sys.stdout
    sys.stdout = f

    print("something")

if(PrintToOutputFile):
    sys.stdout.flush()
    f.close
    sys.stdout = old_target
'''
if (Error_Calc):
    # --- CALC again truth with LDRCal0 to reset all 0-values
    GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0,
                                                                                           RotO0, RetO0, DiO0, RotC0,
                                                                                           RetC0, DiC0, TP0, TS0, RP0,
                                                                                           RS0, ERaT0, RotaT0, RetT0,
                                                                                           ERaR0, RotaR0, RetR0,
                                                                                           LDRCal0)

    # --- Start Errors calculation with variable parameters ------------------------------------------------------------------

    iN = -1
    N = ((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("N = ", N, " ", end="")

    if N > 1e6:
        if user_yes_no_query('Warning: processing ' + str(
            N) + ' samples will take very long. Do you want to proceed?') == 0: sys.exit()
    if N > 5e6:
        if user_yes_no_query('Warning: the memory required for ' + str(N) + ' samples might be ' + '{0:5.1f}'.format(
                    N / 4e6) + ' GB. Do you anyway want to proceed?') == 0: sys.exit()

    # if user_yes_no_query('Warning: processing' + str(N) + ' samples will take very long. Do you want to proceed?') == 0: sys.exit()

    # --- Arrays for plotting ------
    LDRmin = np.zeros(5)
    LDRmax = np.zeros(5)
    F11min = np.zeros(5)
    F11max = np.zeros(5)

    LDRrange = np.zeros(5)
    LDRrange = 0.004, 0.02, 0.1, 0.3, 0.45
    # aLDRsimx = np.zeros(N)
    # aLDRsimx2 = np.zeros(N)
    # aLDRcorr = np.zeros(N)
    # aLDRcorr2 = np.zeros(N)
    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)
    aA = np.zeros((5, N))
    aX = np.zeros((5, N))
    aF11corr = np.zeros((5, N))

    atime = clock()
    dtime = clock()

    # --- Calc Error signals
    # GT, HT, GR, HR, K, Eta, LDRsim = Calc(RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC, TP, TS)
    # ---- Do the calculations of bra-ket vectors
    h = -1. if TypeC == 2 else 1

    # from input file: measured LDRm and true LDRtrue, LDRtrue2  =>
    ameas = (1. - LDRmeas) / (1 + LDRmeas)
    atrue = (1. - LDRtrue) / (1 + LDRtrue)
    atrue2 = (1. - LDRtrue2) / (1 + LDRtrue2)

    for iLDRCal in range(-nLDRCal, nLDRCal + 1):
        # from input file:  assumed LDRCal for calibration measurements
        LDRCal = LDRCal0
        if nLDRCal > 0: LDRCal = LDRCal0 + iLDRCal * dLDRCal / nLDRCal

        GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0,
                                                                                               DiE0, RotO0, RetO0, DiO0,
                                                                                               RotC0, RetC0, DiC0, TP0,
                                                                                               TS0, RP0, RS0, ERaT0,
                                                                                               RotaT0, RetT0, ERaR0,
                                                                                               RotaR0, RetR0, LDRCal)
        aCal = (1. - LDRCal) / (1 + LDRCal)
        for iRotL, iRotE, iRetE, iDiE \
                in [(iRotL, iRotE, iRetE, iDiE)
                    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 nRotL > 0: RotL = RotL0 + iRotL * dRotL / nRotL
            if nRotE > 0: RotE = RotE0 + iRotE * dRotE / nRotE
            if nRetE > 0: RetE = RetE0 + iRetE * dRetE / nRetE
            if nDiE > 0:  DiE = DiE0 + iDiE * dDiE / nDiE

            # angles of emitter and laser and calibrator and receiver optics
            # RotL = alpha, RotE = beta, RotO = gamma, RotC = epsilon
            S2a = np.sin(2 * np.deg2rad(RotL))
            C2a = np.cos(2 * np.deg2rad(RotL))
            S2b = np.sin(2 * np.deg2rad(RotE))
            C2b = np.cos(2 * np.deg2rad(RotE))
            S2ab = np.sin(np.deg2rad(2 * RotL - 2 * RotE))
            C2ab = np.cos(np.deg2rad(2 * RotL - 2 * RotE))

            # Laser with Degree of linear polarization DOLP = bL
            IinL = 1.
            QinL = bL
            UinL = 0.
            VinL = (1. - bL ** 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)]:

                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

                if nRotO > 0: RotO = RotO0 + iRotO * dRotO / nRotO
                if nRetO > 0: RetO = RetO0 + iRetO * dRetO / nRetO
                if nDiO > 0:  DiO = DiO0 + iDiO * dDiO / nDiO
                if nRotC > 0: RotC = RotC0 + iRotC * dRotC / nRotC
                if nRetC > 0: RetC = RetC0 + iRetC * dRetC / nRetC
                if nDiC > 0:  DiC = DiC0 + iDiC * dDiC / nDiC
                if nTP > 0:   TP = TP0 + iTP * dTP / nTP
                if nTS > 0:   TS = TS0 + iTS * dTS / nTS
                if nRP > 0:   RP = RP0 + iRP * dRP / nRP
                if nRS > 0:   RS = RS0 + iRS * dRS / nRS
                if nERaT > 0: ERaT = ERaT0 + iERaT * dERaT / nERaT
                if nRotaT > 0: RotaT = RotaT0 + iRotaT * dRotaT / nRotaT
                if nRetT > 0: RetT = RetT0 + iRetT * dRetT / nRetT
                if nERaR > 0: ERaR = ERaR0 + iERaR * dERaR / nERaR
                if nRotaR > 0: RotaR = RotaR0 + iRotaR * dRotaR / nRotaR
                if nRetR > 0: RetR = RetR0 + iRetR * dRetR / nRetR

                # print("{0:5.2f}, {1:5.2f}, {2:5.2f}, {3:10d}".format(RotL, RotE, RotO, iN))

                # receiver optics
                CosO = np.cos(np.deg2rad(RetO))
                SinO = np.sin(np.deg2rad(RetO))
                ZiO = (1. - DiO ** 2) ** 0.5
                WiO = (1. - ZiO * CosO)
                S2g = np.sin(np.deg2rad(2 * RotO))
                C2g = np.cos(np.deg2rad(2 * RotO))
                # calibrator
                CosC = np.cos(np.deg2rad(RetC))
                SinC = np.sin(np.deg2rad(RetC))
                ZiC = (1. - DiC ** 2) ** 0.5
                WiC = (1. - ZiC * CosC)

                # analyser
                # For POLLY_XTs
                if (RS_RP_depend_on_TS_TP):
                    RS = 1 - TS
                    RP = 1 - TP
                TiT = 0.5 * (TP + TS)
                DiT = (TP - TS) / (TP + TS)
                ZiT = (1. - DiT ** 2) ** 0.5
                TiR = 0.5 * (RP + RS)
                DiR = (RP - RS) / (RP + RS)
                ZiR = (1. - DiR ** 2) ** 0.5
                CosT = np.cos(np.deg2rad(RetT))
                SinT = np.sin(np.deg2rad(RetT))
                CosR = np.cos(np.deg2rad(RetR))
                SinR = np.sin(np.deg2rad(RetR))

                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))

                # Aanalyzer As before the PBS Eq. D.5
                ATP1 = (1 + C2aT * DaT * DiT)
                ATP2 = Y * (DiT + C2aT * DaT)
                ATP3 = Y * S2aT * DaT * ZiT * CosT
                ATP4 = S2aT * DaT * ZiT * SinT
                ATP = np.array([ATP1, ATP2, ATP3, ATP4])

                ARP1 = (1 + C2aR * DaR * DiR)
                ARP2 = Y * (DiR + C2aR * DaR)
                ARP3 = Y * S2aR * DaR * ZiR * CosR
                ARP4 = S2aR * DaR * ZiR * SinR
                ARP = np.array([ARP1, ARP2, ARP3, ARP4])

                TTa = TiT * TaT  # *ATP1
                TRa = TiR * TaR  # *ARP1

                # ---- 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 paremeters 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 paremeters 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 paremeters 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 paremeters 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 paremeters 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 paremeters 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 paremeters 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 paremeters 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 paremeters 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()

                # Calibration signals with aCal => Determination of the correction K of the real calibration factor
                IoutTp = TaT * TiT * TiO * TiE * (AT + BT)
                IoutTm = TaT * TiT * TiO * TiE * (AT - BT)
                IoutRp = TaR * TiR * TiO * TiE * (AR + BR)
                IoutRm = TaR * TiR * TiO * TiE * (AR - BR)
                # --- Results and Corrections; electronic etaR and etaT are assumed to be 1
                # Eta = TiR/TiT   # Eta = Eta*/K  Eq. 84
                Etapx = IoutRp / IoutTp
                Etamx = IoutRm / IoutTm
                Etax = (Etapx * Etamx) ** 0.5
                K = Etax / Eta0
                # 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

                # *************************************************************************
                iLDR = -1
                for LDRTrue in LDRrange:
                    iLDR = iLDR + 1
                    atrue = (1 - LDRTrue) / (1 + LDRTrue)
                    # ----- Forward simulated signals and LDRsim with atrue; from input file
                    It = TaT * TiT * TiO * TiE * (GT + atrue * HT)  # TaT*TiT*TiC*TiO*IinL*(GT+atrue*HT)
                    Ir = TaR * TiR * TiO * TiE * (GR + atrue * HR)  # TaR*TiR*TiC*TiO*IinL*(GR+atrue*HR)

                    # LDRsim = 1/Eta*Ir/It  # simulated LDR* with Y from input file
                    LDRsim = Ir / It  # simulated uncorrected LDR with Y from input file
                    '''
                    if Y == 1.:
                        LDRsimx = LDRsim
                        LDRsimx2 = LDRsim2
                    else:
                        LDRsimx = 1./LDRsim
                        LDRsimx2 = 1./LDRsim2
                    '''
                    # ----- Backward correction
                    # Corrected LDRCorr from forward simulated LDRsim (atrue) with assumed true G0,H0,K0
                    LDRCorr = (LDRsim * K0 / Etax * (GT0 + HT0) - (GR0 + HR0)) / (
                    (GR0 - HR0) - LDRsim * K0 / Etax * (GT0 - HT0))

                    # -- 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)

                    # 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
                    aA[iLDR, iN] = LDRCorr

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

                    aLDRCal[iN] = iLDRCal
                    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

    # --- END loop
    btime = clock()
    print("\r done in      ", "{0:5.0f}".format(btime - atime), "sec")  # , end="\r")

    # --- Plot -----------------------------------------------------------------
    if (sns_loaded):
        sns.set_style("whitegrid")
        sns.set_palette("bright", 6)

    '''
    fig2 = plt.figure()
    plt.plot(aA[2,:],'b.')
    plt.plot(aA[3,:],'r.')
    plt.plot(aA[4,:],'g.')
    #plt.plot(aA[6,:],'c.')
    plt.show
    '''


    # Plot LDR
    def PlotSubHist(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

            LDRmin[iLDR] = np.min(aA[iLDR, :])
            LDRmax[iLDR] = np.max(aA[iLDR, :])
            Rmin = LDRmin[iLDR] * 0.995  # np.min(aA[iLDR,:])    * 0.995
            Rmax = LDRmax[iLDR] * 1.005  # np.max(aA[iLDR,:])    * 1.005

            # plt.subplot(5,2,iLDR+1)
            plt.subplot(1, 5, iLDR + 1)
            (n, bins, patches) = plt.hist(aA[iLDR, :],
                                          bins=100, log=False,
                                          range=[Rmin, Rmax],
                                          alpha=0.5, normed=False, color='0.5', histtype='stepfilled')

            for iaX in range(-naX, naX + 1):
                plt.hist(aA[iLDR, aX == iaX],
                         range=[Rmin, Rmax],
                         bins=100, log=False, alpha=0.3, normed=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 + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar, fontsize=14, y=1.05)
        # plt.show()
        # fig.savefig(LID + '_' + aVar + '.png', dpi=150, bbox_inches='tight', pad_inches=0)
        # plt.close
        return


    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)

    plt.show()
    plt.close
    '''
    print()
    #print("IT(LDRtrue) devided by IT(LDRtrue = 0.004)")
    print(" ############################################################################## ")
    print(Text1)
    print()

    iLDR = 5
    for LDRTrue in LDRrange:
        iLDR = iLDR - 1
        aF11corr[iLDR,:] = aF11corr[iLDR,:] / aF11corr[0,:] - 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(aA[iLDR,:])    * 0.995
            #Rmax = F11max[iLDR] * 1.005 #  np.max(aA[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, normed=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, normed=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 (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)

    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(aA[iLDR,:])
        LDRmax[iLDR] = np.max(aA[iLDR,:])
        Rmin = np.min(aA[iLDR,:])    * 0.999
        Rmax = np.max(aA[iLDR,:])    * 1.001
        plt.subplot(5,2,iLDR+1)
        (n, bins, patches) = plt.hist(aA[iLDR,:],
                 range=[Rmin, Rmax],
                 bins=200, log=False, alpha=0.2, normed=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
    '''

    # --- Plot LDRmin, LDRmax
    fig2 = plt.figure()
    plt.plot(LDRrange, LDRmax - LDRrange, linewidth=2.0, color='b')
    plt.plot(LDRrange, LDRmin - LDRrange, 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\LDR_min_max_' + LID + '.dat', 'w') as f:
        with redirect_stdout(f):
            print(LID)
            print("LDRtrue, LDRmin, LDRmax")
            for i in range(len(LDRrange)):
                print("{0:7.4f},{1:7.4f},{2:7.4f}".format(LDRrange[i], LDRmin[i], LDRmax[i]))

    '''
    # --- Plot K over LDRCal
    fig3 = plt.figure()
    plt.plot(LDRCal0+aLDRCal*dLDRCal/nLDRCal,aX[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')
#
'''