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