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