Lidar polarisation GHK-parameter and error calculation: comparison lidar_correction

-:b66dbbe6c4fa
+:40d55af749b6
 Equation reference: http://www.atmos-meas-tech-discuss.net/amt-2015-338/amt-2015-338.pdf
 With equations code from Appendix C
 Python 3.4.2
 """
-#!/usr/bin/env python3
+# !/usr/bin/env python3
 # Comment:  The code works with Python 2.7  with the help of following line, which enables Python2 to correctly interpret the Python 3 print statements.
 from __future__ import print_function
-#import math
+import os
+import sys
 import numpy as np
-import sys
-import os
 # Comment: the seaborn library makes nicer plots, but the code works also without it.
 try:
 import seaborn as sns
 sns_loaded = True
 except ImportError:
 sns_loaded = False
 import matplotlib.pyplot as plt
 from time import clock
-#from matplotlib.backends.backend_pdf import PdfPages
+# from matplotlib.backends.backend_pdf import PdfPages
-#pdffile = '{}.pdf'.format('path')
+# pdffile = '{}.pdf'.format('path')
-#pp = PdfPages(pdffile)
+# pp = PdfPages(pdffile)
 ## pp.savefig can be called multiple times to save to multiple pages
-#pp.savefig()
+# pp.savefig()
-#pp.close()
+# pp.close()
 from contextlib import contextmanager
 @contextmanager
 def redirect_stdout(new_target):
-old_target, sys.stdout = sys.stdout, new_target # replace sys.stdout
+old_target, sys.stdout = sys.stdout, new_target  # replace sys.stdout
 try:
-yield new_target # run some code with the replaced stdout
+yield new_target  # run some code with the replaced stdout
 finally:
 sys.stdout.flush()
-sys.stdout = old_target # restore to the previous value
+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()
+# 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
+# PrintToOutputFile = True
-sqr05 = 0.5**0.5
+sqr05 = 0.5 ** 0.5
 # Do you want to calculate the errors? If not, just the GHK-parameters are determined.
 Error_Calc = True
 # ---- Initial definition of variables; the actual values will be read in with exec(open('./optic_input.py').read()) below
 LID = "internal"
 EID = "internal"
 # --- IL Laser IL and +-Uncertainty
-bL = 1.    #degree of linear polarization; default 1
+bL = 1.  # degree of linear polarization; default 1
-RotL, dRotL, nRotL     = 0.0, 0.0,     1    #alpha; rotation of laser polarization in degrees; default 0
+RotL, dRotL, nRotL = 0.0, 0.0, 1  # alpha; rotation of laser polarization in degrees; default 0
 # --- ME Emitter and +-Uncertainty
-DiE, dDiE, nDiE     = 0., 0.00,     1    # Diattenuation
+DiE, dDiE, nDiE = 0., 0.00, 1  # Diattenuation
-TiE         = 1.        # Unpolarized transmittance
+TiE = 1.  # Unpolarized transmittance
-RetE, dRetE, nRetE     = 0., 180.0,     0    # Retardance in degrees
+RetE, dRetE, nRetE = 0., 180.0, 0  # Retardance in degrees
-RotE, dRotE, nRotE     = 0., 0.0,     0    # beta: Rotation of optical element 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
+DiO, dDiO, nDiO = -0.055, 0.003, 1
-TiO                 = 0.9
+TiO = 0.9
-RetO, dRetO, nRetO     = 0., 180.0,     2
+RetO, dRetO, nRetO = 0., 180.0, 2
-RotO, dRotO, nRotO     = 0., 0.1,     1    #gamma
+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
+TP, dTP, nTP = 0.98, 0.02, 1
-TS,   dTS, nTS     = 0.001, 0.001,    1
+TS, dTS, nTS = 0.001, 0.001, 1
 TiT = 0.5 * (TP + TS)
-DiT = (TP-TS)/(TP+TS)
+DiT = (TP - TS) / (TP + TS)
 # PolFilter
-RetT, dRetT, nRetT  = 0.,     180.,    0
+RetT, dRetT, nRetT = 0., 180., 0
-ERaT, dERaT, nERaT  = 0.001, 0.001,    1
+ERaT, dERaT, nERaT = 0.001, 0.001, 1
-RotaT, dRotaT, nRotaT = 0.,     3., 1
+RotaT, dRotaT, nRotaT = 0., 3., 1
-DaT = (1-ERaT)/(1+ERaT)
+DaT = (1 - ERaT) / (1 + ERaT)
-TaT = 0.5*(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):
+if (RS_RP_depend_on_TS_TP):
-RP, dRP, nRP        = 1-TP,  0.00, 0
+RP, dRP, nRP = 1 - TP, 0.00, 0
-RS, dRS, nRS        = 1-TS,  0.00, 0
+RS, dRS, nRS = 1 - TS, 0.00, 0
 else:
-RP, dRP, nRP        = 0.05,  0.01, 1
+RP, dRP, nRP = 0.05, 0.01, 1
-RS, dRS, nRS        = 0.98,  0.01, 1
+RS, dRS, nRS = 0.98, 0.01, 1
 TiR = 0.5 * (RP + RS)
-DiR = (RP-RS)/(RP+RS)
+DiR = (RP - RS) / (RP + RS)
 # PolFilter
-RetR, dRetR, nRetR  = 0.,     180.,     0
+RetR, dRetR, nRetR = 0., 180., 0
-ERaR, dERaR, nERaR  = 0.001, 0.001, 1
+ERaR, dERaR, nERaR = 0.001, 0.001, 1
-RotaR,dRotaR,nRotaR = 90.,     3.,        1
+RotaR, dRotaR, nRotaR = 90., 3., 1
-DaR = (1-ERaR)/(1+ERaR)
+DaR = (1 - ERaR) / (1 + ERaR)
-TaR = 0.5*(1+ERaR)
+TaR = 0.5 * (1 + ERaR)
 # Parellel signal detected in the transmitted channel => Y = 1, or in the reflected channel => Y = -1
 Y = -1.
 # Calibrator =  type defined by matrix values
-LocC = 4     # location of calibrator: behind laser = 1; behind emitter = 2; before receiver = 3; before PBS = 4
+LocC = 4  # location of calibrator: behind laser = 1; behind emitter = 2; before receiver = 3; before PBS = 4
-TypeC = 3   # linear polarizer calibrator
+TypeC = 3  # linear polarizer calibrator
 # example with extinction ratio 0.001
-DiC, dDiC, nDiC     = 1.0,     0.,     0    # ideal 1.0
+DiC, dDiC, nDiC = 1.0, 0., 0  # ideal 1.0
-TiC = 0.5    # ideal 0.5
+TiC = 0.5  # ideal 0.5
-RetC, dRetC, nRetC     = 0.,     0.,     0
+RetC, dRetC, nRetC = 0., 0., 0
-RotC, dRotC, nRotC     = 0.0,     0.1,     0    #constant calibrator offset epsilon
+RotC, dRotC, nRotC = 0.0, 0.1, 0  # constant calibrator offset epsilon
-RotationErrorEpsilonForNormalMeasurements = False    #     is in general False for TypeC == 3 calibrator
+RotationErrorEpsilonForNormalMeasurements = False  # is in general False for TypeC == 3 calibrator
 # Rotation error without calibrator: if False, then epsilon = 0 for normal measurements
 RotationErrorEpsilonForNormalMeasurements = True
 # LDRCal assumed atmospheric linear depolarization ratio during the calibration measurements (first guess)
-LDRCal0,dLDRCal,nLDRCal= 0.25, 0.04, 1
+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
 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°']
+Type = ['', 'mechanical rotator', 'hwp rotator', 'linear polarizer', 'qwp rotator', 'circular polarizer',
+'real HWP +-22.5°']
 dY = ['reflected channel', '', 'transmitted channel']
 #  end of initial definition of variables
 # *******************************************************************************************************************************
 # --- Read actual lidar system parameters from ./optic_input.py  (must be in the same directory)
 InputFile = 'optic_input_example_lidar.py'
-#InputFile = 'optic_input_ver6e_POLIS_355.py'
+# InputFile = 'optic_input_ver6e_POLIS_355.py'
-#InputFile = 'optic_input_ver6e_POLIS_355_JA.py'
+# InputFile = 'optic_input_ver6e_POLIS_355_JA.py'
-#InputFile = 'optic_input_ver6c_POLIS_532.py'
+# InputFile = 'optic_input_ver6c_POLIS_532.py'
-#InputFile = 'optic_input_ver6e_POLIS_532.py'
+# InputFile = 'optic_input_ver6e_POLIS_532.py'
-#InputFile = 'optic_input_ver8c_POLIS_532.py'
+# InputFile = 'optic_input_ver8c_POLIS_532.py'
-#InputFile = 'optic_input_ver6e_MUSA.py'
+# InputFile = 'optic_input_ver6e_MUSA.py'
-#InputFile = 'optic_input_ver6e_MUSA_JA.py'
+# InputFile = 'optic_input_ver6e_MUSA_JA.py'
-#InputFile = 'optic_input_ver6e_PollyXTSea.py'
+# InputFile = 'optic_input_ver6e_PollyXTSea.py'
-#InputFile = 'optic_input_ver6e_PollyXTSea_JA.py'
+# InputFile = 'optic_input_ver6e_PollyXTSea_JA.py'
-#InputFile = 'optic_input_ver6e_PollyXT_RALPH.py'
+# InputFile = 'optic_input_ver6e_PollyXT_RALPH.py'
-#InputFile = 'optic_input_ver8c_PollyXT_RALPH.py'
+# InputFile = 'optic_input_ver8c_PollyXT_RALPH.py'
-#InputFile = 'optic_input_ver8c_PollyXT_RALPH_2.py'
+# InputFile = 'optic_input_ver8c_PollyXT_RALPH_2.py'
-#InputFile = 'optic_input_ver8c_PollyXT_RALPH_3.py'
+# InputFile = 'optic_input_ver8c_PollyXT_RALPH_3.py'
-#InputFile = 'optic_input_ver8c_PollyXT_RALPH_4.py'
+# InputFile = 'optic_input_ver8c_PollyXT_RALPH_4.py'
-#InputFile = 'optic_input_ver8c_PollyXT_RALPH_5.py'
+# InputFile = 'optic_input_ver8c_PollyXT_RALPH_5.py'
-#InputFile = 'optic_input_ver8c_PollyXT_RALPH_6.py'
+# InputFile = 'optic_input_ver8c_PollyXT_RALPH_6.py'
-#InputFile = 'optic_input_ver8c_PollyXT_RALPH_7.py'
+# InputFile = 'optic_input_ver8c_PollyXT_RALPH_7.py'
-#InputFile = 'optic_input_ver8a_MOHP_DPL_355.py'
+# InputFile = 'optic_input_ver8a_MOHP_DPL_355.py'
-#InputFile = 'optic_input_ver9_MOHP_DPL_355.py'
+# InputFile = 'optic_input_ver9_MOHP_DPL_355.py'
-#InputFile = 'optic_input_ver6e_RALI.py'
+# InputFile = 'optic_input_ver6e_RALI.py'
-#InputFile = 'optic_input_ver6e_RALI_JA.py'
+# InputFile = 'optic_input_ver6e_RALI_JA.py'
-#InputFile = 'optic_input_ver6e_RALI_new.py'
+# InputFile = 'optic_input_ver6e_RALI_new.py'
-#InputFile = 'optic_input_ver6e_RALI_act.py'
+# InputFile = 'optic_input_ver6e_RALI_act.py'
-#InputFile = 'optic_input_ver6e_MULHACEN.py'
+# InputFile = 'optic_input_ver6e_MULHACEN.py'
-#InputFile = 'optic_input_ver6e_MULHACEN_JA.py'
+# InputFile = 'optic_input_ver6e_MULHACEN_JA.py'
-#InputFile = 'optic_input_ver6e-IPRAL.py'
+# InputFile = 'optic_input_ver6e-IPRAL.py'
-#InputFile = 'optic_input_ver6e-IPRAL_JA.py'
+# InputFile = 'optic_input_ver6e-IPRAL_JA.py'
-#InputFile = 'optic_input_ver6e-LB21.py'
+# InputFile = 'optic_input_ver6e-LB21.py'
-#InputFile = 'optic_input_ver6e-LB21_JA.py'
+# InputFile = 'optic_input_ver6e-LB21_JA.py'
-#InputFile = 'optic_input_ver6e_Bertha_b_355.py'
+# InputFile = 'optic_input_ver6e_Bertha_b_355.py'
-#InputFile = 'optic_input_ver6e_Bertha_b_532.py'
+# InputFile = 'optic_input_ver6e_Bertha_b_532.py'
-#InputFile = 'optic_input_ver6e_Bertha_b_1064.py'
+# InputFile = 'optic_input_ver6e_Bertha_b_1064.py'
 '''
 print("From ", dname)
 print("Running ", fname)
 print("Reading input file ", InputFile, " for")
 '''
 input_path = os.path.join('.', 'system_settings', InputFile)
 # this works with Python 2 - and 3?
-exec(open(input_path).read(), globals())
+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.
+# DiO = 0.
-#LDRtrue = 0.45
+# LDRtrue = 0.45
-#LDRtrue2 = 0.004
+# LDRtrue2 = 0.004
-#Y = -1
+# Y = -1
-#LocC = 4 #location of calibrator: 1 = behind laser; 2 = behind emitter; 3 = before receiver; 4 = before PBS
+# 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
+# RotationErrorEpsilonForNormalMeasurements = True
-#LDRCal = 0.25
+# LDRCal = 0.25
-#bL = 0.8
+# bL = 0.8
 ## --- Errors
 RotL0, dRotL, nRotL = RotL, dRotL, nRotL
-DiE0,  dDiE,  nDiE  = DiE,  dDiE,  nDiE
+DiE0, dDiE, nDiE = DiE, dDiE, nDiE
 RetE0, dRetE, nRetE = RetE, dRetE, nRetE
 RotE0, dRotE, nRotE = RotE, dRotE, nRotE
-DiO0,  dDiO,  nDiO  = DiO,  dDiO,  nDiO
+DiO0, dDiO, nDiO = DiO, dDiO, nDiO
 RetO0, dRetO, nRetO = RetO, dRetO, nRetO
 RotO0, dRotO, nRotO = RotO, dRotO, nRotO
-DiC0,  dDiC,  nDiC  = DiC,  dDiC,  nDiC
+DiC0, dDiC, nDiC = DiC, dDiC, nDiC
 RetC0, dRetC, nRetC = RetC, dRetC, nRetC
 RotC0, dRotC, nRotC = RotC, dRotC, nRotC
-TP0,   dTP,   nTP   = TP,   dTP,     nTP
+TP0, dTP, nTP = TP, dTP, nTP
-TS0,   dTS,   nTS   = TS,   dTS,     nTS
+TS0, dTS, nTS = TS, dTS, nTS
 RetT0, dRetT, nRetT = RetT, dRetT, nRetT
 ERaT0, dERaT, nERaT = ERaT, dERaT, nERaT
-RotaT0,dRotaT,nRotaT= RotaT,dRotaT,nRotaT
+RotaT0, dRotaT, nRotaT = RotaT, dRotaT, nRotaT
-RP0,   dRP,   nRP   = RP,   dRP,   nRP
+RP0, dRP, nRP = RP, dRP, nRP
-RS0,   dRS,   nRS   = RS,   dRS,   nRS
+RS0, dRS, nRS = RS, dRS, nRS
 RetR0, dRetR, nRetR = RetR, dRetR, nRetR
 ERaR0, dERaR, nERaR = ERaR, dERaR, nERaR
-RotaR0,dRotaR,nRotaR= RotaR,dRotaR,nRotaR
+RotaR0, dRotaR, nRotaR = RotaR, dRotaR, nRotaR
-LDRCal0,dLDRCal,nLDRCal=LDRCal,dLDRCal,nLDRCal
+LDRCal0, dLDRCal, nLDRCal = LDRCal, dLDRCal, nLDRCal
-#LDRCal0,dLDRCal,nLDRCal=LDRCal,dLDRCal,0
+# LDRCal0,dLDRCal,nLDRCal=LDRCal,dLDRCal,0
 # ---------- End of manual parameter change
 RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC = RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0
-TP, TS, RP, RS, ERaT, RotaT, RetT, ERaR, RotaR, RetR = TP0, TS0, RP0, RS0 , ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0
+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
+DTa0, TTa0, DRa0, TRa0, LDRsimx, LDRCorr = 0, 0, 0, 0, 0, 0
 TiT = 0.5 * (TP + TS)
-DiT = (TP-TS)/(TP+TS)
+DiT = (TP - TS) / (TP + TS)
-ZiT = (1. - DiT**2)**0.5
+ZiT = (1. - DiT ** 2) ** 0.5
 TiR = 0.5 * (RP + RS)
-DiR = (RP-RS)/(RP+RS)
+DiR = (RP - RS) / (RP + RS)
-ZiR = (1. - DiR**2)**0.5
+ZiR = (1. - DiR ** 2) ** 0.5
 # --- this subroutine is for the calculation with certain fixed parameters -----------------------------------------------------
-def Calc(RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC, TP, TS, RP, RS, ERaT, RotaT, RetT, ERaR, RotaR, RetR, LDRCal):
+def Calc(RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC, TP, TS, RP, RS, ERaT, RotaT, RetT, ERaR, RotaR, RetR,
+LDRCal):
 # ---- Do the calculations of bra-ket vectors
 h = -1. if TypeC == 2 else 1
 # from input file:  assumed LDRCal for calibration measurements
-aCal = (1.-LDRCal)/(1+LDRCal)
+aCal = (1. - LDRCal) / (1 + LDRCal)
 # from input file: measured LDRm and true LDRtrue, LDRtrue2  =>
-#ameas = (1.-LDRmeas)/(1+LDRmeas)
+# ameas = (1.-LDRmeas)/(1+LDRmeas)
-atrue = (1.-LDRtrue)/(1+LDRtrue)
+atrue = (1. - LDRtrue) / (1 + LDRtrue)
-#atrue2 = (1.-LDRtrue2)/(1+LDRtrue2)
+# 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))
+S2a = np.sin(2 * np.deg2rad(RotL))
-C2a = np.cos(2*np.deg2rad(RotL))
+C2a = np.cos(2 * np.deg2rad(RotL))
-S2b = np.sin(2*np.deg2rad(RotE))
+S2b = np.sin(2 * np.deg2rad(RotE))
-C2b = np.cos(2*np.deg2rad(RotE))
+C2b = np.cos(2 * np.deg2rad(RotE))
-S2ab = np.sin(np.deg2rad(2*RotL-2*RotE))
+S2ab = np.sin(np.deg2rad(2 * RotL - 2 * RotE))
-C2ab = np.cos(np.deg2rad(2*RotL-2*RotE))
+C2ab = np.cos(np.deg2rad(2 * RotL - 2 * RotE))
-S2g = np.sin(np.deg2rad(2*RotO))
+S2g = np.sin(np.deg2rad(2 * RotO))
-C2g = np.cos(np.deg2rad(2*RotO))
+C2g = np.cos(np.deg2rad(2 * RotO))
 # Laser with Degree of linear polarization DOLP = bL
 IinL = 1.
 QinL = bL
 UinL = 0.
-VinL = (1. - bL**2)**0.5
+VinL = (1. - bL ** 2) ** 0.5
 # Stokes Input Vector rotation Eq. E.4
-A = C2a*QinL - S2a*UinL
+A = C2a * QinL - S2a * UinL
-B = S2a*QinL + C2a*UinL
+B = S2a * QinL + C2a * UinL
 # Stokes Input Vector rotation Eq. E.9
-C = C2ab*QinL - S2ab*UinL
+C = C2ab * QinL - S2ab * UinL
-D = S2ab*QinL + C2ab*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
+ZiE = (1. - DiE ** 2) ** 0.5
-WiE = (1. - ZiE*CosE)
+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)
+IinE = (IinL + DiE * C)
-QinE = (C2b*DiE*IinL + A + S2b*(WiE*D - ZiE*SinE*VinL))
+QinE = (C2b * DiE * IinL + A + S2b * (WiE * D - ZiE * SinE * VinL))
-UinE = (S2b*DiE*IinL + B - C2b*(WiE*D - ZiE*SinE*VinL))
+UinE = (S2b * DiE * IinL + B - C2b * (WiE * D - ZiE * SinE * VinL))
-VinE = (-ZiE*SinE*D + ZiE*CosE*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
+QinF = aCal * QinE
-UinF = -aCal*UinE
+UinF = -aCal * UinE
-VinF = (1.-2.*aCal)*VinE
+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
+ZiO = (1. - DiO ** 2) ** 0.5
-WiO = (1. - ZiO*CosO)
+WiO = (1. - ZiO * CosO)
 # calibrator
 CosC = np.cos(np.deg2rad(RetC))
 SinC = np.sin(np.deg2rad(RetC))
-ZiC = (1. - DiC**2)**0.5
+ZiC = (1. - DiC ** 2) ** 0.5
-WiC = (1. - ZiC*CosC)
+WiC = (1. - ZiC * CosC)
 # Stokes Input Vector before the polarising beam splitter Eq. E.31
-A = C2g*QinE - S2g*UinE
+A = C2g * QinE - S2g * UinE
-B = S2g*QinE + C2g*UinE
+B = S2g * QinE + C2g * UinE
-IinP = (IinE + DiO*aCal*A)
+IinP = (IinE + DiO * aCal * A)
-QinP = (C2g*DiO*IinE + aCal*QinE - S2g*(WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE))
+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))
+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)
+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
+# F11sim = TiO*(IinE + DiO*atrue*A)/IinE
-#-------------------------
+# -------------------------
 # analyser
-#For POLLY_XTs
+# For POLLY_XTs
-if(RS_RP_depend_on_TS_TP):
+if (RS_RP_depend_on_TS_TP):
 RS = 1 - TS
 RP = 1 - TP
 TiT = 0.5 * (TP + TS)
-DiT = (TP-TS)/(TP+TS)
+DiT = (TP - TS) / (TP + TS)
-ZiT = (1. - DiT**2)**0.5
+ZiT = (1. - DiT ** 2) ** 0.5
 TiR = 0.5 * (RP + RS)
-DiR = (RP-RS)/(RP+RS)
+DiR = (RP - RS) / (RP + RS)
-ZiR = (1. - DiR**2)**0.5
+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)
+DaT = (1 - ERaT) / (1 + ERaT)
-DaR = (1-ERaR)/(1+ERaR)
+DaR = (1 - ERaR) / (1 + ERaR)
-TaT = 0.5*(1+ERaT)
+TaT = 0.5 * (1 + ERaT)
-TaR = 0.5*(1+ERaR)
+TaR = 0.5 * (1 + ERaR)
-S2aT = np.sin(np.deg2rad(h*2*RotaT))
+S2aT = np.sin(np.deg2rad(h * 2 * RotaT))
-C2aT = np.cos(np.deg2rad(2*RotaT))
+C2aT = np.cos(np.deg2rad(2 * RotaT))
-S2aR = np.sin(np.deg2rad(h*2*RotaR))
+S2aR = np.sin(np.deg2rad(h * 2 * RotaR))
-C2aR = np.cos(np.deg2rad(2*RotaR))
+C2aR = np.cos(np.deg2rad(2 * RotaR))
 # Aanalyzer As before the PBS Eq. D.5
-ATP1 = (1+C2aT*DaT*DiT)
+ATP1 = (1 + C2aT * DaT * DiT)
-ATP2 = Y*(DiT+C2aT*DaT)
+ATP2 = Y * (DiT + C2aT * DaT)
-ATP3 = Y*S2aT*DaT*ZiT*CosT
+ATP3 = Y * S2aT * DaT * ZiT * CosT
-ATP4 = S2aT*DaT*ZiT*SinT
+ATP4 = S2aT * DaT * ZiT * SinT
-ATP = np.array([ATP1,ATP2,ATP3,ATP4])
+ATP = np.array([ATP1, ATP2, ATP3, ATP4])
-ARP1 = (1+C2aR*DaR*DiR)
+ARP1 = (1 + C2aR * DaR * DiR)
-ARP2 = Y*(DiR+C2aR*DaR)
+ARP2 = Y * (DiR + C2aR * DaR)
-ARP3 = Y*S2aR*DaR*ZiR*CosR
+ARP3 = Y * S2aR * DaR * ZiR * CosR
-ARP4 = S2aR*DaR*ZiR*SinR
+ARP4 = S2aR * DaR * ZiR * SinR
-ARP = np.array([ARP1,ARP2,ARP3,ARP4])
+ARP = np.array([ARP1, ARP2, ARP3, ARP4])
-DTa = ATP2*Y/ATP1
+DTa = ATP2 * Y / ATP1
-DRa = ARP2*Y/ARP1
+DRa = ARP2 * Y / ARP1
 # ---- Calculate signals and correction parameters for diffeent locations and calibrators
 if LocC == 4:  # Calibrator before the PBS
-#print("Calibrator location not implemented yet")
+# print("Calibrator location not implemented yet")
-#S2ge = np.sin(np.deg2rad(2*RotO + h*2*RotC))
+# S2ge = np.sin(np.deg2rad(2*RotO + h*2*RotC))
-#C2ge = np.cos(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))
+S2e = np.sin(np.deg2rad(h * 2 * RotC))
-C2e = np.cos(np.deg2rad(2*RotC))
+C2e = np.cos(np.deg2rad(2 * RotC))
 # rotated AinP by epsilon Eq. C.3
-ATP2e = C2e*ATP2 + S2e*ATP3
+ATP2e = C2e * ATP2 + S2e * ATP3
-ATP3e = C2e*ATP3 - S2e*ATP2
+ATP3e = C2e * ATP3 - S2e * ATP2
-ARP2e = C2e*ARP2 + S2e*ARP3
+ARP2e = C2e * ARP2 + S2e * ARP3
-ARP3e = C2e*ARP3 - S2e*ARP2
+ARP3e = C2e * ARP3 - S2e * ARP2
-ATPe = np.array([ATP1,ATP2e,ATP3e,ATP4])
+ATPe = np.array([ATP1, ATP2e, ATP3e, ATP4])
-ARPe = np.array([ARP1,ARP2e,ARP3e,ARP4])
+ARPe = np.array([ARP1, ARP2e, ARP3e, ARP4])
 # Stokes Input Vector before the polarising beam splitter Eq. E.31
-A = C2g*QinE - S2g*UinE
+A = C2g * QinE - S2g * UinE
-B = S2g*QinE + C2g*UinE
+B = S2g * QinE + C2g * UinE
-#C = (WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE)
+# C = (WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE)
-Co = ZiO*SinO*VinE
+Co = ZiO * SinO * VinE
-Ca = (WiO*B - 2*ZiO*SinO*VinE)
+Ca = (WiO * B - 2 * ZiO * SinO * VinE)
-#C = Co + aCal*Ca
+# C = Co + aCal*Ca
-#IinP = (IinE + DiO*aCal*A)
+# IinP = (IinE + DiO*aCal*A)
-#QinP = (C2g*DiO*IinE + aCal*QinE - S2g*C)
+# QinP = (C2g*DiO*IinE + aCal*QinE - S2g*C)
-#UinP = (S2g*DiO*IinE - aCal*UinE + C2g*C)
+# UinP = (S2g*DiO*IinE - aCal*UinE + C2g*C)
-#VinP = (ZiO*SinO*aCal*B + ZiO*CosO*(1-2*aCal)*VinE)
+# VinP = (ZiO*SinO*aCal*B + ZiO*CosO*(1-2*aCal)*VinE)
 IinPo = IinE
-QinPo = (C2g*DiO*IinE - S2g*Co)
+QinPo = (C2g * DiO * IinE - S2g * Co)
-UinPo = (S2g*DiO*IinE + C2g*Co)
+UinPo = (S2g * DiO * IinE + C2g * Co)
-VinPo = ZiO*CosO*VinE
+VinPo = ZiO * CosO * VinE
-IinPa = DiO*A
+IinPa = DiO * A
-QinPa = QinE - S2g*Ca
+QinPa = QinE - S2g * Ca
-UinPa = -UinE + C2g*Ca
+UinPa = -UinE + C2g * Ca
-VinPa = ZiO*(SinO*B - 2*CosO*VinE)
+VinPa = ZiO * (SinO * B - 2 * CosO * VinE)
-IinP = IinPo + aCal*IinPa
+IinP = IinPo + aCal * IinPa
-QinP = QinPo + aCal*QinPa
+QinP = QinPo + aCal * QinPa
-UinP = UinPo + aCal*UinPa
+UinP = UinPo + aCal * UinPa
-VinP = VinPo + aCal*VinPa
+VinP = VinPo + aCal * VinPa
 # Stokes Input Vector before the polarising beam splitter rotated by epsilon Eq. C.3
-#QinPe = C2e*QinP + S2e*UinP
+# QinPe = C2e*QinP + S2e*UinP
-#UinPe = C2e*UinP - S2e*QinP
+# UinPe = C2e*UinP - S2e*QinP
-QinPoe = C2e*QinPo + S2e*UinPo
+QinPoe = C2e * QinPo + S2e * UinPo
-UinPoe = C2e*UinPo - S2e*QinPo
+UinPoe = C2e * UinPo - S2e * QinPo
-QinPae = C2e*QinPa + S2e*UinPa
+QinPae = C2e * QinPa + S2e * UinPa
-UinPae = C2e*UinPa - S2e*QinPa
+UinPae = C2e * UinPa - S2e * QinPa
-QinPe = C2e*QinP + S2e*UinP
+QinPe = C2e * QinP + S2e * UinP
-UinPe = C2e*UinP - S2e*QinP
+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
+AT = ATP1 * IinP + h * ATP4 * VinP
-BT = ATP3e*QinP - h*ATP2e*UinP
+BT = ATP3e * QinP - h * ATP2e * UinP
-AR = ARP1*IinP + h*ARP4*VinP
+AR = ARP1 * IinP + h * ARP4 * VinP
-BR = ARP3e*QinP - h*ARP2e*UinP
+BR = ARP3e * QinP - h * ARP2e * UinP
 # Correction paremeters for normal measurements; they are independent of LDR
-if (not RotationErrorEpsilonForNormalMeasurements):   # calibrator taken out
+if (not RotationErrorEpsilonForNormalMeasurements):  # calibrator taken out
-IS1 = np.array([IinPo,QinPo,UinPo,VinPo])
+IS1 = np.array([IinPo, QinPo, UinPo, VinPo])
-IS2 = np.array([IinPa,QinPa,UinPa,VinPa])
+IS2 = np.array([IinPa, QinPa, UinPa, VinPa])
-GT = np.dot(ATP,IS1)
+GT = np.dot(ATP, IS1)
-GR = np.dot(ARP,IS1)
+GR = np.dot(ARP, IS1)
-HT = np.dot(ATP,IS2)
+HT = np.dot(ATP, IS2)
-HR = np.dot(ARP,IS2)
+HR = np.dot(ARP, IS2)
 else:
-IS1 = np.array([IinPo,QinPo,UinPo,VinPo])
+IS1 = np.array([IinPo, QinPo, UinPo, VinPo])
-IS2 = np.array([IinPa,QinPa,UinPa,VinPa])
+IS2 = np.array([IinPa, QinPa, UinPa, VinPa])
-GT = np.dot(ATPe,IS1)
+GT = np.dot(ATPe, IS1)
-GR = np.dot(ARPe,IS1)
+GR = np.dot(ARPe, IS1)
-HT = np.dot(ATPe,IS2)
+HT = np.dot(ATPe, IS2)
-HR = np.dot(ARPe,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)
+AT = ATP1 * IinP + ATP3e * UinPe + ZiC * CosC * (ATP2e * QinPe + ATP4 * VinP)
-BT = DiC*(ATP1*UinPe + ATP3e*IinP) - ZiC*SinC*(ATP2e*VinP - ATP4*QinPe)
+BT = DiC * (ATP1 * UinPe + ATP3e * IinP) - ZiC * SinC * (ATP2e * VinP - ATP4 * QinPe)
-AR = ARP1*IinP + ARP3e*UinPe + ZiC*CosC*(ARP2e*QinPe + ARP4*VinP)
+AR = ARP1 * IinP + ARP3e * UinPe + ZiC * CosC * (ARP2e * QinPe + ARP4 * VinP)
-BR = DiC*(ARP1*UinPe + ARP3e*IinP) - ZiC*SinC*(ARP2e*VinP - ARP4*QinPe)
+BR = DiC * (ARP1 * UinPe + ARP3e * IinP) - ZiC * SinC * (ARP2e * VinP - ARP4 * QinPe)
 # Correction paremeters for normal measurements; they are independent of LDR
-if (not RotationErrorEpsilonForNormalMeasurements):   # calibrator taken out
+if (not RotationErrorEpsilonForNormalMeasurements):  # calibrator taken out
-IS1 = np.array([IinPo,QinPo,UinPo,VinPo])
+IS1 = np.array([IinPo, QinPo, UinPo, VinPo])
-IS2 = np.array([IinPa,QinPa,UinPa,VinPa])
+IS2 = np.array([IinPa, QinPa, UinPa, VinPa])
-GT = np.dot(ATP,IS1)
+GT = np.dot(ATP, IS1)
-GR = np.dot(ARP,IS1)
+GR = np.dot(ARP, IS1)
-HT = np.dot(ATP,IS2)
+HT = np.dot(ATP, IS2)
-HR = np.dot(ARP,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)])
+IS1e = np.array([IinPo + DiC * QinPoe, DiC * IinPo + QinPoe, ZiC * (CosC * UinPoe + SinC * VinPo),
-IS2e = np.array([IinPa+DiC*QinPae,DiC*IinPa+QinPae,ZiC*(CosC*UinPae+SinC*VinPa),-ZiC*(SinC*UinPae-CosC*VinPa)])
+-ZiC * (SinC * UinPoe - CosC * VinPo)])
-GT = np.dot(ATPe,IS1e)
+IS2e = np.array([IinPa + DiC * QinPae, DiC * IinPa + QinPae, ZiC * (CosC * UinPae + SinC * VinPa),
-GR = np.dot(ARPe,IS1e)
+-ZiC * (SinC * UinPae - CosC * VinPa)])
-HT = np.dot(ATPe,IS2e)
+GT = np.dot(ATPe, IS1e)
-HR = np.dot(ARPe,IS2e)
+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)
+AT = ATP1 * IinP + sqr05 * DiC * (ATP1 * QinPe + ATP2e * IinP) + (1 - 0.5 * WiC) * (
-BT = sqr05*DiC*(ATP1*UinPe + ATP3e*IinP) + 0.5*WiC*(ATP2e*UinPe + ATP3e*QinPe) - sqr05*ZiC*SinC*(ATP2e*VinP - ATP4*QinPe)
+ATP2e * QinPe + ATP3e * UinPe) + ZiC * (sqr05 * SinC * (ATP3e * VinP - ATP4 * UinPe) + ATP4 * CosC * VinP)
-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)
+BT = sqr05 * DiC * (ATP1 * UinPe + ATP3e * IinP) + 0.5 * WiC * (
-BR = sqr05*DiC*(ARP1*UinPe + ARP3e*IinP) + 0.5*WiC*(ARP2e*UinPe + ARP3e*QinPe) - sqr05*ZiC*SinC*(ARP2e*VinP - ARP4*QinPe)
+ATP2e * UinPe + ATP3e * QinPe) - sqr05 * ZiC * SinC * (ATP2e * VinP - ATP4 * QinPe)
+AR = ARP1 * IinP + sqr05 * DiC * (ARP1 * QinPe + ARP2e * IinP) + (1 - 0.5 * WiC) * (
+ARP2e * QinPe + ARP3e * UinPe) + ZiC * (sqr05 * SinC * (ARP3e * VinP - ARP4 * UinPe) + ARP4 * CosC * VinP)
+BR = sqr05 * DiC * (ARP1 * UinPe + ARP3e * IinP) + 0.5 * WiC * (
+ARP2e * UinPe + ARP3e * QinPe) - sqr05 * ZiC * SinC * (ARP2e * VinP - ARP4 * QinPe)
 # Correction paremeters for normal measurements; they are independent of LDR
-if (not RotationErrorEpsilonForNormalMeasurements):   # calibrator taken out
+if (not RotationErrorEpsilonForNormalMeasurements):  # calibrator taken out
-IS1 = np.array([IinPo,QinPo,UinPo,VinPo])
+IS1 = np.array([IinPo, QinPo, UinPo, VinPo])
-IS2 = np.array([IinPa,QinPa,UinPa,VinPa])
+IS2 = np.array([IinPa, QinPa, UinPa, VinPa])
-GT = np.dot(ATP,IS1)
+GT = np.dot(ATP, IS1)
-GR = np.dot(ARP,IS1)
+GR = np.dot(ARP, IS1)
-HT = np.dot(ATP,IS2)
+HT = np.dot(ATP, IS2)
-HR = np.dot(ARP,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)])
+IS1e = np.array([IinPo + DiC * QinPoe, DiC * IinPo + QinPoe, ZiC * (CosC * UinPoe + SinC * VinPo),
-IS2e = np.array([IinPa+DiC*QinPae,DiC*IinPa+QinPae,ZiC*(CosC*UinPae+SinC*VinPa),-ZiC*(SinC*UinPae-CosC*VinPa)])
+-ZiC * (SinC * UinPoe - CosC * VinPo)])
-GT = np.dot(ATPe,IS1e)
+IS2e = np.array([IinPa + DiC * QinPae, DiC * IinPa + QinPae, ZiC * (CosC * UinPae + SinC * VinPa),
-GR = np.dot(ARPe,IS1e)
+-ZiC * (SinC * UinPae - CosC * VinPa)])
-HT = np.dot(ATPe,IS2e)
+GT = np.dot(ATPe, IS1e)
-HR = np.dot(ARPe,IS2e)
+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))
+# S2ge = np.sin(np.deg2rad(2*RotO - 2*RotC))
-#C2ge = np.cos(np.deg2rad(2*RotO - 2*RotC))
+# C2ge = np.cos(np.deg2rad(2*RotO - 2*RotC))
-S2e = np.sin(np.deg2rad(2*RotC))
+S2e = np.sin(np.deg2rad(2 * RotC))
-C2e = np.cos(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])
+AF1 = np.array([1, C2g * DiO, S2g * DiO, 0])
-AF2 = np.array([C2g*DiO,1-S2g**2*WiO,S2g*C2g*WiO,-S2g*ZiO*SinO])
+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])
+AF3 = np.array([S2g * DiO, S2g * C2g * WiO, 1 - C2g ** 2 * WiO, C2g * ZiO * SinO])
-AF4 = np.array([0,S2g*SinO,-C2g*SinO,CosO])
+AF4 = np.array([0, S2g * SinO, -C2g * SinO, CosO])
-ATF = (ATP1*AF1+ATP2*AF2+ATP3*AF3+ATP4*AF4)
+ATF = (ATP1 * AF1 + ATP2 * AF2 + ATP3 * AF3 + ATP4 * AF4)
-ARF = (ARP1*AF1+ARP2*AF2+ARP3*AF3+ARP4*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]
+ATF2e = C2e * ATF[1] + S2e * ATF[2]
-ATF3e = C2e*ATF[2] - S2e*ATF[1]
+ATF3e = C2e * ATF[2] - S2e * ATF[1]
 ARF1 = ARF[0]
 ARF4 = ARF[3]
-ARF2e = C2e*ARF[1] + S2e*ARF[2]
+ARF2e = C2e * ARF[1] + S2e * ARF[2]
-ARF3e = C2e*ARF[2] - S2e*ARF[1]
+ARF3e = C2e * ARF[2] - S2e * ARF[1]
-ATFe = np.array([ATF1,ATF2e,ATF3e,ATF4])
+ATFe = np.array([ATF1, ATF2e, ATF3e, ATF4])
-ARFe = np.array([ARF1,ARF2e,ARF3e,ARF4])
+ARFe = np.array([ARF1, ARF2e, ARF3e, ARF4])
-QinEe = C2e*QinE + S2e*UinE
+QinEe = C2e * QinE + S2e * UinE
-UinEe = C2e*UinE - S2e*QinE
+UinEe = C2e * UinE - S2e * QinE
 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F)
 IinF = IinE
-QinF = aCal*QinE
+QinF = aCal * QinE
-UinF = -aCal*UinE
+UinF = -aCal * UinE
-VinF = (1.-2.*aCal)*VinE
+VinF = (1. - 2. * aCal) * VinE
 IinFo = IinE
 QinFo = 0.
 UinFo = 0.
 VinFo = VinE
 IinFa = 0.
 QinFa = QinE
 UinFa = -UinE
-VinFa = -2.*VinE
+VinFa = -2. * VinE
 # Stokes Input Vector before receiver optics rotated by epsilon Eq. C.3
-QinFe = C2e*QinF + S2e*UinF
+QinFe = C2e * QinF + S2e * UinF
-UinFe = C2e*UinF - S2e*QinF
+UinFe = C2e * UinF - S2e * QinF
-QinFoe = C2e*QinFo + S2e*UinFo
+QinFoe = C2e * QinFo + S2e * UinFo
-UinFoe = C2e*UinFo - S2e*QinFo
+UinFoe = C2e * UinFo - S2e * QinFo
-QinFae = C2e*QinFa + S2e*UinFa
+QinFae = C2e * QinFa + S2e * UinFa
-UinFae = C2e*UinFa - S2e*QinFa
+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
+if (TypeC == 2) or (TypeC == 1):  # rotator calibration Eq. C.4
 # parameters for calibration with aCal
-AT = ATF1*IinF + ATF4*h*VinF
+AT = ATF1 * IinF + ATF4 * h * VinF
-BT = ATF3e*QinF - ATF2e*h*UinF
+BT = ATF3e * QinF - ATF2e * h * UinF
-AR = ARF1*IinF + ARF4*h*VinF
+AR = ARF1 * IinF + ARF4 * h * VinF
-BR = ARF3e*QinF - ARF2e*h*UinF
+BR = ARF3e * QinF - ARF2e * h * UinF
 # Correction paremeters for normal measurements; they are independent of LDR
 if (not RotationErrorEpsilonForNormalMeasurements):
-GT = ATF1*IinE + ATF4*VinE
+GT = ATF1 * IinE + ATF4 * VinE
-GR = ARF1*IinE + ARF4*VinE
+GR = ARF1 * IinE + ARF4 * VinE
-HT = ATF2*QinE - ATF3*UinE - ATF4*2*VinE
+HT = ATF2 * QinE - ATF3 * UinE - ATF4 * 2 * VinE
-HR = ARF2*QinE - ARF3*UinE - ARF4*2*VinE
+HR = ARF2 * QinE - ARF3 * UinE - ARF4 * 2 * VinE
 else:
-GT = ATF1*IinE + ATF4*h*VinE
+GT = ATF1 * IinE + ATF4 * h * VinE
-GR = ARF1*IinE + ARF4*h*VinE
+GR = ARF1 * IinE + ARF4 * h * VinE
-HT = ATF2e*QinE - ATF3e*h*UinE - ATF4*h*2*VinE
+HT = ATF2e * QinE - ATF3e * h * UinE - ATF4 * h * 2 * VinE
-HR = ARF2e*QinE - ARF3e*h*UinE - ARF4*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])
+IF1e = np.array([IinF, ZiC * CosC * QinFe, UinFe, ZiC * CosC * VinF])
-IF2e = np.array([DiC*UinFe, -ZiC*SinC*VinF, DiC*IinF, ZiC*SinC*QinFe])
+IF2e = np.array([DiC * UinFe, -ZiC * SinC * VinF, DiC * IinF, ZiC * SinC * QinFe])
-AT = np.dot(ATFe,IF1e)
+AT = np.dot(ATFe, IF1e)
-AR = np.dot(ARFe,IF1e)
+AR = np.dot(ARFe, IF1e)
-BT = np.dot(ATFe,IF2e)
+BT = np.dot(ATFe, IF2e)
-BR = np.dot(ARFe,IF2e)
+BR = np.dot(ARFe, IF2e)
 # Correction paremeters for normal measurements; they are independent of LDR  --- the same as for TypeC = 6
-if (not RotationErrorEpsilonForNormalMeasurements):   # calibrator taken out
+if (not RotationErrorEpsilonForNormalMeasurements):  # calibrator taken out
-IS1 = np.array([IinE,0,0,VinE])
+IS1 = np.array([IinE, 0, 0, VinE])
-IS2 = np.array([0,QinE,-UinE,-2*VinE])
+IS2 = np.array([0, QinE, -UinE, -2 * VinE])
-GT = np.dot(ATF,IS1)
+GT = np.dot(ATF, IS1)
-GR = np.dot(ARF,IS1)
+GR = np.dot(ARF, IS1)
-HT = np.dot(ATF,IS2)
+HT = np.dot(ATF, IS2)
-HR = np.dot(ARF,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)])
+IS1e = np.array([IinFo + DiC * QinFoe, DiC * IinFo + QinFoe, ZiC * (CosC * UinFoe + SinC * VinFo),
-IS2e = np.array([IinFa+DiC*QinFae,DiC*IinFa+QinFae,ZiC*(CosC*UinFae+SinC*VinFa),-ZiC*(SinC*UinFae-CosC*VinFa)])
+-ZiC * (SinC * UinFoe - CosC * VinFo)])
-GT = np.dot(ATFe,IS1e)
+IS2e = np.array([IinFa + DiC * QinFae, DiC * IinFa + QinFae, ZiC * (CosC * UinFae + SinC * VinFa),
-GR = np.dot(ARFe,IS1e)
+-ZiC * (SinC * UinFae - CosC * VinFa)])
-HT = np.dot(ATFe,IS2e)
+GT = np.dot(ATFe, IS1e)
-HR = np.dot(ARFe,IS2e)
+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])
+IF1e = np.array([IinF + sqr05 * DiC * QinFe, sqr05 * DiC * IinF + (1 - 0.5 * WiC) * QinFe,
-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])
+(1 - 0.5 * WiC) * UinFe + sqr05 * ZiC * SinC * VinF,
-AT = np.dot(ATFe,IF1e)
+-sqr05 * ZiC * SinC * UinFe + ZiC * CosC * VinF])
-AR = np.dot(ARFe,IF1e)
+IF2e = np.array([sqr05 * DiC * UinFe, 0.5 * WiC * UinFe - sqr05 * ZiC * SinC * VinF,
-BT = np.dot(ATFe,IF2e)
+sqr05 * DiC * IinF + 0.5 * WiC * QinFe, sqr05 * ZiC * SinC * QinFe])
-BR = np.dot(ARFe,IF2e)
+AT = np.dot(ATFe, IF1e)
+AR = np.dot(ARFe, IF1e)
+BT = np.dot(ATFe, IF2e)
+BR = np.dot(ARFe, IF2e)
 # Correction paremeters for normal measurements; they are independent of LDR
-if (not RotationErrorEpsilonForNormalMeasurements):   # calibrator taken out
+if (not RotationErrorEpsilonForNormalMeasurements):  # calibrator taken out
-#IS1 = np.array([IinE,0,0,VinE])
+# IS1 = np.array([IinE,0,0,VinE])
-#IS2 = np.array([0,QinE,-UinE,-2*VinE])
+# IS2 = np.array([0,QinE,-UinE,-2*VinE])
-IS1 = np.array([IinFo,0,0,VinFo])
+IS1 = np.array([IinFo, 0, 0, VinFo])
-IS2 = np.array([0,QinFa,UinFa,VinFa])
+IS2 = np.array([0, QinFa, UinFa, VinFa])
-GT = np.dot(ATF,IS1)
+GT = np.dot(ATF, IS1)
-GR = np.dot(ARF,IS1)
+GR = np.dot(ARF, IS1)
-HT = np.dot(ATF,IS2)
+HT = np.dot(ATF, IS2)
-HR = np.dot(ARF,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)])
+IS1e = np.array([IinFo + DiC * QinFoe, DiC * IinFo + QinFoe, ZiC * (CosC * UinFoe + SinC * VinFo),
-IS2e = np.array([IinFa+DiC*QinFae,DiC*IinFa+QinFae,ZiC*(CosC*UinFae+SinC*VinFa),-ZiC*(SinC*UinFae-CosC*VinFa)])
+-ZiC * (SinC * UinFoe - CosC * VinFo)])
-#IS1e = np.array([IinFo,0,0,VinFo])
+IS2e = np.array([IinFa + DiC * QinFae, DiC * IinFa + QinFae, ZiC * (CosC * UinFae + SinC * VinFa),
-#IS2e = np.array([0,QinFae,UinFae,VinFa])
+-ZiC * (SinC * UinFae - CosC * VinFa)])
-GT = np.dot(ATFe,IS1e)
+# IS1e = np.array([IinFo,0,0,VinFo])
-GR = np.dot(ARFe,IS1e)
+# IS2e = np.array([0,QinFae,UinFae,VinFa])
-HT = np.dot(ATFe,IS2e)
+GT = np.dot(ATFe, IS1e)
-HR = np.dot(ARFe,IS2e)
+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")
+# print("Calibrator location not implemented yet")
-S2e = np.sin(np.deg2rad(2*RotC))
+S2e = np.sin(np.deg2rad(2 * RotC))
-C2e = np.cos(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])
+AF1 = np.array([1, C2g * DiO, S2g * DiO, 0])
-AF2 = np.array([C2g*DiO,1-S2g**2*WiO,S2g*C2g*WiO,-S2g*ZiO*SinO])
+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])
+AF3 = np.array([S2g * DiO, S2g * C2g * WiO, 1 - C2g ** 2 * WiO, C2g * ZiO * SinO])
-AF4 = np.array([0, S2g*SinO, -C2g*SinO, CosO])
+AF4 = np.array([0, S2g * SinO, -C2g * SinO, CosO])
-ATF = (ATP1*AF1+ATP2*AF2+ATP3*AF3+ATP4*AF4)
+ATF = (ATP1 * AF1 + ATP2 * AF2 + ATP3 * AF3 + ATP4 * AF4)
-ARF = (ARP1*AF1+ARP2*AF2+ARP3*AF3+ARP4*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]
 ATE1o, ARE1o = ATF1, ARF1
 ATE2o, ARE2o = 0., 0.
 ATE3o, ARE3o = 0., 0.
 ATE4o, ARE4o = ATF4, ARF4
 # terms with aCal
-ATE1a, ARE1a = 0. , 0.
+ATE1a, ARE1a = 0., 0.
 ATE2a, ARE2a = ATF2, ARF2
 ATE3a, ARE3a = -ATF3, -ARF3
-ATE4a, ARE4a = -2*ATF4, -2*ARF4
+ATE4a, ARE4a = -2 * ATF4, -2 * ARF4
 # rotated AinEa by epsilon
-ATE2ae =  C2e*ATF2 + S2e*ATF3
+ATE2ae = C2e * ATF2 + S2e * ATF3
-ATE3ae = -S2e*ATF2 - C2e*ATF3
+ATE3ae = -S2e * ATF2 - C2e * ATF3
-ARE2ae =  C2e*ARF2 + S2e*ARF3
+ARE2ae = C2e * ARF2 + S2e * ARF3
-ARE3ae = -S2e*ARF2 - C2e*ARF3
+ARE3ae = -S2e * ARF2 - C2e * ARF3
 ATE1 = ATE1o
-ATE2e = aCal*ATE2ae
+ATE2e = aCal * ATE2ae
-ATE3e = aCal*ATE3ae
+ATE3e = aCal * ATE3ae
-ATE4 = (1-2*aCal)*ATF4
+ATE4 = (1 - 2 * aCal) * ATF4
 ARE1 = ARE1o
-ARE2e = aCal*ARE2ae
+ARE2e = aCal * ARE2ae
-ARE3e = aCal*ARE3ae
+ARE3e = aCal * ARE3ae
-ARE4 = (1-2*aCal)*ARF4
+ARE4 = (1 - 2 * aCal) * ARF4
 # rotated IinE
-QinEe = C2e*QinE + S2e*UinE
+QinEe = C2e * QinE + S2e * UinE
-UinEe = C2e*UinE - S2e*QinE
+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
+if (TypeC == 2) or (TypeC == 1):  # +++++++++ rotator calibration Eq. C.4
-AT = ATE1o*IinE + (ATE4o+aCal*ATE4a)*h*VinE
+AT = ATE1o * IinE + (ATE4o + aCal * ATE4a) * h * VinE
-BT = aCal * (ATE3ae*QinEe - ATE2ae*h*UinEe)
+BT = aCal * (ATE3ae * QinEe - ATE2ae * h * UinEe)
-AR = ARE1o*IinE + (ARE4o+aCal*ARE4a)*h*VinE
+AR = ARE1o * IinE + (ARE4o + aCal * ARE4a) * h * VinE
-BR = aCal * (ARE3ae*QinEe - ARE2ae*h*UinEe)
+BR = aCal * (ARE3ae * QinEe - ARE2ae * h * UinEe)
 # Correction paremeters for normal measurements; they are independent of LDR
 if (not RotationErrorEpsilonForNormalMeasurements):
 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F)
-GT = ATE1o*IinE + ATE4o*h*VinE
+GT = ATE1o * IinE + ATE4o * h * VinE
-GR = ARE1o*IinE + ARE4o*h*VinE
+GR = ARE1o * IinE + ARE4o * h * VinE
-HT = ATE2a*QinE + ATE3a*h*UinEe + ATE4a*h*VinE
+HT = ATE2a * QinE + ATE3a * h * UinEe + ATE4a * h * VinE
-HR = ARE2a*QinE + ARE3a*h*UinEe + ARE4a*h*VinE
+HR = ARE2a * QinE + ARE3a * h * UinEe + ARE4a * h * VinE
 else:
-GT = ATE1o*IinE + ATE4o*h*VinE
+GT = ATE1o * IinE + ATE4o * h * VinE
-GR = ARE1o*IinE + ARE4o*h*VinE
+GR = ARE1o * IinE + ARE4o * h * VinE
-HT = ATE2ae*QinE + ATE3ae*h*UinEe + ATE4a*h*VinE
+HT = ATE2ae * QinE + ATE3ae * h * UinEe + ATE4a * h * VinE
-HR = ARE2ae*QinE + ARE3ae*h*UinEe + ARE4a*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
+AT = ATE1 * IinE + ZiC * CosC * (ATE2e * QinEe + ATE4 * VinE) + ATE3e * UinEe
-BT = DiC*(ATE1*UinEe + ATE3e*IinE) + ZiC*SinC*(ATE4*QinEe - ATE2e*VinE)
+BT = DiC * (ATE1 * UinEe + ATE3e * IinE) + ZiC * SinC * (ATE4 * QinEe - ATE2e * VinE)
-AR = ARE1*IinE + ZiC*CosC*(ARE2e*QinEe + ARE4*VinE) + ARE3e*UinEe
+AR = ARE1 * IinE + ZiC * CosC * (ARE2e * QinEe + ARE4 * VinE) + ARE3e * UinEe
-BR = DiC*(ARE1*UinEe + ARE3e*IinE) + ZiC*SinC*(ARE4*QinEe - ARE2e*VinE)
+BR = DiC * (ARE1 * UinEe + ARE3e * IinE) + ZiC * SinC * (ARE4 * QinEe - ARE2e * VinE)
 # Correction paremeters for normal measurements; they are independent of LDR
 if (not RotationErrorEpsilonForNormalMeasurements):
 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F)
-GT = ATE1o*IinE + ATE4o*VinE
+GT = ATE1o * IinE + ATE4o * VinE
-GR = ARE1o*IinE + ARE4o*VinE
+GR = ARE1o * IinE + ARE4o * VinE
-HT = ATE2a*QinE + ATE3a*UinE + ATE4a*VinE
+HT = ATE2a * QinE + ATE3a * UinE + ATE4a * VinE
-HR = ARE2a*QinE + ARE3a*UinE + ARE4a*VinE
+HR = ARE2a * QinE + ARE3a * UinE + ARE4a * VinE
 else:
-D = IinE + DiC*QinEe
+D = IinE + DiC * QinEe
-A = DiC*IinE + QinEe
+A = DiC * IinE + QinEe
-B = ZiC*(CosC*UinEe + SinC*VinE)
+B = ZiC * (CosC * UinEe + SinC * VinE)
-C = -ZiC*(SinC*UinEe - CosC*VinE)
+C = -ZiC * (SinC * UinEe - CosC * VinE)
-GT = ATE1o*D + ATE4o*C
+GT = ATE1o * D + ATE4o * C
-GR = ARE1o*D + ARE4o*C
+GR = ARE1o * D + ARE4o * C
-HT = ATE2a*A + ATE3a*B + ATE4a*C
+HT = ATE2a * A + ATE3a * B + ATE4a * C
-HR = ARE2a*A + ARE3a*B + ARE4a*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])
+IE1e = np.array([IinE + sqr05 * DiC * QinEe, sqr05 * DiC * IinE + (1 - 0.5 * WiC) * QinEe,
-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])
+(1 - 0.5 * WiC) * UinEe + sqr05 * ZiC * SinC * VinE,
-ATEe = np.array([ATE1,ATE2e,ATE3e,ATE4])
+-sqr05 * ZiC * SinC * UinEe + ZiC * CosC * VinE])
-AREe = np.array([ARE1,ARE2e,ARE3e,ARE4])
+IE2e = np.array([sqr05 * DiC * UinEe, 0.5 * WiC * UinEe - sqr05 * ZiC * SinC * VinE,
-AT = np.dot(ATEe,IE1e)
+sqr05 * DiC * IinE + 0.5 * WiC * QinEe, sqr05 * ZiC * SinC * QinEe])
-AR = np.dot(AREe,IE1e)
+ATEe = np.array([ATE1, ATE2e, ATE3e, ATE4])
-BT = np.dot(ATEe,IE2e)
+AREe = np.array([ARE1, ARE2e, ARE3e, ARE4])
-BR = np.dot(AREe,IE2e)
+AT = np.dot(ATEe, IE1e)
+AR = np.dot(AREe, IE1e)
+BT = np.dot(ATEe, IE2e)
+BR = np.dot(AREe, IE2e)
 # Correction paremeters for normal measurements; they are independent of LDR
-if (not RotationErrorEpsilonForNormalMeasurements):   # calibrator taken out
+if (not RotationErrorEpsilonForNormalMeasurements):  # calibrator taken out
-GT = ATE1o*IinE + ATE4o*VinE
+GT = ATE1o * IinE + ATE4o * VinE
-GR = ARE1o*IinE + ARE4o*VinE
+GR = ARE1o * IinE + ARE4o * VinE
-HT = ATE2a*QinE + ATE3a*UinE + ATE4a*VinE
+HT = ATE2a * QinE + ATE3a * UinE + ATE4a * VinE
-HR = ARE2a*QinE + ARE3a*UinE + ARE4a*VinE
+HR = ARE2a * QinE + ARE3a * UinE + ARE4a * VinE
 else:
-D = IinE + DiC*QinEe
+D = IinE + DiC * QinEe
-A = DiC*IinE + QinEe
+A = DiC * IinE + QinEe
-B = ZiC*(CosC*UinEe + SinC*VinE)
+B = ZiC * (CosC * UinEe + SinC * VinE)
-C = -ZiC*(SinC*UinEe - CosC*VinE)
+C = -ZiC * (SinC * UinEe - CosC * VinE)
-GT = ATE1o*D + ATE4o*C
+GT = ATE1o * D + ATE4o * C
-GR = ARE1o*D + ARE4o*C
+GR = ARE1o * D + ARE4o * C
-HT = ATE2a*A + ATE3a*B + ATE4a*C
+HT = ATE2a * A + ATE3a * B + ATE4a * C
-HR = ARE2a*A + ARE3a*B + ARE4a*C
+HR = ARE2a * A + ARE3a * B + ARE4a * C
 else:
 print('Calibrator not implemented yet')
 sys.exit()
 else:
 print("Calibrator location not implemented yet")
 sys.exit()
 # Determination of the correction K of the calibration factor
-IoutTp = TaT*TiT*TiO*TiE*(AT + BT)
+IoutTp = TaT * TiT * TiO * TiE * (AT + BT)
-IoutTm = TaT*TiT*TiO*TiE*(AT - BT)
+IoutTm = TaT * TiT * TiO * TiE * (AT - BT)
-IoutRp = TaR*TiR*TiO*TiE*(AR + BR)
+IoutRp = TaR * TiR * TiO * TiE * (AR + BR)
-IoutRm = TaR*TiR*TiO*TiE*(AR - BR)
+IoutRm = TaR * TiR * TiO * TiE * (AR - BR)
 # --- Results and Corrections; electronic etaR and etaT are assumed to be 1
-Etapx = IoutRp/IoutTp
+Etapx = IoutRp / IoutTp
-Etamx = IoutRm/IoutTm
+Etamx = IoutRm / IoutTm
-Etax = (Etapx*Etamx)**0.5
+Etax = (Etapx * Etamx) ** 0.5
-Eta = (TaR*TiR)/(TaT*TiT)   # Eta = Eta*/K  Eq. 84
+Eta = (TaR * TiR) / (TaT * TiT)  # Eta = Eta*/K  Eq. 84
 K = Etax / Eta
 #  For comparison with Volkers Libreoffice Müller Matrix spreadsheet
-#Eta_test_p = (IoutRp/IoutTp)
+# Eta_test_p = (IoutRp/IoutTp)
-#Eta_test_m = (IoutRm/IoutTm)
+# Eta_test_m = (IoutRm/IoutTm)
-#Eta_test = (Eta_test_p*Eta_test_m)**0.5
+# Eta_test = (Eta_test_p*Eta_test_m)**0.5
 # ----- Forward simulated signals and LDRsim with atrue; from input file
-It = TaT*TiT*TiO*TiE*(GT+atrue*HT)
+It = TaT * TiT * TiO * TiE * (GT + atrue * HT)
-Ir = TaR*TiR*TiO*TiE*(GR+atrue*HR)
+Ir = TaR * TiR * TiO * TiE * (GR + atrue * HR)
 # LDRsim = 1/Eta*Ir/It  # simulated LDR* with Y from input file
-LDRsim = Ir/It  # simulated uncorrected LDR with Y from input file
+LDRsim = Ir / It  # simulated uncorrected LDR with Y from input file
 # Corrected LDRsimCorr from forward simulated LDRsim (atrue)
 # LDRsimCorr = (1./Eta*LDRsim*(GT+HT)-(GR+HR))/((GR-HR)-1./Eta*LDRsim*(GT-HT))
 if Y == -1.:
-LDRsimx = 1./LDRsim
+LDRsimx = 1. / LDRsim
 else:
 LDRsimx = LDRsim
 # The following is correct without doubt
-#LDRCorr = (LDRsim*K/Etax*(GT+HT)-(GR+HR))/((GR-HR)-LDRsim*K/Etax*(GT-HT))
+# LDRCorr = (LDRsim*K/Etax*(GT+HT)-(GR+HR))/((GR-HR)-LDRsim*K/Etax*(GT-HT))
 # The following is a test whether the equations for calibration Etax and normal  signal (GHK, LDRsim) are consistent
-LDRCorr = (LDRsim/Eta*(GT+HT)-(GR+HR))/((GR-HR)-LDRsim*K/Etax*(GT-HT))
+LDRCorr = (LDRsim / Eta * (GT + HT) - (GR + HR)) / ((GR - HR) - LDRsim * K / Etax * (GT - HT))
-TTa = TiT*TaT #*ATP1
+TTa = TiT * TaT  # *ATP1
-TRa = TiR*TaR #*ARP1
+TRa = TiR * TaR  # *ARP1
-F11sim = 1/(TiO*TiE)*((HR*Etax/K*It/TTa-HT*Ir/TRa)/(HR*GT-HT*GR))    # IL = 1, Etat = Etar = 1
+F11sim = 1 / (TiO * TiE) * (
+(HR * Etax / K * It / TTa - HT * Ir / TRa) / (HR * GT - HT * GR))  # IL = 1, Etat = Etar = 1
 return (GT, HT, GR, HR, K, Eta, LDRsimx, LDRCorr, DTa, DRa, TTa, TRa, F11sim)
 # *******************************************************************************************************************************
 # --- CALC truth
-GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0)
+GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0,
+RetO0, DiO0, RotC0, RetC0, DiC0,
+TP0, TS0, RP0, RS0, ERaT0,
+RotaT0, RetT0, ERaR0, RotaR0,
+RetR0, LDRCal0)
 # --- Print parameters to console and output file
 with open('output_files\output_' + LID + '.dat', 'w') as f:
 with redirect_stdout(f):
 print("From ", dname)
 print("Running ", fname)
-print("Reading input file ", InputFile) #, "  for Lidar system :", EID, ", ", LID)
+print("Reading input file ", InputFile)  # , "  for Lidar system :", EID, ", ", LID)
 print("for Lidar system: ", EID, ", ", LID)
 # --- Print iput information*********************************
 print(" --- Input parameters: value ±error / ±steps  ----------------------")
-print("{0:8} {1:8} {2:8.5f}; {3:8} {4:7.4f}±{5:7.4f}/{6:2d}".format("Laser: ", "DOLP = ", bL, "        rotation alpha = ", RotL0, dRotL, nRotL))
+print("{0:8} {1:8} {2:8.5f}; {3:8} {4:7.4f}±{5:7.4f}/{6:2d}".format("Laser: ", "DOLP = ", bL,
+"        rotation alpha = ", RotL0, dRotL,
+nRotL))
 print("              Diatt.,             Tunpol,   Retard.,   Rotation (deg)")
-print("{0:12} {1:7.4f}±{2:7.4f}/{8:2d}, {3:7.4f}, {4:3.0f}±{5:3.0f}/{9:2d}, {6:7.4f}±{7:7.4f}/{10:2d}".format("Emitter    ", DiE0, dDiE, TiE, RetE0, dRetE, RotE0, dRotE, nDiE, nRetE, nRotE))
+print("{0:12} {1:7.4f}±{2:7.4f}/{8:2d}, {3:7.4f}, {4:3.0f}±{5:3.0f}/{9:2d}, {6:7.4f}±{7:7.4f}/{10:2d}".format(
-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))
+"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("Calibrator ", DiC0, dDiC, TiC, RetC0, dRetC, RotC0, dRotC, nDiC, nRetC, nRotC))
+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(
-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))
+"{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("TP,TS           :", TP0, dTP, nTP, TS0,
-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))
+dTS, nTS))
+print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format("RP,RS           :", RP0, dRP, nRP, RS0,
+dRS, nRS))
 print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}, {5:1.0f}".format("DT,TT,DR,TR,Y   :", DiT, TiT, DiR, TiR, Y))
 print("{0:12}".format(" --- Combined PBS + Pol.-filter ---"))
 print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}".format("DT,TT,DR,TR     :", DTa0, TTa0, DRa0, TRa0))
 print()
 print("Rotation Error Epsilon For Normal Measurements = ", RotationErrorEpsilonForNormalMeasurements)
 print(Type[TypeC], Loc[LocC])
-print("Parallel signal detected in", dY[int(Y+1)])
+print("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()
-#print(" --- LDRCal during calibration | simulated and corrected LDRs -------------")
+# 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} |{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.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: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: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("{0:8}       |{1:8}->{2:8}->{3:8}".format(" LDRCal"," LDRtrue", " LDRsimx", " LDRcorr"))
-#print(" --- LDRCal during calibration")
+# print(" --- LDRCal during calibration")
-print("{0:26}: {1:6.3f}±{2:5.3f}/{3:2d}".format("LDRCal during calibration in calibration range", LDRCal0, dLDRCal, nLDRCal))
+print("{0:26}: {1:6.3f}±{2:5.3f}/{3:2d}".format("LDRCal during calibration in calibration range", LDRCal0,
+dLDRCal, nLDRCal))
-#print("{0:8}={1:8.5f};{2:8}={3:8.5f}".format(" IinP",IinP," F11sim",F11sim))
+# print("{0:8}={1:8.5f};{2:8}={3:8.5f}".format(" IinP",IinP," F11sim",F11sim))
 print()
 K0List = np.zeros(3)
 LDRsimxList = np.zeros(3)
 LDRCalList = 0.004, 0.2, 0.45
-for i,LDRCal in enumerate(LDRCalList):
+for i, LDRCal in enumerate(LDRCalList):
-GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal)
+GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0,
+DiE0, RotO0, RetO0,
+DiO0, RotC0, RetC0,
+DiC0, TP0, TS0, RP0,
+RS0, ERaT0, RotaT0,
+RetT0, ERaR0, RotaR0,
+RetR0, LDRCal)
 K0List[i] = K0
 LDRsimxList[i] = LDRsimx
 print('========================================================================')
-print("{0:8},{1:8},{2:8},{3:8},{4:9},{5:8},{6:9}".format(" GR", " GT", " HR", " HT", "  K(0.004)", " K(0.2)", "  K(0.45)"))
+print("{0:8},{1:8},{2:8},{3:8},{4:9},{5:8},{6:9}".format(" GR", " GT", " HR", " HT", "  K(0.004)", " K(0.2)",
-print("{0:8.5f},{1:8.5f},{2:8.5f},{3:8.5f},{4:9.5f},{5:9.5f},{6:9.5f}".format(GR0, GT0, HR0, HT0, K0List[0], K0List[1], K0List[2]))
+"  K(0.45)"))
+print("{0:8.5f},{1:8.5f},{2:8.5f},{3:8.5f},{4:9.5f},{5:9.5f},{6:9.5f}".format(GR0, GT0, HR0, HT0, K0List[0],
+K0List[1], K0List[2]))
 print('========================================================================')
 print("{0:9},{1:9},{2:9}".format("  LDRtrue", "  LDRsimx", "  LDRCorr"))
 LDRtrueList = 0.004, 0.02, 0.2, 0.45
-for i,LDRtrue in enumerate(LDRtrueList):
+for i, LDRtrue in enumerate(LDRtrueList):
-GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0)
+GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0,
+DiE0, RotO0, RetO0,
+DiO0, RotC0, RetC0,
+DiC0, TP0, TS0, RP0,
+RS0, ERaT0, RotaT0,
+RetT0, ERaR0, RotaR0,
+RetR0, LDRCal0)
 print("{0:9.5f},{1:9.5f},{2:9.5f}".format(LDRtrue, LDRsimx, LDRCorr))
 file = open('output_files\output_' + LID + '.dat', 'r')
-print (file.read())
+print(file.read())
 file.close()
 '''
 if(PrintToOutputFile):
 f = open('output_ver7.dat', 'w')
 if(PrintToOutputFile):
 sys.stdout.flush()
 f.close
 sys.stdout = old_target
 '''
-if(Error_Calc):
+if (Error_Calc):
 # --- CALC again truth with LDRCal0 to reset all 0-values
-GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0)
+GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0,
+RotO0, RetO0, DiO0, RotC0,
+RetC0, DiC0, TP0, TS0, RP0,
+RS0, ERaT0, RotaT0, RetT0,
+ERaR0, RotaR0, RetR0,
+LDRCal0)
 # --- Start Errors calculation with variable parameters ------------------------------------------------------------------
 iN = -1
-N = ((nRotL*2+1)*
+N = ((nRotL * 2 + 1) *
-(nRotE*2+1)*(nRetE*2+1)*(nDiE*2+1)*
+(nRotE * 2 + 1) * (nRetE * 2 + 1) * (nDiE * 2 + 1) *
-(nRotO*2+1)*(nRetO*2+1)*(nDiO*2+1)*
+(nRotO * 2 + 1) * (nRetO * 2 + 1) * (nDiO * 2 + 1) *
-(nRotC*2+1)*(nRetC*2+1)*(nDiC*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)*
+(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))
+(nRotaT * 2 + 1) * (nRotaR * 2 + 1) * (nRetT * 2 + 1) * (nRetR * 2 + 1) * (nLDRCal * 2 + 1))
-print("N = ",N ," ", end="")
+print("N = ", N, " ", end="")
 if N > 1e6:
-if user_yes_no_query('Warning: processing ' + str(N) + ' samples will take very long. Do you want to proceed?') == 0: sys.exit()
+if 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: 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()
+# if user_yes_no_query('Warning: processing' + str(N) + ' samples will take very long. Do you want to proceed?') == 0: sys.exit()
 # --- Arrays for plotting ------
 LDRmin = np.zeros(5)
 LDRmax = np.zeros(5)
 F11min = np.zeros(5)
 F11max = np.zeros(5)
 LDRrange = np.zeros(5)
 LDRrange = 0.004, 0.02, 0.1, 0.3, 0.45
-#aLDRsimx = np.zeros(N)
+# aLDRsimx = np.zeros(N)
-#aLDRsimx2 = np.zeros(N)
+# aLDRsimx2 = np.zeros(N)
-#aLDRcorr = np.zeros(N)
+# aLDRcorr = np.zeros(N)
-#aLDRcorr2 = np.zeros(N)
+# aLDRcorr2 = np.zeros(N)
 aERaT = np.zeros(N)
 aERaR = np.zeros(N)
 aRotaT = np.zeros(N)
 aRotaR = np.zeros(N)
 aRetT = np.zeros(N)
 aRetE = np.zeros(N)
 aRotE = np.zeros(N)
 aRetO = np.zeros(N)
 aRotO = np.zeros(N)
 aLDRCal = np.zeros(N)
-aA = np.zeros((5,N))
+aA = np.zeros((5, N))
-aX = np.zeros((5,N))
+aX = np.zeros((5, N))
-aF11corr = np.zeros((5,N))
+aF11corr = np.zeros((5, N))
 atime = clock()
 dtime = clock()
 # --- Calc Error signals
-#GT, HT, GR, HR, K, Eta, LDRsim = Calc(RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC, TP, TS)
+# GT, HT, GR, HR, K, Eta, LDRsim = Calc(RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC, TP, TS)
 # ---- Do the calculations of bra-ket vectors
 h = -1. if TypeC == 2 else 1
 # from input file: measured LDRm and true LDRtrue, LDRtrue2  =>
-ameas = (1.-LDRmeas)/(1+LDRmeas)
+ameas = (1. - LDRmeas) / (1 + LDRmeas)
-atrue = (1.-LDRtrue)/(1+LDRtrue)
+atrue = (1. - LDRtrue) / (1 + LDRtrue)
-atrue2 = (1.-LDRtrue2)/(1+LDRtrue2)
+atrue2 = (1. - LDRtrue2) / (1 + LDRtrue2)
-for iLDRCal in range(-nLDRCal,nLDRCal+1):
+for iLDRCal in range(-nLDRCal, nLDRCal + 1):
 # from input file:  assumed LDRCal for calibration measurements
 LDRCal = LDRCal0
-if nLDRCal > 0: LDRCal = LDRCal0 + iLDRCal*dLDRCal/nLDRCal
+if nLDRCal > 0: LDRCal = LDRCal0 + iLDRCal * dLDRCal / nLDRCal
-GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal)
+GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0,
-aCal = (1.-LDRCal)/(1+LDRCal)
+DiE0, RotO0, RetO0, DiO0,
+RotC0, RetC0, DiC0, TP0,
+TS0, RP0, RS0, ERaT0,
+RotaT0, RetT0, ERaR0,
+RotaR0, RetR0, LDRCal)
+aCal = (1. - LDRCal) / (1 + LDRCal)
 for iRotL, iRotE, iRetE, iDiE \
-in [(iRotL,iRotE,iRetE,iDiE)
+in [(iRotL, iRotE, iRetE, iDiE)
-for iRotL in range(-nRotL,nRotL+1)
+for iRotL in range(-nRotL, nRotL + 1)
-for iRotE in range(-nRotE,nRotE+1)
+for iRotE in range(-nRotE, nRotE + 1)
-for iRetE in range(-nRetE,nRetE+1)
+for iRetE in range(-nRetE, nRetE + 1)
-for iDiE in range(-nDiE,nDiE+1)]:
+for iDiE in range(-nDiE, nDiE + 1)]:
-if nRotL > 0: RotL = RotL0 + iRotL*dRotL/nRotL
+if nRotL > 0: RotL = RotL0 + iRotL * dRotL / nRotL
-if nRotE > 0: RotE = RotE0 + iRotE*dRotE/nRotE
+if nRotE > 0: RotE = RotE0 + iRotE * dRotE / nRotE
-if nRetE > 0: RetE = RetE0 + iRetE*dRetE/nRetE
+if nRetE > 0: RetE = RetE0 + iRetE * dRetE / nRetE
-if nDiE > 0:  DiE  = DiE0  + iDiE*dDiE/nDiE
+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))
+S2a = np.sin(2 * np.deg2rad(RotL))
-C2a = np.cos(2*np.deg2rad(RotL))
+C2a = np.cos(2 * np.deg2rad(RotL))
-S2b = np.sin(2*np.deg2rad(RotE))
+S2b = np.sin(2 * np.deg2rad(RotE))
-C2b = np.cos(2*np.deg2rad(RotE))
+C2b = np.cos(2 * np.deg2rad(RotE))
-S2ab = np.sin(np.deg2rad(2*RotL-2*RotE))
+S2ab = np.sin(np.deg2rad(2 * RotL - 2 * RotE))
-C2ab = np.cos(np.deg2rad(2*RotL-2*RotE))
+C2ab = np.cos(np.deg2rad(2 * RotL - 2 * RotE))
 # Laser with Degree of linear polarization DOLP = bL
 IinL = 1.
 QinL = bL
 UinL = 0.
-VinL = (1. - bL**2)**0.5
+VinL = (1. - bL ** 2) ** 0.5
 # Stokes Input Vector rotation Eq. E.4
-A = C2a*QinL - S2a*UinL
+A = C2a * QinL - S2a * UinL
-B = S2a*QinL + C2a*UinL
+B = S2a * QinL + C2a * UinL
 # Stokes Input Vector rotation Eq. E.9
-C = C2ab*QinL - S2ab*UinL
+C = C2ab * QinL - S2ab * UinL
-D = S2ab*QinL + C2ab*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
+ZiE = (1. - DiE ** 2) ** 0.5
-WiE = (1. - ZiE*CosE)
+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)
+IinE = (IinL + DiE * C)
-QinE = (C2b*DiE*IinL + A + S2b*(WiE*D - ZiE*SinE*VinL))
+QinE = (C2b * DiE * IinL + A + S2b * (WiE * D - ZiE * SinE * VinL))
-UinE = (S2b*DiE*IinL + B - C2b*(WiE*D - ZiE*SinE*VinL))
+UinE = (S2b * DiE * IinL + B - C2b * (WiE * D - ZiE * SinE * VinL))
-VinE = (-ZiE*SinE*D + ZiE*CosE*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
+# 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 )
+in [
-for iRotO in range(-nRotO,nRotO+1)
+(iRotO, iRetO, iDiO, iRotC, iRetC, iDiC, iTP, iTS, iRP, iRS, iERaT, iRotaT, iRetT, iERaR, iRotaR, iRetR)
-for iRetO in range(-nRetO,nRetO+1)
+for iRotO in range(-nRotO, nRotO + 1)
-for iDiO in range(-nDiO,nDiO+1)
+for iRetO in range(-nRetO, nRetO + 1)
-for iRotC in range(-nRotC,nRotC+1)
+for iDiO in range(-nDiO, nDiO + 1)
-for iRetC in range(-nRetC,nRetC+1)
+for iRotC in range(-nRotC, nRotC + 1)
-for iDiC in range(-nDiC,nDiC+1)
+for iRetC in range(-nRetC, nRetC + 1)
-for iTP in range(-nTP,nTP+1)
+for iDiC in range(-nDiC, nDiC + 1)
-for iTS in range(-nTS,nTS+1)
+for iTP in range(-nTP, nTP + 1)
-for iRP in range(-nRP,nRP+1)
+for iTS in range(-nTS, nTS + 1)
-for iRS in range(-nRS,nRS+1)
+for iRP in range(-nRP, nRP + 1)
-for iERaT in range(-nERaT,nERaT+1)
+for iRS in range(-nRS, nRS + 1)
-for iRotaT in range(-nRotaT,nRotaT+1)
+for iERaT in range(-nERaT, nERaT + 1)
-for iRetT in range(-nRetT,nRetT+1)
+for iRotaT in range(-nRotaT, nRotaT + 1)
-for iERaR in range(-nERaR,nERaR+1)
+for iRetT in range(-nRetT, nRetT + 1)
-for iRotaR in range(-nRotaR,nRotaR+1)
+for iERaR in range(-nERaR, nERaR + 1)
-for iRetR in range(-nRetR,nRetR+1)]:
+for iRotaR in range(-nRotaR, nRotaR + 1)
+for iRetR in range(-nRetR, nRetR + 1)]:
 iN = iN + 1
 if (iN == 10001):
 ctime = clock()
-print(" estimated time ", "{0:4.2f}".format(N/10000 * (ctime-atime)), "sec ") #, end="")
+print(" 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")
+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")
+print("\r elapsed time ", "{0:5.0f}".format((ctime - atime)), "sec ", end="\r")
 dtime = ctime
-if nRotO > 0: RotO = RotO0 + iRotO*dRotO/nRotO
+if nRotO > 0: RotO = RotO0 + iRotO * dRotO / nRotO
-if nRetO > 0: RetO = RetO0 + iRetO*dRetO/nRetO
+if nRetO > 0: RetO = RetO0 + iRetO * dRetO / nRetO
-if nDiO > 0:  DiO  = DiO0  + iDiO*dDiO/nDiO
+if nDiO > 0:  DiO = DiO0 + iDiO * dDiO / nDiO
-if nRotC > 0: RotC = RotC0 + iRotC*dRotC/nRotC
+if nRotC > 0: RotC = RotC0 + iRotC * dRotC / nRotC
-if nRetC > 0: RetC = RetC0 + iRetC*dRetC/nRetC
+if nRetC > 0: RetC = RetC0 + iRetC * dRetC / nRetC
-if nDiC > 0:  DiC  = DiC0  + iDiC*dDiC/nDiC
+if nDiC > 0:  DiC = DiC0 + iDiC * dDiC / nDiC
-if nTP > 0:   TP   = TP0   + iTP*dTP/nTP
+if nTP > 0:   TP = TP0 + iTP * dTP / nTP
-if nTS > 0:   TS   = TS0   + iTS*dTS/nTS
+if nTS > 0:   TS = TS0 + iTS * dTS / nTS
-if nRP > 0:   RP   = RP0   + iRP*dRP/nRP
+if nRP > 0:   RP = RP0 + iRP * dRP / nRP
-if nRS > 0:   RS   = RS0   + iRS*dRS/nRS
+if nRS > 0:   RS = RS0 + iRS * dRS / nRS
-if nERaT > 0: ERaT = ERaT0 + iERaT*dERaT/nERaT
+if nERaT > 0: ERaT = ERaT0 + iERaT * dERaT / nERaT
-if nRotaT > 0:RotaT= RotaT0+ iRotaT*dRotaT/nRotaT
+if nRotaT > 0: RotaT = RotaT0 + iRotaT * dRotaT / nRotaT
-if nRetT > 0: RetT = RetT0 + iRetT*dRetT/nRetT
+if nRetT > 0: RetT = RetT0 + iRetT * dRetT / nRetT
-if nERaR > 0: ERaR = ERaR0 + iERaR*dERaR/nERaR
+if nERaR > 0: ERaR = ERaR0 + iERaR * dERaR / nERaR
-if nRotaR > 0:RotaR= RotaR0+ iRotaR*dRotaR/nRotaR
+if nRotaR > 0: RotaR = RotaR0 + iRotaR * dRotaR / nRotaR
-if nRetR > 0: RetR = RetR0 + iRetR*dRetR/nRetR
+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))
+# 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
+ZiO = (1. - DiO ** 2) ** 0.5
-WiO = (1. - ZiO*CosO)
+WiO = (1. - ZiO * CosO)
-S2g = np.sin(np.deg2rad(2*RotO))
+S2g = np.sin(np.deg2rad(2 * RotO))
-C2g = np.cos(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
+ZiC = (1. - DiC ** 2) ** 0.5
-WiC = (1. - ZiC*CosC)
+WiC = (1. - ZiC * CosC)
 # analyser
-#For POLLY_XTs
+# For POLLY_XTs
-if(RS_RP_depend_on_TS_TP):
+if (RS_RP_depend_on_TS_TP):
 RS = 1 - TS
 RP = 1 - TP
 TiT = 0.5 * (TP + TS)
-DiT = (TP-TS)/(TP+TS)
+DiT = (TP - TS) / (TP + TS)
-ZiT = (1. - DiT**2)**0.5
+ZiT = (1. - DiT ** 2) ** 0.5
 TiR = 0.5 * (RP + RS)
-DiR = (RP-RS)/(RP+RS)
+DiR = (RP - RS) / (RP + RS)
-ZiR = (1. - DiR**2)**0.5
+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)
+DaT = (1 - ERaT) / (1 + ERaT)
-DaR = (1-ERaR)/(1+ERaR)
+DaR = (1 - ERaR) / (1 + ERaR)
-TaT = 0.5*(1+ERaT)
+TaT = 0.5 * (1 + ERaT)
-TaR = 0.5*(1+ERaR)
+TaR = 0.5 * (1 + ERaR)
-S2aT = np.sin(np.deg2rad(h*2*RotaT))
+S2aT = np.sin(np.deg2rad(h * 2 * RotaT))
-C2aT = np.cos(np.deg2rad(2*RotaT))
+C2aT = np.cos(np.deg2rad(2 * RotaT))
-S2aR = np.sin(np.deg2rad(h*2*RotaR))
+S2aR = np.sin(np.deg2rad(h * 2 * RotaR))
-C2aR = np.cos(np.deg2rad(2*RotaR))
+C2aR = np.cos(np.deg2rad(2 * RotaR))
 # Aanalyzer As before the PBS Eq. D.5
-ATP1 = (1+C2aT*DaT*DiT)
+ATP1 = (1 + C2aT * DaT * DiT)
-ATP2 = Y*(DiT+C2aT*DaT)
+ATP2 = Y * (DiT + C2aT * DaT)
-ATP3 = Y*S2aT*DaT*ZiT*CosT
+ATP3 = Y * S2aT * DaT * ZiT * CosT
-ATP4 = S2aT*DaT*ZiT*SinT
+ATP4 = S2aT * DaT * ZiT * SinT
-ATP = np.array([ATP1,ATP2,ATP3,ATP4])
+ATP = np.array([ATP1, ATP2, ATP3, ATP4])
-ARP1 = (1+C2aR*DaR*DiR)
+ARP1 = (1 + C2aR * DaR * DiR)
-ARP2 = Y*(DiR+C2aR*DaR)
+ARP2 = Y * (DiR + C2aR * DaR)
-ARP3 = Y*S2aR*DaR*ZiR*CosR
+ARP3 = Y * S2aR * DaR * ZiR * CosR
-ARP4 = S2aR*DaR*ZiR*SinR
+ARP4 = S2aR * DaR * ZiR * SinR
-ARP = np.array([ARP1,ARP2,ARP3,ARP4])
+ARP = np.array([ARP1, ARP2, ARP3, ARP4])
-TTa = TiT*TaT #*ATP1
+TTa = TiT * TaT  # *ATP1
-TRa = TiR*TaR #*ARP1
+TRa = TiR * TaR  # *ARP1
 # ---- Calculate signals and correction parameters for diffeent locations and calibrators
 if LocC == 4:  # Calibrator before the PBS
-#print("Calibrator location not implemented yet")
+# print("Calibrator location not implemented yet")
-#S2ge = np.sin(np.deg2rad(2*RotO + h*2*RotC))
+# S2ge = np.sin(np.deg2rad(2*RotO + h*2*RotC))
-#C2ge = np.cos(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))
+S2e = np.sin(np.deg2rad(h * 2 * RotC))
-C2e = np.cos(np.deg2rad(2*RotC))
+C2e = np.cos(np.deg2rad(2 * RotC))
 # rotated AinP by epsilon Eq. C.3
-ATP2e = C2e*ATP2 + S2e*ATP3
+ATP2e = C2e * ATP2 + S2e * ATP3
-ATP3e = C2e*ATP3 - S2e*ATP2
+ATP3e = C2e * ATP3 - S2e * ATP2
-ARP2e = C2e*ARP2 + S2e*ARP3
+ARP2e = C2e * ARP2 + S2e * ARP3
-ARP3e = C2e*ARP3 - S2e*ARP2
+ARP3e = C2e * ARP3 - S2e * ARP2
-ATPe = np.array([ATP1,ATP2e,ATP3e,ATP4])
+ATPe = np.array([ATP1, ATP2e, ATP3e, ATP4])
-ARPe = np.array([ARP1,ARP2e,ARP3e,ARP4])
+ARPe = np.array([ARP1, ARP2e, ARP3e, ARP4])
 # Stokes Input Vector before the polarising beam splitter Eq. E.31
-A = C2g*QinE - S2g*UinE
+A = C2g * QinE - S2g * UinE
-B = S2g*QinE + C2g*UinE
+B = S2g * QinE + C2g * UinE
-#C = (WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE)
+# C = (WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE)
-Co = ZiO*SinO*VinE
+Co = ZiO * SinO * VinE
-Ca = (WiO*B - 2*ZiO*SinO*VinE)
+Ca = (WiO * B - 2 * ZiO * SinO * VinE)
-#C = Co + aCal*Ca
+# C = Co + aCal*Ca
-#IinP = (IinE + DiO*aCal*A)
+# IinP = (IinE + DiO*aCal*A)
-#QinP = (C2g*DiO*IinE + aCal*QinE - S2g*C)
+# QinP = (C2g*DiO*IinE + aCal*QinE - S2g*C)
-#UinP = (S2g*DiO*IinE - aCal*UinE + C2g*C)
+# UinP = (S2g*DiO*IinE - aCal*UinE + C2g*C)
-#VinP = (ZiO*SinO*aCal*B + ZiO*CosO*(1-2*aCal)*VinE)
+# VinP = (ZiO*SinO*aCal*B + ZiO*CosO*(1-2*aCal)*VinE)
 IinPo = IinE
-QinPo = (C2g*DiO*IinE - S2g*Co)
+QinPo = (C2g * DiO * IinE - S2g * Co)
-UinPo = (S2g*DiO*IinE + C2g*Co)
+UinPo = (S2g * DiO * IinE + C2g * Co)
-VinPo = ZiO*CosO*VinE
+VinPo = ZiO * CosO * VinE
-IinPa = DiO*A
+IinPa = DiO * A
-QinPa = QinE - S2g*Ca
+QinPa = QinE - S2g * Ca
-UinPa = -UinE + C2g*Ca
+UinPa = -UinE + C2g * Ca
-VinPa = ZiO*(SinO*B - 2*CosO*VinE)
+VinPa = ZiO * (SinO * B - 2 * CosO * VinE)
-IinP = IinPo + aCal*IinPa
+IinP = IinPo + aCal * IinPa
-QinP = QinPo + aCal*QinPa
+QinP = QinPo + aCal * QinPa
-UinP = UinPo + aCal*UinPa
+UinP = UinPo + aCal * UinPa
-VinP = VinPo + aCal*VinPa
+VinP = VinPo + aCal * VinPa
 # Stokes Input Vector before the polarising beam splitter rotated by epsilon Eq. C.3
-#QinPe = C2e*QinP + S2e*UinP
+# QinPe = C2e*QinP + S2e*UinP
-#UinPe = C2e*UinP - S2e*QinP
+# UinPe = C2e*UinP - S2e*QinP
-QinPoe = C2e*QinPo + S2e*UinPo
+QinPoe = C2e * QinPo + S2e * UinPo
-UinPoe = C2e*UinPo - S2e*QinPo
+UinPoe = C2e * UinPo - S2e * QinPo
-QinPae = C2e*QinPa + S2e*UinPa
+QinPae = C2e * QinPa + S2e * UinPa
-UinPae = C2e*UinPa - S2e*QinPa
+UinPae = C2e * UinPa - S2e * QinPa
-QinPe = C2e*QinP + S2e*UinP
+QinPe = C2e * QinP + S2e * UinP
-UinPe = C2e*UinP - S2e*QinP
+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
+AT = ATP1 * IinP + h * ATP4 * VinP
-BT = ATP3e*QinP - h*ATP2e*UinP
+BT = ATP3e * QinP - h * ATP2e * UinP
-AR = ARP1*IinP + h*ARP4*VinP
+AR = ARP1 * IinP + h * ARP4 * VinP
-BR = ARP3e*QinP - h*ARP2e*UinP
+BR = ARP3e * QinP - h * ARP2e * UinP
 # Correction paremeters for normal measurements; they are independent of LDR
-if (not RotationErrorEpsilonForNormalMeasurements):   # calibrator taken out
+if (not RotationErrorEpsilonForNormalMeasurements):  # calibrator taken out
-IS1 = np.array([IinPo,QinPo,UinPo,VinPo])
+IS1 = np.array([IinPo, QinPo, UinPo, VinPo])
-IS2 = np.array([IinPa,QinPa,UinPa,VinPa])
+IS2 = np.array([IinPa, QinPa, UinPa, VinPa])
-GT = np.dot(ATP,IS1)
+GT = np.dot(ATP, IS1)
-GR = np.dot(ARP,IS1)
+GR = np.dot(ARP, IS1)
-HT = np.dot(ATP,IS2)
+HT = np.dot(ATP, IS2)
-HR = np.dot(ARP,IS2)
+HR = np.dot(ARP, IS2)
 else:
-IS1 = np.array([IinPo,QinPo,UinPo,VinPo])
+IS1 = np.array([IinPo, QinPo, UinPo, VinPo])
-IS2 = np.array([IinPa,QinPa,UinPa,VinPa])
+IS2 = np.array([IinPa, QinPa, UinPa, VinPa])
-GT = np.dot(ATPe,IS1)
+GT = np.dot(ATPe, IS1)
-GR = np.dot(ARPe,IS1)
+GR = np.dot(ARPe, IS1)
-HT = np.dot(ATPe,IS2)
+HT = np.dot(ATPe, IS2)
-HR = np.dot(ARPe,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)
+AT = ATP1 * IinP + ATP3e * UinPe + ZiC * CosC * (ATP2e * QinPe + ATP4 * VinP)
-BT = DiC*(ATP1*UinPe + ATP3e*IinP) - ZiC*SinC*(ATP2e*VinP - ATP4*QinPe)
+BT = DiC * (ATP1 * UinPe + ATP3e * IinP) - ZiC * SinC * (ATP2e * VinP - ATP4 * QinPe)
-AR = ARP1*IinP + ARP3e*UinPe + ZiC*CosC*(ARP2e*QinPe + ARP4*VinP)
+AR = ARP1 * IinP + ARP3e * UinPe + ZiC * CosC * (ARP2e * QinPe + ARP4 * VinP)
-BR = DiC*(ARP1*UinPe + ARP3e*IinP) - ZiC*SinC*(ARP2e*VinP - ARP4*QinPe)
+BR = DiC * (ARP1 * UinPe + ARP3e * IinP) - ZiC * SinC * (ARP2e * VinP - ARP4 * QinPe)
 # Correction paremeters for normal measurements; they are independent of LDR
-if (not RotationErrorEpsilonForNormalMeasurements):   # calibrator taken out
+if (not RotationErrorEpsilonForNormalMeasurements):  # calibrator taken out
-IS1 = np.array([IinPo,QinPo,UinPo,VinPo])
+IS1 = np.array([IinPo, QinPo, UinPo, VinPo])
-IS2 = np.array([IinPa,QinPa,UinPa,VinPa])
+IS2 = np.array([IinPa, QinPa, UinPa, VinPa])
-GT = np.dot(ATP,IS1)
+GT = np.dot(ATP, IS1)
-GR = np.dot(ARP,IS1)
+GR = np.dot(ARP, IS1)
-HT = np.dot(ATP,IS2)
+HT = np.dot(ATP, IS2)
-HR = np.dot(ARP,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)])
+IS1e = np.array(
-IS2e = np.array([IinPa+DiC*QinPae,DiC*IinPa+QinPae,ZiC*(CosC*UinPae+SinC*VinPa),-ZiC*(SinC*UinPae-CosC*VinPa)])
+[IinPo + DiC * QinPoe, DiC * IinPo + QinPoe, ZiC * (CosC * UinPoe + SinC * VinPo),
-GT = np.dot(ATPe,IS1e)
+-ZiC * (SinC * UinPoe - CosC * VinPo)])
-GR = np.dot(ARPe,IS1e)
+IS2e = np.array(
-HT = np.dot(ATPe,IS2e)
+[IinPa + DiC * QinPae, DiC * IinPa + QinPae, ZiC * (CosC * UinPae + SinC * VinPa),
-HR = np.dot(ARPe,IS2e)
+-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)
+AT = ATP1 * IinP + sqr05 * DiC * (ATP1 * QinPe + ATP2e * IinP) + (1 - 0.5 * WiC) * (
-BT = sqr05*DiC*(ATP1*UinPe + ATP3e*IinP) + 0.5*WiC*(ATP2e*UinPe + ATP3e*QinPe) - sqr05*ZiC*SinC*(ATP2e*VinP - ATP4*QinPe)
+ATP2e * QinPe + ATP3e * UinPe) + ZiC * (
-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)
+sqr05 * SinC * (ATP3e * VinP - ATP4 * UinPe) + ATP4 * CosC * VinP)
-BR = sqr05*DiC*(ARP1*UinPe + ARP3e*IinP) + 0.5*WiC*(ARP2e*UinPe + ARP3e*QinPe) - sqr05*ZiC*SinC*(ARP2e*VinP - ARP4*QinPe)
+BT = sqr05 * DiC * (ATP1 * UinPe + ATP3e * IinP) + 0.5 * WiC * (
+ATP2e * UinPe + ATP3e * QinPe) - sqr05 * ZiC * SinC * (ATP2e * VinP - ATP4 * QinPe)
+AR = ARP1 * IinP + sqr05 * DiC * (ARP1 * QinPe + ARP2e * IinP) + (1 - 0.5 * WiC) * (
+ARP2e * QinPe + ARP3e * UinPe) + ZiC * (
+sqr05 * SinC * (ARP3e * VinP - ARP4 * UinPe) + ARP4 * CosC * VinP)
+BR = sqr05 * DiC * (ARP1 * UinPe + ARP3e * IinP) + 0.5 * WiC * (
+ARP2e * UinPe + ARP3e * QinPe) - sqr05 * ZiC * SinC * (ARP2e * VinP - ARP4 * QinPe)
 # Correction paremeters for normal measurements; they are independent of LDR
-if (not RotationErrorEpsilonForNormalMeasurements):   # calibrator taken out
+if (not RotationErrorEpsilonForNormalMeasurements):  # calibrator taken out
-IS1 = np.array([IinPo,QinPo,UinPo,VinPo])
+IS1 = np.array([IinPo, QinPo, UinPo, VinPo])
-IS2 = np.array([IinPa,QinPa,UinPa,VinPa])
+IS2 = np.array([IinPa, QinPa, UinPa, VinPa])
-GT = np.dot(ATP,IS1)
+GT = np.dot(ATP, IS1)
-GR = np.dot(ARP,IS1)
+GR = np.dot(ARP, IS1)
-HT = np.dot(ATP,IS2)
+HT = np.dot(ATP, IS2)
-HR = np.dot(ARP,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)])
+IS1e = np.array(
-IS2e = np.array([IinPa+DiC*QinPae,DiC*IinPa+QinPae,ZiC*(CosC*UinPae+SinC*VinPa),-ZiC*(SinC*UinPae-CosC*VinPa)])
+[IinPo + DiC * QinPoe, DiC * IinPo + QinPoe, ZiC * (CosC * UinPoe + SinC * VinPo),
-GT = np.dot(ATPe,IS1e)
+-ZiC * (SinC * UinPoe - CosC * VinPo)])
-GR = np.dot(ARPe,IS1e)
+IS2e = np.array(
-HT = np.dot(ATPe,IS2e)
+[IinPa + DiC * QinPae, DiC * IinPa + QinPae, ZiC * (CosC * UinPae + SinC * VinPa),
-HR = np.dot(ARPe,IS2e)
+-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))
+# S2ge = np.sin(np.deg2rad(2*RotO - 2*RotC))
-#C2ge = np.cos(np.deg2rad(2*RotO - 2*RotC))
+# C2ge = np.cos(np.deg2rad(2*RotO - 2*RotC))
-S2e = np.sin(np.deg2rad(2*RotC))
+S2e = np.sin(np.deg2rad(2 * RotC))
-C2e = np.cos(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])
+AF1 = np.array([1, C2g * DiO, S2g * DiO, 0])
-AF2 = np.array([C2g*DiO,1-S2g**2*WiO,S2g*C2g*WiO,-S2g*ZiO*SinO])
+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])
+AF3 = np.array([S2g * DiO, S2g * C2g * WiO, 1 - C2g ** 2 * WiO, C2g * ZiO * SinO])
-AF4 = np.array([0, S2g*SinO, -C2g*SinO, CosO])
+AF4 = np.array([0, S2g * SinO, -C2g * SinO, CosO])
-ATF = (ATP1*AF1+ATP2*AF2+ATP3*AF3+ATP4*AF4)
+ATF = (ATP1 * AF1 + ATP2 * AF2 + ATP3 * AF3 + ATP4 * AF4)
-ARF = (ARP1*AF1+ARP2*AF2+ARP3*AF3+ARP4*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]
+ATF2e = C2e * ATF[1] + S2e * ATF[2]
-ATF3e = C2e*ATF[2] - S2e*ATF[1]
+ATF3e = C2e * ATF[2] - S2e * ATF[1]
-ARF2e = C2e*ARF[1] + S2e*ARF[2]
+ARF2e = C2e * ARF[1] + S2e * ARF[2]
-ARF3e = C2e*ARF[2] - S2e*ARF[1]
+ARF3e = C2e * ARF[2] - S2e * ARF[1]
-ATFe = np.array([ATF1,ATF2e,ATF3e,ATF4])
+ATFe = np.array([ATF1, ATF2e, ATF3e, ATF4])
-ARFe = np.array([ARF1,ARF2e,ARF3e,ARF4])
+ARFe = np.array([ARF1, ARF2e, ARF3e, ARF4])
-QinEe = C2e*QinE + S2e*UinE
+QinEe = C2e * QinE + S2e * UinE
-UinEe = C2e*UinE - S2e*QinE
+UinEe = C2e * UinE - S2e * QinE
 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F)
 IinF = IinE
-QinF = aCal*QinE
+QinF = aCal * QinE
-UinF = -aCal*UinE
+UinF = -aCal * UinE
-VinF = (1.-2.*aCal)*VinE
+VinF = (1. - 2. * aCal) * VinE
 IinFo = IinE
 QinFo = 0.
 UinFo = 0.
 VinFo = VinE
 IinFa = 0.
 QinFa = QinE
 UinFa = -UinE
-VinFa = -2.*VinE
+VinFa = -2. * VinE
 # Stokes Input Vector before receiver optics rotated by epsilon Eq. C.3
-QinFe = C2e*QinF + S2e*UinF
+QinFe = C2e * QinF + S2e * UinF
-UinFe = C2e*UinF - S2e*QinF
+UinFe = C2e * UinF - S2e * QinF
-QinFoe = C2e*QinFo + S2e*UinFo
+QinFoe = C2e * QinFo + S2e * UinFo
-UinFoe = C2e*UinFo - S2e*QinFo
+UinFoe = C2e * UinFo - S2e * QinFo
-QinFae = C2e*QinFa + S2e*UinFa
+QinFae = C2e * QinFa + S2e * UinFa
-UinFae = C2e*UinFa - S2e*QinFa
+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
+if (TypeC == 2) or (TypeC == 1):  # rotator calibration Eq. C.4
-AT = ATF1*IinF + ATF4*h*VinF
+AT = ATF1 * IinF + ATF4 * h * VinF
-BT = ATF3e*QinF - ATF2e*h*UinF
+BT = ATF3e * QinF - ATF2e * h * UinF
-AR = ARF1*IinF + ARF4*h*VinF
+AR = ARF1 * IinF + ARF4 * h * VinF
-BR = ARF3e*QinF - ARF2e*h*UinF
+BR = ARF3e * QinF - ARF2e * h * UinF
 # Correction paremeters for normal measurements; they are independent of LDR
 if (not RotationErrorEpsilonForNormalMeasurements):
-GT = ATF1*IinE + ATF4*VinE
+GT = ATF1 * IinE + ATF4 * VinE
-GR = ARF1*IinE + ARF4*VinE
+GR = ARF1 * IinE + ARF4 * VinE
-HT = ATF2*QinE - ATF3*UinE - ATF4*2*VinE
+HT = ATF2 * QinE - ATF3 * UinE - ATF4 * 2 * VinE
-HR = ARF2*QinE - ARF3*UinE - ARF4*2*VinE
+HR = ARF2 * QinE - ARF3 * UinE - ARF4 * 2 * VinE
 else:
-GT = ATF1*IinE + ATF4*h*VinE
+GT = ATF1 * IinE + ATF4 * h * VinE
-GR = ARF1*IinE + ARF4*h*VinE
+GR = ARF1 * IinE + ARF4 * h * VinE
-HT = ATF2e*QinE - ATF3e*h*UinE - ATF4*h*2*VinE
+HT = ATF2e * QinE - ATF3e * h * UinE - ATF4 * h * 2 * VinE
-HR = ARF2e*QinE - ARF3e*h*UinE - ARF4*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])
+IF1e = np.array([IinF, ZiC * CosC * QinFe, UinFe, ZiC * CosC * VinF])
-IF2e = np.array([DiC*UinFe, -ZiC*SinC*VinF, DiC*IinF, ZiC*SinC*QinFe])
+IF2e = np.array([DiC * UinFe, -ZiC * SinC * VinF, DiC * IinF, ZiC * SinC * QinFe])
-AT = np.dot(ATFe,IF1e)
+AT = np.dot(ATFe, IF1e)
-AR = np.dot(ARFe,IF1e)
+AR = np.dot(ARFe, IF1e)
-BT = np.dot(ATFe,IF2e)
+BT = np.dot(ATFe, IF2e)
-BR = np.dot(ARFe,IF2e)
+BR = np.dot(ARFe, IF2e)
 # Correction paremeters for normal measurements; they are independent of LDR  --- the same as for TypeC = 6
-if (not RotationErrorEpsilonForNormalMeasurements):   # calibrator taken out
+if (not RotationErrorEpsilonForNormalMeasurements):  # calibrator taken out
-IS1 = np.array([IinE,0,0,VinE])
+IS1 = np.array([IinE, 0, 0, VinE])
-IS2 = np.array([0,QinE,-UinE,-2*VinE])
+IS2 = np.array([0, QinE, -UinE, -2 * VinE])
-GT = np.dot(ATF,IS1)
+GT = np.dot(ATF, IS1)
-GR = np.dot(ARF,IS1)
+GR = np.dot(ARF, IS1)
-HT = np.dot(ATF,IS2)
+HT = np.dot(ATF, IS2)
-HR = np.dot(ARF,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)])
+IS1e = np.array(
-IS2e = np.array([IinFa+DiC*QinFae,DiC*IinFa+QinFae,ZiC*(CosC*UinFae+SinC*VinFa),-ZiC*(SinC*UinFae-CosC*VinFa)])
+[IinFo + DiC * QinFoe, DiC * IinFo + QinFoe, ZiC * (CosC * UinFoe + SinC * VinFo),
-GT = np.dot(ATFe,IS1e)
+-ZiC * (SinC * UinFoe - CosC * VinFo)])
-GR = np.dot(ARFe,IS1e)
+IS2e = np.array(
-HT = np.dot(ATFe,IS2e)
+[IinFa + DiC * QinFae, DiC * IinFa + QinFae, ZiC * (CosC * UinFae + SinC * VinFa),
-HR = np.dot(ARFe,IS2e)
+-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])
+IF1e = np.array([IinF + sqr05 * DiC * QinFe, sqr05 * DiC * IinF + (1 - 0.5 * WiC) * QinFe,
-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])
+(1 - 0.5 * WiC) * UinFe + sqr05 * ZiC * SinC * VinF,
+-sqr05 * ZiC * SinC * UinFe + ZiC * CosC * VinF])
-AT = np.dot(ATFe,IF1e)
+IF2e = np.array([sqr05 * DiC * UinFe, 0.5 * WiC * UinFe - sqr05 * ZiC * SinC * VinF,
-AR = np.dot(ARFe,IF1e)
+sqr05 * DiC * IinF + 0.5 * WiC * QinFe, sqr05 * ZiC * SinC * QinFe])
-BT = np.dot(ATFe,IF2e)
-BR = np.dot(ARFe,IF2e)
+AT = np.dot(ATFe, IF1e)
+AR = np.dot(ARFe, IF1e)
+BT = np.dot(ATFe, IF2e)
+BR = np.dot(ARFe, IF2e)
 # Correction paremeters for normal measurements; they are independent of LDR
-if (not RotationErrorEpsilonForNormalMeasurements):   # calibrator taken out
+if (not RotationErrorEpsilonForNormalMeasurements):  # calibrator taken out
-#IS1 = np.array([IinE,0,0,VinE])
+# IS1 = np.array([IinE,0,0,VinE])
-#IS2 = np.array([0,QinE,-UinE,-2*VinE])
+# IS2 = np.array([0,QinE,-UinE,-2*VinE])
-IS1 = np.array([IinFo,0,0,VinFo])
+IS1 = np.array([IinFo, 0, 0, VinFo])
-IS2 = np.array([0,QinFa,UinFa,VinFa])
+IS2 = np.array([0, QinFa, UinFa, VinFa])
-GT = np.dot(ATF,IS1)
+GT = np.dot(ATF, IS1)
-GR = np.dot(ARF,IS1)
+GR = np.dot(ARF, IS1)
-HT = np.dot(ATF,IS2)
+HT = np.dot(ATF, IS2)
-HR = np.dot(ARF,IS2)
+HR = np.dot(ARF, IS2)
 else:
-#IS1e = np.array([IinE,DiC*IinE,ZiC*SinC*VinE,ZiC*CosC*VinE])
+# 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)])
+# 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)])
+IS1e = np.array(
-IS2e = np.array([IinFa+DiC*QinFae,DiC*IinFa+QinFae,ZiC*(CosC*UinFae+SinC*VinFa),-ZiC*(SinC*UinFae-CosC*VinFa)])
+[IinFo + DiC * QinFoe, DiC * IinFo + QinFoe, ZiC * (CosC * UinFoe + SinC * VinFo),
-GT = np.dot(ATFe,IS1e)
+-ZiC * (SinC * UinFoe - CosC * VinFo)])
-GR = np.dot(ARFe,IS1e)
+IS2e = np.array(
-HT = np.dot(ATFe,IS2e)
+[IinFa + DiC * QinFae, DiC * IinFa + QinFae, ZiC * (CosC * UinFae + SinC * VinFa),
-HR = np.dot(ARFe,IS2e)
+-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")
+# print("Calibrator location not implemented yet")
-S2e = np.sin(np.deg2rad(2*RotC))
+S2e = np.sin(np.deg2rad(2 * RotC))
-C2e = np.cos(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])
+AF1 = np.array([1, C2g * DiO, S2g * DiO, 0])
-AF2 = np.array([C2g*DiO,1-S2g**2*WiO,S2g*C2g*WiO,-S2g*ZiO*SinO])
+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])
+AF3 = np.array([S2g * DiO, S2g * C2g * WiO, 1 - C2g ** 2 * WiO, C2g * ZiO * SinO])
-AF4 = np.array([0, S2g*SinO, -C2g*SinO, CosO])
+AF4 = np.array([0, S2g * SinO, -C2g * SinO, CosO])
-ATF = (ATP1*AF1+ATP2*AF2+ATP3*AF3+ATP4*AF4)
+ATF = (ATP1 * AF1 + ATP2 * AF2 + ATP3 * AF3 + ATP4 * AF4)
-ARF = (ARP1*AF1+ARP2*AF2+ARP3*AF3+ARP4*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]
 ATE1o, ARE1o = ATF1, ARF1
 ATE2o, ARE2o = 0., 0.
 ATE3o, ARE3o = 0., 0.
 ATE4o, ARE4o = ATF4, ARF4
 # terms with aCal
-ATE1a, ARE1a = 0. , 0.
+ATE1a, ARE1a = 0., 0.
 ATE2a, ARE2a = ATF2, ARF2
 ATE3a, ARE3a = -ATF3, -ARF3
-ATE4a, ARE4a = -2*ATF4, -2*ARF4
+ATE4a, ARE4a = -2 * ATF4, -2 * ARF4
 # rotated AinEa by epsilon
-ATE2ae =  C2e*ATF2 + S2e*ATF3
+ATE2ae = C2e * ATF2 + S2e * ATF3
-ATE3ae = -S2e*ATF2 - C2e*ATF3
+ATE3ae = -S2e * ATF2 - C2e * ATF3
-ARE2ae =  C2e*ARF2 + S2e*ARF3
+ARE2ae = C2e * ARF2 + S2e * ARF3
-ARE3ae = -S2e*ARF2 - C2e*ARF3
+ARE3ae = -S2e * ARF2 - C2e * ARF3
 ATE1 = ATE1o
-ATE2e = aCal*ATE2ae
+ATE2e = aCal * ATE2ae
-ATE3e = aCal*ATE3ae
+ATE3e = aCal * ATE3ae
-ATE4 = (1-2*aCal)*ATF4
+ATE4 = (1 - 2 * aCal) * ATF4
 ARE1 = ARE1o
-ARE2e = aCal*ARE2ae
+ARE2e = aCal * ARE2ae
-ARE3e = aCal*ARE3ae
+ARE3e = aCal * ARE3ae
-ARE4 = (1-2*aCal)*ARF4
+ARE4 = (1 - 2 * aCal) * ARF4
 # rotated IinE
-QinEe = C2e*QinE + S2e*UinE
+QinEe = C2e * QinE + S2e * UinE
-UinEe = C2e*UinE - S2e*QinE
+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
+if (TypeC == 2) or (TypeC == 1):  # +++++++++ rotator calibration Eq. C.4
-AT = ATE1o*IinE + (ATE4o+aCal*ATE4a)*h*VinE
+AT = ATE1o * IinE + (ATE4o + aCal * ATE4a) * h * VinE
-BT = aCal * (ATE3ae*QinEe - ATE2ae*h*UinEe)
+BT = aCal * (ATE3ae * QinEe - ATE2ae * h * UinEe)
-AR = ARE1o*IinE + (ARE4o+aCal*ARE4a)*h*VinE
+AR = ARE1o * IinE + (ARE4o + aCal * ARE4a) * h * VinE
-BR = aCal * (ARE3ae*QinEe - ARE2ae*h*UinEe)
+BR = aCal * (ARE3ae * QinEe - ARE2ae * h * UinEe)
 # Correction paremeters for normal measurements; they are independent of LDR
 if (not RotationErrorEpsilonForNormalMeasurements):
 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F)
-GT = ATE1o*IinE + ATE4o*h*VinE
+GT = ATE1o * IinE + ATE4o * h * VinE
-GR = ARE1o*IinE + ARE4o*h*VinE
+GR = ARE1o * IinE + ARE4o * h * VinE
-HT = ATE2a*QinE + ATE3a*h*UinEe + ATE4a*h*VinE
+HT = ATE2a * QinE + ATE3a * h * UinEe + ATE4a * h * VinE
-HR = ARE2a*QinE + ARE3a*h*UinEe + ARE4a*h*VinE
+HR = ARE2a * QinE + ARE3a * h * UinEe + ARE4a * h * VinE
 else:
-GT = ATE1o*IinE + ATE4o*h*VinE
+GT = ATE1o * IinE + ATE4o * h * VinE
-GR = ARE1o*IinE + ARE4o*h*VinE
+GR = ARE1o * IinE + ARE4o * h * VinE
-HT = ATE2ae*QinE + ATE3ae*h*UinEe + ATE4a*h*VinE
+HT = ATE2ae * QinE + ATE3ae * h * UinEe + ATE4a * h * VinE
-HR = ARE2ae*QinE + ARE3ae*h*UinEe + ARE4a*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
+AT = ATE1 * IinE + ZiC * CosC * (ATE2e * QinEe + ATE4 * VinE) + ATE3e * UinEe
-BT = DiC*(ATE1*UinEe + ATE3e*IinE) + ZiC*SinC*(ATE4*QinEe - ATE2e*VinE)
+BT = DiC * (ATE1 * UinEe + ATE3e * IinE) + ZiC * SinC * (ATE4 * QinEe - ATE2e * VinE)
-AR = ARE1*IinE + ZiC*CosC*(ARE2e*QinEe + ARE4*VinE) + ARE3e*UinEe
+AR = ARE1 * IinE + ZiC * CosC * (ARE2e * QinEe + ARE4 * VinE) + ARE3e * UinEe
-BR = DiC*(ARE1*UinEe + ARE3e*IinE) + ZiC*SinC*(ARE4*QinEe - ARE2e*VinE)
+BR = DiC * (ARE1 * UinEe + ARE3e * IinE) + ZiC * SinC * (ARE4 * QinEe - ARE2e * VinE)
 # Correction paremeters for normal measurements; they are independent of LDR
 if (not RotationErrorEpsilonForNormalMeasurements):
 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F)
-GT = ATE1o*IinE + ATE4o*VinE
+GT = ATE1o * IinE + ATE4o * VinE
-GR = ARE1o*IinE + ARE4o*VinE
+GR = ARE1o * IinE + ARE4o * VinE
-HT = ATE2a*QinE + ATE3a*UinE + ATE4a*VinE
+HT = ATE2a * QinE + ATE3a * UinE + ATE4a * VinE
-HR = ARE2a*QinE + ARE3a*UinE + ARE4a*VinE
+HR = ARE2a * QinE + ARE3a * UinE + ARE4a * VinE
 else:
-D = IinE + DiC*QinEe
+D = IinE + DiC * QinEe
-A = DiC*IinE + QinEe
+A = DiC * IinE + QinEe
-B = ZiC*(CosC*UinEe + SinC*VinE)
+B = ZiC * (CosC * UinEe + SinC * VinE)
-C = -ZiC*(SinC*UinEe - CosC*VinE)
+C = -ZiC * (SinC * UinEe - CosC * VinE)
-GT = ATE1o*D + ATE4o*C
+GT = ATE1o * D + ATE4o * C
-GR = ARE1o*D + ARE4o*C
+GR = ARE1o * D + ARE4o * C
-HT = ATE2a*A + ATE3a*B + ATE4a*C
+HT = ATE2a * A + ATE3a * B + ATE4a * C
-HR = ARE2a*A + ARE3a*B + ARE4a*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])
+IE1e = np.array([IinE + sqr05 * DiC * QinEe, sqr05 * DiC * IinE + (1 - 0.5 * WiC) * QinEe,
-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])
+(1 - 0.5 * WiC) * UinEe + sqr05 * ZiC * SinC * VinE,
-ATEe = np.array([ATE1,ATE2e,ATE3e,ATE4])
+-sqr05 * ZiC * SinC * UinEe + ZiC * CosC * VinE])
-AREe = np.array([ARE1,ARE2e,ARE3e,ARE4])
+IE2e = np.array([sqr05 * DiC * UinEe, 0.5 * WiC * UinEe - sqr05 * ZiC * SinC * VinE,
-AT = np.dot(ATEe,IE1e)
+sqr05 * DiC * IinE + 0.5 * WiC * QinEe, sqr05 * ZiC * SinC * QinEe])
-AR = np.dot(AREe,IE1e)
+ATEe = np.array([ATE1, ATE2e, ATE3e, ATE4])
-BT = np.dot(ATEe,IE2e)
+AREe = np.array([ARE1, ARE2e, ARE3e, ARE4])
-BR = np.dot(AREe,IE2e)
+AT = np.dot(ATEe, IE1e)
+AR = np.dot(AREe, IE1e)
+BT = np.dot(ATEe, IE2e)
+BR = np.dot(AREe, IE2e)
 # Correction paremeters for normal measurements; they are independent of LDR
-if (not RotationErrorEpsilonForNormalMeasurements):   # calibrator taken out
+if (not RotationErrorEpsilonForNormalMeasurements):  # calibrator taken out
-GT = ATE1o*IinE + ATE4o*VinE
+GT = ATE1o * IinE + ATE4o * VinE
-GR = ARE1o*IinE + ARE4o*VinE
+GR = ARE1o * IinE + ARE4o * VinE
-HT = ATE2a*QinE + ATE3a*UinE + ATE4a*VinE
+HT = ATE2a * QinE + ATE3a * UinE + ATE4a * VinE
-HR = ARE2a*QinE + ARE3a*UinE + ARE4a*VinE
+HR = ARE2a * QinE + ARE3a * UinE + ARE4a * VinE
 else:
-D = IinE + DiC*QinEe
+D = IinE + DiC * QinEe
-A = DiC*IinE + QinEe
+A = DiC * IinE + QinEe
-B = ZiC*(CosC*UinEe + SinC*VinE)
+B = ZiC * (CosC * UinEe + SinC * VinE)
-C = -ZiC*(SinC*UinEe - CosC*VinE)
+C = -ZiC * (SinC * UinEe - CosC * VinE)
-GT = ATE1o*D + ATE4o*C
+GT = ATE1o * D + ATE4o * C
-GR = ARE1o*D + ARE4o*C
+GR = ARE1o * D + ARE4o * C
-HT = ATE2a*A + ATE3a*B + ATE4a*C
+HT = ATE2a * A + ATE3a * B + ATE4a * C
-HR = ARE2a*A + ARE3a*B + ARE4a*C
+HR = ARE2a * A + ARE3a * B + ARE4a * C
 else:
 print('Calibrator not implemented yet')
 sys.exit()
 # Calibration signals with aCal => Determination of the correction K of the real calibration factor
-IoutTp = TaT*TiT*TiO*TiE*(AT + BT)
+IoutTp = TaT * TiT * TiO * TiE * (AT + BT)
-IoutTm = TaT*TiT*TiO*TiE*(AT - BT)
+IoutTm = TaT * TiT * TiO * TiE * (AT - BT)
-IoutRp = TaR*TiR*TiO*TiE*(AR + BR)
+IoutRp = TaR * TiR * TiO * TiE * (AR + BR)
-IoutRm = TaR*TiR*TiO*TiE*(AR - BR)
+IoutRm = TaR * TiR * TiO * TiE * (AR - BR)
 # --- Results and Corrections; electronic etaR and etaT are assumed to be 1
-#Eta = TiR/TiT   # Eta = Eta*/K  Eq. 84
+# Eta = TiR/TiT   # Eta = Eta*/K  Eq. 84
-Etapx = IoutRp/IoutTp
+Etapx = IoutRp / IoutTp
-Etamx = IoutRm/IoutTm
+Etamx = IoutRm / IoutTm
-Etax = (Etapx*Etamx)**0.5
+Etax = (Etapx * Etamx) ** 0.5
 K = Etax / Eta0
-#print("{0:6.3f},{1:6.3f},{2:6.3f},{3:6.3f},{4:6.3f},{5:6.3f},{6:6.3f},{7:6.3f},{8:6.3f},{9:6.3f},{10:6.3f}".format(AT, BT, AR, BR, DiC, ZiC, RetO, TP, TS, Kp, Km))
+# print("{0:6.3f},{1:6.3f},{2:6.3f},{3:6.3f},{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))
+# 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_p = (IoutRp/IoutTp)
-#Eta_test_m = (IoutRm/IoutTm)
+# Eta_test_m = (IoutRm/IoutTm)
-#Eta_test = (Eta_test_p*Eta_test_m)**0.5
+# Eta_test = (Eta_test_p*Eta_test_m)**0.5
 # *************************************************************************
 iLDR = -1
 for LDRTrue in LDRrange:
 iLDR = iLDR + 1
-atrue = (1-LDRTrue)/(1+LDRTrue)
+atrue = (1 - LDRTrue) / (1 + LDRTrue)
 # ----- Forward simulated signals and LDRsim with atrue; from input file
-It = TaT*TiT*TiO*TiE*(GT+atrue*HT) #  TaT*TiT*TiC*TiO*IinL*(GT+atrue*HT)
+It = TaT * TiT * TiO * TiE * (GT + atrue * HT)  # TaT*TiT*TiC*TiO*IinL*(GT+atrue*HT)
-Ir = TaR*TiR*TiO*TiE*(GR+atrue*HR) #  TaR*TiR*TiC*TiO*IinL*(GR+atrue*HR)
+Ir = TaR * TiR * TiO * TiE * (GR + atrue * HR)  # TaR*TiR*TiC*TiO*IinL*(GR+atrue*HR)
 # LDRsim = 1/Eta*Ir/It  # simulated LDR* with Y from input file
-LDRsim = Ir/It  # simulated uncorrected LDR with Y from input file
+LDRsim = Ir / It  # simulated uncorrected LDR with Y from input file
 '''
 if Y == 1.:
 LDRsimx = LDRsim
 LDRsimx2 = LDRsim2
 else:
 LDRsimx = 1./LDRsim
 LDRsimx2 = 1./LDRsim2
 '''
 # ----- Backward correction
 # Corrected LDRCorr from forward simulated LDRsim (atrue) with assumed true G0,H0,K0
-LDRCorr = (LDRsim*K0/Etax*(GT0+HT0)-(GR0+HR0))/((GR0-HR0)-LDRsim*K0/Etax*(GT0-HT0))
+LDRCorr = (LDRsim * K0 / Etax * (GT0 + HT0) - (GR0 + HR0)) / (
+(GR0 - HR0) - LDRsim * K0 / Etax * (GT0 - HT0))
 # -- F11corr from It and Ir and calibration EtaX
 Text1 = "!!! EXPERIMENTAL !!!  F11corr from It and Ir with calibration EtaX: x-axis: F11corr(LDRtrue) / F11corr(LDRtrue = 0.004) - 1"
-F11corr = 1/(TiO*TiE)*((HR0*Etax/K0*It/TTa-HT0*Ir/TRa)/(HR0*GT0-HT0*GR0))    # IL = 1  Eq.(64)
+F11corr = 1 / (TiO * TiE) * (
+(HR0 * Etax / K0 * It / TTa - HT0 * Ir / TRa) / (HR0 * GT0 - HT0 * GR0))  # IL = 1  Eq.(64)
-#Text1 = "F11corr from It and Ir without corrections but with calibration EtaX: x-axis: F11corr(LDRtrue) devided by F11corr(LDRtrue = 0.004)"
-#F11corr = 0.5/(TiO*TiE)*(Etax*It/TTa+Ir/TRa)    # IL = 1  Eq.(64)
+# 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"
+# Text1 = " x-axis: IT(LDRtrue) / IT(LDRtrue = 0.004) - 1"
-#F11corr = It/(TaT*TiT*TiO*TiE)   #/(TaT*TiT*TiO*TiE*(GT0+atrue*HT0))
+# F11corr = It/(TaT*TiT*TiO*TiE)   #/(TaT*TiT*TiO*TiE*(GT0+atrue*HT0))
 # !!! see below line 1673ff
-aF11corr[iLDR,iN] = F11corr
+aF11corr[iLDR, iN] = F11corr
-aA[iLDR,iN] = LDRCorr
+aA[iLDR, iN] = LDRCorr
-aX[0,iN] = GR
+aX[0, iN] = GR
-aX[1,iN] = GT
+aX[1, iN] = GT
-aX[2,iN] = HR
+aX[2, iN] = HR
-aX[3,iN] = HT
+aX[3, iN] = HT
-aX[4,iN] = K
+aX[4, iN] = K
 aLDRCal[iN] = iLDRCal
 aERaT[iN] = iERaT
 aERaR[iN] = iERaR
 aRotaT[iN] = iRotaT
 aRP[iN] = iRP
 aRS[iN] = iRS
 # --- END loop
 btime = clock()
-print("\r done in      ", "{0:5.0f}".format(btime-atime), "sec") #, end="\r")
+print("\r done in      ", "{0:5.0f}".format(btime - atime), "sec")  # , end="\r")
 # --- Plot -----------------------------------------------------------------
 if (sns_loaded):
 sns.set_style("whitegrid")
 sns.set_palette("bright", 6)
 plt.plot(aA[3,:],'r.')
 plt.plot(aA[4,:],'g.')
 #plt.plot(aA[6,:],'c.')
 plt.show
 '''
 # Plot LDR
 def PlotSubHist(aVar, aX, X0, daX, iaX, naX):
 fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2))
 iLDR = -1
 for LDRTrue in LDRrange:
 iLDR = iLDR + 1
-LDRmin[iLDR] = np.min(aA[iLDR,:])
+LDRmin[iLDR] = np.min(aA[iLDR, :])
-LDRmax[iLDR] = np.max(aA[iLDR,:])
+LDRmax[iLDR] = np.max(aA[iLDR, :])
-Rmin = LDRmin[iLDR] * 0.995 #  np.min(aA[iLDR,:])    * 0.995
+Rmin = LDRmin[iLDR] * 0.995  # np.min(aA[iLDR,:])    * 0.995
-Rmax = LDRmax[iLDR] * 1.005 #  np.max(aA[iLDR,:])    * 1.005
+Rmax = LDRmax[iLDR] * 1.005  # np.max(aA[iLDR,:])    * 1.005
-#plt.subplot(5,2,iLDR+1)
+# plt.subplot(5,2,iLDR+1)
-plt.subplot(1,5,iLDR+1)
+plt.subplot(1, 5, iLDR + 1)
-(n, bins, patches) = plt.hist(aA[iLDR,:],
+(n, bins, patches) = plt.hist(aA[iLDR, :],
 bins=100, log=False,
 range=[Rmin, Rmax],
-alpha=0.5, normed=False, color = '0.5', histtype='stepfilled')
+alpha=0.5, normed=False, color='0.5', histtype='stepfilled')
-for iaX in range(-naX,naX+1):
+for iaX in range(-naX, naX + 1):
-plt.hist(aA[iLDR,aX == iaX],
+plt.hist(aA[iLDR, aX == iaX],
 range=[Rmin, Rmax],
-bins=100, log=False, alpha=0.3, normed=False, histtype='stepfilled', label = str(round(X0 + iaX*daX/naX,5)))
+bins=100, log=False, alpha=0.3, normed=False, histtype='stepfilled',
+label=str(round(X0 + iaX * daX / naX, 5)))
 if (iLDR == 2): plt.legend()
 plt.tick_params(axis='both', labelsize=9)
 plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2)
-#plt.title(LID + '  ' + aVar, fontsize=18)
+# plt.title(LID + '  ' + aVar, fontsize=18)
-#plt.ylabel('frequency', fontsize=10)
+# plt.ylabel('frequency', fontsize=10)
-#plt.xlabel('LDRcorr', fontsize=10)
+# plt.xlabel('LDRcorr', fontsize=10)
-#fig.tight_layout()
+# fig.tight_layout()
-fig.suptitle(LID + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC])  + ' - ' + aVar, fontsize=14, y=1.05)
+fig.suptitle(LID + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar, fontsize=14, y=1.05)
-#plt.show()
+# plt.show()
-#fig.savefig(LID + '_' + aVar + '.png', dpi=150, bbox_inches='tight', pad_inches=0)
+# fig.savefig(LID + '_' + aVar + '.png', dpi=150, bbox_inches='tight', pad_inches=0)
-#plt.close
+# plt.close
 return
 if (nRotL > 0): PlotSubHist("RotL", aRotL, RotL0, dRotL, iRotL, nRotL)
 if (nRetE > 0): PlotSubHist("RetE", aRetE, RetE0, dRetE, iRetE, nRetE)
 if (nRotE > 0): PlotSubHist("RotE", aRotE, RotE0, dRotE, iRotE, nRotE)
 if (nDiE > 0): PlotSubHist("DiE", aDiE, DiE0, dDiE, iDiE, nDiE)
 plt.close
 '''
 # --- Plot LDRmin, LDRmax
 fig2 = plt.figure()
-plt.plot(LDRrange,LDRmax-LDRrange, linewidth=2.0, color='b')
+plt.plot(LDRrange, LDRmax - LDRrange, linewidth=2.0, color='b')
-plt.plot(LDRrange,LDRmin-LDRrange, linewidth=2.0, color='g')
+plt.plot(LDRrange, LDRmin - LDRrange, linewidth=2.0, color='g')
 plt.xlabel('LDRtrue', fontsize=18)
 plt.ylabel('LDRTrue-LDRmin, LDRTrue-LDRmax', fontsize=14)
 plt.title(LID + ' ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]), fontsize=18)
-#plt.ylimit(-0.07, 0.07)
+# plt.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

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