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@7: Python 3.4.2 ,, ulalume3@0: """ ulalume3@0: #!/usr/bin/env python3 ulalume3@0: from __future__ import print_function ulalume3@0: #import math ulalume3@0: import numpy as np ulalume3@0: import sys ulalume3@0: import os ulalume3@0: ulalume3@0: #import seaborn as sns ulalume3@0: import matplotlib.pyplot as plt ulalume3@0: from time import clock ulalume3@0: ulalume3@0: #from matplotlib.backends.backend_pdf import PdfPages ulalume3@0: #pdffile = '{}.pdf'.format('path') ulalume3@0: #pp = PdfPages(pdffile) ulalume3@0: ## pp.savefig can be called multiple times to save to multiple pages ulalume3@0: #pp.savefig() ulalume3@0: #pp.close() ulalume3@0: ulalume3@0: from contextlib import contextmanager ulalume3@0: @contextmanager ulalume3@0: def redirect_stdout(new_target): ulalume3@0: old_target, sys.stdout = sys.stdout, new_target # replace sys.stdout ulalume3@0: try: ulalume3@0: yield new_target # run some code with the replaced stdout ulalume3@0: finally: ulalume3@0: sys.stdout.flush() ulalume3@0: sys.stdout = old_target # restore to the previous value ulalume3@0: ''' ulalume3@0: real_raw_input = vars(__builtins__).get('raw_input',input) ulalume3@0: ''' ulalume3@0: try: ulalume3@0: import __builtin__ ulalume3@0: input = getattr(__builtin__, 'raw_input') ulalume3@0: except (ImportError, AttributeError): ulalume3@0: pass ulalume3@0: ulalume3@0: from distutils.util import strtobool 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: ulalume3@0: #if user_yes_no_query('want to exit?') == 1: sys.exit() ulalume3@0: ulalume3@0: ''' ulalume3@0: ## {{{ http://code.activestate.com/recipes/577058/ (r2) ulalume3@0: def query_yes_no(question, default="yes"): ulalume3@0: valid = {"yes":"yes", "y":"yes", "ye":"yes", ulalume3@0: "no":"no", "n":"no"} ulalume3@0: if default == None: ulalume3@0: prompt = " [y/n] " ulalume3@0: elif default == "yes": ulalume3@0: prompt = " [Y/n] " ulalume3@0: elif default == "no": ulalume3@0: prompt = " [y/N] " ulalume3@0: else: ulalume3@0: raise ValueError("invalid default answer: '%s'" % default) ulalume3@0: ulalume3@0: while 1: ulalume3@0: sys.stdout.write(question + prompt) ulalume3@0: choice = input().lower() ulalume3@0: if default is not None and choice == '': ulalume3@0: return default ulalume3@0: elif choice in valid.keys(): ulalume3@0: return valid[choice] ulalume3@0: else: ulalume3@0: sys.stdout.write("Please respond with 'yes' or 'no' "\ ulalume3@0: "(or 'y' or 'n').\n") ulalume3@0: ## end of http://code.activestate.com/recipes/577058/ }}} 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: ulalume3@0: #PrintToOutputFile = True ulalume3@0: ulalume3@0: sqr05 = 0.5**0.5 ulalume3@0: 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 ulalume3@0: bL = 1. #degree of linear polarization; default 1 ulalume3@0: RotL, dRotL, nRotL = 0.0, 0.0, 1 #alpha; rotation of laser polarization in degrees; default 0 ulalume3@0: # --- ME Emitter and +-Uncertainty ulalume3@0: DiE, dDiE, nDiE = 0., 0.00, 1 # Diattenuation ulalume3@0: TiE = 1. # Unpolarized transmittance ulalume3@0: RetE, dRetE, nRetE = 0., 180.0, 0 # Retardance in degrees ulalume3@0: RotE, dRotE, nRotE = 0., 0.0, 0 # beta: Rotation of optical element in degrees ulalume3@0: # --- MO Receiver Optics including telescope ulalume3@0: DiO, dDiO, nDiO = -0.055, 0.003, 1 ulalume3@0: TiO = 0.9 ulalume3@0: RetO, dRetO, nRetO = 0., 180.0, 2 ulalume3@0: RotO, dRotO, nRotO = 0., 0.1, 1 #gamma ulalume3@0: # --- PBS MT transmitting path defined with (TS,TP); and +-Uncertainty ulalume3@0: TP, dTP, nTP = 0.98, 0.02, 1 ulalume3@0: TS, dTS, nTS = 0.001, 0.001, 1 ulalume3@0: TiT = 0.5 * (TP + TS) ulalume3@0: DiT = (TP-TS)/(TP+TS) ulalume3@0: # PolFilter ulalume3@0: RetT, dRetT, nRetT = 0., 180., 0 ulalume3@0: ERaT, dERaT, nERaT = 0.001, 0.001, 1 ulalume3@0: RotaT, dRotaT, nRotaT = 0., 3., 1 ulalume3@0: DaT = (1-ERaT)/(1+ERaT) ulalume3@0: TaT = 0.5*(1+ERaT) ulalume3@0: # --- PBS MR reflecting path defined with (RS,RP); and +-Uncertainty ulalume3@0: RS, dRS, nRS = 1 - TS, 0., 0 ulalume3@0: RP, dRP, nRP = 1 - TP, 0., 0 ulalume3@0: TiR = 0.5 * (RP + RS) ulalume3@0: DiR = (RP-RS)/(RP+RS) ulalume3@0: # PolFilter ulalume3@0: RetR, dRetR, nRetR = 0., 180., 0 ulalume3@0: ERaR, dERaR, nERaR = 0.001, 0.001, 1 ulalume3@0: RotaR,dRotaR,nRotaR = 90., 3., 1 ulalume3@0: DaR = (1-ERaR)/(1+ERaR) ulalume3@0: 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 ulalume3@0: LocC = 4 # location of calibrator: behind laser = 1; behind emitter = 2; before receiver = 3; before PBS = 4 ulalume3@0: ulalume3@0: TypeC = 3 # linear polarizer calibrator ulalume3@0: # example with extinction ratio 0.001 ulalume3@0: DiC, dDiC, nDiC = 1.0, 0., 0 # ideal 1.0 ulalume3@0: TiC = 0.5 # ideal 0.5 ulalume3@0: RetC, dRetC, nRetC = 0., 0., 0 ulalume3@0: RotC, dRotC, nRotC = 0.0, 0.1, 0 #constant calibrator offset epsilon ulalume3@0: 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) ulalume3@0: 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'] ulalume3@0: Type = ['', 'mechanical rotator', 'hwp rotator', 'linear polarizer', 'qwp rotator', 'circular polarizer', '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) ulalume3@0: ulalume3@0: #InputFile = 'optic_input_ver6e_POLIS_355.py' ulalume3@0: #InputFile = 'optic_input_ver6e_POLIS_355_JA.py' ulalume3@0: #InputFile = 'optic_input_ver6c_POLIS_532.py' ulalume3@0: #InputFile = 'optic_input_ver6e_POLIS_532.py' ulalume3@0: #InputFile = 'optic_input_ver8c_POLIS_532.py' ulalume3@0: #InputFile = 'optic_input_ver6e_MUSA.py' ulalume3@0: #InputFile = 'optic_input_ver6e_MUSA_JA.py' ulalume3@0: #InputFile = 'optic_input_ver6e_PollyXTSea.py' ulalume3@0: #InputFile = 'optic_input_ver6e_PollyXTSea_JA.py' ulalume3@0: #InputFile = 'optic_input_ver6e_PollyXT_RALPH.py' ulalume3@0: #InputFile = 'optic_input_ver8c_PollyXT_RALPH.py' ulalume3@0: #InputFile = 'optic_input_ver8c_PollyXT_RALPH_2.py' ulalume3@0: #InputFile = 'optic_input_ver8c_PollyXT_RALPH_3.py' ulalume3@0: #InputFile = 'optic_input_ver8c_PollyXT_RALPH_4.py' ulalume3@0: #InputFile = 'optic_input_ver8c_PollyXT_RALPH_5.py' ulalume3@0: #InputFile = 'optic_input_ver8c_PollyXT_RALPH_6.py' ulalume3@0: InputFile = 'optic_input_ver8c_PollyXT_RALPH_7.py' ulalume3@0: #InputFile = 'optic_input_ver8a_MOHP_DPL_355.py' ulalume3@0: #InputFile = 'optic_input_ver9_MOHP_DPL_355.py' ulalume3@0: #InputFile = 'optic_input_ver6e_RALI.py' ulalume3@0: #InputFile = 'optic_input_ver6e_RALI_JA.py' ulalume3@0: #InputFile = 'optic_input_ver6e_RALI_new.py' ulalume3@0: #InputFile = 'optic_input_ver6e_RALI_act.py' ulalume3@0: #InputFile = 'optic_input_ver6e_MULHACEN.py' ulalume3@0: #InputFile = 'optic_input_ver6e_MULHACEN_JA.py' ulalume3@0: #InputFile = 'optic_input_ver6e-IPRAL.py' ulalume3@0: #InputFile = 'optic_input_ver6e-IPRAL_JA.py' ulalume3@0: #InputFile = 'optic_input_ver6e-LB21.py' ulalume3@0: #InputFile = 'optic_input_ver6e-LB21_JA.py' ulalume3@0: #InputFile = 'optic_input_ver6e_Bertha_b_355.py' ulalume3@0: #InputFile = 'optic_input_ver6e_Bertha_b_532.py' ulalume3@0: #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? ulalume3@0: 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 ) ulalume3@0: #DiO = 0. ulalume3@0: #LDRtrue = 0.45 ulalume3@0: #LDRtrue2 = 0.004 ulalume3@0: #Y = -1 ulalume3@0: #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 ulalume3@0: #RotationErrorEpsilonForNormalMeasurements = True ulalume3@0: #LDRCal = 0.25 ulalume3@0: #bL = 0.8 ulalume3@0: ## --- Errors ulalume3@0: RotL0, dRotL, nRotL = RotL, dRotL, nRotL ulalume3@0: ulalume3@0: DiE0, dDiE, nDiE = DiE, dDiE, nDiE ulalume3@0: RetE0, dRetE, nRetE = RetE, dRetE, nRetE ulalume3@0: RotE0, dRotE, nRotE = RotE, dRotE, nRotE ulalume3@0: ulalume3@0: DiO0, dDiO, nDiO = DiO, dDiO, nDiO ulalume3@0: RetO0, dRetO, nRetO = RetO, dRetO, nRetO ulalume3@0: RotO0, dRotO, nRotO = RotO, dRotO, nRotO ulalume3@0: ulalume3@0: DiC0, dDiC, nDiC = DiC, dDiC, nDiC ulalume3@0: RetC0, dRetC, nRetC = RetC, dRetC, nRetC ulalume3@0: RotC0, dRotC, nRotC = RotC, dRotC, nRotC ulalume3@0: ulalume3@0: TP0, dTP, nTP = TP, dTP, nTP ulalume3@0: TS0, dTS, nTS = TS, dTS, nTS ulalume3@0: RetT0, dRetT, nRetT = RetT, dRetT, nRetT ulalume3@0: ulalume3@0: ERaT0, dERaT, nERaT = ERaT, dERaT, nERaT ulalume3@0: RotaT0,dRotaT,nRotaT= RotaT,dRotaT,nRotaT ulalume3@0: ulalume3@0: RP0, dRP, nRP = RP, dRP, nRP ulalume3@0: RS0, dRS, nRS = RS, dRS, nRS ulalume3@0: RetR0, dRetR, nRetR = RetR, dRetR, nRetR ulalume3@0: ulalume3@0: ERaR0, dERaR, nERaR = ERaR, dERaR, nERaR ulalume3@0: RotaR0,dRotaR,nRotaR= RotaR,dRotaR,nRotaR ulalume3@0: ulalume3@0: LDRCal0,dLDRCal,nLDRCal=LDRCal,dLDRCal,nLDRCal ulalume3@0: #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 ulalume3@0: TP, TS, RP, RS, ERaT, RotaT, RetT, ERaR, RotaR, RetR = TP0, TS0, RP0, RS0 , ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0 ulalume3@0: LDRCal = LDRCal0 ulalume3@0: DTa0, TTa0, DRa0, TRa0, LDRsimx, LDRCorr = 0,0,0,0,0,0 ulalume3@0: ulalume3@0: TiT = 0.5 * (TP + TS) ulalume3@0: DiT = (TP-TS)/(TP+TS) ulalume3@0: ZiT = (1. - DiT**2)**0.5 ulalume3@0: TiR = 0.5 * (RP + RS) ulalume3@0: DiR = (RP-RS)/(RP+RS) ulalume3@0: ZiR = (1. - DiR**2)**0.5 ulalume3@0: ulalume3@0: # -------------------------------------------------------- ulalume3@0: def Calc(RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC, TP, TS, RP, RS, ERaT, RotaT, RetT, ERaR, RotaR, RetR, 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 ulalume3@0: aCal = (1.-LDRCal)/(1+LDRCal) ulalume3@0: # from input file: measured LDRm and true LDRtrue, LDRtrue2 => ulalume3@0: #ameas = (1.-LDRmeas)/(1+LDRmeas) ulalume3@0: atrue = (1.-LDRtrue)/(1+LDRtrue) ulalume3@0: #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 ulalume3@0: S2a = np.sin(2*np.deg2rad(RotL)) ulalume3@0: C2a = np.cos(2*np.deg2rad(RotL)) ulalume3@0: S2b = np.sin(2*np.deg2rad(RotE)) ulalume3@0: C2b = np.cos(2*np.deg2rad(RotE)) ulalume3@0: S2ab = np.sin(np.deg2rad(2*RotL-2*RotE)) ulalume3@0: C2ab = np.cos(np.deg2rad(2*RotL-2*RotE)) ulalume3@0: S2g = np.sin(np.deg2rad(2*RotO)) ulalume3@0: 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. ulalume3@0: VinL = (1. - bL**2)**0.5 ulalume3@0: ulalume3@0: # Stokes Input Vector rotation Eq. E.4 ulalume3@0: A = C2a*QinL - S2a*UinL ulalume3@0: B = S2a*QinL + C2a*UinL ulalume3@0: # Stokes Input Vector rotation Eq. E.9 ulalume3@0: C = C2ab*QinL - S2ab*UinL ulalume3@0: 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)) ulalume3@0: ZiE = (1. - DiE**2)**0.5 ulalume3@0: 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 ulalume3@0: IinE = (IinL + DiE*C) ulalume3@0: QinE = (C2b*DiE*IinL + A + S2b*(WiE*D - ZiE*SinE*VinL)) ulalume3@0: UinE = (S2b*DiE*IinL + B - C2b*(WiE*D - ZiE*SinE*VinL)) ulalume3@0: 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 ulalume3@0: QinF = aCal*QinE ulalume3@0: UinF = -aCal*UinE ulalume3@0: 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)) ulalume3@0: ZiO = (1. - DiO**2)**0.5 ulalume3@0: 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)) ulalume3@0: ZiC = (1. - DiC**2)**0.5 ulalume3@0: WiC = (1. - ZiC*CosC) ulalume3@0: ulalume3@0: # Stokes Input Vector before the polarising beam splitter Eq. E.31 ulalume3@0: A = C2g*QinE - S2g*UinE ulalume3@0: B = S2g*QinE + C2g*UinE ulalume3@0: ulalume3@0: IinP = (IinE + DiO*aCal*A) ulalume3@0: QinP = (C2g*DiO*IinE + aCal*QinE - S2g*(WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE)) ulalume3@0: UinP = (S2g*DiO*IinE - aCal*UinE + C2g*(WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE)) ulalume3@0: VinP = (ZiO*SinO*aCal*B + ZiO*CosO*(1-2*aCal)*VinE) ulalume3@0: ulalume3@0: #------------------------- ulalume3@0: # F11 assuemd to be = 1 => measured: F11m = IinP / IinE with atrue ulalume3@0: #F11sim = TiO*(IinE + DiO*atrue*A)/IinE ulalume3@0: #------------------------- ulalume3@0: ulalume3@0: # For PollyXT ulalume3@0: # analyser ulalume3@0: #RS = 1 - TS ulalume3@0: #RP = 1 - TP ulalume3@0: ulalume3@0: TiT = 0.5 * (TP + TS) ulalume3@0: DiT = (TP-TS)/(TP+TS) ulalume3@0: ZiT = (1. - DiT**2)**0.5 ulalume3@0: TiR = 0.5 * (RP + RS) ulalume3@0: DiR = (RP-RS)/(RP+RS) ulalume3@0: 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: ulalume3@0: DaT = (1-ERaT)/(1+ERaT) ulalume3@0: DaR = (1-ERaR)/(1+ERaR) ulalume3@0: TaT = 0.5*(1+ERaT) ulalume3@0: TaR = 0.5*(1+ERaR) ulalume3@0: ulalume3@0: S2aT = np.sin(np.deg2rad(h*2*RotaT)) ulalume3@0: C2aT = np.cos(np.deg2rad(2*RotaT)) ulalume3@0: S2aR = np.sin(np.deg2rad(h*2*RotaR)) ulalume3@0: C2aR = np.cos(np.deg2rad(2*RotaR)) ulalume3@0: ulalume3@0: # Aanalyzer As before the PBS Eq. D.5 ulalume3@0: ATP1 = (1+C2aT*DaT*DiT) ulalume3@0: ATP2 = Y*(DiT+C2aT*DaT) ulalume3@0: ATP3 = Y*S2aT*DaT*ZiT*CosT ulalume3@0: ATP4 = S2aT*DaT*ZiT*SinT ulalume3@0: ATP = np.array([ATP1,ATP2,ATP3,ATP4]) ulalume3@0: ulalume3@0: ARP1 = (1+C2aR*DaR*DiR) ulalume3@0: ARP2 = Y*(DiR+C2aR*DaR) ulalume3@0: ARP3 = Y*S2aR*DaR*ZiR*CosR ulalume3@0: ARP4 = S2aR*DaR*ZiR*SinR ulalume3@0: ARP = np.array([ARP1,ARP2,ARP3,ARP4]) ulalume3@0: ulalume3@0: DTa = ATP2*Y/ATP1 ulalume3@0: 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 ulalume3@0: #print("Calibrator location not implemented yet") ulalume3@0: ulalume3@0: #S2ge = np.sin(np.deg2rad(2*RotO + h*2*RotC)) ulalume3@0: #C2ge = np.cos(np.deg2rad(2*RotO + h*2*RotC)) ulalume3@0: S2e = np.sin(np.deg2rad(h*2*RotC)) ulalume3@0: C2e = np.cos(np.deg2rad(2*RotC)) ulalume3@0: # rotated AinP by epsilon Eq. C.3 ulalume3@0: ATP2e = C2e*ATP2 + S2e*ATP3 ulalume3@0: ATP3e = C2e*ATP3 - S2e*ATP2 ulalume3@0: ARP2e = C2e*ARP2 + S2e*ARP3 ulalume3@0: ARP3e = C2e*ARP3 - S2e*ARP2 ulalume3@0: ATPe = np.array([ATP1,ATP2e,ATP3e,ATP4]) ulalume3@0: ARPe = np.array([ARP1,ARP2e,ARP3e,ARP4]) ulalume3@0: # Stokes Input Vector before the polarising beam splitter Eq. E.31 ulalume3@0: A = C2g*QinE - S2g*UinE ulalume3@0: B = S2g*QinE + C2g*UinE ulalume3@0: #C = (WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE) ulalume3@0: Co = ZiO*SinO*VinE ulalume3@0: Ca = (WiO*B - 2*ZiO*SinO*VinE) ulalume3@0: #C = Co + aCal*Ca ulalume3@0: #IinP = (IinE + DiO*aCal*A) ulalume3@0: #QinP = (C2g*DiO*IinE + aCal*QinE - S2g*C) ulalume3@0: #UinP = (S2g*DiO*IinE - aCal*UinE + C2g*C) ulalume3@0: #VinP = (ZiO*SinO*aCal*B + ZiO*CosO*(1-2*aCal)*VinE) ulalume3@0: IinPo = IinE ulalume3@0: QinPo = (C2g*DiO*IinE - S2g*Co) ulalume3@0: UinPo = (S2g*DiO*IinE + C2g*Co) ulalume3@0: VinPo = ZiO*CosO*VinE ulalume3@0: ulalume3@0: IinPa = DiO*A ulalume3@0: QinPa = QinE - S2g*Ca ulalume3@0: UinPa = -UinE + C2g*Ca ulalume3@0: VinPa = ZiO*(SinO*B - 2*CosO*VinE) ulalume3@0: ulalume3@0: IinP = IinPo + aCal*IinPa ulalume3@0: QinP = QinPo + aCal*QinPa ulalume3@0: UinP = UinPo + aCal*UinPa ulalume3@0: VinP = VinPo + aCal*VinPa ulalume3@0: # Stokes Input Vector before the polarising beam splitter rotated by epsilon Eq. C.3 ulalume3@0: #QinPe = C2e*QinP + S2e*UinP ulalume3@0: #UinPe = C2e*UinP - S2e*QinP ulalume3@0: QinPoe = C2e*QinPo + S2e*UinPo ulalume3@0: UinPoe = C2e*UinPo - S2e*QinPo ulalume3@0: QinPae = C2e*QinPa + S2e*UinPa ulalume3@0: UinPae = C2e*UinPa - S2e*QinPa ulalume3@0: QinPe = C2e*QinP + S2e*UinP ulalume3@0: 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 ulalume3@0: AT = ATP1*IinP + h*ATP4*VinP ulalume3@0: BT = ATP3e*QinP - h*ATP2e*UinP ulalume3@0: AR = ARP1*IinP + h*ARP4*VinP ulalume3@0: BR = ARP3e*QinP - h*ARP2e*UinP ulalume3@0: # Correction paremeters for normal measurements; they are independent of LDR ulalume3@0: if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out ulalume3@0: IS1 = np.array([IinPo,QinPo,UinPo,VinPo]) ulalume3@0: IS2 = np.array([IinPa,QinPa,UinPa,VinPa]) ulalume3@0: GT = np.dot(ATP,IS1) ulalume3@0: GR = np.dot(ARP,IS1) ulalume3@0: HT = np.dot(ATP,IS2) ulalume3@0: HR = np.dot(ARP,IS2) ulalume3@0: else: ulalume3@0: IS1 = np.array([IinPo,QinPo,UinPo,VinPo]) ulalume3@0: IS2 = np.array([IinPa,QinPa,UinPa,VinPa]) ulalume3@0: GT = np.dot(ATPe,IS1) ulalume3@0: GR = np.dot(ARPe,IS1) ulalume3@0: HT = np.dot(ATPe,IS2) ulalume3@0: 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 ulalume3@0: AT = ATP1*IinP + ATP3e*UinPe + ZiC*CosC*(ATP2e*QinPe + ATP4*VinP) ulalume3@0: BT = DiC*(ATP1*UinPe + ATP3e*IinP) - ZiC*SinC*(ATP2e*VinP - ATP4*QinPe) ulalume3@0: AR = ARP1*IinP + ARP3e*UinPe + ZiC*CosC*(ARP2e*QinPe + ARP4*VinP) ulalume3@0: BR = DiC*(ARP1*UinPe + ARP3e*IinP) - ZiC*SinC*(ARP2e*VinP - ARP4*QinPe) ulalume3@0: # Correction paremeters for normal measurements; they are independent of LDR ulalume3@0: if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out ulalume3@0: IS1 = np.array([IinPo,QinPo,UinPo,VinPo]) ulalume3@0: IS2 = np.array([IinPa,QinPa,UinPa,VinPa]) ulalume3@0: GT = np.dot(ATP,IS1) ulalume3@0: GR = np.dot(ARP,IS1) ulalume3@0: HT = np.dot(ATP,IS2) ulalume3@0: HR = np.dot(ARP,IS2) ulalume3@0: else: ulalume3@0: IS1e = np.array([IinPo+DiC*QinPoe,DiC*IinPo+QinPoe,ZiC*(CosC*UinPoe+SinC*VinPo),-ZiC*(SinC*UinPoe-CosC*VinPo)]) ulalume3@0: IS2e = np.array([IinPa+DiC*QinPae,DiC*IinPa+QinPae,ZiC*(CosC*UinPae+SinC*VinPa),-ZiC*(SinC*UinPae-CosC*VinPa)]) ulalume3@0: GT = np.dot(ATPe,IS1e) ulalume3@0: GR = np.dot(ARPe,IS1e) ulalume3@0: HT = np.dot(ATPe,IS2e) ulalume3@0: 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 ulalume3@0: AT = ATP1*IinP + sqr05*DiC*(ATP1*QinPe + ATP2e*IinP) + (1-0.5*WiC)*(ATP2e*QinPe + ATP3e*UinPe) + ZiC*(sqr05*SinC*(ATP3e*VinP-ATP4*UinPe) + ATP4*CosC*VinP) ulalume3@0: BT = sqr05*DiC*(ATP1*UinPe + ATP3e*IinP) + 0.5*WiC*(ATP2e*UinPe + ATP3e*QinPe) - sqr05*ZiC*SinC*(ATP2e*VinP - ATP4*QinPe) ulalume3@0: AR = ARP1*IinP + sqr05*DiC*(ARP1*QinPe + ARP2e*IinP) + (1-0.5*WiC)*(ARP2e*QinPe + ARP3e*UinPe) + ZiC*(sqr05*SinC*(ARP3e*VinP-ARP4*UinPe) + ARP4*CosC*VinP) ulalume3@0: BR = sqr05*DiC*(ARP1*UinPe + ARP3e*IinP) + 0.5*WiC*(ARP2e*UinPe + ARP3e*QinPe) - sqr05*ZiC*SinC*(ARP2e*VinP - ARP4*QinPe) ulalume3@0: # Correction paremeters for normal measurements; they are independent of LDR ulalume3@0: if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out ulalume3@0: IS1 = np.array([IinPo,QinPo,UinPo,VinPo]) ulalume3@0: IS2 = np.array([IinPa,QinPa,UinPa,VinPa]) ulalume3@0: GT = np.dot(ATP,IS1) ulalume3@0: GR = np.dot(ARP,IS1) ulalume3@0: HT = np.dot(ATP,IS2) ulalume3@0: HR = np.dot(ARP,IS2) ulalume3@0: else: ulalume3@0: IS1e = np.array([IinPo+DiC*QinPoe,DiC*IinPo+QinPoe,ZiC*(CosC*UinPoe+SinC*VinPo),-ZiC*(SinC*UinPoe-CosC*VinPo)]) ulalume3@0: IS2e = np.array([IinPa+DiC*QinPae,DiC*IinPa+QinPae,ZiC*(CosC*UinPae+SinC*VinPa),-ZiC*(SinC*UinPae-CosC*VinPa)]) ulalume3@0: GT = np.dot(ATPe,IS1e) ulalume3@0: GR = np.dot(ARPe,IS1e) ulalume3@0: HT = np.dot(ATPe,IS2e) ulalume3@0: 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: ulalume3@0: #S2ge = np.sin(np.deg2rad(2*RotO - 2*RotC)) ulalume3@0: #C2ge = np.cos(np.deg2rad(2*RotO - 2*RotC)) ulalume3@0: S2e = np.sin(np.deg2rad(2*RotC)) ulalume3@0: C2e = np.cos(np.deg2rad(2*RotC)) ulalume3@0: ulalume3@0: # As with C before the receiver optics (rotated_diattenuator_X22x5deg.odt) ulalume3@0: AF1 = np.array([1,C2g*DiO,S2g*DiO,0]) ulalume3@0: AF2 = np.array([C2g*DiO,1-S2g**2*WiO,S2g*C2g*WiO,-S2g*ZiO*SinO]) ulalume3@0: AF3 = np.array([S2g*DiO,S2g*C2g*WiO,1-C2g**2*WiO,C2g*ZiO*SinO]) ulalume3@0: AF4 = np.array([0,S2g*SinO,-C2g*SinO,CosO]) ulalume3@0: ulalume3@0: ATF = (ATP1*AF1+ATP2*AF2+ATP3*AF3+ATP4*AF4) ulalume3@0: 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] ulalume3@0: ATF2e = C2e*ATF[1] + S2e*ATF[2] ulalume3@0: ATF3e = C2e*ATF[2] - S2e*ATF[1] ulalume3@0: ARF1 = ARF[0] ulalume3@0: ARF4 = ARF[3] ulalume3@0: ARF2e = C2e*ARF[1] + S2e*ARF[2] ulalume3@0: ARF3e = C2e*ARF[2] - S2e*ARF[1] ulalume3@0: ulalume3@0: ATFe = np.array([ATF1,ATF2e,ATF3e,ATF4]) ulalume3@0: ARFe = np.array([ARF1,ARF2e,ARF3e,ARF4]) ulalume3@0: ulalume3@0: QinEe = C2e*QinE + S2e*UinE ulalume3@0: 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 ulalume3@0: QinF = aCal*QinE ulalume3@0: UinF = -aCal*UinE ulalume3@0: 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 ulalume3@0: VinFa = -2.*VinE ulalume3@0: ulalume3@0: # Stokes Input Vector before receiver optics rotated by epsilon Eq. C.3 ulalume3@0: QinFe = C2e*QinF + S2e*UinF ulalume3@0: UinFe = C2e*UinF - S2e*QinF ulalume3@0: QinFoe = C2e*QinFo + S2e*UinFo ulalume3@0: UinFoe = C2e*UinFo - S2e*QinFo ulalume3@0: QinFae = C2e*QinFa + S2e*UinFa ulalume3@0: UinFae = C2e*UinFa - S2e*QinFa 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 ulalume3@0: AT = ATF1*IinF + ATF4*h*VinF ulalume3@0: BT = ATF3e*QinF - ATF2e*h*UinF ulalume3@0: AR = ARF1*IinF + ARF4*h*VinF ulalume3@0: BR = ARF3e*QinF - ARF2e*h*UinF ulalume3@0: # Correction paremeters for normal measurements; they are independent of LDR ulalume3@0: if (not RotationErrorEpsilonForNormalMeasurements): ulalume3@0: GT = ATF1*IinE + ATF4*VinE ulalume3@0: GR = ARF1*IinE + ARF4*VinE ulalume3@0: HT = ATF2*QinE - ATF3*UinE - ATF4*2*VinE ulalume3@0: HR = ARF2*QinE - ARF3*UinE - ARF4*2*VinE ulalume3@0: else: ulalume3@0: GT = ATF1*IinE + ATF4*h*VinE ulalume3@0: GR = ARF1*IinE + ARF4*h*VinE ulalume3@0: HT = ATF2e*QinE - ATF3e*h*UinE - ATF4*h*2*VinE ulalume3@0: 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° ulalume3@0: IF1e = np.array([IinF, ZiC*CosC*QinFe, UinFe, ZiC*CosC*VinF]) ulalume3@0: IF2e = np.array([DiC*UinFe, -ZiC*SinC*VinF, DiC*IinF, ZiC*SinC*QinFe]) ulalume3@0: AT = np.dot(ATFe,IF1e) ulalume3@0: AR = np.dot(ARFe,IF1e) ulalume3@0: BT = np.dot(ATFe,IF2e) ulalume3@0: 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 ulalume3@0: if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out ulalume3@0: IS1 = np.array([IinE,0,0,VinE]) ulalume3@0: IS2 = np.array([0,QinE,-UinE,-2*VinE]) ulalume3@0: GT = np.dot(ATF,IS1) ulalume3@0: GR = np.dot(ARF,IS1) ulalume3@0: HT = np.dot(ATF,IS2) ulalume3@0: HR = np.dot(ARF,IS2) ulalume3@0: else: ulalume3@0: IS1e = np.array([IinFo+DiC*QinFoe,DiC*IinFo+QinFoe,ZiC*(CosC*UinFoe+SinC*VinFo),-ZiC*(SinC*UinFoe-CosC*VinFo)]) ulalume3@0: IS2e = np.array([IinFa+DiC*QinFae,DiC*IinFa+QinFae,ZiC*(CosC*UinFae+SinC*VinFa),-ZiC*(SinC*UinFae-CosC*VinFa)]) ulalume3@0: GT = np.dot(ATFe,IS1e) ulalume3@0: GR = np.dot(ARFe,IS1e) ulalume3@0: HT = np.dot(ATFe,IS2e) ulalume3@0: 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 ulalume3@0: IF1e = np.array([IinF+sqr05*DiC*QinFe, sqr05*DiC*IinF+(1-0.5*WiC)*QinFe, (1-0.5*WiC)*UinFe+sqr05*ZiC*SinC*VinF, -sqr05*ZiC*SinC*UinFe+ZiC*CosC*VinF]) ulalume3@0: IF2e = np.array([sqr05*DiC*UinFe, 0.5*WiC*UinFe-sqr05*ZiC*SinC*VinF, sqr05*DiC*IinF+0.5*WiC*QinFe, sqr05*ZiC*SinC*QinFe]) ulalume3@0: AT = np.dot(ATFe,IF1e) ulalume3@0: AR = np.dot(ARFe,IF1e) ulalume3@0: BT = np.dot(ATFe,IF2e) ulalume3@0: BR = np.dot(ARFe,IF2e) ulalume3@0: ulalume3@0: # Correction paremeters for normal measurements; they are independent of LDR ulalume3@0: if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out ulalume3@0: #IS1 = np.array([IinE,0,0,VinE]) ulalume3@0: #IS2 = np.array([0,QinE,-UinE,-2*VinE]) ulalume3@0: IS1 = np.array([IinFo,0,0,VinFo]) ulalume3@0: IS2 = np.array([0,QinFa,UinFa,VinFa]) ulalume3@0: GT = np.dot(ATF,IS1) ulalume3@0: GR = np.dot(ARF,IS1) ulalume3@0: HT = np.dot(ATF,IS2) ulalume3@0: HR = np.dot(ARF,IS2) ulalume3@0: else: ulalume3@0: IS1e = np.array([IinFo+DiC*QinFoe,DiC*IinFo+QinFoe,ZiC*(CosC*UinFoe+SinC*VinFo),-ZiC*(SinC*UinFoe-CosC*VinFo)]) ulalume3@0: IS2e = np.array([IinFa+DiC*QinFae,DiC*IinFa+QinFae,ZiC*(CosC*UinFae+SinC*VinFa),-ZiC*(SinC*UinFae-CosC*VinFa)]) ulalume3@0: #IS1e = np.array([IinFo,0,0,VinFo]) ulalume3@0: #IS2e = np.array([0,QinFae,UinFae,VinFa]) ulalume3@0: GT = np.dot(ATFe,IS1e) ulalume3@0: GR = np.dot(ARFe,IS1e) ulalume3@0: HT = np.dot(ATFe,IS2e) ulalume3@0: 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 ------------------------------------------------------- ulalume3@0: #print("Calibrator location not implemented yet") ulalume3@0: S2e = np.sin(np.deg2rad(2*RotC)) ulalume3@0: C2e = np.cos(np.deg2rad(2*RotC)) ulalume3@0: ulalume3@0: # AS with C before the receiver optics (see document rotated_diattenuator_X22x5deg.odt) ulalume3@0: AF1 = np.array([1,C2g*DiO,S2g*DiO,0]) ulalume3@0: AF2 = np.array([C2g*DiO,1-S2g**2*WiO,S2g*C2g*WiO,-S2g*ZiO*SinO]) ulalume3@0: AF3 = np.array([S2g*DiO, S2g*C2g*WiO, 1-C2g**2*WiO, C2g*ZiO*SinO]) ulalume3@0: AF4 = np.array([0, S2g*SinO, -C2g*SinO, CosO]) ulalume3@0: ulalume3@0: ATF = (ATP1*AF1+ATP2*AF2+ATP3*AF3+ATP4*AF4) ulalume3@0: 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 ulalume3@0: ATE1a, ARE1a = 0. , 0. ulalume3@0: ATE2a, ARE2a = ATF2, ARF2 ulalume3@0: ATE3a, ARE3a = -ATF3, -ARF3 ulalume3@0: ATE4a, ARE4a = -2*ATF4, -2*ARF4 ulalume3@0: # rotated AinEa by epsilon ulalume3@0: ATE2ae = C2e*ATF2 + S2e*ATF3 ulalume3@0: ATE3ae = -S2e*ATF2 - C2e*ATF3 ulalume3@0: ARE2ae = C2e*ARF2 + S2e*ARF3 ulalume3@0: ARE3ae = -S2e*ARF2 - C2e*ARF3 ulalume3@0: ulalume3@0: ATE1 = ATE1o ulalume3@0: ATE2e = aCal*ATE2ae ulalume3@0: ATE3e = aCal*ATE3ae ulalume3@0: ATE4 = (1-2*aCal)*ATF4 ulalume3@0: ARE1 = ARE1o ulalume3@0: ARE2e = aCal*ARE2ae ulalume3@0: ARE3e = aCal*ARE3ae ulalume3@0: ARE4 = (1-2*aCal)*ARF4 ulalume3@0: ulalume3@0: # rotated IinE ulalume3@0: QinEe = C2e*QinE + S2e*UinE ulalume3@0: UinEe = C2e*UinE - S2e*QinE 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: AT = ATE1o*IinE + (ATE4o+aCal*ATE4a)*h*VinE ulalume3@0: BT = aCal * (ATE3ae*QinEe - ATE2ae*h*UinEe) ulalume3@0: AR = ARE1o*IinE + (ARE4o+aCal*ARE4a)*h*VinE ulalume3@0: 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) ulalume3@0: GT = ATE1o*IinE + ATE4o*h*VinE ulalume3@0: GR = ARE1o*IinE + ARE4o*h*VinE ulalume3@0: HT = ATE2a*QinE + ATE3a*h*UinEe + ATE4a*h*VinE ulalume3@0: HR = ARE2a*QinE + ARE3a*h*UinEe + ARE4a*h*VinE ulalume3@0: else: ulalume3@0: GT = ATE1o*IinE + ATE4o*h*VinE ulalume3@0: GR = ARE1o*IinE + ARE4o*h*VinE ulalume3@0: HT = ATE2ae*QinE + ATE3ae*h*UinEe + ATE4a*h*VinE ulalume3@0: 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° ulalume3@0: AT = ATE1*IinE + ZiC*CosC*(ATE2e*QinEe + ATE4*VinE) + ATE3e*UinEe ulalume3@0: BT = DiC*(ATE1*UinEe + ATE3e*IinE) + ZiC*SinC*(ATE4*QinEe - ATE2e*VinE) ulalume3@0: AR = ARE1*IinE + ZiC*CosC*(ARE2e*QinEe + ARE4*VinE) + ARE3e*UinEe ulalume3@0: 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) ulalume3@0: GT = ATE1o*IinE + ATE4o*VinE ulalume3@0: GR = ARE1o*IinE + ARE4o*VinE ulalume3@0: HT = ATE2a*QinE + ATE3a*UinE + ATE4a*VinE ulalume3@0: HR = ARE2a*QinE + ARE3a*UinE + ARE4a*VinE ulalume3@0: else: ulalume3@0: D = IinE + DiC*QinEe ulalume3@0: A = DiC*IinE + QinEe ulalume3@0: B = ZiC*(CosC*UinEe + SinC*VinE) ulalume3@0: C = -ZiC*(SinC*UinEe - CosC*VinE) ulalume3@0: GT = ATE1o*D + ATE4o*C ulalume3@0: GR = ARE1o*D + ARE4o*C ulalume3@0: HT = ATE2a*A + ATE3a*B + ATE4a*C ulalume3@0: 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° ulalume3@0: IE1e = np.array([IinE+sqr05*DiC*QinEe, sqr05*DiC*IinE+(1-0.5*WiC)*QinEe, (1-0.5*WiC)*UinEe+sqr05*ZiC*SinC*VinE, -sqr05*ZiC*SinC*UinEe+ZiC*CosC*VinE]) ulalume3@0: IE2e = np.array([sqr05*DiC*UinEe, 0.5*WiC*UinEe-sqr05*ZiC*SinC*VinE, sqr05*DiC*IinE+0.5*WiC*QinEe, sqr05*ZiC*SinC*QinEe]) ulalume3@0: ATEe = np.array([ATE1,ATE2e,ATE3e,ATE4]) ulalume3@0: AREe = np.array([ARE1,ARE2e,ARE3e,ARE4]) ulalume3@0: AT = np.dot(ATEe,IE1e) ulalume3@0: AR = np.dot(AREe,IE1e) ulalume3@0: BT = np.dot(ATEe,IE2e) ulalume3@0: BR = np.dot(AREe,IE2e) ulalume3@0: ulalume3@0: # Correction paremeters for normal measurements; they are independent of LDR ulalume3@0: if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out ulalume3@0: GT = ATE1o*IinE + ATE4o*VinE ulalume3@0: GR = ARE1o*IinE + ARE4o*VinE ulalume3@0: HT = ATE2a*QinE + ATE3a*UinE + ATE4a*VinE ulalume3@0: HR = ARE2a*QinE + ARE3a*UinE + ARE4a*VinE ulalume3@0: else: ulalume3@0: D = IinE + DiC*QinEe ulalume3@0: A = DiC*IinE + QinEe ulalume3@0: B = ZiC*(CosC*UinEe + SinC*VinE) ulalume3@0: C = -ZiC*(SinC*UinEe - CosC*VinE) ulalume3@0: GT = ATE1o*D + ATE4o*C ulalume3@0: GR = ARE1o*D + ARE4o*C ulalume3@0: HT = ATE2a*A + ATE3a*B + ATE4a*C ulalume3@0: 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 ulalume3@0: IoutTp = TaT*TiT*TiO*TiE*(AT + BT) ulalume3@0: IoutTm = TaT*TiT*TiO*TiE*(AT - BT) ulalume3@0: IoutRp = TaR*TiR*TiO*TiE*(AR + BR) ulalume3@0: IoutRm = TaR*TiR*TiO*TiE*(AR - BR) ulalume3@0: ulalume3@0: # --- Results and Corrections; electronic etaR and etaT are assumed to be 1 ulalume3@0: Etapx = IoutRp/IoutTp ulalume3@0: Etamx = IoutRm/IoutTm ulalume3@0: Etax = (Etapx*Etamx)**0.5 ulalume3@0: ulalume3@0: 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 ulalume3@0: #Eta_test_p = (IoutRp/IoutTp) ulalume3@0: #Eta_test_m = (IoutRm/IoutTm) ulalume3@0: #Eta_test = (Eta_test_p*Eta_test_m)**0.5 ulalume3@0: ulalume3@0: # ----- Forward simulated signals and LDRsim with atrue; from input file ulalume3@0: It = TaT*TiT*TiO*TiE*(GT+atrue*HT) ulalume3@0: Ir = TaR*TiR*TiO*TiE*(GR+atrue*HR) ulalume3@0: # LDRsim = 1/Eta*Ir/It # simulated LDR* with Y from input file ulalume3@0: 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.: ulalume3@0: LDRsimx = 1./LDRsim ulalume3@0: else: ulalume3@0: LDRsimx = LDRsim ulalume3@0: ulalume3@0: # The following is correct without doubt ulalume3@0: #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 ulalume3@0: LDRCorr = (LDRsim/Eta*(GT+HT)-(GR+HR))/((GR-HR)-LDRsim*K/Etax*(GT-HT)) ulalume3@0: ulalume3@0: TTa = TiT*TaT #*ATP1 ulalume3@0: TRa = TiR*TaR #*ARP1 ulalume3@0: ulalume3@0: F11sim = 1/(TiO*TiE)*((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) ulalume3@0: # ******************************************************************************************************************************* ulalume3@0: ulalume3@0: # --- CALC truth ulalume3@0: GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0) ulalume3@0: ulalume3@0: # -------------------------------------------------------- ulalume3@0: with open('output_' + LID + '.dat', 'w') as f: ulalume3@0: with redirect_stdout(f): ulalume3@0: print("From ", dname) ulalume3@0: print("Running ", fname) ulalume3@0: 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 ----------------------") ulalume3@0: print("{0:8} {1:8} {2:8.5f}; {3:8} {4:7.4f}±{5:7.4f}/{6:2d}".format("Laser: ", "DOLP = ", bL, " rotation alpha = ", RotL0, dRotL, nRotL)) ulalume3@0: print(" Diatt., Tunpol, Retard., Rotation (deg)") ulalume3@0: print("{0:12} {1:7.4f}±{2:7.4f}/{8:2d}, {3:7.4f}, {4:3.0f}±{5:3.0f}/{9:2d}, {6:7.4f}±{7:7.4f}/{10:2d}".format("Emitter ", DiE0, dDiE, TiE, RetE0, dRetE, RotE0, dRotE, nDiE, nRetE, nRotE)) ulalume3@0: print("{0:12} {1:7.4f}±{2:7.4f}/{8:2d}, {3:7.4f}, {4:3.0f}±{5:3.0f}/{9:2d}, {6:7.4f}±{7:7.4f}/{10:2d}".format("Receiver ", DiO0, dDiO, TiO, RetO0, dRetO, RotO0, dRotO, nDiO, nRetO, nRotO)) ulalume3@0: print("{0:12} {1:7.4f}±{2:7.4f}/{8:2d}, {3:7.4f}, {4:3.0f}±{5:3.0f}/{9:2d}, {6:7.4f}±{7:7.4f}/{10:2d}".format("Calibrator ", DiC0, dDiC, TiC, RetC0, dRetC, RotC0, dRotC, nDiC, nRetC, nRotC)) ulalume3@0: print("{0:12}".format(" --- Pol.-filter ---")) ulalume3@0: print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format("ERT, ERR :", ERaT0, dERaT, nERaT, ERaR0, dERaR, nERaR)) ulalume3@0: print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format("RotaT , RotaR :", RotaT0, dRotaT, nRotaT, RotaR0,dRotaR,nRotaR)) ulalume3@0: print("{0:12}".format(" --- PBS ---")) ulalume3@0: print("{0:12}{1:7.4f}±{2:7.4f}/{9:2d}, {3:7.4f}±{4:7.4f}/{10:2d}, {5:7.4f}±{6:7.4f}/{11:2d},{7:7.4f}±{8:7.4f}/{12:2d}".format("TP,TS,RP,RS :", TP0, dTP, TS0, dTS, RP0, dRP, RS0, dRS, nTP, nTS, nRP, 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 ---")) ulalume3@0: print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}".format("DTa,TTa,DRa,TRa: ", DTa0, TTa0, DRa0, TRa0)) ulalume3@0: print() ulalume3@0: print("Rotation Error Epsilon For Normal Measurements = ", RotationErrorEpsilonForNormalMeasurements) ulalume3@0: #print ('LocC = ', LocC, Loc[LocC], '; TypeC = ',TypeC, Type[TypeC]) ulalume3@0: print(Type[TypeC], Loc[LocC], "; Parallel signal detected in", dY[int(Y+1)]) ulalume3@0: # end of print actual system parameters ulalume3@0: # ****************************************************************************** ulalume3@0: ulalume3@0: #print() ulalume3@0: #print(" --- LDRCal during calibration | simulated and corrected LDRs -------------") ulalume3@0: #print("{0:8} |{1:8}->{2:8},{3:9}->{4:9} |{5:8}->{6:8}".format(" LDRCal"," LDRtrue", " LDRsim"," LDRtrue2", " LDRsim2", " LDRmeas", " LDRcorr")) ulalume3@0: #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)) ulalume3@0: #print("{0:8} |{1:8}->{2:8}->{3:8}".format(" LDRCal"," LDRtrue", " LDRsimx", " LDRcorr")) ulalume3@0: #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)) ulalume3@0: #print("{0:8} |{1:8}->{2:8}->{3:8}".format(" LDRCal"," LDRtrue", " LDRsimx", " LDRcorr")) ulalume3@0: #print(" --- LDRCal during calibration") ulalume3@0: print("{0:26}: {1:6.3f}±{2:5.3f}/{3:2d}".format("LDRCal during calibration", LDRCal0, dLDRCal, nLDRCal)) ulalume3@0: ulalume3@0: #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 ulalume3@0: for i,LDRCal in enumerate(LDRCalList): ulalume3@0: GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal) ulalume3@0: K0List[i] = K0 ulalume3@0: LDRsimxList[i] = LDRsimx ulalume3@0: ulalume3@0: print("{0:8},{1:8},{2:8},{3:8},{4:9},{5:9},{6:9}".format(" GR", " GT", " HR", " HT", " K(0.004)", " K(0.2)", " K(0.45)")) ulalume3@0: print("{0:8.5f},{1:8.5f},{2:8.5f},{3:8.5f},{4:9.5f},{5:9.5f},{6:9.5f}".format(GR0, GT0, HR0, HT0, K0List[0], K0List[1], K0List[2])) 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 ulalume3@0: for i,LDRtrue in enumerate(LDRtrueList): ulalume3@0: GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0) ulalume3@0: print("{0:9.5f},{1:9.5f},{2:9.5f}".format(LDRtrue, LDRsimx, LDRCorr)) ulalume3@0: ulalume3@0: ulalume3@0: file = open('output_' + LID + '.dat', 'r') ulalume3@0: 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: ''' ulalume3@0: # --- CALC again truth with LDRCal0 to reset all 0-values ulalume3@0: GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0) ulalume3@0: ulalume3@0: # --- Start Errors calculation ulalume3@0: ulalume3@0: iN = -1 ulalume3@0: N = ((nRotL*2+1)* ulalume3@0: (nRotE*2+1)*(nRetE*2+1)*(nDiE*2+1)* ulalume3@0: (nRotO*2+1)*(nRetO*2+1)*(nDiO*2+1)* ulalume3@0: (nRotC*2+1)*(nRetC*2+1)*(nDiC*2+1)* ulalume3@0: (nTP*2+1)*(nTS*2+1)*(nRP*2+1)*(nRS*2+1)*(nERaT*2+1)*(nERaR*2+1)* ulalume3@0: (nRotaT*2+1)*(nRotaR*2+1)*(nRetT*2+1)*(nRetR*2+1)*(nLDRCal*2+1)) ulalume3@0: print("N = ",N ," ", end="") ulalume3@0: ulalume3@0: if N > 1e6: ulalume3@0: if user_yes_no_query('Warning: processing ' + str(N) + ' samples will take very long. Do you want to proceed?') == 0: sys.exit() ulalume3@0: if N > 5e6: ulalume3@0: if user_yes_no_query('Warning: the memory required for ' + str(N) + ' samples might be ' + '{0:5.1f}'.format(N/4e6) + ' GB. Do you anyway want to proceed?') == 0: sys.exit() ulalume3@0: ulalume3@0: #if user_yes_no_query('Warning: processing' + str(N) + ' samples will take very long. Do you want to proceed?') == 0: sys.exit() ulalume3@0: ulalume3@0: # --- Arrays for plotting ------ ulalume3@0: LDRmin = np.zeros(5) ulalume3@0: LDRmax = np.zeros(5) ulalume3@0: F11min = np.zeros(5) ulalume3@0: F11max = np.zeros(5) ulalume3@0: ulalume3@0: LDRrange = np.zeros(5) ulalume3@0: LDRrange = 0.004, 0.02, 0.1, 0.3, 0.45 ulalume3@0: #aLDRsimx = np.zeros(N) ulalume3@0: #aLDRsimx2 = np.zeros(N) ulalume3@0: #aLDRcorr = np.zeros(N) ulalume3@0: #aLDRcorr2 = np.zeros(N) ulalume3@0: aERaT = np.zeros(N) ulalume3@0: aERaR = np.zeros(N) ulalume3@0: aRotaT = np.zeros(N) ulalume3@0: aRotaR = np.zeros(N) ulalume3@0: aRetT = np.zeros(N) ulalume3@0: aRetR = np.zeros(N) ulalume3@0: aTP = np.zeros(N) ulalume3@0: aTS = np.zeros(N) ulalume3@0: aRP = np.zeros(N) ulalume3@0: aRS = np.zeros(N) ulalume3@0: aDiE = np.zeros(N) ulalume3@0: aDiO = np.zeros(N) ulalume3@0: aDiC = np.zeros(N) ulalume3@0: aRotC = np.zeros(N) ulalume3@0: aRetC = np.zeros(N) ulalume3@0: aRotL = np.zeros(N) ulalume3@0: aRetE = np.zeros(N) ulalume3@0: aRotE = np.zeros(N) ulalume3@0: aRetO = np.zeros(N) ulalume3@0: aRotO = np.zeros(N) ulalume3@0: aLDRCal = np.zeros(N) ulalume3@0: aA = np.zeros((5,N)) ulalume3@0: aX = np.zeros((5,N)) ulalume3@0: aF11corr = np.zeros((5,N)) ulalume3@0: ulalume3@0: atime = clock() ulalume3@0: dtime = clock() ulalume3@0: ulalume3@0: # --- Calc Error signals ulalume3@0: #GT, HT, GR, HR, K, Eta, LDRsim = Calc(RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC, TP, TS) ulalume3@0: # ---- Do the calculations of bra-ket vectors ulalume3@0: h = -1. if TypeC == 2 else 1 ulalume3@0: ulalume3@0: # from input file: measured LDRm and true LDRtrue, LDRtrue2 => ulalume3@0: ameas = (1.-LDRmeas)/(1+LDRmeas) ulalume3@0: atrue = (1.-LDRtrue)/(1+LDRtrue) ulalume3@0: atrue2 = (1.-LDRtrue2)/(1+LDRtrue2) ulalume3@0: ulalume3@0: for iLDRCal in range(-nLDRCal,nLDRCal+1): ulalume3@0: # from input file: assumed LDRCal for calibration measurements ulalume3@0: LDRCal = LDRCal0 ulalume3@0: if nLDRCal > 0: LDRCal = LDRCal0 + iLDRCal*dLDRCal/nLDRCal ulalume3@0: ulalume3@0: GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal) ulalume3@0: aCal = (1.-LDRCal)/(1+LDRCal) ulalume3@0: for iRotL, iRotE, iRetE, iDiE \ ulalume3@0: in [(iRotL,iRotE,iRetE,iDiE) ulalume3@0: for iRotL in range(-nRotL,nRotL+1) ulalume3@0: for iRotE in range(-nRotE,nRotE+1) ulalume3@0: for iRetE in range(-nRetE,nRetE+1) ulalume3@0: for iDiE in range(-nDiE,nDiE+1)]: ulalume3@0: ulalume3@0: if nRotL > 0: RotL = RotL0 + iRotL*dRotL/nRotL ulalume3@0: if nRotE > 0: RotE = RotE0 + iRotE*dRotE/nRotE ulalume3@0: if nRetE > 0: RetE = RetE0 + iRetE*dRetE/nRetE ulalume3@0: if nDiE > 0: DiE = DiE0 + iDiE*dDiE/nDiE ulalume3@0: ulalume3@0: # angles of emitter and laser and calibrator and receiver optics ulalume3@0: # RotL = alpha, RotE = beta, RotO = gamma, RotC = epsilon ulalume3@0: S2a = np.sin(2*np.deg2rad(RotL)) ulalume3@0: C2a = np.cos(2*np.deg2rad(RotL)) ulalume3@0: S2b = np.sin(2*np.deg2rad(RotE)) ulalume3@0: C2b = np.cos(2*np.deg2rad(RotE)) ulalume3@0: S2ab = np.sin(np.deg2rad(2*RotL-2*RotE)) ulalume3@0: C2ab = np.cos(np.deg2rad(2*RotL-2*RotE)) ulalume3@0: ulalume3@0: # Laser with Degree of linear polarization DOLP = bL ulalume3@0: IinL = 1. ulalume3@0: QinL = bL ulalume3@0: UinL = 0. ulalume3@0: VinL = (1. - bL**2)**0.5 ulalume3@0: ulalume3@0: # Stokes Input Vector rotation Eq. E.4 ulalume3@0: A = C2a*QinL - S2a*UinL ulalume3@0: B = S2a*QinL + C2a*UinL ulalume3@0: # Stokes Input Vector rotation Eq. E.9 ulalume3@0: C = C2ab*QinL - S2ab*UinL ulalume3@0: 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)) ulalume3@0: ZiE = (1. - DiE**2)**0.5 ulalume3@0: 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 ulalume3@0: IinE = (IinL + DiE*C) ulalume3@0: QinE = (C2b*DiE*IinL + A + S2b*(WiE*D - ZiE*SinE*VinL)) ulalume3@0: UinE = (S2b*DiE*IinL + B - C2b*(WiE*D - ZiE*SinE*VinL)) ulalume3@0: VinE = (-ZiE*SinE*D + ZiE*CosE*VinL) ulalume3@0: ulalume3@0: #------------------------- ulalume3@0: # F11 assuemd to be = 1 => measured: F11m = IinP / IinE with atrue ulalume3@0: #F11sim = (IinE + DiO*atrue*(C2g*QinE - S2g*UinE))/IinE ulalume3@0: #------------------------- ulalume3@0: ulalume3@0: for iRotO, iRetO, iDiO, iRotC, iRetC, iDiC, iTP, iTS, iRP, iRS, iERaT, iRotaT, iRetT, iERaR, iRotaR, iRetR \ ulalume3@0: in [(iRotO,iRetO,iDiO,iRotC,iRetC,iDiC,iTP,iTS,iRP,iRS,iERaT,iRotaT,iRetT,iERaR,iRotaR,iRetR ) ulalume3@0: for iRotO in range(-nRotO,nRotO+1) ulalume3@0: for iRetO in range(-nRetO,nRetO+1) ulalume3@0: for iDiO in range(-nDiO,nDiO+1) ulalume3@0: for iRotC in range(-nRotC,nRotC+1) ulalume3@0: for iRetC in range(-nRetC,nRetC+1) ulalume3@0: for iDiC in range(-nDiC,nDiC+1) ulalume3@0: for iTP in range(-nTP,nTP+1) ulalume3@0: for iTS in range(-nTS,nTS+1) ulalume3@0: for iRP in range(-nRP,nRP+1) ulalume3@0: for iRS in range(-nRS,nRS+1) ulalume3@0: for iERaT in range(-nERaT,nERaT+1) ulalume3@0: for iRotaT in range(-nRotaT,nRotaT+1) ulalume3@0: for iRetT in range(-nRetT,nRetT+1) ulalume3@0: for iERaR in range(-nERaR,nERaR+1) ulalume3@0: for iRotaR in range(-nRotaR,nRotaR+1) ulalume3@0: for iRetR in range(-nRetR,nRetR+1)]: ulalume3@0: ulalume3@0: iN = iN + 1 ulalume3@0: if (iN == 10001): ulalume3@0: ctime = clock() ulalume3@0: print(" estimated time ", "{0:4.2f}".format(N/10000 * (ctime-atime)), "sec ") #, end="") ulalume3@0: print("\r elapsed time ", "{0:5.0f}".format((ctime-atime)), "sec ", end="\r") ulalume3@0: ctime = clock() ulalume3@0: if ((ctime - dtime) > 10): ulalume3@0: print("\r elapsed time ", "{0:5.0f}".format((ctime-atime)), "sec ", end="\r") ulalume3@0: dtime = ctime ulalume3@0: ulalume3@0: if nRotO > 0: RotO = RotO0 + iRotO*dRotO/nRotO ulalume3@0: if nRetO > 0: RetO = RetO0 + iRetO*dRetO/nRetO ulalume3@0: if nDiO > 0: DiO = DiO0 + iDiO*dDiO/nDiO ulalume3@0: if nRotC > 0: RotC = RotC0 + iRotC*dRotC/nRotC ulalume3@0: if nRetC > 0: RetC = RetC0 + iRetC*dRetC/nRetC ulalume3@0: if nDiC > 0: DiC = DiC0 + iDiC*dDiC/nDiC ulalume3@0: if nTP > 0: TP = TP0 + iTP*dTP/nTP ulalume3@0: if nTS > 0: TS = TS0 + iTS*dTS/nTS ulalume3@0: if nRP > 0: RP = RP0 + iRP*dRP/nRP ulalume3@0: if nRS > 0: RS = RS0 + iRS*dRS/nRS ulalume3@0: if nERaT > 0: ERaT = ERaT0 + iERaT*dERaT/nERaT ulalume3@0: if nRotaT > 0:RotaT= RotaT0+ iRotaT*dRotaT/nRotaT ulalume3@0: if nRetT > 0: RetT = RetT0 + iRetT*dRetT/nRetT ulalume3@0: if nERaR > 0: ERaR = ERaR0 + iERaR*dERaR/nERaR ulalume3@0: if nRotaR > 0:RotaR= RotaR0+ iRotaR*dRotaR/nRotaR ulalume3@0: if nRetR > 0: RetR = RetR0 + iRetR*dRetR/nRetR ulalume3@0: ulalume3@0: #print("{0:5.2f}, {1:5.2f}, {2:5.2f}, {3:10d}".format(RotL, RotE, RotO, iN)) ulalume3@0: ulalume3@0: # receiver optics ulalume3@0: CosO = np.cos(np.deg2rad(RetO)) ulalume3@0: SinO = np.sin(np.deg2rad(RetO)) ulalume3@0: ZiO = (1. - DiO**2)**0.5 ulalume3@0: WiO = (1. - ZiO*CosO) ulalume3@0: S2g = np.sin(np.deg2rad(2*RotO)) ulalume3@0: C2g = np.cos(np.deg2rad(2*RotO)) ulalume3@0: # calibrator ulalume3@0: CosC = np.cos(np.deg2rad(RetC)) ulalume3@0: SinC = np.sin(np.deg2rad(RetC)) ulalume3@0: ZiC = (1. - DiC**2)**0.5 ulalume3@0: WiC = (1. - ZiC*CosC) ulalume3@0: ulalume3@0: #For POLLY_XT ulalume3@0: # analyser ulalume3@0: #RS = 1 - TS ulalume3@0: #RP = 1 - TP ulalume3@0: TiT = 0.5 * (TP + TS) ulalume3@0: DiT = (TP-TS)/(TP+TS) ulalume3@0: ZiT = (1. - DiT**2)**0.5 ulalume3@0: TiR = 0.5 * (RP + RS) ulalume3@0: DiR = (RP-RS)/(RP+RS) ulalume3@0: 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: ulalume3@0: DaT = (1-ERaT)/(1+ERaT) ulalume3@0: DaR = (1-ERaR)/(1+ERaR) ulalume3@0: TaT = 0.5*(1+ERaT) ulalume3@0: TaR = 0.5*(1+ERaR) ulalume3@0: ulalume3@0: S2aT = np.sin(np.deg2rad(h*2*RotaT)) ulalume3@0: C2aT = np.cos(np.deg2rad(2*RotaT)) ulalume3@0: S2aR = np.sin(np.deg2rad(h*2*RotaR)) ulalume3@0: C2aR = np.cos(np.deg2rad(2*RotaR)) ulalume3@0: ulalume3@0: # Aanalyzer As before the PBS Eq. D.5 ulalume3@0: ATP1 = (1+C2aT*DaT*DiT) ulalume3@0: ATP2 = Y*(DiT+C2aT*DaT) ulalume3@0: ATP3 = Y*S2aT*DaT*ZiT*CosT ulalume3@0: ATP4 = S2aT*DaT*ZiT*SinT ulalume3@0: ATP = np.array([ATP1,ATP2,ATP3,ATP4]) ulalume3@0: ulalume3@0: ARP1 = (1+C2aR*DaR*DiR) ulalume3@0: ARP2 = Y*(DiR+C2aR*DaR) ulalume3@0: ARP3 = Y*S2aR*DaR*ZiR*CosR ulalume3@0: ARP4 = S2aR*DaR*ZiR*SinR ulalume3@0: ARP = np.array([ARP1,ARP2,ARP3,ARP4]) ulalume3@0: ulalume3@0: TTa = TiT*TaT #*ATP1 ulalume3@0: TRa = TiR*TaR #*ARP1 ulalume3@0: ulalume3@0: # ---- Calculate signals and correction parameters for diffeent locations and calibrators ulalume3@0: if LocC == 4: # Calibrator before the PBS ulalume3@0: #print("Calibrator location not implemented yet") ulalume3@0: ulalume3@0: #S2ge = np.sin(np.deg2rad(2*RotO + h*2*RotC)) ulalume3@0: #C2ge = np.cos(np.deg2rad(2*RotO + h*2*RotC)) ulalume3@0: S2e = np.sin(np.deg2rad(h*2*RotC)) ulalume3@0: C2e = np.cos(np.deg2rad(2*RotC)) ulalume3@0: # rotated AinP by epsilon Eq. C.3 ulalume3@0: ATP2e = C2e*ATP2 + S2e*ATP3 ulalume3@0: ATP3e = C2e*ATP3 - S2e*ATP2 ulalume3@0: ARP2e = C2e*ARP2 + S2e*ARP3 ulalume3@0: ARP3e = C2e*ARP3 - S2e*ARP2 ulalume3@0: ATPe = np.array([ATP1,ATP2e,ATP3e,ATP4]) ulalume3@0: ARPe = np.array([ARP1,ARP2e,ARP3e,ARP4]) ulalume3@0: # Stokes Input Vector before the polarising beam splitter Eq. E.31 ulalume3@0: A = C2g*QinE - S2g*UinE ulalume3@0: B = S2g*QinE + C2g*UinE ulalume3@0: #C = (WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE) ulalume3@0: Co = ZiO*SinO*VinE ulalume3@0: Ca = (WiO*B - 2*ZiO*SinO*VinE) ulalume3@0: #C = Co + aCal*Ca ulalume3@0: #IinP = (IinE + DiO*aCal*A) ulalume3@0: #QinP = (C2g*DiO*IinE + aCal*QinE - S2g*C) ulalume3@0: #UinP = (S2g*DiO*IinE - aCal*UinE + C2g*C) ulalume3@0: #VinP = (ZiO*SinO*aCal*B + ZiO*CosO*(1-2*aCal)*VinE) ulalume3@0: IinPo = IinE ulalume3@0: QinPo = (C2g*DiO*IinE - S2g*Co) ulalume3@0: UinPo = (S2g*DiO*IinE + C2g*Co) ulalume3@0: VinPo = ZiO*CosO*VinE ulalume3@0: ulalume3@0: IinPa = DiO*A ulalume3@0: QinPa = QinE - S2g*Ca ulalume3@0: UinPa = -UinE + C2g*Ca ulalume3@0: VinPa = ZiO*(SinO*B - 2*CosO*VinE) ulalume3@0: ulalume3@0: IinP = IinPo + aCal*IinPa ulalume3@0: QinP = QinPo + aCal*QinPa ulalume3@0: UinP = UinPo + aCal*UinPa ulalume3@0: VinP = VinPo + aCal*VinPa ulalume3@0: # Stokes Input Vector before the polarising beam splitter rotated by epsilon Eq. C.3 ulalume3@0: #QinPe = C2e*QinP + S2e*UinP ulalume3@0: #UinPe = C2e*UinP - S2e*QinP ulalume3@0: QinPoe = C2e*QinPo + S2e*UinPo ulalume3@0: UinPoe = C2e*UinPo - S2e*QinPo ulalume3@0: QinPae = C2e*QinPa + S2e*UinPa ulalume3@0: UinPae = C2e*UinPa - S2e*QinPa ulalume3@0: QinPe = C2e*QinP + S2e*UinP ulalume3@0: 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 ulalume3@0: AT = ATP1*IinP + h*ATP4*VinP ulalume3@0: BT = ATP3e*QinP - h*ATP2e*UinP ulalume3@0: AR = ARP1*IinP + h*ARP4*VinP ulalume3@0: BR = ARP3e*QinP - h*ARP2e*UinP ulalume3@0: # Correction paremeters for normal measurements; they are independent of LDR ulalume3@0: if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out ulalume3@0: IS1 = np.array([IinPo,QinPo,UinPo,VinPo]) ulalume3@0: IS2 = np.array([IinPa,QinPa,UinPa,VinPa]) ulalume3@0: GT = np.dot(ATP,IS1) ulalume3@0: GR = np.dot(ARP,IS1) ulalume3@0: HT = np.dot(ATP,IS2) ulalume3@0: HR = np.dot(ARP,IS2) ulalume3@0: else: ulalume3@0: IS1 = np.array([IinPo,QinPo,UinPo,VinPo]) ulalume3@0: IS2 = np.array([IinPa,QinPa,UinPa,VinPa]) ulalume3@0: GT = np.dot(ATPe,IS1) ulalume3@0: GR = np.dot(ARPe,IS1) ulalume3@0: HT = np.dot(ATPe,IS2) ulalume3@0: 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 ulalume3@0: AT = ATP1*IinP + ATP3e*UinPe + ZiC*CosC*(ATP2e*QinPe + ATP4*VinP) ulalume3@0: BT = DiC*(ATP1*UinPe + ATP3e*IinP) - ZiC*SinC*(ATP2e*VinP - ATP4*QinPe) ulalume3@0: AR = ARP1*IinP + ARP3e*UinPe + ZiC*CosC*(ARP2e*QinPe + ARP4*VinP) ulalume3@0: BR = DiC*(ARP1*UinPe + ARP3e*IinP) - ZiC*SinC*(ARP2e*VinP - ARP4*QinPe) ulalume3@0: # Correction paremeters for normal measurements; they are independent of LDR ulalume3@0: if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out ulalume3@0: IS1 = np.array([IinPo,QinPo,UinPo,VinPo]) ulalume3@0: IS2 = np.array([IinPa,QinPa,UinPa,VinPa]) ulalume3@0: GT = np.dot(ATP,IS1) ulalume3@0: GR = np.dot(ARP,IS1) ulalume3@0: HT = np.dot(ATP,IS2) ulalume3@0: HR = np.dot(ARP,IS2) ulalume3@0: else: ulalume3@0: IS1e = np.array([IinPo+DiC*QinPoe,DiC*IinPo+QinPoe,ZiC*(CosC*UinPoe+SinC*VinPo),-ZiC*(SinC*UinPoe-CosC*VinPo)]) ulalume3@0: IS2e = np.array([IinPa+DiC*QinPae,DiC*IinPa+QinPae,ZiC*(CosC*UinPae+SinC*VinPa),-ZiC*(SinC*UinPae-CosC*VinPa)]) ulalume3@0: GT = np.dot(ATPe,IS1e) ulalume3@0: GR = np.dot(ARPe,IS1e) ulalume3@0: HT = np.dot(ATPe,IS2e) ulalume3@0: 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 ulalume3@0: AT = ATP1*IinP + sqr05*DiC*(ATP1*QinPe + ATP2e*IinP) + (1-0.5*WiC)*(ATP2e*QinPe + ATP3e*UinPe) + ZiC*(sqr05*SinC*(ATP3e*VinP-ATP4*UinPe) + ATP4*CosC*VinP) ulalume3@0: BT = sqr05*DiC*(ATP1*UinPe + ATP3e*IinP) + 0.5*WiC*(ATP2e*UinPe + ATP3e*QinPe) - sqr05*ZiC*SinC*(ATP2e*VinP - ATP4*QinPe) ulalume3@0: AR = ARP1*IinP + sqr05*DiC*(ARP1*QinPe + ARP2e*IinP) + (1-0.5*WiC)*(ARP2e*QinPe + ARP3e*UinPe) + ZiC*(sqr05*SinC*(ARP3e*VinP-ARP4*UinPe) + ARP4*CosC*VinP) ulalume3@0: BR = sqr05*DiC*(ARP1*UinPe + ARP3e*IinP) + 0.5*WiC*(ARP2e*UinPe + ARP3e*QinPe) - sqr05*ZiC*SinC*(ARP2e*VinP - ARP4*QinPe) ulalume3@0: # Correction paremeters for normal measurements; they are independent of LDR ulalume3@0: if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out ulalume3@0: IS1 = np.array([IinPo,QinPo,UinPo,VinPo]) ulalume3@0: IS2 = np.array([IinPa,QinPa,UinPa,VinPa]) ulalume3@0: GT = np.dot(ATP,IS1) ulalume3@0: GR = np.dot(ARP,IS1) ulalume3@0: HT = np.dot(ATP,IS2) ulalume3@0: HR = np.dot(ARP,IS2) ulalume3@0: else: ulalume3@0: IS1e = np.array([IinPo+DiC*QinPoe,DiC*IinPo+QinPoe,ZiC*(CosC*UinPoe+SinC*VinPo),-ZiC*(SinC*UinPoe-CosC*VinPo)]) ulalume3@0: IS2e = np.array([IinPa+DiC*QinPae,DiC*IinPa+QinPae,ZiC*(CosC*UinPae+SinC*VinPa),-ZiC*(SinC*UinPae-CosC*VinPa)]) ulalume3@0: GT = np.dot(ATPe,IS1e) ulalume3@0: GR = np.dot(ARPe,IS1e) ulalume3@0: HT = np.dot(ATPe,IS2e) ulalume3@0: 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: ulalume3@0: #S2ge = np.sin(np.deg2rad(2*RotO - 2*RotC)) ulalume3@0: #C2ge = np.cos(np.deg2rad(2*RotO - 2*RotC)) ulalume3@0: S2e = np.sin(np.deg2rad(2*RotC)) ulalume3@0: C2e = np.cos(np.deg2rad(2*RotC)) ulalume3@0: ulalume3@0: # AS with C before the receiver optics (see document rotated_diattenuator_X22x5deg.odt) ulalume3@0: AF1 = np.array([1,C2g*DiO,S2g*DiO,0]) ulalume3@0: AF2 = np.array([C2g*DiO,1-S2g**2*WiO,S2g*C2g*WiO,-S2g*ZiO*SinO]) ulalume3@0: AF3 = np.array([S2g*DiO, S2g*C2g*WiO, 1-C2g**2*WiO, C2g*ZiO*SinO]) ulalume3@0: AF4 = np.array([0, S2g*SinO, -C2g*SinO, CosO]) ulalume3@0: ulalume3@0: ATF = (ATP1*AF1+ATP2*AF2+ATP3*AF3+ATP4*AF4) ulalume3@0: 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: # rotated AinF by epsilon ulalume3@0: ATF2e = C2e*ATF[1] + S2e*ATF[2] ulalume3@0: ATF3e = C2e*ATF[2] - S2e*ATF[1] ulalume3@0: ARF2e = C2e*ARF[1] + S2e*ARF[2] ulalume3@0: ARF3e = C2e*ARF[2] - S2e*ARF[1] ulalume3@0: ulalume3@0: ATFe = np.array([ATF1,ATF2e,ATF3e,ATF4]) ulalume3@0: ARFe = np.array([ARF1,ARF2e,ARF3e,ARF4]) ulalume3@0: ulalume3@0: QinEe = C2e*QinE + S2e*UinE ulalume3@0: 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 ulalume3@0: QinF = aCal*QinE ulalume3@0: UinF = -aCal*UinE ulalume3@0: 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 ulalume3@0: VinFa = -2.*VinE ulalume3@0: ulalume3@0: # Stokes Input Vector before receiver optics rotated by epsilon Eq. C.3 ulalume3@0: QinFe = C2e*QinF + S2e*UinF ulalume3@0: UinFe = C2e*UinF - S2e*QinF ulalume3@0: QinFoe = C2e*QinFo + S2e*UinFo ulalume3@0: UinFoe = C2e*UinFo - S2e*QinFo ulalume3@0: QinFae = C2e*QinFa + S2e*UinFa ulalume3@0: UinFae = C2e*UinFa - S2e*QinFa 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: AT = ATF1*IinF + ATF4*h*VinF ulalume3@0: BT = ATF3e*QinF - ATF2e*h*UinF ulalume3@0: AR = ARF1*IinF + ARF4*h*VinF ulalume3@0: BR = ARF3e*QinF - ARF2e*h*UinF ulalume3@0: ulalume3@0: # Correction paremeters for normal measurements; they are independent of LDR ulalume3@0: if (not RotationErrorEpsilonForNormalMeasurements): ulalume3@0: GT = ATF1*IinE + ATF4*VinE ulalume3@0: GR = ARF1*IinE + ARF4*VinE ulalume3@0: HT = ATF2*QinE - ATF3*UinE - ATF4*2*VinE ulalume3@0: HR = ARF2*QinE - ARF3*UinE - ARF4*2*VinE ulalume3@0: else: ulalume3@0: GT = ATF1*IinE + ATF4*h*VinE ulalume3@0: GR = ARF1*IinE + ARF4*h*VinE ulalume3@0: HT = ATF2e*QinE - ATF3e*h*UinE - ATF4*h*2*VinE ulalume3@0: HR = ARF2e*QinE - ARF3e*h*UinE - ARF4*h*2*VinE ulalume3@0: ulalume3@0: elif (TypeC == 3) or (TypeC == 4): # linear polariser calibration Eq. C.5 ulalume3@0: # p = +45°, m = -45° ulalume3@0: IF1e = np.array([IinF, ZiC*CosC*QinFe, UinFe, ZiC*CosC*VinF]) ulalume3@0: IF2e = np.array([DiC*UinFe, -ZiC*SinC*VinF, DiC*IinF, ZiC*SinC*QinFe]) ulalume3@0: ulalume3@0: AT = np.dot(ATFe,IF1e) ulalume3@0: AR = np.dot(ARFe,IF1e) ulalume3@0: BT = np.dot(ATFe,IF2e) ulalume3@0: 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 ulalume3@0: if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out ulalume3@0: IS1 = np.array([IinE,0,0,VinE]) ulalume3@0: IS2 = np.array([0,QinE,-UinE,-2*VinE]) ulalume3@0: ulalume3@0: GT = np.dot(ATF,IS1) ulalume3@0: GR = np.dot(ARF,IS1) ulalume3@0: HT = np.dot(ATF,IS2) ulalume3@0: HR = np.dot(ARF,IS2) ulalume3@0: else: ulalume3@0: IS1e = np.array([IinFo+DiC*QinFoe,DiC*IinFo+QinFoe,ZiC*(CosC*UinFoe+SinC*VinFo),-ZiC*(SinC*UinFoe-CosC*VinFo)]) ulalume3@0: IS2e = np.array([IinFa+DiC*QinFae,DiC*IinFa+QinFae,ZiC*(CosC*UinFae+SinC*VinFa),-ZiC*(SinC*UinFae-CosC*VinFa)]) ulalume3@0: GT = np.dot(ATFe,IS1e) ulalume3@0: GR = np.dot(ARFe,IS1e) ulalume3@0: HT = np.dot(ATFe,IS2e) ulalume3@0: HR = np.dot(ARFe,IS2e) ulalume3@0: ulalume3@0: elif (TypeC == 6): # diattenuator calibration +-22.5° rotated_diattenuator_X22x5deg.odt ulalume3@0: # p = +22.5°, m = -22.5° ulalume3@0: IF1e = np.array([IinF+sqr05*DiC*QinFe, sqr05*DiC*IinF+(1-0.5*WiC)*QinFe, (1-0.5*WiC)*UinFe+sqr05*ZiC*SinC*VinF, -sqr05*ZiC*SinC*UinFe+ZiC*CosC*VinF]) ulalume3@0: IF2e = np.array([sqr05*DiC*UinFe, 0.5*WiC*UinFe-sqr05*ZiC*SinC*VinF, sqr05*DiC*IinF+0.5*WiC*QinFe, sqr05*ZiC*SinC*QinFe]) ulalume3@0: ulalume3@0: AT = np.dot(ATFe,IF1e) ulalume3@0: AR = np.dot(ARFe,IF1e) ulalume3@0: BT = np.dot(ATFe,IF2e) ulalume3@0: BR = np.dot(ARFe,IF2e) ulalume3@0: ulalume3@0: # Correction paremeters for normal measurements; they are independent of LDR ulalume3@0: if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out ulalume3@0: #IS1 = np.array([IinE,0,0,VinE]) ulalume3@0: #IS2 = np.array([0,QinE,-UinE,-2*VinE]) ulalume3@0: IS1 = np.array([IinFo,0,0,VinFo]) ulalume3@0: IS2 = np.array([0,QinFa,UinFa,VinFa]) ulalume3@0: GT = np.dot(ATF,IS1) ulalume3@0: GR = np.dot(ARF,IS1) ulalume3@0: HT = np.dot(ATF,IS2) ulalume3@0: HR = np.dot(ARF,IS2) ulalume3@0: else: ulalume3@0: #IS1e = np.array([IinE,DiC*IinE,ZiC*SinC*VinE,ZiC*CosC*VinE]) ulalume3@0: #IS2e = np.array([DiC*QinEe,QinEe,-ZiC*(CosC*UinEe+2*SinC*VinE),ZiC*(SinC*UinEe-2*CosC*VinE)]) ulalume3@0: IS1e = np.array([IinFo+DiC*QinFoe,DiC*IinFo+QinFoe,ZiC*(CosC*UinFoe+SinC*VinFo),-ZiC*(SinC*UinFoe-CosC*VinFo)]) ulalume3@0: IS2e = np.array([IinFa+DiC*QinFae,DiC*IinFa+QinFae,ZiC*(CosC*UinFae+SinC*VinFa),-ZiC*(SinC*UinFae-CosC*VinFa)]) ulalume3@0: GT = np.dot(ATFe,IS1e) ulalume3@0: GR = np.dot(ARFe,IS1e) ulalume3@0: HT = np.dot(ATFe,IS2e) ulalume3@0: HR = np.dot(ARFe,IS2e) ulalume3@0: 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 ulalume3@0: #print("Calibrator location not implemented yet") ulalume3@0: S2e = np.sin(np.deg2rad(2*RotC)) ulalume3@0: C2e = np.cos(np.deg2rad(2*RotC)) ulalume3@0: ulalume3@0: # AS with C before the receiver optics (see document rotated_diattenuator_X22x5deg.odt) ulalume3@0: AF1 = np.array([1,C2g*DiO,S2g*DiO,0]) ulalume3@0: AF2 = np.array([C2g*DiO,1-S2g**2*WiO,S2g*C2g*WiO,-S2g*ZiO*SinO]) ulalume3@0: AF3 = np.array([S2g*DiO, S2g*C2g*WiO, 1-C2g**2*WiO, C2g*ZiO*SinO]) ulalume3@0: AF4 = np.array([0, S2g*SinO, -C2g*SinO, CosO]) ulalume3@0: ulalume3@0: ATF = (ATP1*AF1+ATP2*AF2+ATP3*AF3+ATP4*AF4) ulalume3@0: 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 ulalume3@0: ATE1a, ARE1a = 0. , 0. ulalume3@0: ATE2a, ARE2a = ATF2, ARF2 ulalume3@0: ATE3a, ARE3a = -ATF3, -ARF3 ulalume3@0: ATE4a, ARE4a = -2*ATF4, -2*ARF4 ulalume3@0: # rotated AinEa by epsilon ulalume3@0: ATE2ae = C2e*ATF2 + S2e*ATF3 ulalume3@0: ATE3ae = -S2e*ATF2 - C2e*ATF3 ulalume3@0: ARE2ae = C2e*ARF2 + S2e*ARF3 ulalume3@0: ARE3ae = -S2e*ARF2 - C2e*ARF3 ulalume3@0: ulalume3@0: ATE1 = ATE1o ulalume3@0: ATE2e = aCal*ATE2ae ulalume3@0: ATE3e = aCal*ATE3ae ulalume3@0: ATE4 = (1-2*aCal)*ATF4 ulalume3@0: ARE1 = ARE1o ulalume3@0: ARE2e = aCal*ARE2ae ulalume3@0: ARE3e = aCal*ARE3ae ulalume3@0: ARE4 = (1-2*aCal)*ARF4 ulalume3@0: ulalume3@0: # rotated IinE ulalume3@0: QinEe = C2e*QinE + S2e*UinE ulalume3@0: UinEe = C2e*UinE - S2e*QinE 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: AT = ATE1o*IinE + (ATE4o+aCal*ATE4a)*h*VinE ulalume3@0: BT = aCal * (ATE3ae*QinEe - ATE2ae*h*UinEe) ulalume3@0: AR = ARE1o*IinE + (ARE4o+aCal*ARE4a)*h*VinE ulalume3@0: 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) ulalume3@0: GT = ATE1o*IinE + ATE4o*h*VinE ulalume3@0: GR = ARE1o*IinE + ARE4o*h*VinE ulalume3@0: HT = ATE2a*QinE + ATE3a*h*UinEe + ATE4a*h*VinE ulalume3@0: HR = ARE2a*QinE + ARE3a*h*UinEe + ARE4a*h*VinE ulalume3@0: else: ulalume3@0: GT = ATE1o*IinE + ATE4o*h*VinE ulalume3@0: GR = ARE1o*IinE + ARE4o*h*VinE ulalume3@0: HT = ATE2ae*QinE + ATE3ae*h*UinEe + ATE4a*h*VinE ulalume3@0: 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° ulalume3@0: AT = ATE1*IinE + ZiC*CosC*(ATE2e*QinEe + ATE4*VinE) + ATE3e*UinEe ulalume3@0: BT = DiC*(ATE1*UinEe + ATE3e*IinE) + ZiC*SinC*(ATE4*QinEe - ATE2e*VinE) ulalume3@0: AR = ARE1*IinE + ZiC*CosC*(ARE2e*QinEe + ARE4*VinE) + ARE3e*UinEe ulalume3@0: 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) ulalume3@0: GT = ATE1o*IinE + ATE4o*VinE ulalume3@0: GR = ARE1o*IinE + ARE4o*VinE ulalume3@0: HT = ATE2a*QinE + ATE3a*UinE + ATE4a*VinE ulalume3@0: HR = ARE2a*QinE + ARE3a*UinE + ARE4a*VinE ulalume3@0: else: ulalume3@0: D = IinE + DiC*QinEe ulalume3@0: A = DiC*IinE + QinEe ulalume3@0: B = ZiC*(CosC*UinEe + SinC*VinE) ulalume3@0: C = -ZiC*(SinC*UinEe - CosC*VinE) ulalume3@0: GT = ATE1o*D + ATE4o*C ulalume3@0: GR = ARE1o*D + ARE4o*C ulalume3@0: HT = ATE2a*A + ATE3a*B + ATE4a*C ulalume3@0: 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° ulalume3@0: IE1e = np.array([IinE+sqr05*DiC*QinEe, sqr05*DiC*IinE+(1-0.5*WiC)*QinEe, (1-0.5*WiC)*UinEe+sqr05*ZiC*SinC*VinE, -sqr05*ZiC*SinC*UinEe+ZiC*CosC*VinE]) ulalume3@0: IE2e = np.array([sqr05*DiC*UinEe, 0.5*WiC*UinEe-sqr05*ZiC*SinC*VinE, sqr05*DiC*IinE+0.5*WiC*QinEe, sqr05*ZiC*SinC*QinEe]) ulalume3@0: ATEe = np.array([ATE1,ATE2e,ATE3e,ATE4]) ulalume3@0: AREe = np.array([ARE1,ARE2e,ARE3e,ARE4]) ulalume3@0: AT = np.dot(ATEe,IE1e) ulalume3@0: AR = np.dot(AREe,IE1e) ulalume3@0: BT = np.dot(ATEe,IE2e) ulalume3@0: BR = np.dot(AREe,IE2e) ulalume3@0: ulalume3@0: # Correction paremeters for normal measurements; they are independent of LDR ulalume3@0: if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out ulalume3@0: GT = ATE1o*IinE + ATE4o*VinE ulalume3@0: GR = ARE1o*IinE + ARE4o*VinE ulalume3@0: HT = ATE2a*QinE + ATE3a*UinE + ATE4a*VinE ulalume3@0: HR = ARE2a*QinE + ARE3a*UinE + ARE4a*VinE ulalume3@0: else: ulalume3@0: D = IinE + DiC*QinEe ulalume3@0: A = DiC*IinE + QinEe ulalume3@0: B = ZiC*(CosC*UinEe + SinC*VinE) ulalume3@0: C = -ZiC*(SinC*UinEe - CosC*VinE) ulalume3@0: GT = ATE1o*D + ATE4o*C ulalume3@0: GR = ARE1o*D + ARE4o*C ulalume3@0: HT = ATE2a*A + ATE3a*B + ATE4a*C ulalume3@0: 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: # Calibration signals with aCal => Determination of the correction K of the real calibration factor ulalume3@0: IoutTp = TaT*TiT*TiO*TiE*(AT + BT) ulalume3@0: IoutTm = TaT*TiT*TiO*TiE*(AT - BT) ulalume3@0: IoutRp = TaR*TiR*TiO*TiE*(AR + BR) ulalume3@0: IoutRm = TaR*TiR*TiO*TiE*(AR - BR) ulalume3@0: # --- Results and Corrections; electronic etaR and etaT are assumed to be 1 ulalume3@0: #Eta = TiR/TiT # Eta = Eta*/K Eq. 84 ulalume3@0: Etapx = IoutRp/IoutTp ulalume3@0: Etamx = IoutRm/IoutTm ulalume3@0: Etax = (Etapx*Etamx)**0.5 ulalume3@0: K = Etax / Eta0 ulalume3@0: #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)) ulalume3@0: #print("{0:6.3f},{1:6.3f},{2:6.3f},{3:6.3f}".format(DiC, ZiC, Kp, Km)) ulalume3@0: ulalume3@0: # For comparison with Volkers Libreoffice Müller Matrix spreadsheet ulalume3@0: #Eta_test_p = (IoutRp/IoutTp) ulalume3@0: #Eta_test_m = (IoutRm/IoutTm) ulalume3@0: #Eta_test = (Eta_test_p*Eta_test_m)**0.5 ulalume3@0: ulalume3@0: # ************************************************************************* ulalume3@0: iLDR = -1 ulalume3@0: for LDRTrue in LDRrange: ulalume3@0: iLDR = iLDR + 1 ulalume3@0: atrue = (1-LDRTrue)/(1+LDRTrue) ulalume3@0: # ----- Forward simulated signals and LDRsim with atrue; from input file ulalume3@0: It = TaT*TiT*TiO*TiE*(GT+atrue*HT) # TaT*TiT*TiC*TiO*IinL*(GT+atrue*HT) ulalume3@0: Ir = TaR*TiR*TiO*TiE*(GR+atrue*HR) # TaR*TiR*TiC*TiO*IinL*(GR+atrue*HR) ulalume3@0: ulalume3@0: # LDRsim = 1/Eta*Ir/It # simulated LDR* with Y from input file ulalume3@0: LDRsim = Ir/It # simulated uncorrected LDR with Y from input file ulalume3@0: ''' ulalume3@0: if Y == 1.: ulalume3@0: LDRsimx = LDRsim ulalume3@0: LDRsimx2 = LDRsim2 ulalume3@0: else: ulalume3@0: LDRsimx = 1./LDRsim ulalume3@0: LDRsimx2 = 1./LDRsim2 ulalume3@0: ''' ulalume3@0: # ----- Backward correction ulalume3@0: # Corrected LDRCorr from forward simulated LDRsim (atrue) with assumed true G0,H0,K0 ulalume3@0: LDRCorr = (LDRsim*K0/Etax*(GT0+HT0)-(GR0+HR0))/((GR0-HR0)-LDRsim*K0/Etax*(GT0-HT0)) ulalume3@0: ulalume3@0: # -- F11corr from It and Ir and calibration EtaX ulalume3@0: Text1 = "F11corr from It and Ir with calibration EtaX: x-axis: F11corr(LDRtrue) / F11corr(LDRtrue = 0.004) - 1" ulalume3@0: F11corr = 1/(TiO*TiE)*((HR0*Etax/K0*It/TTa-HT0*Ir/TRa)/(HR0*GT0-HT0*GR0)) # IL = 1 Eq.(64) ulalume3@0: ulalume3@0: #Text1 = "F11corr from It and Ir without corrections but with calibration EtaX: x-axis: F11corr(LDRtrue) devided by F11corr(LDRtrue = 0.004)" ulalume3@0: #F11corr = 0.5/(TiO*TiE)*(Etax*It/TTa+Ir/TRa) # IL = 1 Eq.(64) ulalume3@0: ulalume3@0: # -- It from It only with atrue without corrections - for BERTHA (and PollyXTs) ulalume3@0: #Text1 = " x-axis: IT(LDRtrue) / IT(LDRtrue = 0.004) - 1" ulalume3@0: #F11corr = It/(TaT*TiT*TiO*TiE) #/(TaT*TiT*TiO*TiE*(GT0+atrue*HT0)) ulalume3@0: # !!! see below line 1673ff ulalume3@0: ulalume3@0: aF11corr[iLDR,iN] = F11corr ulalume3@0: aA[iLDR,iN] = LDRCorr ulalume3@0: ulalume3@0: aX[0,iN] = GR ulalume3@0: aX[1,iN] = GT ulalume3@0: aX[2,iN] = HR ulalume3@0: aX[3,iN] = HT ulalume3@0: aX[4,iN] = K ulalume3@0: ulalume3@0: aLDRCal[iN] = iLDRCal ulalume3@0: aERaT[iN] = iERaT ulalume3@0: aERaR[iN] = iERaR ulalume3@0: aRotaT[iN] = iRotaT ulalume3@0: aRotaR[iN] = iRotaR ulalume3@0: aRetT[iN] = iRetT ulalume3@0: aRetR[iN] = iRetR ulalume3@0: ulalume3@0: aRotL[iN] = iRotL ulalume3@0: aRotE[iN] = iRotE ulalume3@0: aRetE[iN] = iRetE ulalume3@0: aRotO[iN] = iRotO ulalume3@0: aRetO[iN] = iRetO ulalume3@0: aRotC[iN] = iRotC ulalume3@0: aRetC[iN] = iRetC ulalume3@0: aDiO[iN] = iDiO ulalume3@0: aDiE[iN] = iDiE ulalume3@0: aDiC[iN] = iDiC ulalume3@0: aTP[iN] = iTP ulalume3@0: aTS[iN] = iTS ulalume3@0: aRP[iN] = iRP ulalume3@0: aRS[iN] = iRS ulalume3@0: ulalume3@0: # --- END loop ulalume3@0: btime = clock() ulalume3@0: print("\r done in ", "{0:5.0f}".format(btime-atime), "sec") #, end="\r") ulalume3@0: ulalume3@0: # --- Plot ----------------------------------------------------------------- ulalume3@0: #sns.set_style("whitegrid") ulalume3@0: #sns.set_palette("bright", 6) ulalume3@0: ulalume3@0: ''' ulalume3@0: fig2 = plt.figure() ulalume3@0: plt.plot(aA[2,:],'b.') ulalume3@0: plt.plot(aA[3,:],'r.') ulalume3@0: plt.plot(aA[4,:],'g.') ulalume3@0: #plt.plot(aA[6,:],'c.') ulalume3@0: plt.show ulalume3@0: ''' ulalume3@0: # Plot LDR ulalume3@0: def PlotSubHist(aVar, aX, X0, daX, iaX, naX): ulalume3@0: fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) ulalume3@0: iLDR = -1 ulalume3@0: for LDRTrue in LDRrange: ulalume3@0: iLDR = iLDR + 1 ulalume3@0: ulalume3@0: LDRmin[iLDR] = np.min(aA[iLDR,:]) ulalume3@0: LDRmax[iLDR] = np.max(aA[iLDR,:]) ulalume3@0: Rmin = LDRmin[iLDR] * 0.995 # np.min(aA[iLDR,:]) * 0.995 ulalume3@0: Rmax = LDRmax[iLDR] * 1.005 # np.max(aA[iLDR,:]) * 1.005 ulalume3@0: ulalume3@0: #plt.subplot(5,2,iLDR+1) ulalume3@0: plt.subplot(1,5,iLDR+1) ulalume3@0: (n, bins, patches) = plt.hist(aA[iLDR,:], ulalume3@0: bins=100, log=False, ulalume3@0: range=[Rmin, Rmax], ulalume3@0: alpha=0.5, normed=False, color = '0.5', histtype='stepfilled') ulalume3@0: ulalume3@0: for iaX in range(-naX,naX+1): ulalume3@0: plt.hist(aA[iLDR,aX == iaX], ulalume3@0: range=[Rmin, Rmax], ulalume3@0: bins=100, log=False, alpha=0.3, normed=False, histtype='stepfilled', label = str(round(X0 + iaX*daX/naX,5))) ulalume3@0: ulalume3@0: if (iLDR == 2): plt.legend() ulalume3@0: ulalume3@0: plt.tick_params(axis='both', labelsize=9) ulalume3@0: plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2) ulalume3@0: ulalume3@0: #plt.title(LID + ' ' + aVar, fontsize=18) ulalume3@0: #plt.ylabel('frequency', fontsize=10) ulalume3@0: #plt.xlabel('LDRcorr', fontsize=10) ulalume3@0: #fig.tight_layout() ulalume3@0: fig.suptitle(LID + ' ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar, fontsize=14, y=1.05) ulalume3@0: #plt.show() ulalume3@0: #fig.savefig(LID + '_' + aVar + '.png', dpi=150, bbox_inches='tight', pad_inches=0) ulalume3@0: #plt.close ulalume3@0: return ulalume3@0: ulalume3@0: if (nRotL > 0): PlotSubHist("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) ulalume3@0: if (nRetE > 0): PlotSubHist("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) ulalume3@0: if (nRotE > 0): PlotSubHist("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) ulalume3@0: if (nDiE > 0): PlotSubHist("DiE", aDiE, DiE0, dDiE, iDiE, nDiE) ulalume3@0: if (nRetO > 0): PlotSubHist("RetO", aRetO, RetO0, dRetO, iRetO, nRetO) ulalume3@0: if (nRotO > 0): PlotSubHist("RotO", aRotO, RotO0, dRotO, iRotO, nRotO) ulalume3@0: if (nDiO > 0): PlotSubHist("DiO", aDiO, DiO0, dDiO, iDiO, nDiO) ulalume3@0: if (nDiC > 0): PlotSubHist("DiC", aDiC, DiC0, dDiC, iDiC, nDiC) ulalume3@0: if (nRotC > 0): PlotSubHist("RotC", aRotC, RotC0, dRotC, iRotC, nRotC) ulalume3@0: if (nRetC > 0): PlotSubHist("RetC", aRetC, RetC0, dRetC, iRetC, nRetC) ulalume3@0: if (nTP > 0): PlotSubHist("TP", aTP, TP0, dTP, iTP, nTP) ulalume3@0: if (nTS > 0): PlotSubHist("TS", aTS, TS0, dTS, iTS, nTS) ulalume3@0: if (nRP > 0): PlotSubHist("RP", aRP, RP0, dRP, iRP, nRP) ulalume3@0: if (nRS > 0): PlotSubHist("RS", aRS, RS0, dRS, iRS, nRS) ulalume3@0: if (nRetT > 0): PlotSubHist("RetT", aRetT, RetT0, dRetT, iRetT, nRetT) ulalume3@0: if (nRetR > 0): PlotSubHist("RetR", aRetR, RetR0, dRetR, iRetR, nRetR) ulalume3@0: if (nERaT > 0): PlotSubHist("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) ulalume3@0: if (nERaR > 0): PlotSubHist("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) ulalume3@0: if (nRotaT > 0): PlotSubHist("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) ulalume3@0: if (nRotaR > 0): PlotSubHist("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) ulalume3@0: if (nLDRCal > 0): PlotSubHist("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) ulalume3@0: ulalume3@0: plt.show() ulalume3@0: plt.close ulalume3@0: ulalume3@0: print() ulalume3@0: #print("IT(LDRtrue) devided by IT(LDRtrue = 0.004)") ulalume3@0: print(Text1) ulalume3@0: print() ulalume3@0: ulalume3@0: iLDR = 5 ulalume3@0: for LDRTrue in LDRrange: ulalume3@0: iLDR = iLDR - 1 ulalume3@0: aF11corr[iLDR,:] = aF11corr[iLDR,:] / aF11corr[0,:] - 1.0 ulalume3@0: ulalume3@0: # Plot F11 ulalume3@0: def PlotSubHistF11(aVar, aX, X0, daX, iaX, naX): ulalume3@0: fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) ulalume3@0: iLDR = -1 ulalume3@0: for LDRTrue in LDRrange: ulalume3@0: iLDR = iLDR + 1 ulalume3@0: ulalume3@0: ''' ulalume3@0: F11min[iLDR] = np.min(aF11corr[iLDR,:]) ulalume3@0: F11max[iLDR] = np.max(aF11corr[iLDR,:]) ulalume3@0: Rmin = F11min[iLDR] * 0.995 # np.min(aA[iLDR,:]) * 0.995 ulalume3@0: Rmax = F11max[iLDR] * 1.005 # np.max(aA[iLDR,:]) * 1.005 ulalume3@0: ''' ulalume3@0: #Rmin = 0.8 ulalume3@0: #Rmax = 1.2 ulalume3@0: ulalume3@0: #plt.subplot(5,2,iLDR+1) ulalume3@0: plt.subplot(1,5,iLDR+1) ulalume3@0: (n, bins, patches) = plt.hist(aF11corr[iLDR,:], ulalume3@0: bins=100, log=False, ulalume3@0: alpha=0.5, normed=False, color = '0.5', histtype='stepfilled') ulalume3@0: ulalume3@0: for iaX in range(-naX,naX+1): ulalume3@0: plt.hist(aF11corr[iLDR,aX == iaX], ulalume3@0: bins=100, log=False, alpha=0.3, normed=False, histtype='stepfilled', label = str(round(X0 + iaX*daX/naX,5))) ulalume3@0: ulalume3@0: if (iLDR == 2): plt.legend() ulalume3@0: ulalume3@0: plt.tick_params(axis='both', labelsize=9) ulalume3@0: #plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2) ulalume3@0: ulalume3@0: #plt.title(LID + ' ' + aVar, fontsize=18) ulalume3@0: #plt.ylabel('frequency', fontsize=10) ulalume3@0: #plt.xlabel('LDRcorr', fontsize=10) ulalume3@0: #fig.tight_layout() ulalume3@0: fig.suptitle(LID + ' ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar, fontsize=14, y=1.05) ulalume3@0: #plt.show() ulalume3@0: #fig.savefig(LID + '_' + aVar + '.png', dpi=150, bbox_inches='tight', pad_inches=0) ulalume3@0: #plt.close ulalume3@0: return ulalume3@0: ulalume3@0: if (nRotL > 0): PlotSubHistF11("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) ulalume3@0: if (nRetE > 0): PlotSubHistF11("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) ulalume3@0: if (nRotE > 0): PlotSubHistF11("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) ulalume3@0: if (nDiE > 0): PlotSubHistF11("DiE", aDiE, DiE0, dDiE, iDiE, nDiE) ulalume3@0: if (nRetO > 0): PlotSubHistF11("RetO", aRetO, RetO0, dRetO, iRetO, nRetO) ulalume3@0: if (nRotO > 0): PlotSubHistF11("RotO", aRotO, RotO0, dRotO, iRotO, nRotO) ulalume3@0: if (nDiO > 0): PlotSubHistF11("DiO", aDiO, DiO0, dDiO, iDiO, nDiO) ulalume3@0: if (nDiC > 0): PlotSubHistF11("DiC", aDiC, DiC0, dDiC, iDiC, nDiC) ulalume3@0: if (nRotC > 0): PlotSubHistF11("RotC", aRotC, RotC0, dRotC, iRotC, nRotC) ulalume3@0: if (nRetC > 0): PlotSubHistF11("RetC", aRetC, RetC0, dRetC, iRetC, nRetC) ulalume3@0: if (nTP > 0): PlotSubHistF11("TP", aTP, TP0, dTP, iTP, nTP) ulalume3@0: if (nTS > 0): PlotSubHistF11("TS", aTS, TS0, dTS, iTS, nTS) ulalume3@0: if (nRP > 0): PlotSubHistF11("RP", aRP, RP0, dRP, iRP, nRP) ulalume3@0: if (nRS > 0): PlotSubHistF11("RS", aRS, RS0, dRS, iRS, nRS) ulalume3@0: if (nRetT > 0): PlotSubHistF11("RetT", aRetT, RetT0, dRetT, iRetT, nRetT) ulalume3@0: if (nRetR > 0): PlotSubHistF11("RetR", aRetR, RetR0, dRetR, iRetR, nRetR) ulalume3@0: if (nERaT > 0): PlotSubHistF11("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) ulalume3@0: if (nERaR > 0): PlotSubHistF11("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) ulalume3@0: if (nRotaT > 0): PlotSubHistF11("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) ulalume3@0: if (nRotaR > 0): PlotSubHistF11("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) ulalume3@0: if (nLDRCal > 0): PlotSubHistF11("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) ulalume3@0: ulalume3@0: plt.show() ulalume3@0: plt.close ulalume3@0: ''' ulalume3@0: # only histogram ulalume3@0: #print("******************* " + aVar + " *******************") ulalume3@0: 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,:]) ulalume3@0: Rmin = np.min(aA[iLDR,:]) * 0.999 ulalume3@0: Rmax = np.max(aA[iLDR,:]) * 1.001 ulalume3@0: plt.subplot(5,2,iLDR+1) ulalume3@0: (n, bins, patches) = plt.hist(aA[iLDR,:], ulalume3@0: range=[Rmin, Rmax], ulalume3@0: 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) ulalume3@0: plt.show() ulalume3@0: plt.close ulalume3@0: ''' ulalume3@0: ulalume3@0: # --- Plot LDRmin, LDRmax ulalume3@0: fig2 = plt.figure() ulalume3@0: plt.plot(LDRrange,LDRmax-LDRrange, linewidth=2.0, color='b') ulalume3@0: plt.plot(LDRrange,LDRmin-LDRrange, linewidth=2.0, color='g') ulalume3@0: ulalume3@0: plt.xlabel('LDRtrue', fontsize=18) ulalume3@0: plt.ylabel('LDRTrue-LDRmin, LDRTrue-LDRmax', fontsize=14) ulalume3@0: plt.title(LID + ' ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]), fontsize=18) ulalume3@0: #plt.ylimit(-0.07, 0.07) ulalume3@0: plt.show() ulalume3@0: plt.close ulalume3@0: ulalume3@0: # --- Save LDRmin, LDRmax to file ulalume3@0: # http://stackoverflow.com/questions/4675728/redirect-stdout-to-a-file-in-python ulalume3@0: with open('LDR_min_max_ver7_' + LID + '.dat', 'w') as f: ulalume3@0: with redirect_stdout(f): ulalume3@0: print(LID) ulalume3@0: print("LDRtrue, LDRmin, LDRmax") ulalume3@0: for i in range(len(LDRrange)): ulalume3@0: print("{0:7.4f},{1:7.4f},{2:7.4f}".format(LDRrange[i], LDRmin[i], LDRmax[i])) ulalume3@0: ulalume3@0: ''' ulalume3@0: # --- Plot K over LDRCal ulalume3@0: fig3 = plt.figure() ulalume3@0: plt.plot(LDRCal0+aLDRCal*dLDRCal/nLDRCal,aX[4,:], linewidth=2.0, color='b') ulalume3@0: ulalume3@0: plt.xlabel('LDRCal', fontsize=18) ulalume3@0: plt.ylabel('K', fontsize=14) ulalume3@0: plt.title(LID, fontsize=18) ulalume3@0: plt.show() ulalume3@0: plt.close ulalume3@0: ''' 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: '''