Lidar polarisation GHK-parameter and error calculation: comparison lidar_correction

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

comparison: lidar_correction_ghk.py

lidar_correction_ghk.py

Mercurial > public > atmospheric_lidar_ghk / file comparison

comparison: lidar_correction_ghk.py

lidar_correction_ghk.py