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