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