lidar_correction_ghk_pollyxt.py

Thu, 10 Nov 2016 16:35:00 +0200

author
Iannis <ulalume3@yahoo.com>
date
Thu, 10 Nov 2016 16:35:00 +0200
changeset 3
d55a57e732be
parent 1
5dffdb7caec9
child 11
453b23dd7f94
permissions
-rw-r--r--

* Initial readme file.
* Minor doc changes in example system setting file.

ulalume3@0 1 # -*- coding: utf-8 -*-
ulalume3@0 2 """
ulalume3@0 3 Copyright 2016 Volker Freudenthaler
ulalume3@0 4
ulalume3@0 5 Licensed under the EUPL, Version 1.1 only (the "Licence").
ulalume3@0 6
ulalume3@0 7 You may not use this work except in compliance with the Licence.
ulalume3@0 8 A copy of the Licence is distributed with the code. Alternatively, you may obtain
ulalume3@0 9 a copy of the Licence at:
ulalume3@0 10
ulalume3@0 11 https://joinup.ec.europa.eu/community/eupl/og_page/eupl
ulalume3@0 12
ulalume3@0 13 Unless required by applicable law or agreed to in writing, software distributed
ulalume3@0 14 under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR CONDITIONS
ulalume3@0 15 OF ANY KIND, either express or implied. See the Licence for the specific language governing
ulalume3@0 16 permissions and limitations under the Licence.
ulalume3@0 17
ulalume3@0 18 Equation reference: http://www.atmos-meas-tech-discuss.net/amt-2015-338/amt-2015-338.pdf
ulalume3@0 19 With equations code from Appendix C
ulalume3@0 20 Python 3.4.2
ulalume3@0 21 """
ulalume3@0 22 #!/usr/bin/env python3
ulalume3@0 23 from __future__ import print_function
ulalume3@0 24 #import math
ulalume3@0 25 import numpy as np
ulalume3@0 26 import sys
ulalume3@0 27 import os
ulalume3@0 28 #import seaborn as sns
ulalume3@0 29 import matplotlib.pyplot as plt
ulalume3@0 30 from time import clock
ulalume3@0 31
ulalume3@0 32 #from matplotlib.backends.backend_pdf import PdfPages
ulalume3@0 33 #pdffile = '{}.pdf'.format('path')
ulalume3@0 34 #pp = PdfPages(pdffile)
ulalume3@0 35 ## pp.savefig can be called multiple times to save to multiple pages
ulalume3@0 36 #pp.savefig()
ulalume3@0 37 #pp.close()
ulalume3@0 38
ulalume3@0 39 from contextlib import contextmanager
ulalume3@0 40 @contextmanager
ulalume3@0 41 def redirect_stdout(new_target):
ulalume3@0 42 old_target, sys.stdout = sys.stdout, new_target # replace sys.stdout
ulalume3@0 43 try:
ulalume3@0 44 yield new_target # run some code with the replaced stdout
ulalume3@0 45 finally:
ulalume3@0 46 sys.stdout.flush()
ulalume3@0 47 sys.stdout = old_target # restore to the previous value
ulalume3@0 48 '''
ulalume3@0 49 real_raw_input = vars(__builtins__).get('raw_input',input)
ulalume3@0 50 '''
ulalume3@0 51 try:
ulalume3@0 52 import __builtin__
ulalume3@0 53 input = getattr(__builtin__, 'raw_input')
ulalume3@0 54 except (ImportError, AttributeError):
ulalume3@0 55 pass
ulalume3@0 56
ulalume3@0 57 from distutils.util import strtobool
ulalume3@0 58 def user_yes_no_query(question):
ulalume3@0 59 sys.stdout.write('%s [y/n]\n' % question)
ulalume3@0 60 while True:
ulalume3@0 61 try:
ulalume3@0 62 return strtobool(input().lower())
ulalume3@0 63 except ValueError:
ulalume3@0 64 sys.stdout.write('Please respond with \'y\' or \'n\'.\n')
ulalume3@0 65
ulalume3@0 66 #if user_yes_no_query('want to exit?') == 1: sys.exit()
ulalume3@0 67
ulalume3@0 68 '''
ulalume3@0 69 ## {{{ http://code.activestate.com/recipes/577058/ (r2)
ulalume3@0 70 def query_yes_no(question, default="yes"):
ulalume3@0 71 valid = {"yes":"yes", "y":"yes", "ye":"yes",
ulalume3@0 72 "no":"no", "n":"no"}
ulalume3@0 73 if default == None:
ulalume3@0 74 prompt = " [y/n] "
ulalume3@0 75 elif default == "yes":
ulalume3@0 76 prompt = " [Y/n] "
ulalume3@0 77 elif default == "no":
ulalume3@0 78 prompt = " [y/N] "
ulalume3@0 79 else:
ulalume3@0 80 raise ValueError("invalid default answer: '%s'" % default)
ulalume3@0 81
ulalume3@0 82 while 1:
ulalume3@0 83 sys.stdout.write(question + prompt)
ulalume3@0 84 choice = input().lower()
ulalume3@0 85 if default is not None and choice == '':
ulalume3@0 86 return default
ulalume3@0 87 elif choice in valid.keys():
ulalume3@0 88 return valid[choice]
ulalume3@0 89 else:
ulalume3@0 90 sys.stdout.write("Please respond with 'yes' or 'no' "\
ulalume3@0 91 "(or 'y' or 'n').\n")
ulalume3@0 92 ## end of http://code.activestate.com/recipes/577058/ }}}
ulalume3@0 93 '''
ulalume3@0 94 abspath = os.path.abspath(__file__)
ulalume3@0 95 dname = os.path.dirname(abspath)
ulalume3@0 96 fname = os.path.basename(abspath)
ulalume3@0 97 os.chdir(dname)
ulalume3@0 98
ulalume3@0 99 #PrintToOutputFile = True
ulalume3@0 100
ulalume3@0 101 sqr05 = 0.5**0.5
ulalume3@0 102
ulalume3@0 103 # ---- Initial definition of variables; the actual values will be read in with exec(open('./optic_input.py').read()) below
ulalume3@0 104 LID = "internal"
ulalume3@0 105 EID = "internal"
ulalume3@0 106 # --- IL Laser IL and +-Uncertainty
ulalume3@0 107 bL = 1. #degree of linear polarization; default 1
ulalume3@0 108 RotL, dRotL, nRotL = 0.0, 0.0, 1 #alpha; rotation of laser polarization in degrees; default 0
ulalume3@0 109 # --- ME Emitter and +-Uncertainty
ulalume3@0 110 DiE, dDiE, nDiE = 0., 0.00, 1 # Diattenuation
ulalume3@0 111 TiE = 1. # Unpolarized transmittance
ulalume3@0 112 RetE, dRetE, nRetE = 0., 180.0, 0 # Retardance in degrees
ulalume3@0 113 RotE, dRotE, nRotE = 0., 0.0, 0 # beta: Rotation of optical element in degrees
ulalume3@0 114 # --- MO Receiver Optics including telescope
ulalume3@0 115 DiO, dDiO, nDiO = -0.055, 0.003, 1
ulalume3@0 116 TiO = 0.9
ulalume3@0 117 RetO, dRetO, nRetO = 0., 180.0, 2
ulalume3@0 118 RotO, dRotO, nRotO = 0., 0.1, 1 #gamma
ulalume3@0 119 # --- PBS MT transmitting path defined with (TS,TP); and +-Uncertainty
ulalume3@0 120 TP, dTP, nTP = 0.98, 0.02, 1
ulalume3@0 121 TS, dTS, nTS = 0.001, 0.001, 1
ulalume3@0 122 TiT = 0.5 * (TP + TS)
ulalume3@0 123 DiT = (TP-TS)/(TP+TS)
ulalume3@0 124 # PolFilter
ulalume3@0 125 RetT, dRetT, nRetT = 0., 180., 0
ulalume3@0 126 ERaT, dERaT, nERaT = 0.001, 0.001, 1
ulalume3@0 127 RotaT, dRotaT, nRotaT = 0., 3., 1
ulalume3@0 128 DaT = (1-ERaT)/(1+ERaT)
ulalume3@0 129 TaT = 0.5*(1+ERaT)
ulalume3@0 130 # --- PBS MR reflecting path defined with (RS,RP); and +-Uncertainty
ulalume3@0 131 RS, dRS, nRS = 1 - TS, 0., 0
ulalume3@0 132 RP, dRP, nRP = 1 - TP, 0., 0
ulalume3@0 133 TiR = 0.5 * (RP + RS)
ulalume3@0 134 DiR = (RP-RS)/(RP+RS)
ulalume3@0 135 # PolFilter
ulalume3@0 136 RetR, dRetR, nRetR = 0., 180., 0
ulalume3@0 137 ERaR, dERaR, nERaR = 0.001, 0.001, 1
ulalume3@0 138 RotaR,dRotaR,nRotaR = 90., 3., 1
ulalume3@0 139 DaR = (1-ERaR)/(1+ERaR)
ulalume3@0 140 TaR = 0.5*(1+ERaR)
ulalume3@0 141
ulalume3@0 142 # Parellel signal detected in the transmitted channel => Y = 1, or in the reflected channel => Y = -1
ulalume3@0 143 Y = -1.
ulalume3@0 144
ulalume3@0 145 # Calibrator = type defined by matrix values
ulalume3@0 146 LocC = 4 # location of calibrator: behind laser = 1; behind emitter = 2; before receiver = 3; before PBS = 4
ulalume3@0 147
ulalume3@0 148 TypeC = 3 # linear polarizer calibrator
ulalume3@0 149 # example with extinction ratio 0.001
ulalume3@0 150 DiC, dDiC, nDiC = 1.0, 0., 0 # ideal 1.0
ulalume3@0 151 TiC = 0.5 # ideal 0.5
ulalume3@0 152 RetC, dRetC, nRetC = 0., 0., 0
ulalume3@0 153 RotC, dRotC, nRotC = 0.0, 0.1, 0 #constant calibrator offset epsilon
ulalume3@0 154 RotationErrorEpsilonForNormalMeasurements = False # is in general False for TypeC == 3 calibrator
ulalume3@0 155
ulalume3@0 156 # Rotation error without calibrator: if False, then epsilon = 0 for normal measurements
ulalume3@0 157 RotationErrorEpsilonForNormalMeasurements = True
ulalume3@0 158
ulalume3@0 159 # LDRCal assumed atmospheric linear depolarization ratio during the calibration measurements (first guess)
ulalume3@0 160 LDRCal0,dLDRCal,nLDRCal= 0.25, 0.04, 1
ulalume3@0 161 LDRCal = LDRCal0
ulalume3@0 162 # measured LDRm will be corrected with calculated parameters
ulalume3@0 163 LDRmeas = 0.015
ulalume3@0 164 # LDRtrue for simulation of measurement => LDRsim
ulalume3@0 165 LDRtrue = 0.5
ulalume3@0 166 LDRtrue2 = 0.004
ulalume3@0 167
ulalume3@0 168 # Initialize other values to 0
ulalume3@0 169 ER, nER, dER = 0.001, 0, 0.001
ulalume3@0 170 K = 0.
ulalume3@0 171 Km = 0.
ulalume3@0 172 Kp = 0.
ulalume3@0 173 LDRcorr = 0.
ulalume3@0 174 Eta = 0.
ulalume3@0 175 Ir = 0.
ulalume3@0 176 It = 0.
ulalume3@0 177 h = 1.
ulalume3@0 178
ulalume3@0 179 Loc = ['', 'behind laser', 'behind emitter', 'before receiver', 'before PBS']
ulalume3@0 180 Type = ['', 'mechanical rotator', 'hwp rotator', 'linear polarizer', 'qwp rotator', 'circular polarizer', 'real HWP +-22.5°']
ulalume3@0 181 dY = ['reflected channel', '', 'transmitted channel']
ulalume3@0 182
ulalume3@0 183 # end of initial definition of variables
ulalume3@0 184 # *******************************************************************************************************************************
ulalume3@0 185
ulalume3@0 186 # --- Read actual lidar system parameters from ./optic_input.py (must be in the same directory)
ulalume3@0 187
ulalume3@0 188 #InputFile = 'optic_input_ver6e_PollyXTSea.py'
ulalume3@0 189 #InputFile = 'optic_input_ver6e_PollyXTSea_JA.py'
ulalume3@0 190 #InputFile = 'optic_input_ver6e_PollyXT_RALPH.py'
ulalume3@0 191 #InputFile = 'optic_input_ver8c_PollyXT_RALPH.py'
ulalume3@0 192 #InputFile = 'optic_input_ver8c_PollyXT_RALPH_2.py'
ulalume3@0 193 #InputFile = 'optic_input_ver8c_PollyXT_RALPH_3.py'
ulalume3@0 194 #InputFile = 'optic_input_ver8c_PollyXT_RALPH_4.py'
ulalume3@0 195 #InputFile = 'optic_input_ver8c_PollyXT_RALPH_5.py'
ulalume3@0 196 #InputFile = 'optic_input_ver8c_PollyXT_RALPH_6.py'
ulalume3@0 197 InputFile = 'optic_input_ver8c_PollyXT_RALPH_7.py'
ulalume3@0 198 #InputFile = 'optic_input_ver6e_Bertha_b_355.py'
ulalume3@0 199 #InputFile = 'optic_input_ver6e_Bertha_b_532.py'
ulalume3@0 200 #InputFile = 'optic_input_ver6e_Bertha_b_1064.py'
ulalume3@0 201
ulalume3@0 202 '''
ulalume3@0 203 print("From ", dname)
ulalume3@0 204 print("Running ", fname)
ulalume3@0 205 print("Reading input file ", InputFile, " for")
ulalume3@0 206 '''
ulalume3@0 207 # this works with Python 2 - and 3?
ulalume3@0 208 exec(open('./'+InputFile).read(), globals())
ulalume3@0 209 # end of read actual system parameters
ulalume3@0 210
ulalume3@0 211 # --- Manual Parameter Change ---
ulalume3@0 212 # (use for quick parameter changes without changing the input file )
ulalume3@0 213 #DiO = 0.
ulalume3@0 214 #LDRtrue = 0.45
ulalume3@0 215 #LDRtrue2 = 0.004
ulalume3@0 216 #Y = -1
ulalume3@0 217 #LocC = 4 #location of calibrator: 1 = behind laser; 2 = behind emitter; 3 = before receiver; 4 = before PBS
ulalume3@0 218 ##TypeC = 6 Don't change the TypeC here
ulalume3@0 219 #RotationErrorEpsilonForNormalMeasurements = True
ulalume3@0 220 #LDRCal = 0.25
ulalume3@0 221 #bL = 0.8
ulalume3@0 222 ## --- Errors
ulalume3@0 223 RotL0, dRotL, nRotL = RotL, dRotL, nRotL
ulalume3@0 224
ulalume3@0 225 DiE0, dDiE, nDiE = DiE, dDiE, nDiE
ulalume3@0 226 RetE0, dRetE, nRetE = RetE, dRetE, nRetE
ulalume3@0 227 RotE0, dRotE, nRotE = RotE, dRotE, nRotE
ulalume3@0 228
ulalume3@0 229 DiO0, dDiO, nDiO = DiO, dDiO, nDiO
ulalume3@0 230 RetO0, dRetO, nRetO = RetO, dRetO, nRetO
ulalume3@0 231 RotO0, dRotO, nRotO = RotO, dRotO, nRotO
ulalume3@0 232
ulalume3@0 233 DiC0, dDiC, nDiC = DiC, dDiC, nDiC
ulalume3@0 234 RetC0, dRetC, nRetC = RetC, dRetC, nRetC
ulalume3@0 235 RotC0, dRotC, nRotC = RotC, dRotC, nRotC
ulalume3@0 236
ulalume3@0 237 TP0, dTP, nTP = TP, dTP, nTP
ulalume3@0 238 TS0, dTS, nTS = TS, dTS, nTS
ulalume3@0 239 RetT0, dRetT, nRetT = RetT, dRetT, nRetT
ulalume3@0 240
ulalume3@0 241 ERaT0, dERaT, nERaT = ERaT, dERaT, nERaT
ulalume3@0 242 RotaT0,dRotaT,nRotaT= RotaT,dRotaT,nRotaT
ulalume3@0 243
ulalume3@0 244 RP0, dRP, nRP = RP, dRP, nRP
ulalume3@0 245 RS0, dRS, nRS = RS, dRS, nRS
ulalume3@0 246 RetR0, dRetR, nRetR = RetR, dRetR, nRetR
ulalume3@0 247
ulalume3@0 248 ERaR0, dERaR, nERaR = ERaR, dERaR, nERaR
ulalume3@0 249 RotaR0,dRotaR,nRotaR= RotaR,dRotaR,nRotaR
ulalume3@0 250
ulalume3@0 251 LDRCal0,dLDRCal,nLDRCal=LDRCal,dLDRCal,nLDRCal
ulalume3@0 252 #LDRCal0,dLDRCal,nLDRCal=LDRCal,dLDRCal,0
ulalume3@0 253 # ---------- End of manual parameter change
ulalume3@0 254
ulalume3@0 255 RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC = RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0
ulalume3@0 256 TP, TS, RP, RS, ERaT, RotaT, RetT, ERaR, RotaR, RetR = TP0, TS0, RP0, RS0 , ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0
ulalume3@0 257 LDRCal = LDRCal0
ulalume3@0 258 DTa0, TTa0, DRa0, TRa0, LDRsimx, LDRCorr = 0,0,0,0,0,0
ulalume3@0 259
ulalume3@0 260 TiT = 0.5 * (TP + TS)
ulalume3@0 261 DiT = (TP-TS)/(TP+TS)
ulalume3@0 262 ZiT = (1. - DiT**2)**0.5
ulalume3@0 263 TiR = 0.5 * (RP + RS)
ulalume3@0 264 DiR = (RP-RS)/(RP+RS)
ulalume3@0 265 ZiR = (1. - DiR**2)**0.5
ulalume3@0 266
ulalume3@0 267 # --------------------------------------------------------
ulalume3@0 268 def Calc(RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC, TP, TS, RP, RS, ERaT, RotaT, RetT, ERaR, RotaR, RetR, LDRCal):
ulalume3@0 269 # ---- Do the calculations of bra-ket vectors
ulalume3@0 270 h = -1. if TypeC == 2 else 1
ulalume3@0 271 # from input file: assumed LDRCal for calibration measurements
ulalume3@0 272 aCal = (1.-LDRCal)/(1+LDRCal)
ulalume3@0 273 # from input file: measured LDRm and true LDRtrue, LDRtrue2 =>
ulalume3@0 274 #ameas = (1.-LDRmeas)/(1+LDRmeas)
ulalume3@0 275 atrue = (1.-LDRtrue)/(1+LDRtrue)
ulalume3@0 276 #atrue2 = (1.-LDRtrue2)/(1+LDRtrue2)
ulalume3@0 277
ulalume3@0 278 # angles of emitter and laser and calibrator and receiver optics
ulalume3@0 279 # RotL = alpha, RotE = beta, RotO = gamma, RotC = epsilon
ulalume3@0 280 S2a = np.sin(2*np.deg2rad(RotL))
ulalume3@0 281 C2a = np.cos(2*np.deg2rad(RotL))
ulalume3@0 282 S2b = np.sin(2*np.deg2rad(RotE))
ulalume3@0 283 C2b = np.cos(2*np.deg2rad(RotE))
ulalume3@0 284 S2ab = np.sin(np.deg2rad(2*RotL-2*RotE))
ulalume3@0 285 C2ab = np.cos(np.deg2rad(2*RotL-2*RotE))
ulalume3@0 286 S2g = np.sin(np.deg2rad(2*RotO))
ulalume3@0 287 C2g = np.cos(np.deg2rad(2*RotO))
ulalume3@0 288
ulalume3@0 289 # Laser with Degree of linear polarization DOLP = bL
ulalume3@0 290 IinL = 1.
ulalume3@0 291 QinL = bL
ulalume3@0 292 UinL = 0.
ulalume3@0 293 VinL = (1. - bL**2)**0.5
ulalume3@0 294
ulalume3@0 295 # Stokes Input Vector rotation Eq. E.4
ulalume3@0 296 A = C2a*QinL - S2a*UinL
ulalume3@0 297 B = S2a*QinL + C2a*UinL
ulalume3@0 298 # Stokes Input Vector rotation Eq. E.9
ulalume3@0 299 C = C2ab*QinL - S2ab*UinL
ulalume3@0 300 D = S2ab*QinL + C2ab*UinL
ulalume3@0 301
ulalume3@0 302 # emitter optics
ulalume3@0 303 CosE = np.cos(np.deg2rad(RetE))
ulalume3@0 304 SinE = np.sin(np.deg2rad(RetE))
ulalume3@0 305 ZiE = (1. - DiE**2)**0.5
ulalume3@0 306 WiE = (1. - ZiE*CosE)
ulalume3@0 307
ulalume3@0 308 # Stokes Input Vector after emitter optics equivalent to Eq. E.9 with already rotated input vector from Eq. E.4
ulalume3@0 309 # b = beta
ulalume3@0 310 IinE = (IinL + DiE*C)
ulalume3@0 311 QinE = (C2b*DiE*IinL + A + S2b*(WiE*D - ZiE*SinE*VinL))
ulalume3@0 312 UinE = (S2b*DiE*IinL + B - C2b*(WiE*D - ZiE*SinE*VinL))
ulalume3@0 313 VinE = (-ZiE*SinE*D + ZiE*CosE*VinL)
ulalume3@0 314
ulalume3@0 315 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F)
ulalume3@0 316 IinF = IinE
ulalume3@0 317 QinF = aCal*QinE
ulalume3@0 318 UinF = -aCal*UinE
ulalume3@0 319 VinF = (1.-2.*aCal)*VinE
ulalume3@0 320
ulalume3@0 321 # receiver optics
ulalume3@0 322 CosO = np.cos(np.deg2rad(RetO))
ulalume3@0 323 SinO = np.sin(np.deg2rad(RetO))
ulalume3@0 324 ZiO = (1. - DiO**2)**0.5
ulalume3@0 325 WiO = (1. - ZiO*CosO)
ulalume3@0 326
ulalume3@0 327 # calibrator
ulalume3@0 328 CosC = np.cos(np.deg2rad(RetC))
ulalume3@0 329 SinC = np.sin(np.deg2rad(RetC))
ulalume3@0 330 ZiC = (1. - DiC**2)**0.5
ulalume3@0 331 WiC = (1. - ZiC*CosC)
ulalume3@0 332
ulalume3@0 333 # Stokes Input Vector before the polarising beam splitter Eq. E.31
ulalume3@0 334 A = C2g*QinE - S2g*UinE
ulalume3@0 335 B = S2g*QinE + C2g*UinE
ulalume3@0 336
ulalume3@0 337 IinP = (IinE + DiO*aCal*A)
ulalume3@0 338 QinP = (C2g*DiO*IinE + aCal*QinE - S2g*(WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE))
ulalume3@0 339 UinP = (S2g*DiO*IinE - aCal*UinE + C2g*(WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE))
ulalume3@0 340 VinP = (ZiO*SinO*aCal*B + ZiO*CosO*(1-2*aCal)*VinE)
ulalume3@0 341
ulalume3@0 342 #-------------------------
ulalume3@0 343 # F11 assuemd to be = 1 => measured: F11m = IinP / IinE with atrue
ulalume3@0 344 #F11sim = TiO*(IinE + DiO*atrue*A)/IinE
ulalume3@0 345 #-------------------------
ulalume3@0 346
ulalume3@0 347 # For PollyXT
ulalume3@0 348 # analyser
ulalume3@0 349 RS = 1 - TS
ulalume3@0 350 RP = 1 - TP
ulalume3@0 351
ulalume3@0 352 TiT = 0.5 * (TP + TS)
ulalume3@0 353 DiT = (TP-TS)/(TP+TS)
ulalume3@0 354 ZiT = (1. - DiT**2)**0.5
ulalume3@0 355 TiR = 0.5 * (RP + RS)
ulalume3@0 356 DiR = (RP-RS)/(RP+RS)
ulalume3@0 357 ZiR = (1. - DiR**2)**0.5
ulalume3@0 358 CosT = np.cos(np.deg2rad(RetT))
ulalume3@0 359 SinT = np.sin(np.deg2rad(RetT))
ulalume3@0 360 CosR = np.cos(np.deg2rad(RetR))
ulalume3@0 361 SinR = np.sin(np.deg2rad(RetR))
ulalume3@0 362
ulalume3@0 363 DaT = (1-ERaT)/(1+ERaT)
ulalume3@0 364 DaR = (1-ERaR)/(1+ERaR)
ulalume3@0 365 TaT = 0.5*(1+ERaT)
ulalume3@0 366 TaR = 0.5*(1+ERaR)
ulalume3@0 367
ulalume3@0 368 S2aT = np.sin(np.deg2rad(h*2*RotaT))
ulalume3@0 369 C2aT = np.cos(np.deg2rad(2*RotaT))
ulalume3@0 370 S2aR = np.sin(np.deg2rad(h*2*RotaR))
ulalume3@0 371 C2aR = np.cos(np.deg2rad(2*RotaR))
ulalume3@0 372
ulalume3@0 373 # Aanalyzer As before the PBS Eq. D.5
ulalume3@0 374 ATP1 = (1+C2aT*DaT*DiT)
ulalume3@0 375 ATP2 = Y*(DiT+C2aT*DaT)
ulalume3@0 376 ATP3 = Y*S2aT*DaT*ZiT*CosT
ulalume3@0 377 ATP4 = S2aT*DaT*ZiT*SinT
ulalume3@0 378 ATP = np.array([ATP1,ATP2,ATP3,ATP4])
ulalume3@0 379
ulalume3@0 380 ARP1 = (1+C2aR*DaR*DiR)
ulalume3@0 381 ARP2 = Y*(DiR+C2aR*DaR)
ulalume3@0 382 ARP3 = Y*S2aR*DaR*ZiR*CosR
ulalume3@0 383 ARP4 = S2aR*DaR*ZiR*SinR
ulalume3@0 384 ARP = np.array([ARP1,ARP2,ARP3,ARP4])
ulalume3@0 385
ulalume3@0 386 DTa = ATP2*Y/ATP1
ulalume3@0 387 DRa = ARP2*Y/ARP1
ulalume3@0 388
ulalume3@0 389 # ---- Calculate signals and correction parameters for diffeent locations and calibrators
ulalume3@0 390 if LocC == 4: # Calibrator before the PBS
ulalume3@0 391 #print("Calibrator location not implemented yet")
ulalume3@0 392
ulalume3@0 393 #S2ge = np.sin(np.deg2rad(2*RotO + h*2*RotC))
ulalume3@0 394 #C2ge = np.cos(np.deg2rad(2*RotO + h*2*RotC))
ulalume3@0 395 S2e = np.sin(np.deg2rad(h*2*RotC))
ulalume3@0 396 C2e = np.cos(np.deg2rad(2*RotC))
ulalume3@0 397 # rotated AinP by epsilon Eq. C.3
ulalume3@0 398 ATP2e = C2e*ATP2 + S2e*ATP3
ulalume3@0 399 ATP3e = C2e*ATP3 - S2e*ATP2
ulalume3@0 400 ARP2e = C2e*ARP2 + S2e*ARP3
ulalume3@0 401 ARP3e = C2e*ARP3 - S2e*ARP2
ulalume3@0 402 ATPe = np.array([ATP1,ATP2e,ATP3e,ATP4])
ulalume3@0 403 ARPe = np.array([ARP1,ARP2e,ARP3e,ARP4])
ulalume3@0 404 # Stokes Input Vector before the polarising beam splitter Eq. E.31
ulalume3@0 405 A = C2g*QinE - S2g*UinE
ulalume3@0 406 B = S2g*QinE + C2g*UinE
ulalume3@0 407 #C = (WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE)
ulalume3@0 408 Co = ZiO*SinO*VinE
ulalume3@0 409 Ca = (WiO*B - 2*ZiO*SinO*VinE)
ulalume3@0 410 #C = Co + aCal*Ca
ulalume3@0 411 #IinP = (IinE + DiO*aCal*A)
ulalume3@0 412 #QinP = (C2g*DiO*IinE + aCal*QinE - S2g*C)
ulalume3@0 413 #UinP = (S2g*DiO*IinE - aCal*UinE + C2g*C)
ulalume3@0 414 #VinP = (ZiO*SinO*aCal*B + ZiO*CosO*(1-2*aCal)*VinE)
ulalume3@0 415 IinPo = IinE
ulalume3@0 416 QinPo = (C2g*DiO*IinE - S2g*Co)
ulalume3@0 417 UinPo = (S2g*DiO*IinE + C2g*Co)
ulalume3@0 418 VinPo = ZiO*CosO*VinE
ulalume3@0 419
ulalume3@0 420 IinPa = DiO*A
ulalume3@0 421 QinPa = QinE - S2g*Ca
ulalume3@0 422 UinPa = -UinE + C2g*Ca
ulalume3@0 423 VinPa = ZiO*(SinO*B - 2*CosO*VinE)
ulalume3@0 424
ulalume3@0 425 IinP = IinPo + aCal*IinPa
ulalume3@0 426 QinP = QinPo + aCal*QinPa
ulalume3@0 427 UinP = UinPo + aCal*UinPa
ulalume3@0 428 VinP = VinPo + aCal*VinPa
ulalume3@0 429 # Stokes Input Vector before the polarising beam splitter rotated by epsilon Eq. C.3
ulalume3@0 430 #QinPe = C2e*QinP + S2e*UinP
ulalume3@0 431 #UinPe = C2e*UinP - S2e*QinP
ulalume3@0 432 QinPoe = C2e*QinPo + S2e*UinPo
ulalume3@0 433 UinPoe = C2e*UinPo - S2e*QinPo
ulalume3@0 434 QinPae = C2e*QinPa + S2e*UinPa
ulalume3@0 435 UinPae = C2e*UinPa - S2e*QinPa
ulalume3@0 436 QinPe = C2e*QinP + S2e*UinP
ulalume3@0 437 UinPe = C2e*UinP - S2e*QinP
ulalume3@0 438
ulalume3@0 439 # Calibration signals and Calibration correction K from measurements with LDRCal / aCal
ulalume3@0 440 if (TypeC == 2) or (TypeC == 1): # rotator calibration Eq. C.4
ulalume3@0 441 # parameters for calibration with aCal
ulalume3@0 442 AT = ATP1*IinP + h*ATP4*VinP
ulalume3@0 443 BT = ATP3e*QinP - h*ATP2e*UinP
ulalume3@0 444 AR = ARP1*IinP + h*ARP4*VinP
ulalume3@0 445 BR = ARP3e*QinP - h*ARP2e*UinP
ulalume3@0 446 # Correction paremeters for normal measurements; they are independent of LDR
ulalume3@0 447 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out
ulalume3@0 448 IS1 = np.array([IinPo,QinPo,UinPo,VinPo])
ulalume3@0 449 IS2 = np.array([IinPa,QinPa,UinPa,VinPa])
ulalume3@0 450 GT = np.dot(ATP,IS1)
ulalume3@0 451 GR = np.dot(ARP,IS1)
ulalume3@0 452 HT = np.dot(ATP,IS2)
ulalume3@0 453 HR = np.dot(ARP,IS2)
ulalume3@0 454 else:
ulalume3@0 455 IS1 = np.array([IinPo,QinPo,UinPo,VinPo])
ulalume3@0 456 IS2 = np.array([IinPa,QinPa,UinPa,VinPa])
ulalume3@0 457 GT = np.dot(ATPe,IS1)
ulalume3@0 458 GR = np.dot(ARPe,IS1)
ulalume3@0 459 HT = np.dot(ATPe,IS2)
ulalume3@0 460 HR = np.dot(ARPe,IS2)
ulalume3@0 461 elif (TypeC == 3) or (TypeC == 4): # linear polariser calibration Eq. C.5
ulalume3@0 462 # parameters for calibration with aCal
ulalume3@0 463 AT = ATP1*IinP + ATP3e*UinPe + ZiC*CosC*(ATP2e*QinPe + ATP4*VinP)
ulalume3@0 464 BT = DiC*(ATP1*UinPe + ATP3e*IinP) - ZiC*SinC*(ATP2e*VinP - ATP4*QinPe)
ulalume3@0 465 AR = ARP1*IinP + ARP3e*UinPe + ZiC*CosC*(ARP2e*QinPe + ARP4*VinP)
ulalume3@0 466 BR = DiC*(ARP1*UinPe + ARP3e*IinP) - ZiC*SinC*(ARP2e*VinP - ARP4*QinPe)
ulalume3@0 467 # Correction paremeters for normal measurements; they are independent of LDR
ulalume3@0 468 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out
ulalume3@0 469 IS1 = np.array([IinPo,QinPo,UinPo,VinPo])
ulalume3@0 470 IS2 = np.array([IinPa,QinPa,UinPa,VinPa])
ulalume3@0 471 GT = np.dot(ATP,IS1)
ulalume3@0 472 GR = np.dot(ARP,IS1)
ulalume3@0 473 HT = np.dot(ATP,IS2)
ulalume3@0 474 HR = np.dot(ARP,IS2)
ulalume3@0 475 else:
ulalume3@0 476 IS1e = np.array([IinPo+DiC*QinPoe,DiC*IinPo+QinPoe,ZiC*(CosC*UinPoe+SinC*VinPo),-ZiC*(SinC*UinPoe-CosC*VinPo)])
ulalume3@0 477 IS2e = np.array([IinPa+DiC*QinPae,DiC*IinPa+QinPae,ZiC*(CosC*UinPae+SinC*VinPa),-ZiC*(SinC*UinPae-CosC*VinPa)])
ulalume3@0 478 GT = np.dot(ATPe,IS1e)
ulalume3@0 479 GR = np.dot(ARPe,IS1e)
ulalume3@0 480 HT = np.dot(ATPe,IS2e)
ulalume3@0 481 HR = np.dot(ARPe,IS2e)
ulalume3@0 482 elif (TypeC == 6): # diattenuator calibration +-22.5° rotated_diattenuator_X22x5deg.odt
ulalume3@0 483 # parameters for calibration with aCal
ulalume3@0 484 AT = ATP1*IinP + sqr05*DiC*(ATP1*QinPe + ATP2e*IinP) + (1-0.5*WiC)*(ATP2e*QinPe + ATP3e*UinPe) + ZiC*(sqr05*SinC*(ATP3e*VinP-ATP4*UinPe) + ATP4*CosC*VinP)
ulalume3@0 485 BT = sqr05*DiC*(ATP1*UinPe + ATP3e*IinP) + 0.5*WiC*(ATP2e*UinPe + ATP3e*QinPe) - sqr05*ZiC*SinC*(ATP2e*VinP - ATP4*QinPe)
ulalume3@0 486 AR = ARP1*IinP + sqr05*DiC*(ARP1*QinPe + ARP2e*IinP) + (1-0.5*WiC)*(ARP2e*QinPe + ARP3e*UinPe) + ZiC*(sqr05*SinC*(ARP3e*VinP-ARP4*UinPe) + ARP4*CosC*VinP)
ulalume3@0 487 BR = sqr05*DiC*(ARP1*UinPe + ARP3e*IinP) + 0.5*WiC*(ARP2e*UinPe + ARP3e*QinPe) - sqr05*ZiC*SinC*(ARP2e*VinP - ARP4*QinPe)
ulalume3@0 488 # Correction paremeters for normal measurements; they are independent of LDR
ulalume3@0 489 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out
ulalume3@0 490 IS1 = np.array([IinPo,QinPo,UinPo,VinPo])
ulalume3@0 491 IS2 = np.array([IinPa,QinPa,UinPa,VinPa])
ulalume3@0 492 GT = np.dot(ATP,IS1)
ulalume3@0 493 GR = np.dot(ARP,IS1)
ulalume3@0 494 HT = np.dot(ATP,IS2)
ulalume3@0 495 HR = np.dot(ARP,IS2)
ulalume3@0 496 else:
ulalume3@0 497 IS1e = np.array([IinPo+DiC*QinPoe,DiC*IinPo+QinPoe,ZiC*(CosC*UinPoe+SinC*VinPo),-ZiC*(SinC*UinPoe-CosC*VinPo)])
ulalume3@0 498 IS2e = np.array([IinPa+DiC*QinPae,DiC*IinPa+QinPae,ZiC*(CosC*UinPae+SinC*VinPa),-ZiC*(SinC*UinPae-CosC*VinPa)])
ulalume3@0 499 GT = np.dot(ATPe,IS1e)
ulalume3@0 500 GR = np.dot(ARPe,IS1e)
ulalume3@0 501 HT = np.dot(ATPe,IS2e)
ulalume3@0 502 HR = np.dot(ARPe,IS2e)
ulalume3@0 503 else:
ulalume3@0 504 print("Calibrator not implemented yet")
ulalume3@0 505 sys.exit()
ulalume3@0 506
ulalume3@0 507 elif LocC == 3: # C before receiver optics Eq.57
ulalume3@0 508
ulalume3@0 509 #S2ge = np.sin(np.deg2rad(2*RotO - 2*RotC))
ulalume3@0 510 #C2ge = np.cos(np.deg2rad(2*RotO - 2*RotC))
ulalume3@0 511 S2e = np.sin(np.deg2rad(2*RotC))
ulalume3@0 512 C2e = np.cos(np.deg2rad(2*RotC))
ulalume3@0 513
ulalume3@0 514 # As with C before the receiver optics (rotated_diattenuator_X22x5deg.odt)
ulalume3@0 515 AF1 = np.array([1,C2g*DiO,S2g*DiO,0])
ulalume3@0 516 AF2 = np.array([C2g*DiO,1-S2g**2*WiO,S2g*C2g*WiO,-S2g*ZiO*SinO])
ulalume3@0 517 AF3 = np.array([S2g*DiO,S2g*C2g*WiO,1-C2g**2*WiO,C2g*ZiO*SinO])
ulalume3@0 518 AF4 = np.array([0,S2g*SinO,-C2g*SinO,CosO])
ulalume3@0 519
ulalume3@0 520 ATF = (ATP1*AF1+ATP2*AF2+ATP3*AF3+ATP4*AF4)
ulalume3@0 521 ARF = (ARP1*AF1+ARP2*AF2+ARP3*AF3+ARP4*AF4)
ulalume3@0 522 ATF2 = ATF[1]
ulalume3@0 523 ATF3 = ATF[2]
ulalume3@0 524 ARF2 = ARF[1]
ulalume3@0 525 ARF3 = ARF[2]
ulalume3@0 526
ulalume3@0 527 # rotated AinF by epsilon
ulalume3@0 528 ATF1 = ATF[0]
ulalume3@0 529 ATF4 = ATF[3]
ulalume3@0 530 ATF2e = C2e*ATF[1] + S2e*ATF[2]
ulalume3@0 531 ATF3e = C2e*ATF[2] - S2e*ATF[1]
ulalume3@0 532 ARF1 = ARF[0]
ulalume3@0 533 ARF4 = ARF[3]
ulalume3@0 534 ARF2e = C2e*ARF[1] + S2e*ARF[2]
ulalume3@0 535 ARF3e = C2e*ARF[2] - S2e*ARF[1]
ulalume3@0 536
ulalume3@0 537 ATFe = np.array([ATF1,ATF2e,ATF3e,ATF4])
ulalume3@0 538 ARFe = np.array([ARF1,ARF2e,ARF3e,ARF4])
ulalume3@0 539
ulalume3@0 540 QinEe = C2e*QinE + S2e*UinE
ulalume3@0 541 UinEe = C2e*UinE - S2e*QinE
ulalume3@0 542
ulalume3@0 543 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F)
ulalume3@0 544 IinF = IinE
ulalume3@0 545 QinF = aCal*QinE
ulalume3@0 546 UinF = -aCal*UinE
ulalume3@0 547 VinF = (1.-2.*aCal)*VinE
ulalume3@0 548
ulalume3@0 549 IinFo = IinE
ulalume3@0 550 QinFo = 0.
ulalume3@0 551 UinFo = 0.
ulalume3@0 552 VinFo = VinE
ulalume3@0 553
ulalume3@0 554 IinFa = 0.
ulalume3@0 555 QinFa = QinE
ulalume3@0 556 UinFa = -UinE
ulalume3@0 557 VinFa = -2.*VinE
ulalume3@0 558
ulalume3@0 559 # Stokes Input Vector before receiver optics rotated by epsilon Eq. C.3
ulalume3@0 560 QinFe = C2e*QinF + S2e*UinF
ulalume3@0 561 UinFe = C2e*UinF - S2e*QinF
ulalume3@0 562 QinFoe = C2e*QinFo + S2e*UinFo
ulalume3@0 563 UinFoe = C2e*UinFo - S2e*QinFo
ulalume3@0 564 QinFae = C2e*QinFa + S2e*UinFa
ulalume3@0 565 UinFae = C2e*UinFa - S2e*QinFa
ulalume3@0 566
ulalume3@0 567 # Calibration signals and Calibration correction K from measurements with LDRCal / aCal
ulalume3@0 568 if (TypeC == 2) or (TypeC == 1): # rotator calibration Eq. C.4
ulalume3@0 569 # parameters for calibration with aCal
ulalume3@0 570 AT = ATF1*IinF + ATF4*h*VinF
ulalume3@0 571 BT = ATF3e*QinF - ATF2e*h*UinF
ulalume3@0 572 AR = ARF1*IinF + ARF4*h*VinF
ulalume3@0 573 BR = ARF3e*QinF - ARF2e*h*UinF
ulalume3@0 574 # Correction paremeters for normal measurements; they are independent of LDR
ulalume3@0 575 if (not RotationErrorEpsilonForNormalMeasurements):
ulalume3@0 576 GT = ATF1*IinE + ATF4*VinE
ulalume3@0 577 GR = ARF1*IinE + ARF4*VinE
ulalume3@0 578 HT = ATF2*QinE - ATF3*UinE - ATF4*2*VinE
ulalume3@0 579 HR = ARF2*QinE - ARF3*UinE - ARF4*2*VinE
ulalume3@0 580 else:
ulalume3@0 581 GT = ATF1*IinE + ATF4*h*VinE
ulalume3@0 582 GR = ARF1*IinE + ARF4*h*VinE
ulalume3@0 583 HT = ATF2e*QinE - ATF3e*h*UinE - ATF4*h*2*VinE
ulalume3@0 584 HR = ARF2e*QinE - ARF3e*h*UinE - ARF4*h*2*VinE
ulalume3@0 585 elif (TypeC == 3) or (TypeC == 4): # linear polariser calibration Eq. C.5
ulalume3@0 586 # p = +45°, m = -45°
ulalume3@0 587 IF1e = np.array([IinF, ZiC*CosC*QinFe, UinFe, ZiC*CosC*VinF])
ulalume3@0 588 IF2e = np.array([DiC*UinFe, -ZiC*SinC*VinF, DiC*IinF, ZiC*SinC*QinFe])
ulalume3@0 589 AT = np.dot(ATFe,IF1e)
ulalume3@0 590 AR = np.dot(ARFe,IF1e)
ulalume3@0 591 BT = np.dot(ATFe,IF2e)
ulalume3@0 592 BR = np.dot(ARFe,IF2e)
ulalume3@0 593
ulalume3@0 594 # Correction paremeters for normal measurements; they are independent of LDR --- the same as for TypeC = 6
ulalume3@0 595 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out
ulalume3@0 596 IS1 = np.array([IinE,0,0,VinE])
ulalume3@0 597 IS2 = np.array([0,QinE,-UinE,-2*VinE])
ulalume3@0 598 GT = np.dot(ATF,IS1)
ulalume3@0 599 GR = np.dot(ARF,IS1)
ulalume3@0 600 HT = np.dot(ATF,IS2)
ulalume3@0 601 HR = np.dot(ARF,IS2)
ulalume3@0 602 else:
ulalume3@0 603 IS1e = np.array([IinFo+DiC*QinFoe,DiC*IinFo+QinFoe,ZiC*(CosC*UinFoe+SinC*VinFo),-ZiC*(SinC*UinFoe-CosC*VinFo)])
ulalume3@0 604 IS2e = np.array([IinFa+DiC*QinFae,DiC*IinFa+QinFae,ZiC*(CosC*UinFae+SinC*VinFa),-ZiC*(SinC*UinFae-CosC*VinFa)])
ulalume3@0 605 GT = np.dot(ATFe,IS1e)
ulalume3@0 606 GR = np.dot(ARFe,IS1e)
ulalume3@0 607 HT = np.dot(ATFe,IS2e)
ulalume3@0 608 HR = np.dot(ARFe,IS2e)
ulalume3@0 609
ulalume3@0 610 elif (TypeC == 6): # diattenuator calibration +-22.5° rotated_diattenuator_X22x5deg.odt
ulalume3@0 611 # parameters for calibration with aCal
ulalume3@0 612 IF1e = np.array([IinF+sqr05*DiC*QinFe, sqr05*DiC*IinF+(1-0.5*WiC)*QinFe, (1-0.5*WiC)*UinFe+sqr05*ZiC*SinC*VinF, -sqr05*ZiC*SinC*UinFe+ZiC*CosC*VinF])
ulalume3@0 613 IF2e = np.array([sqr05*DiC*UinFe, 0.5*WiC*UinFe-sqr05*ZiC*SinC*VinF, sqr05*DiC*IinF+0.5*WiC*QinFe, sqr05*ZiC*SinC*QinFe])
ulalume3@0 614 AT = np.dot(ATFe,IF1e)
ulalume3@0 615 AR = np.dot(ARFe,IF1e)
ulalume3@0 616 BT = np.dot(ATFe,IF2e)
ulalume3@0 617 BR = np.dot(ARFe,IF2e)
ulalume3@0 618
ulalume3@0 619 # Correction paremeters for normal measurements; they are independent of LDR
ulalume3@0 620 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out
ulalume3@0 621 #IS1 = np.array([IinE,0,0,VinE])
ulalume3@0 622 #IS2 = np.array([0,QinE,-UinE,-2*VinE])
ulalume3@0 623 IS1 = np.array([IinFo,0,0,VinFo])
ulalume3@0 624 IS2 = np.array([0,QinFa,UinFa,VinFa])
ulalume3@0 625 GT = np.dot(ATF,IS1)
ulalume3@0 626 GR = np.dot(ARF,IS1)
ulalume3@0 627 HT = np.dot(ATF,IS2)
ulalume3@0 628 HR = np.dot(ARF,IS2)
ulalume3@0 629 else:
ulalume3@0 630 IS1e = np.array([IinFo+DiC*QinFoe,DiC*IinFo+QinFoe,ZiC*(CosC*UinFoe+SinC*VinFo),-ZiC*(SinC*UinFoe-CosC*VinFo)])
ulalume3@0 631 IS2e = np.array([IinFa+DiC*QinFae,DiC*IinFa+QinFae,ZiC*(CosC*UinFae+SinC*VinFa),-ZiC*(SinC*UinFae-CosC*VinFa)])
ulalume3@0 632 #IS1e = np.array([IinFo,0,0,VinFo])
ulalume3@0 633 #IS2e = np.array([0,QinFae,UinFae,VinFa])
ulalume3@0 634 GT = np.dot(ATFe,IS1e)
ulalume3@0 635 GR = np.dot(ARFe,IS1e)
ulalume3@0 636 HT = np.dot(ATFe,IS2e)
ulalume3@0 637 HR = np.dot(ARFe,IS2e)
ulalume3@0 638
ulalume3@0 639 else:
ulalume3@0 640 print('Calibrator not implemented yet')
ulalume3@0 641 sys.exit()
ulalume3@0 642
ulalume3@0 643 elif LocC == 2: # C behind emitter optics Eq.57 -------------------------------------------------------
ulalume3@0 644 #print("Calibrator location not implemented yet")
ulalume3@0 645 S2e = np.sin(np.deg2rad(2*RotC))
ulalume3@0 646 C2e = np.cos(np.deg2rad(2*RotC))
ulalume3@0 647
ulalume3@0 648 # AS with C before the receiver optics (see document rotated_diattenuator_X22x5deg.odt)
ulalume3@0 649 AF1 = np.array([1,C2g*DiO,S2g*DiO,0])
ulalume3@0 650 AF2 = np.array([C2g*DiO,1-S2g**2*WiO,S2g*C2g*WiO,-S2g*ZiO*SinO])
ulalume3@0 651 AF3 = np.array([S2g*DiO, S2g*C2g*WiO, 1-C2g**2*WiO, C2g*ZiO*SinO])
ulalume3@0 652 AF4 = np.array([0, S2g*SinO, -C2g*SinO, CosO])
ulalume3@0 653
ulalume3@0 654 ATF = (ATP1*AF1+ATP2*AF2+ATP3*AF3+ATP4*AF4)
ulalume3@0 655 ARF = (ARP1*AF1+ARP2*AF2+ARP3*AF3+ARP4*AF4)
ulalume3@0 656 ATF1 = ATF[0]
ulalume3@0 657 ATF2 = ATF[1]
ulalume3@0 658 ATF3 = ATF[2]
ulalume3@0 659 ATF4 = ATF[3]
ulalume3@0 660 ARF1 = ARF[0]
ulalume3@0 661 ARF2 = ARF[1]
ulalume3@0 662 ARF3 = ARF[2]
ulalume3@0 663 ARF4 = ARF[3]
ulalume3@0 664
ulalume3@0 665 # AS with C behind the emitter
ulalume3@0 666 # terms without aCal
ulalume3@0 667 ATE1o, ARE1o = ATF1, ARF1
ulalume3@0 668 ATE2o, ARE2o = 0., 0.
ulalume3@0 669 ATE3o, ARE3o = 0., 0.
ulalume3@0 670 ATE4o, ARE4o = ATF4, ARF4
ulalume3@0 671 # terms with aCal
ulalume3@0 672 ATE1a, ARE1a = 0. , 0.
ulalume3@0 673 ATE2a, ARE2a = ATF2, ARF2
ulalume3@0 674 ATE3a, ARE3a = -ATF3, -ARF3
ulalume3@0 675 ATE4a, ARE4a = -2*ATF4, -2*ARF4
ulalume3@0 676 # rotated AinEa by epsilon
ulalume3@0 677 ATE2ae = C2e*ATF2 + S2e*ATF3
ulalume3@0 678 ATE3ae = -S2e*ATF2 - C2e*ATF3
ulalume3@0 679 ARE2ae = C2e*ARF2 + S2e*ARF3
ulalume3@0 680 ARE3ae = -S2e*ARF2 - C2e*ARF3
ulalume3@0 681
ulalume3@0 682 ATE1 = ATE1o
ulalume3@0 683 ATE2e = aCal*ATE2ae
ulalume3@0 684 ATE3e = aCal*ATE3ae
ulalume3@0 685 ATE4 = (1-2*aCal)*ATF4
ulalume3@0 686 ARE1 = ARE1o
ulalume3@0 687 ARE2e = aCal*ARE2ae
ulalume3@0 688 ARE3e = aCal*ARE3ae
ulalume3@0 689 ARE4 = (1-2*aCal)*ARF4
ulalume3@0 690
ulalume3@0 691 # rotated IinE
ulalume3@0 692 QinEe = C2e*QinE + S2e*UinE
ulalume3@0 693 UinEe = C2e*UinE - S2e*QinE
ulalume3@0 694
ulalume3@0 695 # Calibration signals and Calibration correction K from measurements with LDRCal / aCal
ulalume3@0 696 if (TypeC == 2) or (TypeC == 1): # +++++++++ rotator calibration Eq. C.4
ulalume3@0 697 AT = ATE1o*IinE + (ATE4o+aCal*ATE4a)*h*VinE
ulalume3@0 698 BT = aCal * (ATE3ae*QinEe - ATE2ae*h*UinEe)
ulalume3@0 699 AR = ARE1o*IinE + (ARE4o+aCal*ARE4a)*h*VinE
ulalume3@0 700 BR = aCal * (ARE3ae*QinEe - ARE2ae*h*UinEe)
ulalume3@0 701
ulalume3@0 702 # Correction paremeters for normal measurements; they are independent of LDR
ulalume3@0 703 if (not RotationErrorEpsilonForNormalMeasurements):
ulalume3@0 704 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F)
ulalume3@0 705 GT = ATE1o*IinE + ATE4o*h*VinE
ulalume3@0 706 GR = ARE1o*IinE + ARE4o*h*VinE
ulalume3@0 707 HT = ATE2a*QinE + ATE3a*h*UinEe + ATE4a*h*VinE
ulalume3@0 708 HR = ARE2a*QinE + ARE3a*h*UinEe + ARE4a*h*VinE
ulalume3@0 709 else:
ulalume3@0 710 GT = ATE1o*IinE + ATE4o*h*VinE
ulalume3@0 711 GR = ARE1o*IinE + ARE4o*h*VinE
ulalume3@0 712 HT = ATE2ae*QinE + ATE3ae*h*UinEe + ATE4a*h*VinE
ulalume3@0 713 HR = ARE2ae*QinE + ARE3ae*h*UinEe + ARE4a*h*VinE
ulalume3@0 714
ulalume3@0 715 elif (TypeC == 3) or (TypeC == 4): # +++++++++ linear polariser calibration Eq. C.5
ulalume3@0 716 # p = +45°, m = -45°
ulalume3@0 717 AT = ATE1*IinE + ZiC*CosC*(ATE2e*QinEe + ATE4*VinE) + ATE3e*UinEe
ulalume3@0 718 BT = DiC*(ATE1*UinEe + ATE3e*IinE) + ZiC*SinC*(ATE4*QinEe - ATE2e*VinE)
ulalume3@0 719 AR = ARE1*IinE + ZiC*CosC*(ARE2e*QinEe + ARE4*VinE) + ARE3e*UinEe
ulalume3@0 720 BR = DiC*(ARE1*UinEe + ARE3e*IinE) + ZiC*SinC*(ARE4*QinEe - ARE2e*VinE)
ulalume3@0 721
ulalume3@0 722 # Correction paremeters for normal measurements; they are independent of LDR
ulalume3@0 723 if (not RotationErrorEpsilonForNormalMeasurements):
ulalume3@0 724 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F)
ulalume3@0 725 GT = ATE1o*IinE + ATE4o*VinE
ulalume3@0 726 GR = ARE1o*IinE + ARE4o*VinE
ulalume3@0 727 HT = ATE2a*QinE + ATE3a*UinE + ATE4a*VinE
ulalume3@0 728 HR = ARE2a*QinE + ARE3a*UinE + ARE4a*VinE
ulalume3@0 729 else:
ulalume3@0 730 D = IinE + DiC*QinEe
ulalume3@0 731 A = DiC*IinE + QinEe
ulalume3@0 732 B = ZiC*(CosC*UinEe + SinC*VinE)
ulalume3@0 733 C = -ZiC*(SinC*UinEe - CosC*VinE)
ulalume3@0 734 GT = ATE1o*D + ATE4o*C
ulalume3@0 735 GR = ARE1o*D + ARE4o*C
ulalume3@0 736 HT = ATE2a*A + ATE3a*B + ATE4a*C
ulalume3@0 737 HR = ARE2a*A + ARE3a*B + ARE4a*C
ulalume3@0 738
ulalume3@0 739 elif (TypeC == 6): # real HWP calibration +-22.5° rotated_diattenuator_X22x5deg.odt
ulalume3@0 740 # p = +22.5°, m = -22.5°
ulalume3@0 741 IE1e = np.array([IinE+sqr05*DiC*QinEe, sqr05*DiC*IinE+(1-0.5*WiC)*QinEe, (1-0.5*WiC)*UinEe+sqr05*ZiC*SinC*VinE, -sqr05*ZiC*SinC*UinEe+ZiC*CosC*VinE])
ulalume3@0 742 IE2e = np.array([sqr05*DiC*UinEe, 0.5*WiC*UinEe-sqr05*ZiC*SinC*VinE, sqr05*DiC*IinE+0.5*WiC*QinEe, sqr05*ZiC*SinC*QinEe])
ulalume3@0 743 ATEe = np.array([ATE1,ATE2e,ATE3e,ATE4])
ulalume3@0 744 AREe = np.array([ARE1,ARE2e,ARE3e,ARE4])
ulalume3@0 745 AT = np.dot(ATEe,IE1e)
ulalume3@0 746 AR = np.dot(AREe,IE1e)
ulalume3@0 747 BT = np.dot(ATEe,IE2e)
ulalume3@0 748 BR = np.dot(AREe,IE2e)
ulalume3@0 749
ulalume3@0 750 # Correction paremeters for normal measurements; they are independent of LDR
ulalume3@0 751 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out
ulalume3@0 752 GT = ATE1o*IinE + ATE4o*VinE
ulalume3@0 753 GR = ARE1o*IinE + ARE4o*VinE
ulalume3@0 754 HT = ATE2a*QinE + ATE3a*UinE + ATE4a*VinE
ulalume3@0 755 HR = ARE2a*QinE + ARE3a*UinE + ARE4a*VinE
ulalume3@0 756 else:
ulalume3@0 757 D = IinE + DiC*QinEe
ulalume3@0 758 A = DiC*IinE + QinEe
ulalume3@0 759 B = ZiC*(CosC*UinEe + SinC*VinE)
ulalume3@0 760 C = -ZiC*(SinC*UinEe - CosC*VinE)
ulalume3@0 761 GT = ATE1o*D + ATE4o*C
ulalume3@0 762 GR = ARE1o*D + ARE4o*C
ulalume3@0 763 HT = ATE2a*A + ATE3a*B + ATE4a*C
ulalume3@0 764 HR = ARE2a*A + ARE3a*B + ARE4a*C
ulalume3@0 765
ulalume3@0 766 else:
ulalume3@0 767 print('Calibrator not implemented yet')
ulalume3@0 768 sys.exit()
ulalume3@0 769
ulalume3@0 770 else:
ulalume3@0 771 print("Calibrator location not implemented yet")
ulalume3@0 772 sys.exit()
ulalume3@0 773
ulalume3@0 774 # Determination of the correction K of the calibration factor
ulalume3@0 775 IoutTp = TaT*TiT*TiO*TiE*(AT + BT)
ulalume3@0 776 IoutTm = TaT*TiT*TiO*TiE*(AT - BT)
ulalume3@0 777 IoutRp = TaR*TiR*TiO*TiE*(AR + BR)
ulalume3@0 778 IoutRm = TaR*TiR*TiO*TiE*(AR - BR)
ulalume3@0 779
ulalume3@0 780 # --- Results and Corrections; electronic etaR and etaT are assumed to be 1
ulalume3@0 781 Etapx = IoutRp/IoutTp
ulalume3@0 782 Etamx = IoutRm/IoutTm
ulalume3@0 783 Etax = (Etapx*Etamx)**0.5
ulalume3@0 784
ulalume3@0 785 Eta = (TaR*TiR)/(TaT*TiT) # Eta = Eta*/K Eq. 84
ulalume3@0 786 K = Etax / Eta
ulalume3@0 787
ulalume3@0 788 # For comparison with Volkers Libreoffice Müller Matrix spreadsheet
ulalume3@0 789 #Eta_test_p = (IoutRp/IoutTp)
ulalume3@0 790 #Eta_test_m = (IoutRm/IoutTm)
ulalume3@0 791 #Eta_test = (Eta_test_p*Eta_test_m)**0.5
ulalume3@0 792
ulalume3@0 793 # ----- Forward simulated signals and LDRsim with atrue; from input file
ulalume3@0 794 It = TaT*TiT*TiO*TiE*(GT+atrue*HT)
ulalume3@0 795 Ir = TaR*TiR*TiO*TiE*(GR+atrue*HR)
ulalume3@0 796 # LDRsim = 1/Eta*Ir/It # simulated LDR* with Y from input file
ulalume3@0 797 LDRsim = Ir/It # simulated uncorrected LDR with Y from input file
ulalume3@0 798 # Corrected LDRsimCorr from forward simulated LDRsim (atrue)
ulalume3@0 799 # LDRsimCorr = (1./Eta*LDRsim*(GT+HT)-(GR+HR))/((GR-HR)-1./Eta*LDRsim*(GT-HT))
ulalume3@0 800 if Y == -1.:
ulalume3@0 801 LDRsimx = 1./LDRsim
ulalume3@0 802 else:
ulalume3@0 803 LDRsimx = LDRsim
ulalume3@0 804
ulalume3@0 805 # The following is correct without doubt
ulalume3@0 806 #LDRCorr = (LDRsim*K/Etax*(GT+HT)-(GR+HR))/((GR-HR)-LDRsim*K/Etax*(GT-HT))
ulalume3@0 807
ulalume3@0 808 # The following is a test whether the equations for calibration Etax and normal signal (GHK, LDRsim) are consistent
ulalume3@0 809 LDRCorr = (LDRsim/Eta*(GT+HT)-(GR+HR))/((GR-HR)-LDRsim*K/Etax*(GT-HT))
ulalume3@0 810
ulalume3@0 811 TTa = TiT*TaT #*ATP1
ulalume3@0 812 TRa = TiR*TaR #*ARP1
ulalume3@0 813
ulalume3@0 814 F11sim = 1/(TiO*TiE)*((HR*Etax/K*It/TTa-HT*Ir/TRa)/(HR*GT-HT*GR)) # IL = 1, Etat = Etar = 1
ulalume3@0 815
ulalume3@0 816 return (GT, HT, GR, HR, K, Eta, LDRsimx, LDRCorr, DTa, DRa, TTa, TRa, F11sim)
ulalume3@0 817 # *******************************************************************************************************************************
ulalume3@0 818
ulalume3@0 819 # --- CALC truth
ulalume3@0 820 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0)
ulalume3@0 821
ulalume3@0 822 # --------------------------------------------------------
ulalume3@0 823 with open('output_' + LID + '.dat', 'w') as f:
ulalume3@0 824 with redirect_stdout(f):
ulalume3@0 825 print("From ", dname)
ulalume3@0 826 print("Running ", fname)
ulalume3@0 827 print("Reading input file ", InputFile) #, " for Lidar system :", EID, ", ", LID)
ulalume3@0 828 print("for Lidar system: ", EID, ", ", LID)
ulalume3@0 829 # --- Print iput information*********************************
ulalume3@0 830 print(" --- Input parameters: value ±error / ±steps ----------------------")
ulalume3@0 831 print("{0:8} {1:8} {2:8.5f}; {3:8} {4:7.4f}±{5:7.4f}/{6:2d}".format("Laser: ", "DOLP = ", bL, " rotation alpha = ", RotL0, dRotL, nRotL))
ulalume3@0 832 print(" Diatt., Tunpol, Retard., Rotation (deg)")
ulalume3@0 833 print("{0:12} {1:7.4f}±{2:7.4f}/{8:2d}, {3:7.4f}, {4:3.0f}±{5:3.0f}/{9:2d}, {6:7.4f}±{7:7.4f}/{10:2d}".format("Emitter ", DiE0, dDiE, TiE, RetE0, dRetE, RotE0, dRotE, nDiE, nRetE, nRotE))
ulalume3@0 834 print("{0:12} {1:7.4f}±{2:7.4f}/{8:2d}, {3:7.4f}, {4:3.0f}±{5:3.0f}/{9:2d}, {6:7.4f}±{7:7.4f}/{10:2d}".format("Receiver ", DiO0, dDiO, TiO, RetO0, dRetO, RotO0, dRotO, nDiO, nRetO, nRotO))
ulalume3@0 835 print("{0:12} {1:7.4f}±{2:7.4f}/{8:2d}, {3:7.4f}, {4:3.0f}±{5:3.0f}/{9:2d}, {6:7.4f}±{7:7.4f}/{10:2d}".format("Calibrator ", DiC0, dDiC, TiC, RetC0, dRetC, RotC0, dRotC, nDiC, nRetC, nRotC))
ulalume3@0 836 print("{0:12}".format(" --- Pol.-filter ---"))
ulalume3@0 837 print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format("ERT, ERR :", ERaT0, dERaT, nERaT, ERaR0, dERaR, nERaR))
ulalume3@0 838 print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format("RotaT , RotaR :", RotaT0, dRotaT, nRotaT, RotaR0,dRotaR,nRotaR))
ulalume3@0 839 print("{0:12}".format(" --- PBS ---"))
ulalume3@0 840 print("{0:12}{1:7.4f}±{2:7.4f}/{9:2d}, {3:7.4f}±{4:7.4f}/{10:2d}, {5:7.4f}±{6:7.4f}/{11:2d},{7:7.4f}±{8:7.4f}/{12:2d}".format("TP,TS,RP,RS :", TP0, dTP, TS0, dTS, RP0, dRP, RS0, dRS, nTP, nTS, nRP, nRS))
ulalume3@0 841 print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}, {5:1.0f}".format("DT,TT,DR,TR,Y :", DiT, TiT, DiR, TiR, Y))
ulalume3@0 842 print("{0:12}".format(" --- Combined PBS + Pol.-filter ---"))
ulalume3@0 843 print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}".format("DTa,TTa,DRa,TRa: ", DTa0, TTa0, DRa0, TRa0))
ulalume3@0 844 print()
ulalume3@0 845 print("Rotation Error Epsilon For Normal Measurements = ", RotationErrorEpsilonForNormalMeasurements)
ulalume3@0 846 #print ('LocC = ', LocC, Loc[LocC], '; TypeC = ',TypeC, Type[TypeC])
ulalume3@0 847 print(Type[TypeC], Loc[LocC], "; Parallel signal detected in", dY[int(Y+1)])
ulalume3@0 848 # end of print actual system parameters
ulalume3@0 849 # ******************************************************************************
ulalume3@0 850
ulalume3@0 851 #print()
ulalume3@0 852 #print(" --- LDRCal during calibration | simulated and corrected LDRs -------------")
ulalume3@0 853 #print("{0:8} |{1:8}->{2:8},{3:9}->{4:9} |{5:8}->{6:8}".format(" LDRCal"," LDRtrue", " LDRsim"," LDRtrue2", " LDRsim2", " LDRmeas", " LDRcorr"))
ulalume3@0 854 #print("{0:8.5f} |{1:8.5f}->{2:8.5f},{3:9.5f}->{4:9.5f} |{5:8.5f}->{6:8.5f}".format(LDRCal, LDRtrue, LDRsim, LDRtrue2, LDRsim2, LDRmeas, LDRCorr))
ulalume3@0 855 #print("{0:8} |{1:8}->{2:8}->{3:8}".format(" LDRCal"," LDRtrue", " LDRsimx", " LDRcorr"))
ulalume3@0 856 #print("{0:6.3f}±{1:5.3f}/{2:2d}|{3:8.5f}->{4:8.5f}->{5:8.5f}".format(LDRCal0, dLDRCal, nLDRCal, LDRtrue, LDRsimx, LDRCorr))
ulalume3@0 857 #print("{0:8} |{1:8}->{2:8}->{3:8}".format(" LDRCal"," LDRtrue", " LDRsimx", " LDRcorr"))
ulalume3@0 858 #print(" --- LDRCal during calibration")
ulalume3@0 859 print("{0:26}: {1:6.3f}±{2:5.3f}/{3:2d}".format("LDRCal during calibration", LDRCal0, dLDRCal, nLDRCal))
ulalume3@0 860
ulalume3@0 861 #print("{0:8}={1:8.5f};{2:8}={3:8.5f}".format(" IinP",IinP," F11sim",F11sim))
ulalume3@0 862 print()
ulalume3@0 863
ulalume3@0 864 K0List = np.zeros(3)
ulalume3@0 865 LDRsimxList = np.zeros(3)
ulalume3@0 866 LDRCalList = 0.004, 0.2, 0.45
ulalume3@0 867 for i,LDRCal in enumerate(LDRCalList):
ulalume3@0 868 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal)
ulalume3@0 869 K0List[i] = K0
ulalume3@0 870 LDRsimxList[i] = LDRsimx
ulalume3@0 871
ulalume3@0 872 print("{0:8},{1:8},{2:8},{3:8},{4:9},{5:9},{6:9}".format(" GR", " GT", " HR", " HT", " K(0.004)", " K(0.2)", " K(0.45)"))
ulalume3@0 873 print("{0:8.5f},{1:8.5f},{2:8.5f},{3:8.5f},{4:9.5f},{5:9.5f},{6:9.5f}".format(GR0, GT0, HR0, HT0, K0List[0], K0List[1], K0List[2]))
ulalume3@0 874 print('========================================================================')
ulalume3@0 875
ulalume3@0 876 print("{0:9},{1:9},{2:9}".format(" LDRtrue", " LDRsimx", " LDRCorr"))
ulalume3@0 877 LDRtrueList = 0.004, 0.02, 0.2, 0.45
ulalume3@0 878 for i,LDRtrue in enumerate(LDRtrueList):
ulalume3@0 879 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0)
ulalume3@0 880 print("{0:9.5f},{1:9.5f},{2:9.5f}".format(LDRtrue, LDRsimx, LDRCorr))
ulalume3@0 881
ulalume3@0 882
ulalume3@0 883 file = open('output_' + LID + '.dat', 'r')
ulalume3@0 884 print (file.read())
ulalume3@0 885 file.close()
ulalume3@0 886
ulalume3@0 887 '''
ulalume3@0 888 if(PrintToOutputFile):
ulalume3@0 889 f = open('output_ver7.dat', 'w')
ulalume3@0 890 old_target = sys.stdout
ulalume3@0 891 sys.stdout = f
ulalume3@0 892
ulalume3@0 893 print("something")
ulalume3@0 894
ulalume3@0 895 if(PrintToOutputFile):
ulalume3@0 896 sys.stdout.flush()
ulalume3@0 897 f.close
ulalume3@0 898 sys.stdout = old_target
ulalume3@0 899 '''
ulalume3@0 900 # --- CALC again truth with LDRCal0 to reset all 0-values
ulalume3@0 901 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0)
ulalume3@0 902
ulalume3@0 903 # --- Start Errors calculation
ulalume3@0 904
ulalume3@0 905 iN = -1
ulalume3@0 906 N = ((nRotL*2+1)*
ulalume3@0 907 (nRotE*2+1)*(nRetE*2+1)*(nDiE*2+1)*
ulalume3@0 908 (nRotO*2+1)*(nRetO*2+1)*(nDiO*2+1)*
ulalume3@0 909 (nRotC*2+1)*(nRetC*2+1)*(nDiC*2+1)*
ulalume3@0 910 (nTP*2+1)*(nTS*2+1)*(nRP*2+1)*(nRS*2+1)*(nERaT*2+1)*(nERaR*2+1)*
ulalume3@0 911 (nRotaT*2+1)*(nRotaR*2+1)*(nRetT*2+1)*(nRetR*2+1)*(nLDRCal*2+1))
ulalume3@0 912 print("N = ",N ," ", end="")
ulalume3@0 913
ulalume3@0 914 if N > 1e6:
ulalume3@0 915 if user_yes_no_query('Warning: processing ' + str(N) + ' samples will take very long. Do you want to proceed?') == 0: sys.exit()
ulalume3@0 916 if N > 5e6:
ulalume3@0 917 if user_yes_no_query('Warning: the memory required for ' + str(N) + ' samples might be ' + '{0:5.1f}'.format(N/4e6) + ' GB. Do you anyway want to proceed?') == 0: sys.exit()
ulalume3@0 918
ulalume3@0 919 #if user_yes_no_query('Warning: processing' + str(N) + ' samples will take very long. Do you want to proceed?') == 0: sys.exit()
ulalume3@0 920
ulalume3@0 921 # --- Arrays for plotting ------
ulalume3@0 922 LDRmin = np.zeros(5)
ulalume3@0 923 LDRmax = np.zeros(5)
ulalume3@0 924 F11min = np.zeros(5)
ulalume3@0 925 F11max = np.zeros(5)
ulalume3@0 926
ulalume3@0 927 LDRrange = np.zeros(5)
ulalume3@0 928 LDRrange = 0.004, 0.02, 0.1, 0.3, 0.45
ulalume3@0 929 #aLDRsimx = np.zeros(N)
ulalume3@0 930 #aLDRsimx2 = np.zeros(N)
ulalume3@0 931 #aLDRcorr = np.zeros(N)
ulalume3@0 932 #aLDRcorr2 = np.zeros(N)
ulalume3@0 933 aERaT = np.zeros(N)
ulalume3@0 934 aERaR = np.zeros(N)
ulalume3@0 935 aRotaT = np.zeros(N)
ulalume3@0 936 aRotaR = np.zeros(N)
ulalume3@0 937 aRetT = np.zeros(N)
ulalume3@0 938 aRetR = np.zeros(N)
ulalume3@0 939 aTP = np.zeros(N)
ulalume3@0 940 aTS = np.zeros(N)
ulalume3@0 941 aRP = np.zeros(N)
ulalume3@0 942 aRS = np.zeros(N)
ulalume3@0 943 aDiE = np.zeros(N)
ulalume3@0 944 aDiO = np.zeros(N)
ulalume3@0 945 aDiC = np.zeros(N)
ulalume3@0 946 aRotC = np.zeros(N)
ulalume3@0 947 aRetC = np.zeros(N)
ulalume3@0 948 aRotL = np.zeros(N)
ulalume3@0 949 aRetE = np.zeros(N)
ulalume3@0 950 aRotE = np.zeros(N)
ulalume3@0 951 aRetO = np.zeros(N)
ulalume3@0 952 aRotO = np.zeros(N)
ulalume3@0 953 aLDRCal = np.zeros(N)
ulalume3@0 954 aA = np.zeros((5,N))
ulalume3@0 955 aX = np.zeros((5,N))
ulalume3@0 956 aF11corr = np.zeros((5,N))
ulalume3@0 957
ulalume3@0 958 atime = clock()
ulalume3@0 959 dtime = clock()
ulalume3@0 960
ulalume3@0 961 # --- Calc Error signals
ulalume3@0 962 #GT, HT, GR, HR, K, Eta, LDRsim = Calc(RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC, TP, TS)
ulalume3@0 963 # ---- Do the calculations of bra-ket vectors
ulalume3@0 964 h = -1. if TypeC == 2 else 1
ulalume3@0 965
ulalume3@0 966 # from input file: measured LDRm and true LDRtrue, LDRtrue2 =>
ulalume3@0 967 ameas = (1.-LDRmeas)/(1+LDRmeas)
ulalume3@0 968 atrue = (1.-LDRtrue)/(1+LDRtrue)
ulalume3@0 969 atrue2 = (1.-LDRtrue2)/(1+LDRtrue2)
ulalume3@0 970
ulalume3@0 971 for iLDRCal in range(-nLDRCal,nLDRCal+1):
ulalume3@0 972 # from input file: assumed LDRCal for calibration measurements
ulalume3@0 973 LDRCal = LDRCal0
ulalume3@0 974 if nLDRCal > 0: LDRCal = LDRCal0 + iLDRCal*dLDRCal/nLDRCal
ulalume3@0 975
ulalume3@0 976 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal)
ulalume3@0 977 aCal = (1.-LDRCal)/(1+LDRCal)
ulalume3@0 978 for iRotL, iRotE, iRetE, iDiE \
ulalume3@0 979 in [(iRotL,iRotE,iRetE,iDiE)
ulalume3@0 980 for iRotL in range(-nRotL,nRotL+1)
ulalume3@0 981 for iRotE in range(-nRotE,nRotE+1)
ulalume3@0 982 for iRetE in range(-nRetE,nRetE+1)
ulalume3@0 983 for iDiE in range(-nDiE,nDiE+1)]:
ulalume3@0 984
ulalume3@0 985 if nRotL > 0: RotL = RotL0 + iRotL*dRotL/nRotL
ulalume3@0 986 if nRotE > 0: RotE = RotE0 + iRotE*dRotE/nRotE
ulalume3@0 987 if nRetE > 0: RetE = RetE0 + iRetE*dRetE/nRetE
ulalume3@0 988 if nDiE > 0: DiE = DiE0 + iDiE*dDiE/nDiE
ulalume3@0 989
ulalume3@0 990 # angles of emitter and laser and calibrator and receiver optics
ulalume3@0 991 # RotL = alpha, RotE = beta, RotO = gamma, RotC = epsilon
ulalume3@0 992 S2a = np.sin(2*np.deg2rad(RotL))
ulalume3@0 993 C2a = np.cos(2*np.deg2rad(RotL))
ulalume3@0 994 S2b = np.sin(2*np.deg2rad(RotE))
ulalume3@0 995 C2b = np.cos(2*np.deg2rad(RotE))
ulalume3@0 996 S2ab = np.sin(np.deg2rad(2*RotL-2*RotE))
ulalume3@0 997 C2ab = np.cos(np.deg2rad(2*RotL-2*RotE))
ulalume3@0 998
ulalume3@0 999 # Laser with Degree of linear polarization DOLP = bL
ulalume3@0 1000 IinL = 1.
ulalume3@0 1001 QinL = bL
ulalume3@0 1002 UinL = 0.
ulalume3@0 1003 VinL = (1. - bL**2)**0.5
ulalume3@0 1004
ulalume3@0 1005 # Stokes Input Vector rotation Eq. E.4
ulalume3@0 1006 A = C2a*QinL - S2a*UinL
ulalume3@0 1007 B = S2a*QinL + C2a*UinL
ulalume3@0 1008 # Stokes Input Vector rotation Eq. E.9
ulalume3@0 1009 C = C2ab*QinL - S2ab*UinL
ulalume3@0 1010 D = S2ab*QinL + C2ab*UinL
ulalume3@0 1011
ulalume3@0 1012 # emitter optics
ulalume3@0 1013 CosE = np.cos(np.deg2rad(RetE))
ulalume3@0 1014 SinE = np.sin(np.deg2rad(RetE))
ulalume3@0 1015 ZiE = (1. - DiE**2)**0.5
ulalume3@0 1016 WiE = (1. - ZiE*CosE)
ulalume3@0 1017
ulalume3@0 1018 # Stokes Input Vector after emitter optics equivalent to Eq. E.9 with already rotated input vector from Eq. E.4
ulalume3@0 1019 # b = beta
ulalume3@0 1020 IinE = (IinL + DiE*C)
ulalume3@0 1021 QinE = (C2b*DiE*IinL + A + S2b*(WiE*D - ZiE*SinE*VinL))
ulalume3@0 1022 UinE = (S2b*DiE*IinL + B - C2b*(WiE*D - ZiE*SinE*VinL))
ulalume3@0 1023 VinE = (-ZiE*SinE*D + ZiE*CosE*VinL)
ulalume3@0 1024
ulalume3@0 1025 #-------------------------
ulalume3@0 1026 # F11 assuemd to be = 1 => measured: F11m = IinP / IinE with atrue
ulalume3@0 1027 #F11sim = (IinE + DiO*atrue*(C2g*QinE - S2g*UinE))/IinE
ulalume3@0 1028 #-------------------------
ulalume3@0 1029
ulalume3@0 1030 for iRotO, iRetO, iDiO, iRotC, iRetC, iDiC, iTP, iTS, iRP, iRS, iERaT, iRotaT, iRetT, iERaR, iRotaR, iRetR \
ulalume3@0 1031 in [(iRotO,iRetO,iDiO,iRotC,iRetC,iDiC,iTP,iTS,iRP,iRS,iERaT,iRotaT,iRetT,iERaR,iRotaR,iRetR )
ulalume3@0 1032 for iRotO in range(-nRotO,nRotO+1)
ulalume3@0 1033 for iRetO in range(-nRetO,nRetO+1)
ulalume3@0 1034 for iDiO in range(-nDiO,nDiO+1)
ulalume3@0 1035 for iRotC in range(-nRotC,nRotC+1)
ulalume3@0 1036 for iRetC in range(-nRetC,nRetC+1)
ulalume3@0 1037 for iDiC in range(-nDiC,nDiC+1)
ulalume3@0 1038 for iTP in range(-nTP,nTP+1)
ulalume3@0 1039 for iTS in range(-nTS,nTS+1)
ulalume3@0 1040 for iRP in range(-nRP,nRP+1)
ulalume3@0 1041 for iRS in range(-nRS,nRS+1)
ulalume3@0 1042 for iERaT in range(-nERaT,nERaT+1)
ulalume3@0 1043 for iRotaT in range(-nRotaT,nRotaT+1)
ulalume3@0 1044 for iRetT in range(-nRetT,nRetT+1)
ulalume3@0 1045 for iERaR in range(-nERaR,nERaR+1)
ulalume3@0 1046 for iRotaR in range(-nRotaR,nRotaR+1)
ulalume3@0 1047 for iRetR in range(-nRetR,nRetR+1)]:
ulalume3@0 1048
ulalume3@0 1049 iN = iN + 1
ulalume3@0 1050 if (iN == 10001):
ulalume3@0 1051 ctime = clock()
ulalume3@0 1052 print(" estimated time ", "{0:4.2f}".format(N/10000 * (ctime-atime)), "sec ") #, end="")
ulalume3@0 1053 print("\r elapsed time ", "{0:5.0f}".format((ctime-atime)), "sec ", end="\r")
ulalume3@0 1054 ctime = clock()
ulalume3@0 1055 if ((ctime - dtime) > 10):
ulalume3@0 1056 print("\r elapsed time ", "{0:5.0f}".format((ctime-atime)), "sec ", end="\r")
ulalume3@0 1057 dtime = ctime
ulalume3@0 1058
ulalume3@0 1059 if nRotO > 0: RotO = RotO0 + iRotO*dRotO/nRotO
ulalume3@0 1060 if nRetO > 0: RetO = RetO0 + iRetO*dRetO/nRetO
ulalume3@0 1061 if nDiO > 0: DiO = DiO0 + iDiO*dDiO/nDiO
ulalume3@0 1062 if nRotC > 0: RotC = RotC0 + iRotC*dRotC/nRotC
ulalume3@0 1063 if nRetC > 0: RetC = RetC0 + iRetC*dRetC/nRetC
ulalume3@0 1064 if nDiC > 0: DiC = DiC0 + iDiC*dDiC/nDiC
ulalume3@0 1065 if nTP > 0: TP = TP0 + iTP*dTP/nTP
ulalume3@0 1066 if nTS > 0: TS = TS0 + iTS*dTS/nTS
ulalume3@0 1067 if nRP > 0: RP = RP0 + iRP*dRP/nRP
ulalume3@0 1068 if nRS > 0: RS = RS0 + iRS*dRS/nRS
ulalume3@0 1069 if nERaT > 0: ERaT = ERaT0 + iERaT*dERaT/nERaT
ulalume3@0 1070 if nRotaT > 0:RotaT= RotaT0+ iRotaT*dRotaT/nRotaT
ulalume3@0 1071 if nRetT > 0: RetT = RetT0 + iRetT*dRetT/nRetT
ulalume3@0 1072 if nERaR > 0: ERaR = ERaR0 + iERaR*dERaR/nERaR
ulalume3@0 1073 if nRotaR > 0:RotaR= RotaR0+ iRotaR*dRotaR/nRotaR
ulalume3@0 1074 if nRetR > 0: RetR = RetR0 + iRetR*dRetR/nRetR
ulalume3@0 1075
ulalume3@0 1076 #print("{0:5.2f}, {1:5.2f}, {2:5.2f}, {3:10d}".format(RotL, RotE, RotO, iN))
ulalume3@0 1077
ulalume3@0 1078 # receiver optics
ulalume3@0 1079 CosO = np.cos(np.deg2rad(RetO))
ulalume3@0 1080 SinO = np.sin(np.deg2rad(RetO))
ulalume3@0 1081 ZiO = (1. - DiO**2)**0.5
ulalume3@0 1082 WiO = (1. - ZiO*CosO)
ulalume3@0 1083 S2g = np.sin(np.deg2rad(2*RotO))
ulalume3@0 1084 C2g = np.cos(np.deg2rad(2*RotO))
ulalume3@0 1085 # calibrator
ulalume3@0 1086 CosC = np.cos(np.deg2rad(RetC))
ulalume3@0 1087 SinC = np.sin(np.deg2rad(RetC))
ulalume3@0 1088 ZiC = (1. - DiC**2)**0.5
ulalume3@0 1089 WiC = (1. - ZiC*CosC)
ulalume3@0 1090
ulalume3@0 1091 #For POLLY_XT
ulalume3@0 1092 # analyser
ulalume3@0 1093 RS = 1 - TS
ulalume3@0 1094 RP = 1 - TP
ulalume3@0 1095 TiT = 0.5 * (TP + TS)
ulalume3@0 1096 DiT = (TP-TS)/(TP+TS)
ulalume3@0 1097 ZiT = (1. - DiT**2)**0.5
ulalume3@0 1098 TiR = 0.5 * (RP + RS)
ulalume3@0 1099 DiR = (RP-RS)/(RP+RS)
ulalume3@0 1100 ZiR = (1. - DiR**2)**0.5
ulalume3@0 1101 CosT = np.cos(np.deg2rad(RetT))
ulalume3@0 1102 SinT = np.sin(np.deg2rad(RetT))
ulalume3@0 1103 CosR = np.cos(np.deg2rad(RetR))
ulalume3@0 1104 SinR = np.sin(np.deg2rad(RetR))
ulalume3@0 1105
ulalume3@0 1106 DaT = (1-ERaT)/(1+ERaT)
ulalume3@0 1107 DaR = (1-ERaR)/(1+ERaR)
ulalume3@0 1108 TaT = 0.5*(1+ERaT)
ulalume3@0 1109 TaR = 0.5*(1+ERaR)
ulalume3@0 1110
ulalume3@0 1111 S2aT = np.sin(np.deg2rad(h*2*RotaT))
ulalume3@0 1112 C2aT = np.cos(np.deg2rad(2*RotaT))
ulalume3@0 1113 S2aR = np.sin(np.deg2rad(h*2*RotaR))
ulalume3@0 1114 C2aR = np.cos(np.deg2rad(2*RotaR))
ulalume3@0 1115
ulalume3@0 1116 # Aanalyzer As before the PBS Eq. D.5
ulalume3@0 1117 ATP1 = (1+C2aT*DaT*DiT)
ulalume3@0 1118 ATP2 = Y*(DiT+C2aT*DaT)
ulalume3@0 1119 ATP3 = Y*S2aT*DaT*ZiT*CosT
ulalume3@0 1120 ATP4 = S2aT*DaT*ZiT*SinT
ulalume3@0 1121 ATP = np.array([ATP1,ATP2,ATP3,ATP4])
ulalume3@0 1122
ulalume3@0 1123 ARP1 = (1+C2aR*DaR*DiR)
ulalume3@0 1124 ARP2 = Y*(DiR+C2aR*DaR)
ulalume3@0 1125 ARP3 = Y*S2aR*DaR*ZiR*CosR
ulalume3@0 1126 ARP4 = S2aR*DaR*ZiR*SinR
ulalume3@0 1127 ARP = np.array([ARP1,ARP2,ARP3,ARP4])
ulalume3@0 1128
ulalume3@0 1129 TTa = TiT*TaT #*ATP1
ulalume3@0 1130 TRa = TiR*TaR #*ARP1
ulalume3@0 1131
ulalume3@0 1132 # ---- Calculate signals and correction parameters for diffeent locations and calibrators
ulalume3@0 1133 if LocC == 4: # Calibrator before the PBS
ulalume3@0 1134 #print("Calibrator location not implemented yet")
ulalume3@0 1135
ulalume3@0 1136 #S2ge = np.sin(np.deg2rad(2*RotO + h*2*RotC))
ulalume3@0 1137 #C2ge = np.cos(np.deg2rad(2*RotO + h*2*RotC))
ulalume3@0 1138 S2e = np.sin(np.deg2rad(h*2*RotC))
ulalume3@0 1139 C2e = np.cos(np.deg2rad(2*RotC))
ulalume3@0 1140 # rotated AinP by epsilon Eq. C.3
ulalume3@0 1141 ATP2e = C2e*ATP2 + S2e*ATP3
ulalume3@0 1142 ATP3e = C2e*ATP3 - S2e*ATP2
ulalume3@0 1143 ARP2e = C2e*ARP2 + S2e*ARP3
ulalume3@0 1144 ARP3e = C2e*ARP3 - S2e*ARP2
ulalume3@0 1145 ATPe = np.array([ATP1,ATP2e,ATP3e,ATP4])
ulalume3@0 1146 ARPe = np.array([ARP1,ARP2e,ARP3e,ARP4])
ulalume3@0 1147 # Stokes Input Vector before the polarising beam splitter Eq. E.31
ulalume3@0 1148 A = C2g*QinE - S2g*UinE
ulalume3@0 1149 B = S2g*QinE + C2g*UinE
ulalume3@0 1150 #C = (WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE)
ulalume3@0 1151 Co = ZiO*SinO*VinE
ulalume3@0 1152 Ca = (WiO*B - 2*ZiO*SinO*VinE)
ulalume3@0 1153 #C = Co + aCal*Ca
ulalume3@0 1154 #IinP = (IinE + DiO*aCal*A)
ulalume3@0 1155 #QinP = (C2g*DiO*IinE + aCal*QinE - S2g*C)
ulalume3@0 1156 #UinP = (S2g*DiO*IinE - aCal*UinE + C2g*C)
ulalume3@0 1157 #VinP = (ZiO*SinO*aCal*B + ZiO*CosO*(1-2*aCal)*VinE)
ulalume3@0 1158 IinPo = IinE
ulalume3@0 1159 QinPo = (C2g*DiO*IinE - S2g*Co)
ulalume3@0 1160 UinPo = (S2g*DiO*IinE + C2g*Co)
ulalume3@0 1161 VinPo = ZiO*CosO*VinE
ulalume3@0 1162
ulalume3@0 1163 IinPa = DiO*A
ulalume3@0 1164 QinPa = QinE - S2g*Ca
ulalume3@0 1165 UinPa = -UinE + C2g*Ca
ulalume3@0 1166 VinPa = ZiO*(SinO*B - 2*CosO*VinE)
ulalume3@0 1167
ulalume3@0 1168 IinP = IinPo + aCal*IinPa
ulalume3@0 1169 QinP = QinPo + aCal*QinPa
ulalume3@0 1170 UinP = UinPo + aCal*UinPa
ulalume3@0 1171 VinP = VinPo + aCal*VinPa
ulalume3@0 1172 # Stokes Input Vector before the polarising beam splitter rotated by epsilon Eq. C.3
ulalume3@0 1173 #QinPe = C2e*QinP + S2e*UinP
ulalume3@0 1174 #UinPe = C2e*UinP - S2e*QinP
ulalume3@0 1175 QinPoe = C2e*QinPo + S2e*UinPo
ulalume3@0 1176 UinPoe = C2e*UinPo - S2e*QinPo
ulalume3@0 1177 QinPae = C2e*QinPa + S2e*UinPa
ulalume3@0 1178 UinPae = C2e*UinPa - S2e*QinPa
ulalume3@0 1179 QinPe = C2e*QinP + S2e*UinP
ulalume3@0 1180 UinPe = C2e*UinP - S2e*QinP
ulalume3@0 1181
ulalume3@0 1182 # Calibration signals and Calibration correction K from measurements with LDRCal / aCal
ulalume3@0 1183 if (TypeC == 2) or (TypeC == 1): # rotator calibration Eq. C.4
ulalume3@0 1184 # parameters for calibration with aCal
ulalume3@0 1185 AT = ATP1*IinP + h*ATP4*VinP
ulalume3@0 1186 BT = ATP3e*QinP - h*ATP2e*UinP
ulalume3@0 1187 AR = ARP1*IinP + h*ARP4*VinP
ulalume3@0 1188 BR = ARP3e*QinP - h*ARP2e*UinP
ulalume3@0 1189 # Correction paremeters for normal measurements; they are independent of LDR
ulalume3@0 1190 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out
ulalume3@0 1191 IS1 = np.array([IinPo,QinPo,UinPo,VinPo])
ulalume3@0 1192 IS2 = np.array([IinPa,QinPa,UinPa,VinPa])
ulalume3@0 1193 GT = np.dot(ATP,IS1)
ulalume3@0 1194 GR = np.dot(ARP,IS1)
ulalume3@0 1195 HT = np.dot(ATP,IS2)
ulalume3@0 1196 HR = np.dot(ARP,IS2)
ulalume3@0 1197 else:
ulalume3@0 1198 IS1 = np.array([IinPo,QinPo,UinPo,VinPo])
ulalume3@0 1199 IS2 = np.array([IinPa,QinPa,UinPa,VinPa])
ulalume3@0 1200 GT = np.dot(ATPe,IS1)
ulalume3@0 1201 GR = np.dot(ARPe,IS1)
ulalume3@0 1202 HT = np.dot(ATPe,IS2)
ulalume3@0 1203 HR = np.dot(ARPe,IS2)
ulalume3@0 1204 elif (TypeC == 3) or (TypeC == 4): # linear polariser calibration Eq. C.5
ulalume3@0 1205 # parameters for calibration with aCal
ulalume3@0 1206 AT = ATP1*IinP + ATP3e*UinPe + ZiC*CosC*(ATP2e*QinPe + ATP4*VinP)
ulalume3@0 1207 BT = DiC*(ATP1*UinPe + ATP3e*IinP) - ZiC*SinC*(ATP2e*VinP - ATP4*QinPe)
ulalume3@0 1208 AR = ARP1*IinP + ARP3e*UinPe + ZiC*CosC*(ARP2e*QinPe + ARP4*VinP)
ulalume3@0 1209 BR = DiC*(ARP1*UinPe + ARP3e*IinP) - ZiC*SinC*(ARP2e*VinP - ARP4*QinPe)
ulalume3@0 1210 # Correction paremeters for normal measurements; they are independent of LDR
ulalume3@0 1211 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out
ulalume3@0 1212 IS1 = np.array([IinPo,QinPo,UinPo,VinPo])
ulalume3@0 1213 IS2 = np.array([IinPa,QinPa,UinPa,VinPa])
ulalume3@0 1214 GT = np.dot(ATP,IS1)
ulalume3@0 1215 GR = np.dot(ARP,IS1)
ulalume3@0 1216 HT = np.dot(ATP,IS2)
ulalume3@0 1217 HR = np.dot(ARP,IS2)
ulalume3@0 1218 else:
ulalume3@0 1219 IS1e = np.array([IinPo+DiC*QinPoe,DiC*IinPo+QinPoe,ZiC*(CosC*UinPoe+SinC*VinPo),-ZiC*(SinC*UinPoe-CosC*VinPo)])
ulalume3@0 1220 IS2e = np.array([IinPa+DiC*QinPae,DiC*IinPa+QinPae,ZiC*(CosC*UinPae+SinC*VinPa),-ZiC*(SinC*UinPae-CosC*VinPa)])
ulalume3@0 1221 GT = np.dot(ATPe,IS1e)
ulalume3@0 1222 GR = np.dot(ARPe,IS1e)
ulalume3@0 1223 HT = np.dot(ATPe,IS2e)
ulalume3@0 1224 HR = np.dot(ARPe,IS2e)
ulalume3@0 1225 elif (TypeC == 6): # diattenuator calibration +-22.5° rotated_diattenuator_X22x5deg.odt
ulalume3@0 1226 # parameters for calibration with aCal
ulalume3@0 1227 AT = ATP1*IinP + sqr05*DiC*(ATP1*QinPe + ATP2e*IinP) + (1-0.5*WiC)*(ATP2e*QinPe + ATP3e*UinPe) + ZiC*(sqr05*SinC*(ATP3e*VinP-ATP4*UinPe) + ATP4*CosC*VinP)
ulalume3@0 1228 BT = sqr05*DiC*(ATP1*UinPe + ATP3e*IinP) + 0.5*WiC*(ATP2e*UinPe + ATP3e*QinPe) - sqr05*ZiC*SinC*(ATP2e*VinP - ATP4*QinPe)
ulalume3@0 1229 AR = ARP1*IinP + sqr05*DiC*(ARP1*QinPe + ARP2e*IinP) + (1-0.5*WiC)*(ARP2e*QinPe + ARP3e*UinPe) + ZiC*(sqr05*SinC*(ARP3e*VinP-ARP4*UinPe) + ARP4*CosC*VinP)
ulalume3@0 1230 BR = sqr05*DiC*(ARP1*UinPe + ARP3e*IinP) + 0.5*WiC*(ARP2e*UinPe + ARP3e*QinPe) - sqr05*ZiC*SinC*(ARP2e*VinP - ARP4*QinPe)
ulalume3@0 1231 # Correction paremeters for normal measurements; they are independent of LDR
ulalume3@0 1232 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out
ulalume3@0 1233 IS1 = np.array([IinPo,QinPo,UinPo,VinPo])
ulalume3@0 1234 IS2 = np.array([IinPa,QinPa,UinPa,VinPa])
ulalume3@0 1235 GT = np.dot(ATP,IS1)
ulalume3@0 1236 GR = np.dot(ARP,IS1)
ulalume3@0 1237 HT = np.dot(ATP,IS2)
ulalume3@0 1238 HR = np.dot(ARP,IS2)
ulalume3@0 1239 else:
ulalume3@0 1240 IS1e = np.array([IinPo+DiC*QinPoe,DiC*IinPo+QinPoe,ZiC*(CosC*UinPoe+SinC*VinPo),-ZiC*(SinC*UinPoe-CosC*VinPo)])
ulalume3@0 1241 IS2e = np.array([IinPa+DiC*QinPae,DiC*IinPa+QinPae,ZiC*(CosC*UinPae+SinC*VinPa),-ZiC*(SinC*UinPae-CosC*VinPa)])
ulalume3@0 1242 GT = np.dot(ATPe,IS1e)
ulalume3@0 1243 GR = np.dot(ARPe,IS1e)
ulalume3@0 1244 HT = np.dot(ATPe,IS2e)
ulalume3@0 1245 HR = np.dot(ARPe,IS2e)
ulalume3@0 1246 else:
ulalume3@0 1247 print("Calibrator not implemented yet")
ulalume3@0 1248 sys.exit()
ulalume3@0 1249
ulalume3@0 1250 elif LocC == 3: # C before receiver optics Eq.57
ulalume3@0 1251
ulalume3@0 1252 #S2ge = np.sin(np.deg2rad(2*RotO - 2*RotC))
ulalume3@0 1253 #C2ge = np.cos(np.deg2rad(2*RotO - 2*RotC))
ulalume3@0 1254 S2e = np.sin(np.deg2rad(2*RotC))
ulalume3@0 1255 C2e = np.cos(np.deg2rad(2*RotC))
ulalume3@0 1256
ulalume3@0 1257 # AS with C before the receiver optics (see document rotated_diattenuator_X22x5deg.odt)
ulalume3@0 1258 AF1 = np.array([1,C2g*DiO,S2g*DiO,0])
ulalume3@0 1259 AF2 = np.array([C2g*DiO,1-S2g**2*WiO,S2g*C2g*WiO,-S2g*ZiO*SinO])
ulalume3@0 1260 AF3 = np.array([S2g*DiO, S2g*C2g*WiO, 1-C2g**2*WiO, C2g*ZiO*SinO])
ulalume3@0 1261 AF4 = np.array([0, S2g*SinO, -C2g*SinO, CosO])
ulalume3@0 1262
ulalume3@0 1263 ATF = (ATP1*AF1+ATP2*AF2+ATP3*AF3+ATP4*AF4)
ulalume3@0 1264 ARF = (ARP1*AF1+ARP2*AF2+ARP3*AF3+ARP4*AF4)
ulalume3@0 1265 ATF1 = ATF[0]
ulalume3@0 1266 ATF2 = ATF[1]
ulalume3@0 1267 ATF3 = ATF[2]
ulalume3@0 1268 ATF4 = ATF[3]
ulalume3@0 1269 ARF1 = ARF[0]
ulalume3@0 1270 ARF2 = ARF[1]
ulalume3@0 1271 ARF3 = ARF[2]
ulalume3@0 1272 ARF4 = ARF[3]
ulalume3@0 1273
ulalume3@0 1274 # rotated AinF by epsilon
ulalume3@0 1275 ATF2e = C2e*ATF[1] + S2e*ATF[2]
ulalume3@0 1276 ATF3e = C2e*ATF[2] - S2e*ATF[1]
ulalume3@0 1277 ARF2e = C2e*ARF[1] + S2e*ARF[2]
ulalume3@0 1278 ARF3e = C2e*ARF[2] - S2e*ARF[1]
ulalume3@0 1279
ulalume3@0 1280 ATFe = np.array([ATF1,ATF2e,ATF3e,ATF4])
ulalume3@0 1281 ARFe = np.array([ARF1,ARF2e,ARF3e,ARF4])
ulalume3@0 1282
ulalume3@0 1283 QinEe = C2e*QinE + S2e*UinE
ulalume3@0 1284 UinEe = C2e*UinE - S2e*QinE
ulalume3@0 1285
ulalume3@0 1286 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F)
ulalume3@0 1287 IinF = IinE
ulalume3@0 1288 QinF = aCal*QinE
ulalume3@0 1289 UinF = -aCal*UinE
ulalume3@0 1290 VinF = (1.-2.*aCal)*VinE
ulalume3@0 1291
ulalume3@0 1292 IinFo = IinE
ulalume3@0 1293 QinFo = 0.
ulalume3@0 1294 UinFo = 0.
ulalume3@0 1295 VinFo = VinE
ulalume3@0 1296
ulalume3@0 1297 IinFa = 0.
ulalume3@0 1298 QinFa = QinE
ulalume3@0 1299 UinFa = -UinE
ulalume3@0 1300 VinFa = -2.*VinE
ulalume3@0 1301
ulalume3@0 1302 # Stokes Input Vector before receiver optics rotated by epsilon Eq. C.3
ulalume3@0 1303 QinFe = C2e*QinF + S2e*UinF
ulalume3@0 1304 UinFe = C2e*UinF - S2e*QinF
ulalume3@0 1305 QinFoe = C2e*QinFo + S2e*UinFo
ulalume3@0 1306 UinFoe = C2e*UinFo - S2e*QinFo
ulalume3@0 1307 QinFae = C2e*QinFa + S2e*UinFa
ulalume3@0 1308 UinFae = C2e*UinFa - S2e*QinFa
ulalume3@0 1309
ulalume3@0 1310 # Calibration signals and Calibration correction K from measurements with LDRCal / aCal
ulalume3@0 1311 if (TypeC == 2) or (TypeC == 1): # rotator calibration Eq. C.4
ulalume3@0 1312 AT = ATF1*IinF + ATF4*h*VinF
ulalume3@0 1313 BT = ATF3e*QinF - ATF2e*h*UinF
ulalume3@0 1314 AR = ARF1*IinF + ARF4*h*VinF
ulalume3@0 1315 BR = ARF3e*QinF - ARF2e*h*UinF
ulalume3@0 1316
ulalume3@0 1317 # Correction paremeters for normal measurements; they are independent of LDR
ulalume3@0 1318 if (not RotationErrorEpsilonForNormalMeasurements):
ulalume3@0 1319 GT = ATF1*IinE + ATF4*VinE
ulalume3@0 1320 GR = ARF1*IinE + ARF4*VinE
ulalume3@0 1321 HT = ATF2*QinE - ATF3*UinE - ATF4*2*VinE
ulalume3@0 1322 HR = ARF2*QinE - ARF3*UinE - ARF4*2*VinE
ulalume3@0 1323 else:
ulalume3@0 1324 GT = ATF1*IinE + ATF4*h*VinE
ulalume3@0 1325 GR = ARF1*IinE + ARF4*h*VinE
ulalume3@0 1326 HT = ATF2e*QinE - ATF3e*h*UinE - ATF4*h*2*VinE
ulalume3@0 1327 HR = ARF2e*QinE - ARF3e*h*UinE - ARF4*h*2*VinE
ulalume3@0 1328
ulalume3@0 1329 elif (TypeC == 3) or (TypeC == 4): # linear polariser calibration Eq. C.5
ulalume3@0 1330 # p = +45°, m = -45°
ulalume3@0 1331 IF1e = np.array([IinF, ZiC*CosC*QinFe, UinFe, ZiC*CosC*VinF])
ulalume3@0 1332 IF2e = np.array([DiC*UinFe, -ZiC*SinC*VinF, DiC*IinF, ZiC*SinC*QinFe])
ulalume3@0 1333
ulalume3@0 1334 AT = np.dot(ATFe,IF1e)
ulalume3@0 1335 AR = np.dot(ARFe,IF1e)
ulalume3@0 1336 BT = np.dot(ATFe,IF2e)
ulalume3@0 1337 BR = np.dot(ARFe,IF2e)
ulalume3@0 1338
ulalume3@0 1339 # Correction paremeters for normal measurements; they are independent of LDR --- the same as for TypeC = 6
ulalume3@0 1340 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out
ulalume3@0 1341 IS1 = np.array([IinE,0,0,VinE])
ulalume3@0 1342 IS2 = np.array([0,QinE,-UinE,-2*VinE])
ulalume3@0 1343
ulalume3@0 1344 GT = np.dot(ATF,IS1)
ulalume3@0 1345 GR = np.dot(ARF,IS1)
ulalume3@0 1346 HT = np.dot(ATF,IS2)
ulalume3@0 1347 HR = np.dot(ARF,IS2)
ulalume3@0 1348 else:
ulalume3@0 1349 IS1e = np.array([IinFo+DiC*QinFoe,DiC*IinFo+QinFoe,ZiC*(CosC*UinFoe+SinC*VinFo),-ZiC*(SinC*UinFoe-CosC*VinFo)])
ulalume3@0 1350 IS2e = np.array([IinFa+DiC*QinFae,DiC*IinFa+QinFae,ZiC*(CosC*UinFae+SinC*VinFa),-ZiC*(SinC*UinFae-CosC*VinFa)])
ulalume3@0 1351 GT = np.dot(ATFe,IS1e)
ulalume3@0 1352 GR = np.dot(ARFe,IS1e)
ulalume3@0 1353 HT = np.dot(ATFe,IS2e)
ulalume3@0 1354 HR = np.dot(ARFe,IS2e)
ulalume3@0 1355
ulalume3@0 1356 elif (TypeC == 6): # diattenuator calibration +-22.5° rotated_diattenuator_X22x5deg.odt
ulalume3@0 1357 # p = +22.5°, m = -22.5°
ulalume3@0 1358 IF1e = np.array([IinF+sqr05*DiC*QinFe, sqr05*DiC*IinF+(1-0.5*WiC)*QinFe, (1-0.5*WiC)*UinFe+sqr05*ZiC*SinC*VinF, -sqr05*ZiC*SinC*UinFe+ZiC*CosC*VinF])
ulalume3@0 1359 IF2e = np.array([sqr05*DiC*UinFe, 0.5*WiC*UinFe-sqr05*ZiC*SinC*VinF, sqr05*DiC*IinF+0.5*WiC*QinFe, sqr05*ZiC*SinC*QinFe])
ulalume3@0 1360
ulalume3@0 1361 AT = np.dot(ATFe,IF1e)
ulalume3@0 1362 AR = np.dot(ARFe,IF1e)
ulalume3@0 1363 BT = np.dot(ATFe,IF2e)
ulalume3@0 1364 BR = np.dot(ARFe,IF2e)
ulalume3@0 1365
ulalume3@0 1366 # Correction paremeters for normal measurements; they are independent of LDR
ulalume3@0 1367 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out
ulalume3@0 1368 #IS1 = np.array([IinE,0,0,VinE])
ulalume3@0 1369 #IS2 = np.array([0,QinE,-UinE,-2*VinE])
ulalume3@0 1370 IS1 = np.array([IinFo,0,0,VinFo])
ulalume3@0 1371 IS2 = np.array([0,QinFa,UinFa,VinFa])
ulalume3@0 1372 GT = np.dot(ATF,IS1)
ulalume3@0 1373 GR = np.dot(ARF,IS1)
ulalume3@0 1374 HT = np.dot(ATF,IS2)
ulalume3@0 1375 HR = np.dot(ARF,IS2)
ulalume3@0 1376 else:
ulalume3@0 1377 #IS1e = np.array([IinE,DiC*IinE,ZiC*SinC*VinE,ZiC*CosC*VinE])
ulalume3@0 1378 #IS2e = np.array([DiC*QinEe,QinEe,-ZiC*(CosC*UinEe+2*SinC*VinE),ZiC*(SinC*UinEe-2*CosC*VinE)])
ulalume3@0 1379 IS1e = np.array([IinFo+DiC*QinFoe,DiC*IinFo+QinFoe,ZiC*(CosC*UinFoe+SinC*VinFo),-ZiC*(SinC*UinFoe-CosC*VinFo)])
ulalume3@0 1380 IS2e = np.array([IinFa+DiC*QinFae,DiC*IinFa+QinFae,ZiC*(CosC*UinFae+SinC*VinFa),-ZiC*(SinC*UinFae-CosC*VinFa)])
ulalume3@0 1381 GT = np.dot(ATFe,IS1e)
ulalume3@0 1382 GR = np.dot(ARFe,IS1e)
ulalume3@0 1383 HT = np.dot(ATFe,IS2e)
ulalume3@0 1384 HR = np.dot(ARFe,IS2e)
ulalume3@0 1385
ulalume3@0 1386
ulalume3@0 1387 else:
ulalume3@0 1388 print('Calibrator not implemented yet')
ulalume3@0 1389 sys.exit()
ulalume3@0 1390
ulalume3@0 1391 elif LocC == 2: # C behind emitter optics Eq.57
ulalume3@0 1392 #print("Calibrator location not implemented yet")
ulalume3@0 1393 S2e = np.sin(np.deg2rad(2*RotC))
ulalume3@0 1394 C2e = np.cos(np.deg2rad(2*RotC))
ulalume3@0 1395
ulalume3@0 1396 # AS with C before the receiver optics (see document rotated_diattenuator_X22x5deg.odt)
ulalume3@0 1397 AF1 = np.array([1,C2g*DiO,S2g*DiO,0])
ulalume3@0 1398 AF2 = np.array([C2g*DiO,1-S2g**2*WiO,S2g*C2g*WiO,-S2g*ZiO*SinO])
ulalume3@0 1399 AF3 = np.array([S2g*DiO, S2g*C2g*WiO, 1-C2g**2*WiO, C2g*ZiO*SinO])
ulalume3@0 1400 AF4 = np.array([0, S2g*SinO, -C2g*SinO, CosO])
ulalume3@0 1401
ulalume3@0 1402 ATF = (ATP1*AF1+ATP2*AF2+ATP3*AF3+ATP4*AF4)
ulalume3@0 1403 ARF = (ARP1*AF1+ARP2*AF2+ARP3*AF3+ARP4*AF4)
ulalume3@0 1404 ATF1 = ATF[0]
ulalume3@0 1405 ATF2 = ATF[1]
ulalume3@0 1406 ATF3 = ATF[2]
ulalume3@0 1407 ATF4 = ATF[3]
ulalume3@0 1408 ARF1 = ARF[0]
ulalume3@0 1409 ARF2 = ARF[1]
ulalume3@0 1410 ARF3 = ARF[2]
ulalume3@0 1411 ARF4 = ARF[3]
ulalume3@0 1412
ulalume3@0 1413 # AS with C behind the emitter --------------------------------------------
ulalume3@0 1414 # terms without aCal
ulalume3@0 1415 ATE1o, ARE1o = ATF1, ARF1
ulalume3@0 1416 ATE2o, ARE2o = 0., 0.
ulalume3@0 1417 ATE3o, ARE3o = 0., 0.
ulalume3@0 1418 ATE4o, ARE4o = ATF4, ARF4
ulalume3@0 1419 # terms with aCal
ulalume3@0 1420 ATE1a, ARE1a = 0. , 0.
ulalume3@0 1421 ATE2a, ARE2a = ATF2, ARF2
ulalume3@0 1422 ATE3a, ARE3a = -ATF3, -ARF3
ulalume3@0 1423 ATE4a, ARE4a = -2*ATF4, -2*ARF4
ulalume3@0 1424 # rotated AinEa by epsilon
ulalume3@0 1425 ATE2ae = C2e*ATF2 + S2e*ATF3
ulalume3@0 1426 ATE3ae = -S2e*ATF2 - C2e*ATF3
ulalume3@0 1427 ARE2ae = C2e*ARF2 + S2e*ARF3
ulalume3@0 1428 ARE3ae = -S2e*ARF2 - C2e*ARF3
ulalume3@0 1429
ulalume3@0 1430 ATE1 = ATE1o
ulalume3@0 1431 ATE2e = aCal*ATE2ae
ulalume3@0 1432 ATE3e = aCal*ATE3ae
ulalume3@0 1433 ATE4 = (1-2*aCal)*ATF4
ulalume3@0 1434 ARE1 = ARE1o
ulalume3@0 1435 ARE2e = aCal*ARE2ae
ulalume3@0 1436 ARE3e = aCal*ARE3ae
ulalume3@0 1437 ARE4 = (1-2*aCal)*ARF4
ulalume3@0 1438
ulalume3@0 1439 # rotated IinE
ulalume3@0 1440 QinEe = C2e*QinE + S2e*UinE
ulalume3@0 1441 UinEe = C2e*UinE - S2e*QinE
ulalume3@0 1442
ulalume3@0 1443 # --- Calibration signals and Calibration correction K from measurements with LDRCal / aCal
ulalume3@0 1444 if (TypeC == 2) or (TypeC == 1): # +++++++++ rotator calibration Eq. C.4
ulalume3@0 1445 AT = ATE1o*IinE + (ATE4o+aCal*ATE4a)*h*VinE
ulalume3@0 1446 BT = aCal * (ATE3ae*QinEe - ATE2ae*h*UinEe)
ulalume3@0 1447 AR = ARE1o*IinE + (ARE4o+aCal*ARE4a)*h*VinE
ulalume3@0 1448 BR = aCal * (ARE3ae*QinEe - ARE2ae*h*UinEe)
ulalume3@0 1449
ulalume3@0 1450 # Correction paremeters for normal measurements; they are independent of LDR
ulalume3@0 1451 if (not RotationErrorEpsilonForNormalMeasurements):
ulalume3@0 1452 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F)
ulalume3@0 1453 GT = ATE1o*IinE + ATE4o*h*VinE
ulalume3@0 1454 GR = ARE1o*IinE + ARE4o*h*VinE
ulalume3@0 1455 HT = ATE2a*QinE + ATE3a*h*UinEe + ATE4a*h*VinE
ulalume3@0 1456 HR = ARE2a*QinE + ARE3a*h*UinEe + ARE4a*h*VinE
ulalume3@0 1457 else:
ulalume3@0 1458 GT = ATE1o*IinE + ATE4o*h*VinE
ulalume3@0 1459 GR = ARE1o*IinE + ARE4o*h*VinE
ulalume3@0 1460 HT = ATE2ae*QinE + ATE3ae*h*UinEe + ATE4a*h*VinE
ulalume3@0 1461 HR = ARE2ae*QinE + ARE3ae*h*UinEe + ARE4a*h*VinE
ulalume3@0 1462
ulalume3@0 1463 elif (TypeC == 3) or (TypeC == 4): # +++++++++ linear polariser calibration Eq. C.5
ulalume3@0 1464 # p = +45°, m = -45°
ulalume3@0 1465 AT = ATE1*IinE + ZiC*CosC*(ATE2e*QinEe + ATE4*VinE) + ATE3e*UinEe
ulalume3@0 1466 BT = DiC*(ATE1*UinEe + ATE3e*IinE) + ZiC*SinC*(ATE4*QinEe - ATE2e*VinE)
ulalume3@0 1467 AR = ARE1*IinE + ZiC*CosC*(ARE2e*QinEe + ARE4*VinE) + ARE3e*UinEe
ulalume3@0 1468 BR = DiC*(ARE1*UinEe + ARE3e*IinE) + ZiC*SinC*(ARE4*QinEe - ARE2e*VinE)
ulalume3@0 1469
ulalume3@0 1470 # Correction paremeters for normal measurements; they are independent of LDR
ulalume3@0 1471 if (not RotationErrorEpsilonForNormalMeasurements):
ulalume3@0 1472 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F)
ulalume3@0 1473 GT = ATE1o*IinE + ATE4o*VinE
ulalume3@0 1474 GR = ARE1o*IinE + ARE4o*VinE
ulalume3@0 1475 HT = ATE2a*QinE + ATE3a*UinE + ATE4a*VinE
ulalume3@0 1476 HR = ARE2a*QinE + ARE3a*UinE + ARE4a*VinE
ulalume3@0 1477 else:
ulalume3@0 1478 D = IinE + DiC*QinEe
ulalume3@0 1479 A = DiC*IinE + QinEe
ulalume3@0 1480 B = ZiC*(CosC*UinEe + SinC*VinE)
ulalume3@0 1481 C = -ZiC*(SinC*UinEe - CosC*VinE)
ulalume3@0 1482 GT = ATE1o*D + ATE4o*C
ulalume3@0 1483 GR = ARE1o*D + ARE4o*C
ulalume3@0 1484 HT = ATE2a*A + ATE3a*B + ATE4a*C
ulalume3@0 1485 HR = ARE2a*A + ARE3a*B + ARE4a*C
ulalume3@0 1486
ulalume3@0 1487 elif (TypeC == 6): # real HWP calibration +-22.5° rotated_diattenuator_X22x5deg.odt
ulalume3@0 1488 # p = +22.5°, m = -22.5°
ulalume3@0 1489 IE1e = np.array([IinE+sqr05*DiC*QinEe, sqr05*DiC*IinE+(1-0.5*WiC)*QinEe, (1-0.5*WiC)*UinEe+sqr05*ZiC*SinC*VinE, -sqr05*ZiC*SinC*UinEe+ZiC*CosC*VinE])
ulalume3@0 1490 IE2e = np.array([sqr05*DiC*UinEe, 0.5*WiC*UinEe-sqr05*ZiC*SinC*VinE, sqr05*DiC*IinE+0.5*WiC*QinEe, sqr05*ZiC*SinC*QinEe])
ulalume3@0 1491 ATEe = np.array([ATE1,ATE2e,ATE3e,ATE4])
ulalume3@0 1492 AREe = np.array([ARE1,ARE2e,ARE3e,ARE4])
ulalume3@0 1493 AT = np.dot(ATEe,IE1e)
ulalume3@0 1494 AR = np.dot(AREe,IE1e)
ulalume3@0 1495 BT = np.dot(ATEe,IE2e)
ulalume3@0 1496 BR = np.dot(AREe,IE2e)
ulalume3@0 1497
ulalume3@0 1498 # Correction paremeters for normal measurements; they are independent of LDR
ulalume3@0 1499 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out
ulalume3@0 1500 GT = ATE1o*IinE + ATE4o*VinE
ulalume3@0 1501 GR = ARE1o*IinE + ARE4o*VinE
ulalume3@0 1502 HT = ATE2a*QinE + ATE3a*UinE + ATE4a*VinE
ulalume3@0 1503 HR = ARE2a*QinE + ARE3a*UinE + ARE4a*VinE
ulalume3@0 1504 else:
ulalume3@0 1505 D = IinE + DiC*QinEe
ulalume3@0 1506 A = DiC*IinE + QinEe
ulalume3@0 1507 B = ZiC*(CosC*UinEe + SinC*VinE)
ulalume3@0 1508 C = -ZiC*(SinC*UinEe - CosC*VinE)
ulalume3@0 1509 GT = ATE1o*D + ATE4o*C
ulalume3@0 1510 GR = ARE1o*D + ARE4o*C
ulalume3@0 1511 HT = ATE2a*A + ATE3a*B + ATE4a*C
ulalume3@0 1512 HR = ARE2a*A + ARE3a*B + ARE4a*C
ulalume3@0 1513
ulalume3@0 1514 else:
ulalume3@0 1515 print('Calibrator not implemented yet')
ulalume3@0 1516 sys.exit()
ulalume3@0 1517
ulalume3@0 1518 # Calibration signals with aCal => Determination of the correction K of the real calibration factor
ulalume3@0 1519 IoutTp = TaT*TiT*TiO*TiE*(AT + BT)
ulalume3@0 1520 IoutTm = TaT*TiT*TiO*TiE*(AT - BT)
ulalume3@0 1521 IoutRp = TaR*TiR*TiO*TiE*(AR + BR)
ulalume3@0 1522 IoutRm = TaR*TiR*TiO*TiE*(AR - BR)
ulalume3@0 1523 # --- Results and Corrections; electronic etaR and etaT are assumed to be 1
ulalume3@0 1524 #Eta = TiR/TiT # Eta = Eta*/K Eq. 84
ulalume3@0 1525 Etapx = IoutRp/IoutTp
ulalume3@0 1526 Etamx = IoutRm/IoutTm
ulalume3@0 1527 Etax = (Etapx*Etamx)**0.5
ulalume3@0 1528 K = Etax / Eta0
ulalume3@0 1529 #print("{0:6.3f},{1:6.3f},{2:6.3f},{3:6.3f},{4:6.3f},{5:6.3f},{6:6.3f},{7:6.3f},{8:6.3f},{9:6.3f},{10:6.3f}".format(AT, BT, AR, BR, DiC, ZiC, RetO, TP, TS, Kp, Km))
ulalume3@0 1530 #print("{0:6.3f},{1:6.3f},{2:6.3f},{3:6.3f}".format(DiC, ZiC, Kp, Km))
ulalume3@0 1531
ulalume3@0 1532 # For comparison with Volkers Libreoffice Müller Matrix spreadsheet
ulalume3@0 1533 #Eta_test_p = (IoutRp/IoutTp)
ulalume3@0 1534 #Eta_test_m = (IoutRm/IoutTm)
ulalume3@0 1535 #Eta_test = (Eta_test_p*Eta_test_m)**0.5
ulalume3@0 1536
ulalume3@0 1537 # *************************************************************************
ulalume3@0 1538 iLDR = -1
ulalume3@0 1539 for LDRTrue in LDRrange:
ulalume3@0 1540 iLDR = iLDR + 1
ulalume3@0 1541 atrue = (1-LDRTrue)/(1+LDRTrue)
ulalume3@0 1542 # ----- Forward simulated signals and LDRsim with atrue; from input file
ulalume3@0 1543 It = TaT*TiT*TiO*TiE*(GT+atrue*HT) # TaT*TiT*TiC*TiO*IinL*(GT+atrue*HT)
ulalume3@0 1544 Ir = TaR*TiR*TiO*TiE*(GR+atrue*HR) # TaR*TiR*TiC*TiO*IinL*(GR+atrue*HR)
ulalume3@0 1545
ulalume3@0 1546 # LDRsim = 1/Eta*Ir/It # simulated LDR* with Y from input file
ulalume3@0 1547 LDRsim = Ir/It # simulated uncorrected LDR with Y from input file
ulalume3@0 1548 '''
ulalume3@0 1549 if Y == 1.:
ulalume3@0 1550 LDRsimx = LDRsim
ulalume3@0 1551 LDRsimx2 = LDRsim2
ulalume3@0 1552 else:
ulalume3@0 1553 LDRsimx = 1./LDRsim
ulalume3@0 1554 LDRsimx2 = 1./LDRsim2
ulalume3@0 1555 '''
ulalume3@0 1556 # ----- Backward correction
ulalume3@0 1557 # Corrected LDRCorr from forward simulated LDRsim (atrue) with assumed true G0,H0,K0
ulalume3@0 1558 LDRCorr = (LDRsim*K0/Etax*(GT0+HT0)-(GR0+HR0))/((GR0-HR0)-LDRsim*K0/Etax*(GT0-HT0))
ulalume3@0 1559
ulalume3@0 1560 # -- F11corr from It and Ir and calibration EtaX
ulalume3@0 1561 Text1 = "F11corr from It and Ir with calibration EtaX: x-axis: F11corr(LDRtrue) / F11corr(LDRtrue = 0.004) - 1"
ulalume3@0 1562 F11corr = 1/(TiO*TiE)*((HR0*Etax/K0*It/TTa-HT0*Ir/TRa)/(HR0*GT0-HT0*GR0)) # IL = 1 Eq.(64)
ulalume3@0 1563
ulalume3@0 1564 #Text1 = "F11corr from It and Ir without corrections but with calibration EtaX: x-axis: F11corr(LDRtrue) devided by F11corr(LDRtrue = 0.004)"
ulalume3@0 1565 #F11corr = 0.5/(TiO*TiE)*(Etax*It/TTa+Ir/TRa) # IL = 1 Eq.(64)
ulalume3@0 1566
ulalume3@0 1567 # -- It from It only with atrue without corrections - for BERTHA (and PollyXTs)
ulalume3@0 1568 #Text1 = " x-axis: IT(LDRtrue) / IT(LDRtrue = 0.004) - 1"
ulalume3@0 1569 #F11corr = It/(TaT*TiT*TiO*TiE) #/(TaT*TiT*TiO*TiE*(GT0+atrue*HT0))
ulalume3@0 1570 # !!! see below line 1673ff
ulalume3@0 1571
ulalume3@0 1572 aF11corr[iLDR,iN] = F11corr
ulalume3@0 1573 aA[iLDR,iN] = LDRCorr
ulalume3@0 1574
ulalume3@0 1575 aX[0,iN] = GR
ulalume3@0 1576 aX[1,iN] = GT
ulalume3@0 1577 aX[2,iN] = HR
ulalume3@0 1578 aX[3,iN] = HT
ulalume3@0 1579 aX[4,iN] = K
ulalume3@0 1580
ulalume3@0 1581 aLDRCal[iN] = iLDRCal
ulalume3@0 1582 aERaT[iN] = iERaT
ulalume3@0 1583 aERaR[iN] = iERaR
ulalume3@0 1584 aRotaT[iN] = iRotaT
ulalume3@0 1585 aRotaR[iN] = iRotaR
ulalume3@0 1586 aRetT[iN] = iRetT
ulalume3@0 1587 aRetR[iN] = iRetR
ulalume3@0 1588
ulalume3@0 1589 aRotL[iN] = iRotL
ulalume3@0 1590 aRotE[iN] = iRotE
ulalume3@0 1591 aRetE[iN] = iRetE
ulalume3@0 1592 aRotO[iN] = iRotO
ulalume3@0 1593 aRetO[iN] = iRetO
ulalume3@0 1594 aRotC[iN] = iRotC
ulalume3@0 1595 aRetC[iN] = iRetC
ulalume3@0 1596 aDiO[iN] = iDiO
ulalume3@0 1597 aDiE[iN] = iDiE
ulalume3@0 1598 aDiC[iN] = iDiC
ulalume3@0 1599 aTP[iN] = iTP
ulalume3@0 1600 aTS[iN] = iTS
ulalume3@0 1601 aRP[iN] = iRP
ulalume3@0 1602 aRS[iN] = iRS
ulalume3@0 1603
ulalume3@0 1604 # --- END loop
ulalume3@0 1605 btime = clock()
ulalume3@0 1606 print("\r done in ", "{0:5.0f}".format(btime-atime), "sec") #, end="\r")
ulalume3@0 1607
ulalume3@0 1608 # --- Plot -----------------------------------------------------------------
ulalume3@0 1609 #sns.set_style("whitegrid")
ulalume3@0 1610 #sns.set_palette("bright", 6)
ulalume3@0 1611
ulalume3@0 1612 '''
ulalume3@0 1613 fig2 = plt.figure()
ulalume3@0 1614 plt.plot(aA[2,:],'b.')
ulalume3@0 1615 plt.plot(aA[3,:],'r.')
ulalume3@0 1616 plt.plot(aA[4,:],'g.')
ulalume3@0 1617 #plt.plot(aA[6,:],'c.')
ulalume3@0 1618 plt.show
ulalume3@0 1619 '''
ulalume3@0 1620
ulalume3@0 1621
ulalume3@0 1622 # Plot LDR
ulalume3@0 1623 def PlotSubHist(aVar, aX, X0, daX, iaX, naX):
ulalume3@0 1624
ulalume3@0 1625 fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2))
ulalume3@0 1626 iLDR = -1
ulalume3@0 1627 for LDRTrue in LDRrange:
ulalume3@0 1628 iLDR = iLDR + 1
ulalume3@0 1629
ulalume3@0 1630 LDRmin[iLDR] = np.min(aA[iLDR,:])
ulalume3@0 1631 LDRmax[iLDR] = np.max(aA[iLDR,:])
ulalume3@0 1632 Rmin = LDRmin[iLDR] * 0.995 # np.min(aA[iLDR,:]) * 0.995
ulalume3@0 1633 Rmax = LDRmax[iLDR] * 1.005 # np.max(aA[iLDR,:]) * 1.005
ulalume3@0 1634
ulalume3@0 1635 #plt.subplot(5,2,iLDR+1)
ulalume3@0 1636 plt.subplot(1,5,iLDR+1)
ulalume3@0 1637 (n, bins, patches) = plt.hist(aA[iLDR,:],
ulalume3@0 1638 bins=100, log=False,
ulalume3@0 1639 range=[Rmin, Rmax],
ulalume3@0 1640 alpha=0.5, normed=False, color = '0.5', histtype='stepfilled')
ulalume3@0 1641
ulalume3@0 1642 for iaX in range(-naX,naX+1):
ulalume3@0 1643 plt.hist(aA[iLDR,aX == iaX],
ulalume3@0 1644 range=[Rmin, Rmax],
ulalume3@0 1645 bins=100, log=False, alpha=0.3, normed=False, histtype='stepfilled', label = str(round(X0 + iaX*daX/naX,5)))
ulalume3@0 1646
ulalume3@0 1647 if (iLDR == 2): plt.legend()
ulalume3@0 1648
ulalume3@0 1649 plt.tick_params(axis='both', labelsize=9)
ulalume3@0 1650 plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2)
ulalume3@0 1651
ulalume3@0 1652 #plt.title(LID + ' ' + aVar, fontsize=18)
ulalume3@0 1653 #plt.ylabel('frequency', fontsize=10)
ulalume3@0 1654 #plt.xlabel('LDRcorr', fontsize=10)
ulalume3@0 1655 #fig.tight_layout()
ulalume3@0 1656 fig.suptitle(LID + ' ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar, fontsize=14, y=1.05)
ulalume3@0 1657 #plt.show()
ulalume3@0 1658 #fig.savefig(LID + '_' + aVar + '.png', dpi=150, bbox_inches='tight', pad_inches=0)
ulalume3@0 1659 #plt.close
ulalume3@0 1660 return
ulalume3@0 1661
ulalume3@0 1662 if (nRotL > 0): PlotSubHist("RotL", aRotL, RotL0, dRotL, iRotL, nRotL)
ulalume3@0 1663 if (nRetE > 0): PlotSubHist("RetE", aRetE, RetE0, dRetE, iRetE, nRetE)
ulalume3@0 1664 if (nRotE > 0): PlotSubHist("RotE", aRotE, RotE0, dRotE, iRotE, nRotE)
ulalume3@0 1665 if (nDiE > 0): PlotSubHist("DiE", aDiE, DiE0, dDiE, iDiE, nDiE)
ulalume3@0 1666 if (nRetO > 0): PlotSubHist("RetO", aRetO, RetO0, dRetO, iRetO, nRetO)
ulalume3@0 1667 if (nRotO > 0): PlotSubHist("RotO", aRotO, RotO0, dRotO, iRotO, nRotO)
ulalume3@0 1668 if (nDiO > 0): PlotSubHist("DiO", aDiO, DiO0, dDiO, iDiO, nDiO)
ulalume3@0 1669 if (nDiC > 0): PlotSubHist("DiC", aDiC, DiC0, dDiC, iDiC, nDiC)
ulalume3@0 1670 if (nRotC > 0): PlotSubHist("RotC", aRotC, RotC0, dRotC, iRotC, nRotC)
ulalume3@0 1671 if (nRetC > 0): PlotSubHist("RetC", aRetC, RetC0, dRetC, iRetC, nRetC)
ulalume3@0 1672 if (nTP > 0): PlotSubHist("TP", aTP, TP0, dTP, iTP, nTP)
ulalume3@0 1673 if (nTS > 0): PlotSubHist("TS", aTS, TS0, dTS, iTS, nTS)
ulalume3@0 1674 if (nRP > 0): PlotSubHist("RP", aRP, RP0, dRP, iRP, nRP)
ulalume3@0 1675 if (nRS > 0): PlotSubHist("RS", aRS, RS0, dRS, iRS, nRS)
ulalume3@0 1676 if (nRetT > 0): PlotSubHist("RetT", aRetT, RetT0, dRetT, iRetT, nRetT)
ulalume3@0 1677 if (nRetR > 0): PlotSubHist("RetR", aRetR, RetR0, dRetR, iRetR, nRetR)
ulalume3@0 1678 if (nERaT > 0): PlotSubHist("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT)
ulalume3@0 1679 if (nERaR > 0): PlotSubHist("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR)
ulalume3@0 1680 if (nRotaT > 0): PlotSubHist("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT)
ulalume3@0 1681 if (nRotaR > 0): PlotSubHist("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR)
ulalume3@0 1682 if (nLDRCal > 0): PlotSubHist("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal)
ulalume3@0 1683
ulalume3@0 1684 plt.show()
ulalume3@0 1685 plt.close
ulalume3@0 1686
ulalume3@0 1687 print()
ulalume3@0 1688 #print("IT(LDRtrue) devided by IT(LDRtrue = 0.004)")
ulalume3@0 1689 print(Text1)
ulalume3@0 1690 print()
ulalume3@0 1691
ulalume3@0 1692 iLDR = 5
ulalume3@0 1693 for LDRTrue in LDRrange:
ulalume3@0 1694 iLDR = iLDR - 1
ulalume3@0 1695 aF11corr[iLDR,:] = aF11corr[iLDR,:] / aF11corr[0,:] - 1.0
ulalume3@0 1696
ulalume3@0 1697 # Plot F11
ulalume3@0 1698 def PlotSubHistF11(aVar, aX, X0, daX, iaX, naX):
ulalume3@0 1699
ulalume3@0 1700 fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2))
ulalume3@0 1701 iLDR = -1
ulalume3@0 1702 for LDRTrue in LDRrange:
ulalume3@0 1703 iLDR = iLDR + 1
ulalume3@0 1704
ulalume3@0 1705 '''
ulalume3@0 1706 F11min[iLDR] = np.min(aF11corr[iLDR,:])
ulalume3@0 1707 F11max[iLDR] = np.max(aF11corr[iLDR,:])
ulalume3@0 1708 Rmin = F11min[iLDR] * 0.995 # np.min(aA[iLDR,:]) * 0.995
ulalume3@0 1709 Rmax = F11max[iLDR] * 1.005 # np.max(aA[iLDR,:]) * 1.005
ulalume3@0 1710 '''
ulalume3@0 1711 #Rmin = 0.8
ulalume3@0 1712 #Rmax = 1.2
ulalume3@0 1713
ulalume3@0 1714 #plt.subplot(5,2,iLDR+1)
ulalume3@0 1715 plt.subplot(1,5,iLDR+1)
ulalume3@0 1716 (n, bins, patches) = plt.hist(aF11corr[iLDR,:],
ulalume3@0 1717 bins=100, log=False,
ulalume3@0 1718 alpha=0.5, normed=False, color = '0.5', histtype='stepfilled')
ulalume3@0 1719
ulalume3@0 1720 for iaX in range(-naX,naX+1):
ulalume3@0 1721 plt.hist(aF11corr[iLDR,aX == iaX],
ulalume3@0 1722 bins=100, log=False, alpha=0.3, normed=False, histtype='stepfilled', label = str(round(X0 + iaX*daX/naX,5)))
ulalume3@0 1723
ulalume3@0 1724 if (iLDR == 2): plt.legend()
ulalume3@0 1725
ulalume3@0 1726 plt.tick_params(axis='both', labelsize=9)
ulalume3@0 1727 #plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2)
ulalume3@0 1728
ulalume3@0 1729 #plt.title(LID + ' ' + aVar, fontsize=18)
ulalume3@0 1730 #plt.ylabel('frequency', fontsize=10)
ulalume3@0 1731 #plt.xlabel('LDRcorr', fontsize=10)
ulalume3@0 1732 #fig.tight_layout()
ulalume3@0 1733 fig.suptitle(LID + ' ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar, fontsize=14, y=1.05)
ulalume3@0 1734 #plt.show()
ulalume3@0 1735 #fig.savefig(LID + '_' + aVar + '.png', dpi=150, bbox_inches='tight', pad_inches=0)
ulalume3@0 1736 #plt.close
ulalume3@0 1737 return
ulalume3@0 1738
ulalume3@0 1739 if (nRotL > 0): PlotSubHistF11("RotL", aRotL, RotL0, dRotL, iRotL, nRotL)
ulalume3@0 1740 if (nRetE > 0): PlotSubHistF11("RetE", aRetE, RetE0, dRetE, iRetE, nRetE)
ulalume3@0 1741 if (nRotE > 0): PlotSubHistF11("RotE", aRotE, RotE0, dRotE, iRotE, nRotE)
ulalume3@0 1742 if (nDiE > 0): PlotSubHistF11("DiE", aDiE, DiE0, dDiE, iDiE, nDiE)
ulalume3@0 1743 if (nRetO > 0): PlotSubHistF11("RetO", aRetO, RetO0, dRetO, iRetO, nRetO)
ulalume3@0 1744 if (nRotO > 0): PlotSubHistF11("RotO", aRotO, RotO0, dRotO, iRotO, nRotO)
ulalume3@0 1745 if (nDiO > 0): PlotSubHistF11("DiO", aDiO, DiO0, dDiO, iDiO, nDiO)
ulalume3@0 1746 if (nDiC > 0): PlotSubHistF11("DiC", aDiC, DiC0, dDiC, iDiC, nDiC)
ulalume3@0 1747 if (nRotC > 0): PlotSubHistF11("RotC", aRotC, RotC0, dRotC, iRotC, nRotC)
ulalume3@0 1748 if (nRetC > 0): PlotSubHistF11("RetC", aRetC, RetC0, dRetC, iRetC, nRetC)
ulalume3@0 1749 if (nTP > 0): PlotSubHistF11("TP", aTP, TP0, dTP, iTP, nTP)
ulalume3@0 1750 if (nTS > 0): PlotSubHistF11("TS", aTS, TS0, dTS, iTS, nTS)
ulalume3@0 1751 if (nRP > 0): PlotSubHistF11("RP", aRP, RP0, dRP, iRP, nRP)
ulalume3@0 1752 if (nRS > 0): PlotSubHistF11("RS", aRS, RS0, dRS, iRS, nRS)
ulalume3@0 1753 if (nRetT > 0): PlotSubHistF11("RetT", aRetT, RetT0, dRetT, iRetT, nRetT)
ulalume3@0 1754 if (nRetR > 0): PlotSubHistF11("RetR", aRetR, RetR0, dRetR, iRetR, nRetR)
ulalume3@0 1755 if (nERaT > 0): PlotSubHistF11("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT)
ulalume3@0 1756 if (nERaR > 0): PlotSubHistF11("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR)
ulalume3@0 1757 if (nRotaT > 0): PlotSubHistF11("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT)
ulalume3@0 1758 if (nRotaR > 0): PlotSubHistF11("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR)
ulalume3@0 1759 if (nLDRCal > 0): PlotSubHistF11("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal)
ulalume3@0 1760
ulalume3@0 1761 plt.show()
ulalume3@0 1762 plt.close
ulalume3@0 1763 '''
ulalume3@0 1764 # only histogram
ulalume3@0 1765 #print("******************* " + aVar + " *******************")
ulalume3@0 1766 fig, ax = plt.subplots(nrows=5, ncols=2, sharex=True, sharey=True, figsize=(10, 10))
ulalume3@0 1767 iLDR = -1
ulalume3@0 1768 for LDRTrue in LDRrange:
ulalume3@0 1769 iLDR = iLDR + 1
ulalume3@0 1770 LDRmin[iLDR] = np.min(aA[iLDR,:])
ulalume3@0 1771 LDRmax[iLDR] = np.max(aA[iLDR,:])
ulalume3@0 1772 Rmin = np.min(aA[iLDR,:]) * 0.999
ulalume3@0 1773 Rmax = np.max(aA[iLDR,:]) * 1.001
ulalume3@0 1774 plt.subplot(5,2,iLDR+1)
ulalume3@0 1775 (n, bins, patches) = plt.hist(aA[iLDR,:],
ulalume3@0 1776 range=[Rmin, Rmax],
ulalume3@0 1777 bins=200, log=False, alpha=0.2, normed=False, color = '0.5', histtype='stepfilled')
ulalume3@0 1778 plt.tick_params(axis='both', labelsize=9)
ulalume3@0 1779 plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2)
ulalume3@0 1780 plt.show()
ulalume3@0 1781 plt.close
ulalume3@0 1782 '''
ulalume3@0 1783
ulalume3@0 1784 # --- Plot LDRmin, LDRmax
ulalume3@0 1785 fig2 = plt.figure()
ulalume3@0 1786 plt.plot(LDRrange,LDRmax-LDRrange, linewidth=2.0, color='b')
ulalume3@0 1787 plt.plot(LDRrange,LDRmin-LDRrange, linewidth=2.0, color='g')
ulalume3@0 1788
ulalume3@0 1789 plt.xlabel('LDRtrue', fontsize=18)
ulalume3@0 1790 plt.ylabel('LDRTrue-LDRmin, LDRTrue-LDRmax', fontsize=14)
ulalume3@0 1791 plt.title(LID + ' ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]), fontsize=18)
ulalume3@0 1792 #plt.ylimit(-0.07, 0.07)
ulalume3@0 1793 plt.show()
ulalume3@0 1794 plt.close
ulalume3@0 1795
ulalume3@0 1796 # --- Save LDRmin, LDRmax to file
ulalume3@0 1797 # http://stackoverflow.com/questions/4675728/redirect-stdout-to-a-file-in-python
ulalume3@0 1798 with open('LDR_min_max_ver7_' + LID + '.dat', 'w') as f:
ulalume3@0 1799 with redirect_stdout(f):
ulalume3@0 1800 print(LID)
ulalume3@0 1801 print("LDRtrue, LDRmin, LDRmax")
ulalume3@0 1802 for i in range(len(LDRrange)):
ulalume3@0 1803 print("{0:7.4f},{1:7.4f},{2:7.4f}".format(LDRrange[i], LDRmin[i], LDRmax[i]))
ulalume3@0 1804
ulalume3@0 1805 '''
ulalume3@0 1806 # --- Plot K over LDRCal
ulalume3@0 1807 fig3 = plt.figure()
ulalume3@0 1808 plt.plot(LDRCal0+aLDRCal*dLDRCal/nLDRCal,aX[4,:], linewidth=2.0, color='b')
ulalume3@0 1809
ulalume3@0 1810 plt.xlabel('LDRCal', fontsize=18)
ulalume3@0 1811 plt.ylabel('K', fontsize=14)
ulalume3@0 1812 plt.title(LID, fontsize=18)
ulalume3@0 1813 plt.show()
ulalume3@0 1814 plt.close
ulalume3@0 1815 '''
ulalume3@0 1816
ulalume3@0 1817 # Additional plot routines ======>
ulalume3@0 1818 '''
ulalume3@0 1819 #******************************************************************************
ulalume3@0 1820 # 1. Plot LDRcorrected - LDR(measured Icross/Iparallel)
ulalume3@0 1821 LDRa = np.arange(1.,100.)*0.005
ulalume3@0 1822 LDRCorra = np.arange(1.,100.)
ulalume3@0 1823 if Y == - 1.: LDRa = 1./LDRa
ulalume3@0 1824 LDRCorra = (1./Eta*LDRa*(GT+HT)-(GR+HR))/((GR-HR)-1./Eta*LDRa*(GT-HT))
ulalume3@0 1825 if Y == - 1.: LDRa = 1./LDRa
ulalume3@0 1826 #
ulalume3@0 1827 #fig = plt.figure()
ulalume3@0 1828 plt.plot(LDRa,LDRCorra-LDRa)
ulalume3@0 1829 plt.plot([0.,0.5],[0.,0.5])
ulalume3@0 1830 plt.suptitle('LDRcorrected - LDR(measured Icross/Iparallel)', fontsize=16)
ulalume3@0 1831 plt.xlabel('LDR', fontsize=18)
ulalume3@0 1832 plt.ylabel('LDRCorr - LDR', fontsize=16)
ulalume3@0 1833 #plt.savefig('test.png')
ulalume3@0 1834 #
ulalume3@0 1835 '''
ulalume3@0 1836 '''
ulalume3@0 1837 #******************************************************************************
ulalume3@0 1838 # 2. Plot LDRsim (simulated measurements without corrections = Icross/Iparallel) over LDRtrue
ulalume3@0 1839 LDRa = np.arange(1.,100.)*0.005
ulalume3@0 1840 LDRsima = np.arange(1.,100.)
ulalume3@0 1841
ulalume3@0 1842 atruea = (1.-LDRa)/(1+LDRa)
ulalume3@0 1843 Ita = TiT*TiO*IinL*(GT+atruea*HT)
ulalume3@0 1844 Ira = TiR*TiO*IinL*(GR+atruea*HR)
ulalume3@0 1845 LDRsima = Ira/Ita # simulated uncorrected LDR with Y from input file
ulalume3@0 1846 if Y == -1.: LDRsima = 1./LDRsima
ulalume3@0 1847 #
ulalume3@0 1848 #fig = plt.figure()
ulalume3@0 1849 plt.plot(LDRa,LDRsima)
ulalume3@0 1850 plt.plot([0.,0.5],[0.,0.5])
ulalume3@0 1851 plt.suptitle('LDRsim (simulated measurements without corrections = Icross/Iparallel) over LDRtrue', fontsize=10)
ulalume3@0 1852 plt.xlabel('LDRtrue', fontsize=18)
ulalume3@0 1853 plt.ylabel('LDRsim', fontsize=16)
ulalume3@0 1854 #plt.savefig('test.png')
ulalume3@0 1855 #
ulalume3@0 1856 '''

mercurial