Mercurial > public > atmospheric_lidar_ghk / file revision
lidar_correction_ghk.py@453b23dd7f94
lidar_correction_ghk.py
Mon, 14 Nov 2016 23:39:33 +0100
- author
- Volker Freudenthaler <volker.freudenthaler@lmu.de>
- date
- Mon, 14 Nov 2016 23:39:33 +0100
- changeset 11
- 453b23dd7f94
- parent 10
- a69e0c8ca55c
- child 13
- f08818615e3a
- permissions
- -rw-r--r--
lidar_correction_ghk.py works with Python 2.7 and 3.5, with and without seaborne installed. No need to change something in the code.
# -*- coding: utf-8 -*- """ Copyright 2016 Volker Freudenthaler Licensed under the EUPL, Version 1.1 only (the "Licence"). You may not use this work except in compliance with the Licence. A copy of the licence is distributed with the code. Alternatively, you may obtain a copy of the Licence at: https://joinup.ec.europa.eu/community/eupl/og_page/eupl Unless required by applicable law or agreed to in writing, software distributed under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the Licence for the specific language governing permissions and limitations under the Licence. Equation reference: http://www.atmos-meas-tech-discuss.net/amt-2015-338/amt-2015-338.pdf With equations code from Appendix C Python 3.4.2 """ #!/usr/bin/env python3 # Comment: The code works with Python 2.7 with the help of following line, which enables Python2 to correctly interpret the Python 3 print statements. from __future__ import print_function #import math import numpy as np import sys import os # Comment: the seaborn library makes nicer plots, but the code works also without it. try: import seaborn as sns sns_loaded = True except ImportError: sns_loaded = False import matplotlib.pyplot as plt from time import clock #from matplotlib.backends.backend_pdf import PdfPages #pdffile = '{}.pdf'.format('path') #pp = PdfPages(pdffile) ## pp.savefig can be called multiple times to save to multiple pages #pp.savefig() #pp.close() from contextlib import contextmanager @contextmanager def redirect_stdout(new_target): old_target, sys.stdout = sys.stdout, new_target # replace sys.stdout try: yield new_target # run some code with the replaced stdout finally: sys.stdout.flush() sys.stdout = old_target # restore to the previous value ''' real_raw_input = vars(__builtins__).get('raw_input',input) ''' try: import __builtin__ input = getattr(__builtin__, 'raw_input') except (ImportError, AttributeError): pass from distutils.util import strtobool def user_yes_no_query(question): sys.stdout.write('%s [y/n]\n' % question) while True: try: return strtobool(input().lower()) except ValueError: sys.stdout.write('Please respond with \'y\' or \'n\'.\n') #if user_yes_no_query('want to exit?') == 1: sys.exit() abspath = os.path.abspath(__file__) dname = os.path.dirname(abspath) fname = os.path.basename(abspath) os.chdir(dname) #PrintToOutputFile = True sqr05 = 0.5**0.5 # ---- Initial definition of variables; the actual values will be read in with exec(open('./optic_input.py').read()) below LID = "internal" EID = "internal" # --- IL Laser IL and +-Uncertainty bL = 1. #degree of linear polarization; default 1 RotL, dRotL, nRotL = 0.0, 0.0, 1 #alpha; rotation of laser polarization in degrees; default 0 # --- ME Emitter and +-Uncertainty DiE, dDiE, nDiE = 0., 0.00, 1 # Diattenuation TiE = 1. # Unpolarized transmittance RetE, dRetE, nRetE = 0., 180.0, 0 # Retardance in degrees RotE, dRotE, nRotE = 0., 0.0, 0 # beta: Rotation of optical element in degrees # --- MO Receiver Optics including telescope DiO, dDiO, nDiO = -0.055, 0.003, 1 TiO = 0.9 RetO, dRetO, nRetO = 0., 180.0, 2 RotO, dRotO, nRotO = 0., 0.1, 1 #gamma # --- PBS MT transmitting path defined with (TS,TP); and +-Uncertainty TP, dTP, nTP = 0.98, 0.02, 1 TS, dTS, nTS = 0.001, 0.001, 1 TiT = 0.5 * (TP + TS) DiT = (TP-TS)/(TP+TS) # PolFilter RetT, dRetT, nRetT = 0., 180., 0 ERaT, dERaT, nERaT = 0.001, 0.001, 1 RotaT, dRotaT, nRotaT = 0., 3., 1 DaT = (1-ERaT)/(1+ERaT) TaT = 0.5*(1+ERaT) # --- PBS MR reflecting path defined with (RS,RP); and +-Uncertainty RS, dRS, nRS = 1 - TS, 0., 0 RP, dRP, nRP = 1 - TP, 0., 0 TiR = 0.5 * (RP + RS) DiR = (RP-RS)/(RP+RS) # PolFilter RetR, dRetR, nRetR = 0., 180., 0 ERaR, dERaR, nERaR = 0.001, 0.001, 1 RotaR,dRotaR,nRotaR = 90., 3., 1 DaR = (1-ERaR)/(1+ERaR) TaR = 0.5*(1+ERaR) # Parellel signal detected in the transmitted channel => Y = 1, or in the reflected channel => Y = -1 Y = -1. # Calibrator = type defined by matrix values LocC = 4 # location of calibrator: behind laser = 1; behind emitter = 2; before receiver = 3; before PBS = 4 TypeC = 3 # linear polarizer calibrator # example with extinction ratio 0.001 DiC, dDiC, nDiC = 1.0, 0., 0 # ideal 1.0 TiC = 0.5 # ideal 0.5 RetC, dRetC, nRetC = 0., 0., 0 RotC, dRotC, nRotC = 0.0, 0.1, 0 #constant calibrator offset epsilon RotationErrorEpsilonForNormalMeasurements = False # is in general False for TypeC == 3 calibrator # Rotation error without calibrator: if False, then epsilon = 0 for normal measurements RotationErrorEpsilonForNormalMeasurements = True # LDRCal assumed atmospheric linear depolarization ratio during the calibration measurements (first guess) LDRCal0,dLDRCal,nLDRCal= 0.25, 0.04, 1 LDRCal = LDRCal0 # measured LDRm will be corrected with calculated parameters LDRmeas = 0.015 # LDRtrue for simulation of measurement => LDRsim LDRtrue = 0.5 LDRtrue2 = 0.004 # Initialize other values to 0 ER, nER, dER = 0.001, 0, 0.001 K = 0. Km = 0. Kp = 0. LDRcorr = 0. Eta = 0. Ir = 0. It = 0. h = 1. Loc = ['', 'behind laser', 'behind emitter', 'before receiver', 'before PBS'] Type = ['', 'mechanical rotator', 'hwp rotator', 'linear polarizer', 'qwp rotator', 'circular polarizer', 'real HWP +-22.5°'] dY = ['reflected channel', '', 'transmitted channel'] # end of initial definition of variables # ******************************************************************************************************************************* # --- Read actual lidar system parameters from ./optic_input.py (must be in the same directory) InputFile = 'optic_input_example_lidar.py' #InputFile = 'optic_input_ver6e_POLIS_355.py' #InputFile = 'optic_input_ver6e_POLIS_355_JA.py' #InputFile = 'optic_input_ver6c_POLIS_532.py' #InputFile = 'optic_input_ver6e_POLIS_532.py' #InputFile = 'optic_input_ver8c_POLIS_532.py' #InputFile = 'optic_input_ver6e_MUSA.py' #InputFile = 'optic_input_ver6e_MUSA_JA.py' #InputFile = 'optic_input_ver6e_PollyXTSea.py' #InputFile = 'optic_input_ver6e_PollyXTSea_JA.py' #InputFile = 'optic_input_ver6e_PollyXT_RALPH.py' #InputFile = 'optic_input_ver8c_PollyXT_RALPH.py' #InputFile = 'optic_input_ver8c_PollyXT_RALPH_2.py' #InputFile = 'optic_input_ver8c_PollyXT_RALPH_3.py' #InputFile = 'optic_input_ver8c_PollyXT_RALPH_4.py' #InputFile = 'optic_input_ver8c_PollyXT_RALPH_5.py' #InputFile = 'optic_input_ver8c_PollyXT_RALPH_6.py' #InputFile = 'optic_input_ver8c_PollyXT_RALPH_7.py' #InputFile = 'optic_input_ver8a_MOHP_DPL_355.py' #InputFile = 'optic_input_ver9_MOHP_DPL_355.py' #InputFile = 'optic_input_ver6e_RALI.py' #InputFile = 'optic_input_ver6e_RALI_JA.py' #InputFile = 'optic_input_ver6e_RALI_new.py' #InputFile = 'optic_input_ver6e_RALI_act.py' #InputFile = 'optic_input_ver6e_MULHACEN.py' #InputFile = 'optic_input_ver6e_MULHACEN_JA.py' #InputFile = 'optic_input_ver6e-IPRAL.py' #InputFile = 'optic_input_ver6e-IPRAL_JA.py' #InputFile = 'optic_input_ver6e-LB21.py' #InputFile = 'optic_input_ver6e-LB21_JA.py' #InputFile = 'optic_input_ver6e_Bertha_b_355.py' #InputFile = 'optic_input_ver6e_Bertha_b_532.py' #InputFile = 'optic_input_ver6e_Bertha_b_1064.py' ''' print("From ", dname) print("Running ", fname) print("Reading input file ", InputFile, " for") ''' input_path = os.path.join('.', 'system_settings', InputFile) # this works with Python 2 - and 3? exec(open(input_path).read(), globals()) # end of read actual system parameters # --- Manual Parameter Change --- # (use for quick parameter changes without changing the input file ) #DiO = 0. #LDRtrue = 0.45 #LDRtrue2 = 0.004 #Y = -1 #LocC = 4 #location of calibrator: 1 = behind laser; 2 = behind emitter; 3 = before receiver; 4 = before PBS ##TypeC = 6 Don't change the TypeC here #RotationErrorEpsilonForNormalMeasurements = True #LDRCal = 0.25 #bL = 0.8 ## --- Errors RotL0, dRotL, nRotL = RotL, dRotL, nRotL DiE0, dDiE, nDiE = DiE, dDiE, nDiE RetE0, dRetE, nRetE = RetE, dRetE, nRetE RotE0, dRotE, nRotE = RotE, dRotE, nRotE DiO0, dDiO, nDiO = DiO, dDiO, nDiO RetO0, dRetO, nRetO = RetO, dRetO, nRetO RotO0, dRotO, nRotO = RotO, dRotO, nRotO DiC0, dDiC, nDiC = DiC, dDiC, nDiC RetC0, dRetC, nRetC = RetC, dRetC, nRetC RotC0, dRotC, nRotC = RotC, dRotC, nRotC TP0, dTP, nTP = TP, dTP, nTP TS0, dTS, nTS = TS, dTS, nTS RetT0, dRetT, nRetT = RetT, dRetT, nRetT ERaT0, dERaT, nERaT = ERaT, dERaT, nERaT RotaT0,dRotaT,nRotaT= RotaT,dRotaT,nRotaT RP0, dRP, nRP = RP, dRP, nRP RS0, dRS, nRS = RS, dRS, nRS RetR0, dRetR, nRetR = RetR, dRetR, nRetR ERaR0, dERaR, nERaR = ERaR, dERaR, nERaR RotaR0,dRotaR,nRotaR= RotaR,dRotaR,nRotaR LDRCal0,dLDRCal,nLDRCal=LDRCal,dLDRCal,nLDRCal #LDRCal0,dLDRCal,nLDRCal=LDRCal,dLDRCal,0 # ---------- End of manual parameter change RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC = RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0 TP, TS, RP, RS, ERaT, RotaT, RetT, ERaR, RotaR, RetR = TP0, TS0, RP0, RS0 , ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0 LDRCal = LDRCal0 DTa0, TTa0, DRa0, TRa0, LDRsimx, LDRCorr = 0,0,0,0,0,0 TiT = 0.5 * (TP + TS) DiT = (TP-TS)/(TP+TS) ZiT = (1. - DiT**2)**0.5 TiR = 0.5 * (RP + RS) DiR = (RP-RS)/(RP+RS) ZiR = (1. - DiR**2)**0.5 # -------------------------------------------------------- def Calc(RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC, TP, TS, RP, RS, ERaT, RotaT, RetT, ERaR, RotaR, RetR, LDRCal): # ---- Do the calculations of bra-ket vectors h = -1. if TypeC == 2 else 1 # from input file: assumed LDRCal for calibration measurements aCal = (1.-LDRCal)/(1+LDRCal) # from input file: measured LDRm and true LDRtrue, LDRtrue2 => #ameas = (1.-LDRmeas)/(1+LDRmeas) atrue = (1.-LDRtrue)/(1+LDRtrue) #atrue2 = (1.-LDRtrue2)/(1+LDRtrue2) # angles of emitter and laser and calibrator and receiver optics # RotL = alpha, RotE = beta, RotO = gamma, RotC = epsilon S2a = np.sin(2*np.deg2rad(RotL)) C2a = np.cos(2*np.deg2rad(RotL)) S2b = np.sin(2*np.deg2rad(RotE)) C2b = np.cos(2*np.deg2rad(RotE)) S2ab = np.sin(np.deg2rad(2*RotL-2*RotE)) C2ab = np.cos(np.deg2rad(2*RotL-2*RotE)) S2g = np.sin(np.deg2rad(2*RotO)) C2g = np.cos(np.deg2rad(2*RotO)) # Laser with Degree of linear polarization DOLP = bL IinL = 1. QinL = bL UinL = 0. VinL = (1. - bL**2)**0.5 # Stokes Input Vector rotation Eq. E.4 A = C2a*QinL - S2a*UinL B = S2a*QinL + C2a*UinL # Stokes Input Vector rotation Eq. E.9 C = C2ab*QinL - S2ab*UinL D = S2ab*QinL + C2ab*UinL # emitter optics CosE = np.cos(np.deg2rad(RetE)) SinE = np.sin(np.deg2rad(RetE)) ZiE = (1. - DiE**2)**0.5 WiE = (1. - ZiE*CosE) # Stokes Input Vector after emitter optics equivalent to Eq. E.9 with already rotated input vector from Eq. E.4 # b = beta IinE = (IinL + DiE*C) QinE = (C2b*DiE*IinL + A + S2b*(WiE*D - ZiE*SinE*VinL)) UinE = (S2b*DiE*IinL + B - C2b*(WiE*D - ZiE*SinE*VinL)) VinE = (-ZiE*SinE*D + ZiE*CosE*VinL) # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) IinF = IinE QinF = aCal*QinE UinF = -aCal*UinE VinF = (1.-2.*aCal)*VinE # receiver optics CosO = np.cos(np.deg2rad(RetO)) SinO = np.sin(np.deg2rad(RetO)) ZiO = (1. - DiO**2)**0.5 WiO = (1. - ZiO*CosO) # calibrator CosC = np.cos(np.deg2rad(RetC)) SinC = np.sin(np.deg2rad(RetC)) ZiC = (1. - DiC**2)**0.5 WiC = (1. - ZiC*CosC) # Stokes Input Vector before the polarising beam splitter Eq. E.31 A = C2g*QinE - S2g*UinE B = S2g*QinE + C2g*UinE IinP = (IinE + DiO*aCal*A) QinP = (C2g*DiO*IinE + aCal*QinE - S2g*(WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE)) UinP = (S2g*DiO*IinE - aCal*UinE + C2g*(WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE)) VinP = (ZiO*SinO*aCal*B + ZiO*CosO*(1-2*aCal)*VinE) #------------------------- # F11 assuemd to be = 1 => measured: F11m = IinP / IinE with atrue #F11sim = TiO*(IinE + DiO*atrue*A)/IinE #------------------------- # For PollyXT # analyser #RS = 1 - TS #RP = 1 - TP TiT = 0.5 * (TP + TS) DiT = (TP-TS)/(TP+TS) ZiT = (1. - DiT**2)**0.5 TiR = 0.5 * (RP + RS) DiR = (RP-RS)/(RP+RS) ZiR = (1. - DiR**2)**0.5 CosT = np.cos(np.deg2rad(RetT)) SinT = np.sin(np.deg2rad(RetT)) CosR = np.cos(np.deg2rad(RetR)) SinR = np.sin(np.deg2rad(RetR)) DaT = (1-ERaT)/(1+ERaT) DaR = (1-ERaR)/(1+ERaR) TaT = 0.5*(1+ERaT) TaR = 0.5*(1+ERaR) S2aT = np.sin(np.deg2rad(h*2*RotaT)) C2aT = np.cos(np.deg2rad(2*RotaT)) S2aR = np.sin(np.deg2rad(h*2*RotaR)) C2aR = np.cos(np.deg2rad(2*RotaR)) # Aanalyzer As before the PBS Eq. D.5 ATP1 = (1+C2aT*DaT*DiT) ATP2 = Y*(DiT+C2aT*DaT) ATP3 = Y*S2aT*DaT*ZiT*CosT ATP4 = S2aT*DaT*ZiT*SinT ATP = np.array([ATP1,ATP2,ATP3,ATP4]) ARP1 = (1+C2aR*DaR*DiR) ARP2 = Y*(DiR+C2aR*DaR) ARP3 = Y*S2aR*DaR*ZiR*CosR ARP4 = S2aR*DaR*ZiR*SinR ARP = np.array([ARP1,ARP2,ARP3,ARP4]) DTa = ATP2*Y/ATP1 DRa = ARP2*Y/ARP1 # ---- Calculate signals and correction parameters for diffeent locations and calibrators if LocC == 4: # Calibrator before the PBS #print("Calibrator location not implemented yet") #S2ge = np.sin(np.deg2rad(2*RotO + h*2*RotC)) #C2ge = np.cos(np.deg2rad(2*RotO + h*2*RotC)) S2e = np.sin(np.deg2rad(h*2*RotC)) C2e = np.cos(np.deg2rad(2*RotC)) # rotated AinP by epsilon Eq. C.3 ATP2e = C2e*ATP2 + S2e*ATP3 ATP3e = C2e*ATP3 - S2e*ATP2 ARP2e = C2e*ARP2 + S2e*ARP3 ARP3e = C2e*ARP3 - S2e*ARP2 ATPe = np.array([ATP1,ATP2e,ATP3e,ATP4]) ARPe = np.array([ARP1,ARP2e,ARP3e,ARP4]) # Stokes Input Vector before the polarising beam splitter Eq. E.31 A = C2g*QinE - S2g*UinE B = S2g*QinE + C2g*UinE #C = (WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE) Co = ZiO*SinO*VinE Ca = (WiO*B - 2*ZiO*SinO*VinE) #C = Co + aCal*Ca #IinP = (IinE + DiO*aCal*A) #QinP = (C2g*DiO*IinE + aCal*QinE - S2g*C) #UinP = (S2g*DiO*IinE - aCal*UinE + C2g*C) #VinP = (ZiO*SinO*aCal*B + ZiO*CosO*(1-2*aCal)*VinE) IinPo = IinE QinPo = (C2g*DiO*IinE - S2g*Co) UinPo = (S2g*DiO*IinE + C2g*Co) VinPo = ZiO*CosO*VinE IinPa = DiO*A QinPa = QinE - S2g*Ca UinPa = -UinE + C2g*Ca VinPa = ZiO*(SinO*B - 2*CosO*VinE) IinP = IinPo + aCal*IinPa QinP = QinPo + aCal*QinPa UinP = UinPo + aCal*UinPa VinP = VinPo + aCal*VinPa # Stokes Input Vector before the polarising beam splitter rotated by epsilon Eq. C.3 #QinPe = C2e*QinP + S2e*UinP #UinPe = C2e*UinP - S2e*QinP QinPoe = C2e*QinPo + S2e*UinPo UinPoe = C2e*UinPo - S2e*QinPo QinPae = C2e*QinPa + S2e*UinPa UinPae = C2e*UinPa - S2e*QinPa QinPe = C2e*QinP + S2e*UinP UinPe = C2e*UinP - S2e*QinP # Calibration signals and Calibration correction K from measurements with LDRCal / aCal if (TypeC == 2) or (TypeC == 1): # rotator calibration Eq. C.4 # parameters for calibration with aCal AT = ATP1*IinP + h*ATP4*VinP BT = ATP3e*QinP - h*ATP2e*UinP AR = ARP1*IinP + h*ARP4*VinP BR = ARP3e*QinP - h*ARP2e*UinP # Correction paremeters for normal measurements; they are independent of LDR if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out IS1 = np.array([IinPo,QinPo,UinPo,VinPo]) IS2 = np.array([IinPa,QinPa,UinPa,VinPa]) GT = np.dot(ATP,IS1) GR = np.dot(ARP,IS1) HT = np.dot(ATP,IS2) HR = np.dot(ARP,IS2) else: IS1 = np.array([IinPo,QinPo,UinPo,VinPo]) IS2 = np.array([IinPa,QinPa,UinPa,VinPa]) GT = np.dot(ATPe,IS1) GR = np.dot(ARPe,IS1) HT = np.dot(ATPe,IS2) HR = np.dot(ARPe,IS2) elif (TypeC == 3) or (TypeC == 4): # linear polariser calibration Eq. C.5 # parameters for calibration with aCal AT = ATP1*IinP + ATP3e*UinPe + ZiC*CosC*(ATP2e*QinPe + ATP4*VinP) BT = DiC*(ATP1*UinPe + ATP3e*IinP) - ZiC*SinC*(ATP2e*VinP - ATP4*QinPe) AR = ARP1*IinP + ARP3e*UinPe + ZiC*CosC*(ARP2e*QinPe + ARP4*VinP) BR = DiC*(ARP1*UinPe + ARP3e*IinP) - ZiC*SinC*(ARP2e*VinP - ARP4*QinPe) # Correction paremeters for normal measurements; they are independent of LDR if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out IS1 = np.array([IinPo,QinPo,UinPo,VinPo]) IS2 = np.array([IinPa,QinPa,UinPa,VinPa]) GT = np.dot(ATP,IS1) GR = np.dot(ARP,IS1) HT = np.dot(ATP,IS2) HR = np.dot(ARP,IS2) else: IS1e = np.array([IinPo+DiC*QinPoe,DiC*IinPo+QinPoe,ZiC*(CosC*UinPoe+SinC*VinPo),-ZiC*(SinC*UinPoe-CosC*VinPo)]) IS2e = np.array([IinPa+DiC*QinPae,DiC*IinPa+QinPae,ZiC*(CosC*UinPae+SinC*VinPa),-ZiC*(SinC*UinPae-CosC*VinPa)]) GT = np.dot(ATPe,IS1e) GR = np.dot(ARPe,IS1e) HT = np.dot(ATPe,IS2e) HR = np.dot(ARPe,IS2e) elif (TypeC == 6): # diattenuator calibration +-22.5° rotated_diattenuator_X22x5deg.odt # parameters for calibration with aCal AT = ATP1*IinP + sqr05*DiC*(ATP1*QinPe + ATP2e*IinP) + (1-0.5*WiC)*(ATP2e*QinPe + ATP3e*UinPe) + ZiC*(sqr05*SinC*(ATP3e*VinP-ATP4*UinPe) + ATP4*CosC*VinP) BT = sqr05*DiC*(ATP1*UinPe + ATP3e*IinP) + 0.5*WiC*(ATP2e*UinPe + ATP3e*QinPe) - sqr05*ZiC*SinC*(ATP2e*VinP - ATP4*QinPe) AR = ARP1*IinP + sqr05*DiC*(ARP1*QinPe + ARP2e*IinP) + (1-0.5*WiC)*(ARP2e*QinPe + ARP3e*UinPe) + ZiC*(sqr05*SinC*(ARP3e*VinP-ARP4*UinPe) + ARP4*CosC*VinP) BR = sqr05*DiC*(ARP1*UinPe + ARP3e*IinP) + 0.5*WiC*(ARP2e*UinPe + ARP3e*QinPe) - sqr05*ZiC*SinC*(ARP2e*VinP - ARP4*QinPe) # Correction paremeters for normal measurements; they are independent of LDR if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out IS1 = np.array([IinPo,QinPo,UinPo,VinPo]) IS2 = np.array([IinPa,QinPa,UinPa,VinPa]) GT = np.dot(ATP,IS1) GR = np.dot(ARP,IS1) HT = np.dot(ATP,IS2) HR = np.dot(ARP,IS2) else: IS1e = np.array([IinPo+DiC*QinPoe,DiC*IinPo+QinPoe,ZiC*(CosC*UinPoe+SinC*VinPo),-ZiC*(SinC*UinPoe-CosC*VinPo)]) IS2e = np.array([IinPa+DiC*QinPae,DiC*IinPa+QinPae,ZiC*(CosC*UinPae+SinC*VinPa),-ZiC*(SinC*UinPae-CosC*VinPa)]) GT = np.dot(ATPe,IS1e) GR = np.dot(ARPe,IS1e) HT = np.dot(ATPe,IS2e) HR = np.dot(ARPe,IS2e) else: print("Calibrator not implemented yet") sys.exit() elif LocC == 3: # C before receiver optics Eq.57 #S2ge = np.sin(np.deg2rad(2*RotO - 2*RotC)) #C2ge = np.cos(np.deg2rad(2*RotO - 2*RotC)) S2e = np.sin(np.deg2rad(2*RotC)) C2e = np.cos(np.deg2rad(2*RotC)) # As with C before the receiver optics (rotated_diattenuator_X22x5deg.odt) AF1 = np.array([1,C2g*DiO,S2g*DiO,0]) AF2 = np.array([C2g*DiO,1-S2g**2*WiO,S2g*C2g*WiO,-S2g*ZiO*SinO]) AF3 = np.array([S2g*DiO,S2g*C2g*WiO,1-C2g**2*WiO,C2g*ZiO*SinO]) AF4 = np.array([0,S2g*SinO,-C2g*SinO,CosO]) ATF = (ATP1*AF1+ATP2*AF2+ATP3*AF3+ATP4*AF4) ARF = (ARP1*AF1+ARP2*AF2+ARP3*AF3+ARP4*AF4) ATF2 = ATF[1] ATF3 = ATF[2] ARF2 = ARF[1] ARF3 = ARF[2] # rotated AinF by epsilon ATF1 = ATF[0] ATF4 = ATF[3] ATF2e = C2e*ATF[1] + S2e*ATF[2] ATF3e = C2e*ATF[2] - S2e*ATF[1] ARF1 = ARF[0] ARF4 = ARF[3] ARF2e = C2e*ARF[1] + S2e*ARF[2] ARF3e = C2e*ARF[2] - S2e*ARF[1] ATFe = np.array([ATF1,ATF2e,ATF3e,ATF4]) ARFe = np.array([ARF1,ARF2e,ARF3e,ARF4]) QinEe = C2e*QinE + S2e*UinE UinEe = C2e*UinE - S2e*QinE # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) IinF = IinE QinF = aCal*QinE UinF = -aCal*UinE VinF = (1.-2.*aCal)*VinE IinFo = IinE QinFo = 0. UinFo = 0. VinFo = VinE IinFa = 0. QinFa = QinE UinFa = -UinE VinFa = -2.*VinE # Stokes Input Vector before receiver optics rotated by epsilon Eq. C.3 QinFe = C2e*QinF + S2e*UinF UinFe = C2e*UinF - S2e*QinF QinFoe = C2e*QinFo + S2e*UinFo UinFoe = C2e*UinFo - S2e*QinFo QinFae = C2e*QinFa + S2e*UinFa UinFae = C2e*UinFa - S2e*QinFa # Calibration signals and Calibration correction K from measurements with LDRCal / aCal if (TypeC == 2) or (TypeC == 1): # rotator calibration Eq. C.4 # parameters for calibration with aCal AT = ATF1*IinF + ATF4*h*VinF BT = ATF3e*QinF - ATF2e*h*UinF AR = ARF1*IinF + ARF4*h*VinF BR = ARF3e*QinF - ARF2e*h*UinF # Correction paremeters for normal measurements; they are independent of LDR if (not RotationErrorEpsilonForNormalMeasurements): GT = ATF1*IinE + ATF4*VinE GR = ARF1*IinE + ARF4*VinE HT = ATF2*QinE - ATF3*UinE - ATF4*2*VinE HR = ARF2*QinE - ARF3*UinE - ARF4*2*VinE else: GT = ATF1*IinE + ATF4*h*VinE GR = ARF1*IinE + ARF4*h*VinE HT = ATF2e*QinE - ATF3e*h*UinE - ATF4*h*2*VinE HR = ARF2e*QinE - ARF3e*h*UinE - ARF4*h*2*VinE elif (TypeC == 3) or (TypeC == 4): # linear polariser calibration Eq. C.5 # p = +45°, m = -45° IF1e = np.array([IinF, ZiC*CosC*QinFe, UinFe, ZiC*CosC*VinF]) IF2e = np.array([DiC*UinFe, -ZiC*SinC*VinF, DiC*IinF, ZiC*SinC*QinFe]) AT = np.dot(ATFe,IF1e) AR = np.dot(ARFe,IF1e) BT = np.dot(ATFe,IF2e) BR = np.dot(ARFe,IF2e) # Correction paremeters for normal measurements; they are independent of LDR --- the same as for TypeC = 6 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out IS1 = np.array([IinE,0,0,VinE]) IS2 = np.array([0,QinE,-UinE,-2*VinE]) GT = np.dot(ATF,IS1) GR = np.dot(ARF,IS1) HT = np.dot(ATF,IS2) HR = np.dot(ARF,IS2) else: IS1e = np.array([IinFo+DiC*QinFoe,DiC*IinFo+QinFoe,ZiC*(CosC*UinFoe+SinC*VinFo),-ZiC*(SinC*UinFoe-CosC*VinFo)]) IS2e = np.array([IinFa+DiC*QinFae,DiC*IinFa+QinFae,ZiC*(CosC*UinFae+SinC*VinFa),-ZiC*(SinC*UinFae-CosC*VinFa)]) GT = np.dot(ATFe,IS1e) GR = np.dot(ARFe,IS1e) HT = np.dot(ATFe,IS2e) HR = np.dot(ARFe,IS2e) elif (TypeC == 6): # diattenuator calibration +-22.5° rotated_diattenuator_X22x5deg.odt # parameters for calibration with aCal IF1e = np.array([IinF+sqr05*DiC*QinFe, sqr05*DiC*IinF+(1-0.5*WiC)*QinFe, (1-0.5*WiC)*UinFe+sqr05*ZiC*SinC*VinF, -sqr05*ZiC*SinC*UinFe+ZiC*CosC*VinF]) IF2e = np.array([sqr05*DiC*UinFe, 0.5*WiC*UinFe-sqr05*ZiC*SinC*VinF, sqr05*DiC*IinF+0.5*WiC*QinFe, sqr05*ZiC*SinC*QinFe]) AT = np.dot(ATFe,IF1e) AR = np.dot(ARFe,IF1e) BT = np.dot(ATFe,IF2e) BR = np.dot(ARFe,IF2e) # Correction paremeters for normal measurements; they are independent of LDR if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out #IS1 = np.array([IinE,0,0,VinE]) #IS2 = np.array([0,QinE,-UinE,-2*VinE]) IS1 = np.array([IinFo,0,0,VinFo]) IS2 = np.array([0,QinFa,UinFa,VinFa]) GT = np.dot(ATF,IS1) GR = np.dot(ARF,IS1) HT = np.dot(ATF,IS2) HR = np.dot(ARF,IS2) else: IS1e = np.array([IinFo+DiC*QinFoe,DiC*IinFo+QinFoe,ZiC*(CosC*UinFoe+SinC*VinFo),-ZiC*(SinC*UinFoe-CosC*VinFo)]) IS2e = np.array([IinFa+DiC*QinFae,DiC*IinFa+QinFae,ZiC*(CosC*UinFae+SinC*VinFa),-ZiC*(SinC*UinFae-CosC*VinFa)]) #IS1e = np.array([IinFo,0,0,VinFo]) #IS2e = np.array([0,QinFae,UinFae,VinFa]) GT = np.dot(ATFe,IS1e) GR = np.dot(ARFe,IS1e) HT = np.dot(ATFe,IS2e) HR = np.dot(ARFe,IS2e) else: print('Calibrator not implemented yet') sys.exit() elif LocC == 2: # C behind emitter optics Eq.57 ------------------------------------------------------- #print("Calibrator location not implemented yet") S2e = np.sin(np.deg2rad(2*RotC)) C2e = np.cos(np.deg2rad(2*RotC)) # AS with C before the receiver optics (see document rotated_diattenuator_X22x5deg.odt) AF1 = np.array([1,C2g*DiO,S2g*DiO,0]) AF2 = np.array([C2g*DiO,1-S2g**2*WiO,S2g*C2g*WiO,-S2g*ZiO*SinO]) AF3 = np.array([S2g*DiO, S2g*C2g*WiO, 1-C2g**2*WiO, C2g*ZiO*SinO]) AF4 = np.array([0, S2g*SinO, -C2g*SinO, CosO]) ATF = (ATP1*AF1+ATP2*AF2+ATP3*AF3+ATP4*AF4) ARF = (ARP1*AF1+ARP2*AF2+ARP3*AF3+ARP4*AF4) ATF1 = ATF[0] ATF2 = ATF[1] ATF3 = ATF[2] ATF4 = ATF[3] ARF1 = ARF[0] ARF2 = ARF[1] ARF3 = ARF[2] ARF4 = ARF[3] # AS with C behind the emitter # terms without aCal ATE1o, ARE1o = ATF1, ARF1 ATE2o, ARE2o = 0., 0. ATE3o, ARE3o = 0., 0. ATE4o, ARE4o = ATF4, ARF4 # terms with aCal ATE1a, ARE1a = 0. , 0. ATE2a, ARE2a = ATF2, ARF2 ATE3a, ARE3a = -ATF3, -ARF3 ATE4a, ARE4a = -2*ATF4, -2*ARF4 # rotated AinEa by epsilon ATE2ae = C2e*ATF2 + S2e*ATF3 ATE3ae = -S2e*ATF2 - C2e*ATF3 ARE2ae = C2e*ARF2 + S2e*ARF3 ARE3ae = -S2e*ARF2 - C2e*ARF3 ATE1 = ATE1o ATE2e = aCal*ATE2ae ATE3e = aCal*ATE3ae ATE4 = (1-2*aCal)*ATF4 ARE1 = ARE1o ARE2e = aCal*ARE2ae ARE3e = aCal*ARE3ae ARE4 = (1-2*aCal)*ARF4 # rotated IinE QinEe = C2e*QinE + S2e*UinE UinEe = C2e*UinE - S2e*QinE # Calibration signals and Calibration correction K from measurements with LDRCal / aCal if (TypeC == 2) or (TypeC == 1): # +++++++++ rotator calibration Eq. C.4 AT = ATE1o*IinE + (ATE4o+aCal*ATE4a)*h*VinE BT = aCal * (ATE3ae*QinEe - ATE2ae*h*UinEe) AR = ARE1o*IinE + (ARE4o+aCal*ARE4a)*h*VinE BR = aCal * (ARE3ae*QinEe - ARE2ae*h*UinEe) # Correction paremeters for normal measurements; they are independent of LDR if (not RotationErrorEpsilonForNormalMeasurements): # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) GT = ATE1o*IinE + ATE4o*h*VinE GR = ARE1o*IinE + ARE4o*h*VinE HT = ATE2a*QinE + ATE3a*h*UinEe + ATE4a*h*VinE HR = ARE2a*QinE + ARE3a*h*UinEe + ARE4a*h*VinE else: GT = ATE1o*IinE + ATE4o*h*VinE GR = ARE1o*IinE + ARE4o*h*VinE HT = ATE2ae*QinE + ATE3ae*h*UinEe + ATE4a*h*VinE HR = ARE2ae*QinE + ARE3ae*h*UinEe + ARE4a*h*VinE elif (TypeC == 3) or (TypeC == 4): # +++++++++ linear polariser calibration Eq. C.5 # p = +45°, m = -45° AT = ATE1*IinE + ZiC*CosC*(ATE2e*QinEe + ATE4*VinE) + ATE3e*UinEe BT = DiC*(ATE1*UinEe + ATE3e*IinE) + ZiC*SinC*(ATE4*QinEe - ATE2e*VinE) AR = ARE1*IinE + ZiC*CosC*(ARE2e*QinEe + ARE4*VinE) + ARE3e*UinEe BR = DiC*(ARE1*UinEe + ARE3e*IinE) + ZiC*SinC*(ARE4*QinEe - ARE2e*VinE) # Correction paremeters for normal measurements; they are independent of LDR if (not RotationErrorEpsilonForNormalMeasurements): # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) GT = ATE1o*IinE + ATE4o*VinE GR = ARE1o*IinE + ARE4o*VinE HT = ATE2a*QinE + ATE3a*UinE + ATE4a*VinE HR = ARE2a*QinE + ARE3a*UinE + ARE4a*VinE else: D = IinE + DiC*QinEe A = DiC*IinE + QinEe B = ZiC*(CosC*UinEe + SinC*VinE) C = -ZiC*(SinC*UinEe - CosC*VinE) GT = ATE1o*D + ATE4o*C GR = ARE1o*D + ARE4o*C HT = ATE2a*A + ATE3a*B + ATE4a*C HR = ARE2a*A + ARE3a*B + ARE4a*C elif (TypeC == 6): # real HWP calibration +-22.5° rotated_diattenuator_X22x5deg.odt # p = +22.5°, m = -22.5° IE1e = np.array([IinE+sqr05*DiC*QinEe, sqr05*DiC*IinE+(1-0.5*WiC)*QinEe, (1-0.5*WiC)*UinEe+sqr05*ZiC*SinC*VinE, -sqr05*ZiC*SinC*UinEe+ZiC*CosC*VinE]) IE2e = np.array([sqr05*DiC*UinEe, 0.5*WiC*UinEe-sqr05*ZiC*SinC*VinE, sqr05*DiC*IinE+0.5*WiC*QinEe, sqr05*ZiC*SinC*QinEe]) ATEe = np.array([ATE1,ATE2e,ATE3e,ATE4]) AREe = np.array([ARE1,ARE2e,ARE3e,ARE4]) AT = np.dot(ATEe,IE1e) AR = np.dot(AREe,IE1e) BT = np.dot(ATEe,IE2e) BR = np.dot(AREe,IE2e) # Correction paremeters for normal measurements; they are independent of LDR if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out GT = ATE1o*IinE + ATE4o*VinE GR = ARE1o*IinE + ARE4o*VinE HT = ATE2a*QinE + ATE3a*UinE + ATE4a*VinE HR = ARE2a*QinE + ARE3a*UinE + ARE4a*VinE else: D = IinE + DiC*QinEe A = DiC*IinE + QinEe B = ZiC*(CosC*UinEe + SinC*VinE) C = -ZiC*(SinC*UinEe - CosC*VinE) GT = ATE1o*D + ATE4o*C GR = ARE1o*D + ARE4o*C HT = ATE2a*A + ATE3a*B + ATE4a*C HR = ARE2a*A + ARE3a*B + ARE4a*C else: print('Calibrator not implemented yet') sys.exit() else: print("Calibrator location not implemented yet") sys.exit() # Determination of the correction K of the calibration factor IoutTp = TaT*TiT*TiO*TiE*(AT + BT) IoutTm = TaT*TiT*TiO*TiE*(AT - BT) IoutRp = TaR*TiR*TiO*TiE*(AR + BR) IoutRm = TaR*TiR*TiO*TiE*(AR - BR) # --- Results and Corrections; electronic etaR and etaT are assumed to be 1 Etapx = IoutRp/IoutTp Etamx = IoutRm/IoutTm Etax = (Etapx*Etamx)**0.5 Eta = (TaR*TiR)/(TaT*TiT) # Eta = Eta*/K Eq. 84 K = Etax / Eta # For comparison with Volkers Libreoffice Müller Matrix spreadsheet #Eta_test_p = (IoutRp/IoutTp) #Eta_test_m = (IoutRm/IoutTm) #Eta_test = (Eta_test_p*Eta_test_m)**0.5 # ----- Forward simulated signals and LDRsim with atrue; from input file It = TaT*TiT*TiO*TiE*(GT+atrue*HT) Ir = TaR*TiR*TiO*TiE*(GR+atrue*HR) # LDRsim = 1/Eta*Ir/It # simulated LDR* with Y from input file LDRsim = Ir/It # simulated uncorrected LDR with Y from input file # Corrected LDRsimCorr from forward simulated LDRsim (atrue) # LDRsimCorr = (1./Eta*LDRsim*(GT+HT)-(GR+HR))/((GR-HR)-1./Eta*LDRsim*(GT-HT)) if Y == -1.: LDRsimx = 1./LDRsim else: LDRsimx = LDRsim # The following is correct without doubt #LDRCorr = (LDRsim*K/Etax*(GT+HT)-(GR+HR))/((GR-HR)-LDRsim*K/Etax*(GT-HT)) # The following is a test whether the equations for calibration Etax and normal signal (GHK, LDRsim) are consistent LDRCorr = (LDRsim/Eta*(GT+HT)-(GR+HR))/((GR-HR)-LDRsim*K/Etax*(GT-HT)) TTa = TiT*TaT #*ATP1 TRa = TiR*TaR #*ARP1 F11sim = 1/(TiO*TiE)*((HR*Etax/K*It/TTa-HT*Ir/TRa)/(HR*GT-HT*GR)) # IL = 1, Etat = Etar = 1 return (GT, HT, GR, HR, K, Eta, LDRsimx, LDRCorr, DTa, DRa, TTa, TRa, F11sim) # ******************************************************************************************************************************* # --- CALC truth GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0) # -------------------------------------------------------- with open('output_' + LID + '.dat', 'w') as f: with redirect_stdout(f): print("From ", dname) print("Running ", fname) print("Reading input file ", InputFile) #, " for Lidar system :", EID, ", ", LID) print("for Lidar system: ", EID, ", ", LID) # --- Print iput information********************************* print(" --- Input parameters: value ±error / ±steps ----------------------") print("{0:8} {1:8} {2:8.5f}; {3:8} {4:7.4f}±{5:7.4f}/{6:2d}".format("Laser: ", "DOLP = ", bL, " rotation alpha = ", RotL0, dRotL, nRotL)) print(" Diatt., Tunpol, Retard., Rotation (deg)") print("{0:12} {1:7.4f}±{2:7.4f}/{8:2d}, {3:7.4f}, {4:3.0f}±{5:3.0f}/{9:2d}, {6:7.4f}±{7:7.4f}/{10:2d}".format("Emitter ", DiE0, dDiE, TiE, RetE0, dRetE, RotE0, dRotE, nDiE, nRetE, nRotE)) print("{0:12} {1:7.4f}±{2:7.4f}/{8:2d}, {3:7.4f}, {4:3.0f}±{5:3.0f}/{9:2d}, {6:7.4f}±{7:7.4f}/{10:2d}".format("Receiver ", DiO0, dDiO, TiO, RetO0, dRetO, RotO0, dRotO, nDiO, nRetO, nRotO)) print("{0:12} {1:7.4f}±{2:7.4f}/{8:2d}, {3:7.4f}, {4:3.0f}±{5:3.0f}/{9:2d}, {6:7.4f}±{7:7.4f}/{10:2d}".format("Calibrator ", DiC0, dDiC, TiC, RetC0, dRetC, RotC0, dRotC, nDiC, nRetC, nRotC)) print("{0:12}".format(" --- Pol.-filter ---")) print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format("ERT, ERR :", ERaT0, dERaT, nERaT, ERaR0, dERaR, nERaR)) print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format("RotaT , RotaR :", RotaT0, dRotaT, nRotaT, RotaR0,dRotaR,nRotaR)) print("{0:12}".format(" --- PBS ---")) print("{0:12}{1:7.4f}±{2:7.4f}/{9:2d}, {3:7.4f}±{4:7.4f}/{10:2d}, {5:7.4f}±{6:7.4f}/{11:2d},{7:7.4f}±{8:7.4f}/{12:2d}".format("TP,TS,RP,RS :", TP0, dTP, TS0, dTS, RP0, dRP, RS0, dRS, nTP, nTS, nRP, nRS)) print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}, {5:1.0f}".format("DT,TT,DR,TR,Y :", DiT, TiT, DiR, TiR, Y)) print("{0:12}".format(" --- Combined PBS + Pol.-filter ---")) print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}".format("DTa,TTa,DRa,TRa: ", DTa0, TTa0, DRa0, TRa0)) print() print("Rotation Error Epsilon For Normal Measurements = ", RotationErrorEpsilonForNormalMeasurements) #print ('LocC = ', LocC, Loc[LocC], '; TypeC = ',TypeC, Type[TypeC]) print(Type[TypeC], Loc[LocC], "; Parallel signal detected in", dY[int(Y+1)]) # end of print actual system parameters # ****************************************************************************** #print() #print(" --- LDRCal during calibration | simulated and corrected LDRs -------------") #print("{0:8} |{1:8}->{2:8},{3:9}->{4:9} |{5:8}->{6:8}".format(" LDRCal"," LDRtrue", " LDRsim"," LDRtrue2", " LDRsim2", " LDRmeas", " LDRcorr")) #print("{0:8.5f} |{1:8.5f}->{2:8.5f},{3:9.5f}->{4:9.5f} |{5:8.5f}->{6:8.5f}".format(LDRCal, LDRtrue, LDRsim, LDRtrue2, LDRsim2, LDRmeas, LDRCorr)) #print("{0:8} |{1:8}->{2:8}->{3:8}".format(" LDRCal"," LDRtrue", " LDRsimx", " LDRcorr")) #print("{0:6.3f}±{1:5.3f}/{2:2d}|{3:8.5f}->{4:8.5f}->{5:8.5f}".format(LDRCal0, dLDRCal, nLDRCal, LDRtrue, LDRsimx, LDRCorr)) #print("{0:8} |{1:8}->{2:8}->{3:8}".format(" LDRCal"," LDRtrue", " LDRsimx", " LDRcorr")) #print(" --- LDRCal during calibration") print("{0:26}: {1:6.3f}±{2:5.3f}/{3:2d}".format("LDRCal during calibration", LDRCal0, dLDRCal, nLDRCal)) #print("{0:8}={1:8.5f};{2:8}={3:8.5f}".format(" IinP",IinP," F11sim",F11sim)) print() K0List = np.zeros(3) LDRsimxList = np.zeros(3) LDRCalList = 0.004, 0.2, 0.45 for i,LDRCal in enumerate(LDRCalList): GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal) K0List[i] = K0 LDRsimxList[i] = LDRsimx print("{0:8},{1:8},{2:8},{3:8},{4:9},{5:9},{6:9}".format(" GR", " GT", " HR", " HT", " K(0.004)", " K(0.2)", " K(0.45)")) print("{0:8.5f},{1:8.5f},{2:8.5f},{3:8.5f},{4:9.5f},{5:9.5f},{6:9.5f}".format(GR0, GT0, HR0, HT0, K0List[0], K0List[1], K0List[2])) print('========================================================================') print("{0:9},{1:9},{2:9}".format(" LDRtrue", " LDRsimx", " LDRCorr")) LDRtrueList = 0.004, 0.02, 0.2, 0.45 for i,LDRtrue in enumerate(LDRtrueList): GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0) print("{0:9.5f},{1:9.5f},{2:9.5f}".format(LDRtrue, LDRsimx, LDRCorr)) file = open('output_' + LID + '.dat', 'r') print (file.read()) file.close() ''' if(PrintToOutputFile): f = open('output_ver7.dat', 'w') old_target = sys.stdout sys.stdout = f print("something") if(PrintToOutputFile): sys.stdout.flush() f.close sys.stdout = old_target ''' # --- CALC again truth with LDRCal0 to reset all 0-values GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0) # --- Start Errors calculation iN = -1 N = ((nRotL*2+1)* (nRotE*2+1)*(nRetE*2+1)*(nDiE*2+1)* (nRotO*2+1)*(nRetO*2+1)*(nDiO*2+1)* (nRotC*2+1)*(nRetC*2+1)*(nDiC*2+1)* (nTP*2+1)*(nTS*2+1)*(nRP*2+1)*(nRS*2+1)*(nERaT*2+1)*(nERaR*2+1)* (nRotaT*2+1)*(nRotaR*2+1)*(nRetT*2+1)*(nRetR*2+1)*(nLDRCal*2+1)) print("N = ",N ," ", end="") if N > 1e6: if user_yes_no_query('Warning: processing ' + str(N) + ' samples will take very long. Do you want to proceed?') == 0: sys.exit() if N > 5e6: if user_yes_no_query('Warning: the memory required for ' + str(N) + ' samples might be ' + '{0:5.1f}'.format(N/4e6) + ' GB. Do you anyway want to proceed?') == 0: sys.exit() #if user_yes_no_query('Warning: processing' + str(N) + ' samples will take very long. Do you want to proceed?') == 0: sys.exit() # --- Arrays for plotting ------ LDRmin = np.zeros(5) LDRmax = np.zeros(5) F11min = np.zeros(5) F11max = np.zeros(5) LDRrange = np.zeros(5) LDRrange = 0.004, 0.02, 0.1, 0.3, 0.45 #aLDRsimx = np.zeros(N) #aLDRsimx2 = np.zeros(N) #aLDRcorr = np.zeros(N) #aLDRcorr2 = np.zeros(N) aERaT = np.zeros(N) aERaR = np.zeros(N) aRotaT = np.zeros(N) aRotaR = np.zeros(N) aRetT = np.zeros(N) aRetR = np.zeros(N) aTP = np.zeros(N) aTS = np.zeros(N) aRP = np.zeros(N) aRS = np.zeros(N) aDiE = np.zeros(N) aDiO = np.zeros(N) aDiC = np.zeros(N) aRotC = np.zeros(N) aRetC = np.zeros(N) aRotL = np.zeros(N) aRetE = np.zeros(N) aRotE = np.zeros(N) aRetO = np.zeros(N) aRotO = np.zeros(N) aLDRCal = np.zeros(N) aA = np.zeros((5,N)) aX = np.zeros((5,N)) aF11corr = np.zeros((5,N)) atime = clock() dtime = clock() # --- Calc Error signals #GT, HT, GR, HR, K, Eta, LDRsim = Calc(RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC, TP, TS) # ---- Do the calculations of bra-ket vectors h = -1. if TypeC == 2 else 1 # from input file: measured LDRm and true LDRtrue, LDRtrue2 => ameas = (1.-LDRmeas)/(1+LDRmeas) atrue = (1.-LDRtrue)/(1+LDRtrue) atrue2 = (1.-LDRtrue2)/(1+LDRtrue2) for iLDRCal in range(-nLDRCal,nLDRCal+1): # from input file: assumed LDRCal for calibration measurements LDRCal = LDRCal0 if nLDRCal > 0: LDRCal = LDRCal0 + iLDRCal*dLDRCal/nLDRCal GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal) aCal = (1.-LDRCal)/(1+LDRCal) for iRotL, iRotE, iRetE, iDiE \ in [(iRotL,iRotE,iRetE,iDiE) for iRotL in range(-nRotL,nRotL+1) for iRotE in range(-nRotE,nRotE+1) for iRetE in range(-nRetE,nRetE+1) for iDiE in range(-nDiE,nDiE+1)]: if nRotL > 0: RotL = RotL0 + iRotL*dRotL/nRotL if nRotE > 0: RotE = RotE0 + iRotE*dRotE/nRotE if nRetE > 0: RetE = RetE0 + iRetE*dRetE/nRetE if nDiE > 0: DiE = DiE0 + iDiE*dDiE/nDiE # angles of emitter and laser and calibrator and receiver optics # RotL = alpha, RotE = beta, RotO = gamma, RotC = epsilon S2a = np.sin(2*np.deg2rad(RotL)) C2a = np.cos(2*np.deg2rad(RotL)) S2b = np.sin(2*np.deg2rad(RotE)) C2b = np.cos(2*np.deg2rad(RotE)) S2ab = np.sin(np.deg2rad(2*RotL-2*RotE)) C2ab = np.cos(np.deg2rad(2*RotL-2*RotE)) # Laser with Degree of linear polarization DOLP = bL IinL = 1. QinL = bL UinL = 0. VinL = (1. - bL**2)**0.5 # Stokes Input Vector rotation Eq. E.4 A = C2a*QinL - S2a*UinL B = S2a*QinL + C2a*UinL # Stokes Input Vector rotation Eq. E.9 C = C2ab*QinL - S2ab*UinL D = S2ab*QinL + C2ab*UinL # emitter optics CosE = np.cos(np.deg2rad(RetE)) SinE = np.sin(np.deg2rad(RetE)) ZiE = (1. - DiE**2)**0.5 WiE = (1. - ZiE*CosE) # Stokes Input Vector after emitter optics equivalent to Eq. E.9 with already rotated input vector from Eq. E.4 # b = beta IinE = (IinL + DiE*C) QinE = (C2b*DiE*IinL + A + S2b*(WiE*D - ZiE*SinE*VinL)) UinE = (S2b*DiE*IinL + B - C2b*(WiE*D - ZiE*SinE*VinL)) VinE = (-ZiE*SinE*D + ZiE*CosE*VinL) #------------------------- # F11 assuemd to be = 1 => measured: F11m = IinP / IinE with atrue #F11sim = (IinE + DiO*atrue*(C2g*QinE - S2g*UinE))/IinE #------------------------- for iRotO, iRetO, iDiO, iRotC, iRetC, iDiC, iTP, iTS, iRP, iRS, iERaT, iRotaT, iRetT, iERaR, iRotaR, iRetR \ in [(iRotO,iRetO,iDiO,iRotC,iRetC,iDiC,iTP,iTS,iRP,iRS,iERaT,iRotaT,iRetT,iERaR,iRotaR,iRetR ) for iRotO in range(-nRotO,nRotO+1) for iRetO in range(-nRetO,nRetO+1) for iDiO in range(-nDiO,nDiO+1) for iRotC in range(-nRotC,nRotC+1) for iRetC in range(-nRetC,nRetC+1) for iDiC in range(-nDiC,nDiC+1) for iTP in range(-nTP,nTP+1) for iTS in range(-nTS,nTS+1) for iRP in range(-nRP,nRP+1) for iRS in range(-nRS,nRS+1) for iERaT in range(-nERaT,nERaT+1) for iRotaT in range(-nRotaT,nRotaT+1) for iRetT in range(-nRetT,nRetT+1) for iERaR in range(-nERaR,nERaR+1) for iRotaR in range(-nRotaR,nRotaR+1) for iRetR in range(-nRetR,nRetR+1)]: iN = iN + 1 if (iN == 10001): ctime = clock() print(" estimated time ", "{0:4.2f}".format(N/10000 * (ctime-atime)), "sec ") #, end="") print("\r elapsed time ", "{0:5.0f}".format((ctime-atime)), "sec ", end="\r") ctime = clock() if ((ctime - dtime) > 10): print("\r elapsed time ", "{0:5.0f}".format((ctime-atime)), "sec ", end="\r") dtime = ctime if nRotO > 0: RotO = RotO0 + iRotO*dRotO/nRotO if nRetO > 0: RetO = RetO0 + iRetO*dRetO/nRetO if nDiO > 0: DiO = DiO0 + iDiO*dDiO/nDiO if nRotC > 0: RotC = RotC0 + iRotC*dRotC/nRotC if nRetC > 0: RetC = RetC0 + iRetC*dRetC/nRetC if nDiC > 0: DiC = DiC0 + iDiC*dDiC/nDiC if nTP > 0: TP = TP0 + iTP*dTP/nTP if nTS > 0: TS = TS0 + iTS*dTS/nTS if nRP > 0: RP = RP0 + iRP*dRP/nRP if nRS > 0: RS = RS0 + iRS*dRS/nRS if nERaT > 0: ERaT = ERaT0 + iERaT*dERaT/nERaT if nRotaT > 0:RotaT= RotaT0+ iRotaT*dRotaT/nRotaT if nRetT > 0: RetT = RetT0 + iRetT*dRetT/nRetT if nERaR > 0: ERaR = ERaR0 + iERaR*dERaR/nERaR if nRotaR > 0:RotaR= RotaR0+ iRotaR*dRotaR/nRotaR if nRetR > 0: RetR = RetR0 + iRetR*dRetR/nRetR #print("{0:5.2f}, {1:5.2f}, {2:5.2f}, {3:10d}".format(RotL, RotE, RotO, iN)) # receiver optics CosO = np.cos(np.deg2rad(RetO)) SinO = np.sin(np.deg2rad(RetO)) ZiO = (1. - DiO**2)**0.5 WiO = (1. - ZiO*CosO) S2g = np.sin(np.deg2rad(2*RotO)) C2g = np.cos(np.deg2rad(2*RotO)) # calibrator CosC = np.cos(np.deg2rad(RetC)) SinC = np.sin(np.deg2rad(RetC)) ZiC = (1. - DiC**2)**0.5 WiC = (1. - ZiC*CosC) #For POLLY_XT # analyser #RS = 1 - TS #RP = 1 - TP TiT = 0.5 * (TP + TS) DiT = (TP-TS)/(TP+TS) ZiT = (1. - DiT**2)**0.5 TiR = 0.5 * (RP + RS) DiR = (RP-RS)/(RP+RS) ZiR = (1. - DiR**2)**0.5 CosT = np.cos(np.deg2rad(RetT)) SinT = np.sin(np.deg2rad(RetT)) CosR = np.cos(np.deg2rad(RetR)) SinR = np.sin(np.deg2rad(RetR)) DaT = (1-ERaT)/(1+ERaT) DaR = (1-ERaR)/(1+ERaR) TaT = 0.5*(1+ERaT) TaR = 0.5*(1+ERaR) S2aT = np.sin(np.deg2rad(h*2*RotaT)) C2aT = np.cos(np.deg2rad(2*RotaT)) S2aR = np.sin(np.deg2rad(h*2*RotaR)) C2aR = np.cos(np.deg2rad(2*RotaR)) # Aanalyzer As before the PBS Eq. D.5 ATP1 = (1+C2aT*DaT*DiT) ATP2 = Y*(DiT+C2aT*DaT) ATP3 = Y*S2aT*DaT*ZiT*CosT ATP4 = S2aT*DaT*ZiT*SinT ATP = np.array([ATP1,ATP2,ATP3,ATP4]) ARP1 = (1+C2aR*DaR*DiR) ARP2 = Y*(DiR+C2aR*DaR) ARP3 = Y*S2aR*DaR*ZiR*CosR ARP4 = S2aR*DaR*ZiR*SinR ARP = np.array([ARP1,ARP2,ARP3,ARP4]) TTa = TiT*TaT #*ATP1 TRa = TiR*TaR #*ARP1 # ---- Calculate signals and correction parameters for diffeent locations and calibrators if LocC == 4: # Calibrator before the PBS #print("Calibrator location not implemented yet") #S2ge = np.sin(np.deg2rad(2*RotO + h*2*RotC)) #C2ge = np.cos(np.deg2rad(2*RotO + h*2*RotC)) S2e = np.sin(np.deg2rad(h*2*RotC)) C2e = np.cos(np.deg2rad(2*RotC)) # rotated AinP by epsilon Eq. C.3 ATP2e = C2e*ATP2 + S2e*ATP3 ATP3e = C2e*ATP3 - S2e*ATP2 ARP2e = C2e*ARP2 + S2e*ARP3 ARP3e = C2e*ARP3 - S2e*ARP2 ATPe = np.array([ATP1,ATP2e,ATP3e,ATP4]) ARPe = np.array([ARP1,ARP2e,ARP3e,ARP4]) # Stokes Input Vector before the polarising beam splitter Eq. E.31 A = C2g*QinE - S2g*UinE B = S2g*QinE + C2g*UinE #C = (WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE) Co = ZiO*SinO*VinE Ca = (WiO*B - 2*ZiO*SinO*VinE) #C = Co + aCal*Ca #IinP = (IinE + DiO*aCal*A) #QinP = (C2g*DiO*IinE + aCal*QinE - S2g*C) #UinP = (S2g*DiO*IinE - aCal*UinE + C2g*C) #VinP = (ZiO*SinO*aCal*B + ZiO*CosO*(1-2*aCal)*VinE) IinPo = IinE QinPo = (C2g*DiO*IinE - S2g*Co) UinPo = (S2g*DiO*IinE + C2g*Co) VinPo = ZiO*CosO*VinE IinPa = DiO*A QinPa = QinE - S2g*Ca UinPa = -UinE + C2g*Ca VinPa = ZiO*(SinO*B - 2*CosO*VinE) IinP = IinPo + aCal*IinPa QinP = QinPo + aCal*QinPa UinP = UinPo + aCal*UinPa VinP = VinPo + aCal*VinPa # Stokes Input Vector before the polarising beam splitter rotated by epsilon Eq. C.3 #QinPe = C2e*QinP + S2e*UinP #UinPe = C2e*UinP - S2e*QinP QinPoe = C2e*QinPo + S2e*UinPo UinPoe = C2e*UinPo - S2e*QinPo QinPae = C2e*QinPa + S2e*UinPa UinPae = C2e*UinPa - S2e*QinPa QinPe = C2e*QinP + S2e*UinP UinPe = C2e*UinP - S2e*QinP # Calibration signals and Calibration correction K from measurements with LDRCal / aCal if (TypeC == 2) or (TypeC == 1): # rotator calibration Eq. C.4 # parameters for calibration with aCal AT = ATP1*IinP + h*ATP4*VinP BT = ATP3e*QinP - h*ATP2e*UinP AR = ARP1*IinP + h*ARP4*VinP BR = ARP3e*QinP - h*ARP2e*UinP # Correction paremeters for normal measurements; they are independent of LDR if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out IS1 = np.array([IinPo,QinPo,UinPo,VinPo]) IS2 = np.array([IinPa,QinPa,UinPa,VinPa]) GT = np.dot(ATP,IS1) GR = np.dot(ARP,IS1) HT = np.dot(ATP,IS2) HR = np.dot(ARP,IS2) else: IS1 = np.array([IinPo,QinPo,UinPo,VinPo]) IS2 = np.array([IinPa,QinPa,UinPa,VinPa]) GT = np.dot(ATPe,IS1) GR = np.dot(ARPe,IS1) HT = np.dot(ATPe,IS2) HR = np.dot(ARPe,IS2) elif (TypeC == 3) or (TypeC == 4): # linear polariser calibration Eq. C.5 # parameters for calibration with aCal AT = ATP1*IinP + ATP3e*UinPe + ZiC*CosC*(ATP2e*QinPe + ATP4*VinP) BT = DiC*(ATP1*UinPe + ATP3e*IinP) - ZiC*SinC*(ATP2e*VinP - ATP4*QinPe) AR = ARP1*IinP + ARP3e*UinPe + ZiC*CosC*(ARP2e*QinPe + ARP4*VinP) BR = DiC*(ARP1*UinPe + ARP3e*IinP) - ZiC*SinC*(ARP2e*VinP - ARP4*QinPe) # Correction paremeters for normal measurements; they are independent of LDR if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out IS1 = np.array([IinPo,QinPo,UinPo,VinPo]) IS2 = np.array([IinPa,QinPa,UinPa,VinPa]) GT = np.dot(ATP,IS1) GR = np.dot(ARP,IS1) HT = np.dot(ATP,IS2) HR = np.dot(ARP,IS2) else: IS1e = np.array([IinPo+DiC*QinPoe,DiC*IinPo+QinPoe,ZiC*(CosC*UinPoe+SinC*VinPo),-ZiC*(SinC*UinPoe-CosC*VinPo)]) IS2e = np.array([IinPa+DiC*QinPae,DiC*IinPa+QinPae,ZiC*(CosC*UinPae+SinC*VinPa),-ZiC*(SinC*UinPae-CosC*VinPa)]) GT = np.dot(ATPe,IS1e) GR = np.dot(ARPe,IS1e) HT = np.dot(ATPe,IS2e) HR = np.dot(ARPe,IS2e) elif (TypeC == 6): # diattenuator calibration +-22.5° rotated_diattenuator_X22x5deg.odt # parameters for calibration with aCal AT = ATP1*IinP + sqr05*DiC*(ATP1*QinPe + ATP2e*IinP) + (1-0.5*WiC)*(ATP2e*QinPe + ATP3e*UinPe) + ZiC*(sqr05*SinC*(ATP3e*VinP-ATP4*UinPe) + ATP4*CosC*VinP) BT = sqr05*DiC*(ATP1*UinPe + ATP3e*IinP) + 0.5*WiC*(ATP2e*UinPe + ATP3e*QinPe) - sqr05*ZiC*SinC*(ATP2e*VinP - ATP4*QinPe) AR = ARP1*IinP + sqr05*DiC*(ARP1*QinPe + ARP2e*IinP) + (1-0.5*WiC)*(ARP2e*QinPe + ARP3e*UinPe) + ZiC*(sqr05*SinC*(ARP3e*VinP-ARP4*UinPe) + ARP4*CosC*VinP) BR = sqr05*DiC*(ARP1*UinPe + ARP3e*IinP) + 0.5*WiC*(ARP2e*UinPe + ARP3e*QinPe) - sqr05*ZiC*SinC*(ARP2e*VinP - ARP4*QinPe) # Correction paremeters for normal measurements; they are independent of LDR if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out IS1 = np.array([IinPo,QinPo,UinPo,VinPo]) IS2 = np.array([IinPa,QinPa,UinPa,VinPa]) GT = np.dot(ATP,IS1) GR = np.dot(ARP,IS1) HT = np.dot(ATP,IS2) HR = np.dot(ARP,IS2) else: IS1e = np.array([IinPo+DiC*QinPoe,DiC*IinPo+QinPoe,ZiC*(CosC*UinPoe+SinC*VinPo),-ZiC*(SinC*UinPoe-CosC*VinPo)]) IS2e = np.array([IinPa+DiC*QinPae,DiC*IinPa+QinPae,ZiC*(CosC*UinPae+SinC*VinPa),-ZiC*(SinC*UinPae-CosC*VinPa)]) GT = np.dot(ATPe,IS1e) GR = np.dot(ARPe,IS1e) HT = np.dot(ATPe,IS2e) HR = np.dot(ARPe,IS2e) else: print("Calibrator not implemented yet") sys.exit() elif LocC == 3: # C before receiver optics Eq.57 #S2ge = np.sin(np.deg2rad(2*RotO - 2*RotC)) #C2ge = np.cos(np.deg2rad(2*RotO - 2*RotC)) S2e = np.sin(np.deg2rad(2*RotC)) C2e = np.cos(np.deg2rad(2*RotC)) # AS with C before the receiver optics (see document rotated_diattenuator_X22x5deg.odt) AF1 = np.array([1,C2g*DiO,S2g*DiO,0]) AF2 = np.array([C2g*DiO,1-S2g**2*WiO,S2g*C2g*WiO,-S2g*ZiO*SinO]) AF3 = np.array([S2g*DiO, S2g*C2g*WiO, 1-C2g**2*WiO, C2g*ZiO*SinO]) AF4 = np.array([0, S2g*SinO, -C2g*SinO, CosO]) ATF = (ATP1*AF1+ATP2*AF2+ATP3*AF3+ATP4*AF4) ARF = (ARP1*AF1+ARP2*AF2+ARP3*AF3+ARP4*AF4) ATF1 = ATF[0] ATF2 = ATF[1] ATF3 = ATF[2] ATF4 = ATF[3] ARF1 = ARF[0] ARF2 = ARF[1] ARF3 = ARF[2] ARF4 = ARF[3] # rotated AinF by epsilon ATF2e = C2e*ATF[1] + S2e*ATF[2] ATF3e = C2e*ATF[2] - S2e*ATF[1] ARF2e = C2e*ARF[1] + S2e*ARF[2] ARF3e = C2e*ARF[2] - S2e*ARF[1] ATFe = np.array([ATF1,ATF2e,ATF3e,ATF4]) ARFe = np.array([ARF1,ARF2e,ARF3e,ARF4]) QinEe = C2e*QinE + S2e*UinE UinEe = C2e*UinE - S2e*QinE # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) IinF = IinE QinF = aCal*QinE UinF = -aCal*UinE VinF = (1.-2.*aCal)*VinE IinFo = IinE QinFo = 0. UinFo = 0. VinFo = VinE IinFa = 0. QinFa = QinE UinFa = -UinE VinFa = -2.*VinE # Stokes Input Vector before receiver optics rotated by epsilon Eq. C.3 QinFe = C2e*QinF + S2e*UinF UinFe = C2e*UinF - S2e*QinF QinFoe = C2e*QinFo + S2e*UinFo UinFoe = C2e*UinFo - S2e*QinFo QinFae = C2e*QinFa + S2e*UinFa UinFae = C2e*UinFa - S2e*QinFa # Calibration signals and Calibration correction K from measurements with LDRCal / aCal if (TypeC == 2) or (TypeC == 1): # rotator calibration Eq. C.4 AT = ATF1*IinF + ATF4*h*VinF BT = ATF3e*QinF - ATF2e*h*UinF AR = ARF1*IinF + ARF4*h*VinF BR = ARF3e*QinF - ARF2e*h*UinF # Correction paremeters for normal measurements; they are independent of LDR if (not RotationErrorEpsilonForNormalMeasurements): GT = ATF1*IinE + ATF4*VinE GR = ARF1*IinE + ARF4*VinE HT = ATF2*QinE - ATF3*UinE - ATF4*2*VinE HR = ARF2*QinE - ARF3*UinE - ARF4*2*VinE else: GT = ATF1*IinE + ATF4*h*VinE GR = ARF1*IinE + ARF4*h*VinE HT = ATF2e*QinE - ATF3e*h*UinE - ATF4*h*2*VinE HR = ARF2e*QinE - ARF3e*h*UinE - ARF4*h*2*VinE elif (TypeC == 3) or (TypeC == 4): # linear polariser calibration Eq. C.5 # p = +45°, m = -45° IF1e = np.array([IinF, ZiC*CosC*QinFe, UinFe, ZiC*CosC*VinF]) IF2e = np.array([DiC*UinFe, -ZiC*SinC*VinF, DiC*IinF, ZiC*SinC*QinFe]) AT = np.dot(ATFe,IF1e) AR = np.dot(ARFe,IF1e) BT = np.dot(ATFe,IF2e) BR = np.dot(ARFe,IF2e) # Correction paremeters for normal measurements; they are independent of LDR --- the same as for TypeC = 6 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out IS1 = np.array([IinE,0,0,VinE]) IS2 = np.array([0,QinE,-UinE,-2*VinE]) GT = np.dot(ATF,IS1) GR = np.dot(ARF,IS1) HT = np.dot(ATF,IS2) HR = np.dot(ARF,IS2) else: IS1e = np.array([IinFo+DiC*QinFoe,DiC*IinFo+QinFoe,ZiC*(CosC*UinFoe+SinC*VinFo),-ZiC*(SinC*UinFoe-CosC*VinFo)]) IS2e = np.array([IinFa+DiC*QinFae,DiC*IinFa+QinFae,ZiC*(CosC*UinFae+SinC*VinFa),-ZiC*(SinC*UinFae-CosC*VinFa)]) GT = np.dot(ATFe,IS1e) GR = np.dot(ARFe,IS1e) HT = np.dot(ATFe,IS2e) HR = np.dot(ARFe,IS2e) elif (TypeC == 6): # diattenuator calibration +-22.5° rotated_diattenuator_X22x5deg.odt # p = +22.5°, m = -22.5° IF1e = np.array([IinF+sqr05*DiC*QinFe, sqr05*DiC*IinF+(1-0.5*WiC)*QinFe, (1-0.5*WiC)*UinFe+sqr05*ZiC*SinC*VinF, -sqr05*ZiC*SinC*UinFe+ZiC*CosC*VinF]) IF2e = np.array([sqr05*DiC*UinFe, 0.5*WiC*UinFe-sqr05*ZiC*SinC*VinF, sqr05*DiC*IinF+0.5*WiC*QinFe, sqr05*ZiC*SinC*QinFe]) AT = np.dot(ATFe,IF1e) AR = np.dot(ARFe,IF1e) BT = np.dot(ATFe,IF2e) BR = np.dot(ARFe,IF2e) # Correction paremeters for normal measurements; they are independent of LDR if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out #IS1 = np.array([IinE,0,0,VinE]) #IS2 = np.array([0,QinE,-UinE,-2*VinE]) IS1 = np.array([IinFo,0,0,VinFo]) IS2 = np.array([0,QinFa,UinFa,VinFa]) GT = np.dot(ATF,IS1) GR = np.dot(ARF,IS1) HT = np.dot(ATF,IS2) HR = np.dot(ARF,IS2) else: #IS1e = np.array([IinE,DiC*IinE,ZiC*SinC*VinE,ZiC*CosC*VinE]) #IS2e = np.array([DiC*QinEe,QinEe,-ZiC*(CosC*UinEe+2*SinC*VinE),ZiC*(SinC*UinEe-2*CosC*VinE)]) IS1e = np.array([IinFo+DiC*QinFoe,DiC*IinFo+QinFoe,ZiC*(CosC*UinFoe+SinC*VinFo),-ZiC*(SinC*UinFoe-CosC*VinFo)]) IS2e = np.array([IinFa+DiC*QinFae,DiC*IinFa+QinFae,ZiC*(CosC*UinFae+SinC*VinFa),-ZiC*(SinC*UinFae-CosC*VinFa)]) GT = np.dot(ATFe,IS1e) GR = np.dot(ARFe,IS1e) HT = np.dot(ATFe,IS2e) HR = np.dot(ARFe,IS2e) else: print('Calibrator not implemented yet') sys.exit() elif LocC == 2: # C behind emitter optics Eq.57 #print("Calibrator location not implemented yet") S2e = np.sin(np.deg2rad(2*RotC)) C2e = np.cos(np.deg2rad(2*RotC)) # AS with C before the receiver optics (see document rotated_diattenuator_X22x5deg.odt) AF1 = np.array([1,C2g*DiO,S2g*DiO,0]) AF2 = np.array([C2g*DiO,1-S2g**2*WiO,S2g*C2g*WiO,-S2g*ZiO*SinO]) AF3 = np.array([S2g*DiO, S2g*C2g*WiO, 1-C2g**2*WiO, C2g*ZiO*SinO]) AF4 = np.array([0, S2g*SinO, -C2g*SinO, CosO]) ATF = (ATP1*AF1+ATP2*AF2+ATP3*AF3+ATP4*AF4) ARF = (ARP1*AF1+ARP2*AF2+ARP3*AF3+ARP4*AF4) ATF1 = ATF[0] ATF2 = ATF[1] ATF3 = ATF[2] ATF4 = ATF[3] ARF1 = ARF[0] ARF2 = ARF[1] ARF3 = ARF[2] ARF4 = ARF[3] # AS with C behind the emitter -------------------------------------------- # terms without aCal ATE1o, ARE1o = ATF1, ARF1 ATE2o, ARE2o = 0., 0. ATE3o, ARE3o = 0., 0. ATE4o, ARE4o = ATF4, ARF4 # terms with aCal ATE1a, ARE1a = 0. , 0. ATE2a, ARE2a = ATF2, ARF2 ATE3a, ARE3a = -ATF3, -ARF3 ATE4a, ARE4a = -2*ATF4, -2*ARF4 # rotated AinEa by epsilon ATE2ae = C2e*ATF2 + S2e*ATF3 ATE3ae = -S2e*ATF2 - C2e*ATF3 ARE2ae = C2e*ARF2 + S2e*ARF3 ARE3ae = -S2e*ARF2 - C2e*ARF3 ATE1 = ATE1o ATE2e = aCal*ATE2ae ATE3e = aCal*ATE3ae ATE4 = (1-2*aCal)*ATF4 ARE1 = ARE1o ARE2e = aCal*ARE2ae ARE3e = aCal*ARE3ae ARE4 = (1-2*aCal)*ARF4 # rotated IinE QinEe = C2e*QinE + S2e*UinE UinEe = C2e*UinE - S2e*QinE # --- Calibration signals and Calibration correction K from measurements with LDRCal / aCal if (TypeC == 2) or (TypeC == 1): # +++++++++ rotator calibration Eq. C.4 AT = ATE1o*IinE + (ATE4o+aCal*ATE4a)*h*VinE BT = aCal * (ATE3ae*QinEe - ATE2ae*h*UinEe) AR = ARE1o*IinE + (ARE4o+aCal*ARE4a)*h*VinE BR = aCal * (ARE3ae*QinEe - ARE2ae*h*UinEe) # Correction paremeters for normal measurements; they are independent of LDR if (not RotationErrorEpsilonForNormalMeasurements): # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) GT = ATE1o*IinE + ATE4o*h*VinE GR = ARE1o*IinE + ARE4o*h*VinE HT = ATE2a*QinE + ATE3a*h*UinEe + ATE4a*h*VinE HR = ARE2a*QinE + ARE3a*h*UinEe + ARE4a*h*VinE else: GT = ATE1o*IinE + ATE4o*h*VinE GR = ARE1o*IinE + ARE4o*h*VinE HT = ATE2ae*QinE + ATE3ae*h*UinEe + ATE4a*h*VinE HR = ARE2ae*QinE + ARE3ae*h*UinEe + ARE4a*h*VinE elif (TypeC == 3) or (TypeC == 4): # +++++++++ linear polariser calibration Eq. C.5 # p = +45°, m = -45° AT = ATE1*IinE + ZiC*CosC*(ATE2e*QinEe + ATE4*VinE) + ATE3e*UinEe BT = DiC*(ATE1*UinEe + ATE3e*IinE) + ZiC*SinC*(ATE4*QinEe - ATE2e*VinE) AR = ARE1*IinE + ZiC*CosC*(ARE2e*QinEe + ARE4*VinE) + ARE3e*UinEe BR = DiC*(ARE1*UinEe + ARE3e*IinE) + ZiC*SinC*(ARE4*QinEe - ARE2e*VinE) # Correction paremeters for normal measurements; they are independent of LDR if (not RotationErrorEpsilonForNormalMeasurements): # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) GT = ATE1o*IinE + ATE4o*VinE GR = ARE1o*IinE + ARE4o*VinE HT = ATE2a*QinE + ATE3a*UinE + ATE4a*VinE HR = ARE2a*QinE + ARE3a*UinE + ARE4a*VinE else: D = IinE + DiC*QinEe A = DiC*IinE + QinEe B = ZiC*(CosC*UinEe + SinC*VinE) C = -ZiC*(SinC*UinEe - CosC*VinE) GT = ATE1o*D + ATE4o*C GR = ARE1o*D + ARE4o*C HT = ATE2a*A + ATE3a*B + ATE4a*C HR = ARE2a*A + ARE3a*B + ARE4a*C elif (TypeC == 6): # real HWP calibration +-22.5° rotated_diattenuator_X22x5deg.odt # p = +22.5°, m = -22.5° IE1e = np.array([IinE+sqr05*DiC*QinEe, sqr05*DiC*IinE+(1-0.5*WiC)*QinEe, (1-0.5*WiC)*UinEe+sqr05*ZiC*SinC*VinE, -sqr05*ZiC*SinC*UinEe+ZiC*CosC*VinE]) IE2e = np.array([sqr05*DiC*UinEe, 0.5*WiC*UinEe-sqr05*ZiC*SinC*VinE, sqr05*DiC*IinE+0.5*WiC*QinEe, sqr05*ZiC*SinC*QinEe]) ATEe = np.array([ATE1,ATE2e,ATE3e,ATE4]) AREe = np.array([ARE1,ARE2e,ARE3e,ARE4]) AT = np.dot(ATEe,IE1e) AR = np.dot(AREe,IE1e) BT = np.dot(ATEe,IE2e) BR = np.dot(AREe,IE2e) # Correction paremeters for normal measurements; they are independent of LDR if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out GT = ATE1o*IinE + ATE4o*VinE GR = ARE1o*IinE + ARE4o*VinE HT = ATE2a*QinE + ATE3a*UinE + ATE4a*VinE HR = ARE2a*QinE + ARE3a*UinE + ARE4a*VinE else: D = IinE + DiC*QinEe A = DiC*IinE + QinEe B = ZiC*(CosC*UinEe + SinC*VinE) C = -ZiC*(SinC*UinEe - CosC*VinE) GT = ATE1o*D + ATE4o*C GR = ARE1o*D + ARE4o*C HT = ATE2a*A + ATE3a*B + ATE4a*C HR = ARE2a*A + ARE3a*B + ARE4a*C else: print('Calibrator not implemented yet') sys.exit() # Calibration signals with aCal => Determination of the correction K of the real calibration factor IoutTp = TaT*TiT*TiO*TiE*(AT + BT) IoutTm = TaT*TiT*TiO*TiE*(AT - BT) IoutRp = TaR*TiR*TiO*TiE*(AR + BR) IoutRm = TaR*TiR*TiO*TiE*(AR - BR) # --- Results and Corrections; electronic etaR and etaT are assumed to be 1 #Eta = TiR/TiT # Eta = Eta*/K Eq. 84 Etapx = IoutRp/IoutTp Etamx = IoutRm/IoutTm Etax = (Etapx*Etamx)**0.5 K = Etax / Eta0 #print("{0:6.3f},{1:6.3f},{2:6.3f},{3:6.3f},{4:6.3f},{5:6.3f},{6:6.3f},{7:6.3f},{8:6.3f},{9:6.3f},{10:6.3f}".format(AT, BT, AR, BR, DiC, ZiC, RetO, TP, TS, Kp, Km)) #print("{0:6.3f},{1:6.3f},{2:6.3f},{3:6.3f}".format(DiC, ZiC, Kp, Km)) # For comparison with Volkers Libreoffice Müller Matrix spreadsheet #Eta_test_p = (IoutRp/IoutTp) #Eta_test_m = (IoutRm/IoutTm) #Eta_test = (Eta_test_p*Eta_test_m)**0.5 # ************************************************************************* iLDR = -1 for LDRTrue in LDRrange: iLDR = iLDR + 1 atrue = (1-LDRTrue)/(1+LDRTrue) # ----- Forward simulated signals and LDRsim with atrue; from input file It = TaT*TiT*TiO*TiE*(GT+atrue*HT) # TaT*TiT*TiC*TiO*IinL*(GT+atrue*HT) Ir = TaR*TiR*TiO*TiE*(GR+atrue*HR) # TaR*TiR*TiC*TiO*IinL*(GR+atrue*HR) # LDRsim = 1/Eta*Ir/It # simulated LDR* with Y from input file LDRsim = Ir/It # simulated uncorrected LDR with Y from input file ''' if Y == 1.: LDRsimx = LDRsim LDRsimx2 = LDRsim2 else: LDRsimx = 1./LDRsim LDRsimx2 = 1./LDRsim2 ''' # ----- Backward correction # Corrected LDRCorr from forward simulated LDRsim (atrue) with assumed true G0,H0,K0 LDRCorr = (LDRsim*K0/Etax*(GT0+HT0)-(GR0+HR0))/((GR0-HR0)-LDRsim*K0/Etax*(GT0-HT0)) # -- F11corr from It and Ir and calibration EtaX Text1 = "!!! EXPERIMENTAL !!! F11corr from It and Ir with calibration EtaX: x-axis: F11corr(LDRtrue) / F11corr(LDRtrue = 0.004) - 1" F11corr = 1/(TiO*TiE)*((HR0*Etax/K0*It/TTa-HT0*Ir/TRa)/(HR0*GT0-HT0*GR0)) # IL = 1 Eq.(64) #Text1 = "F11corr from It and Ir without corrections but with calibration EtaX: x-axis: F11corr(LDRtrue) devided by F11corr(LDRtrue = 0.004)" #F11corr = 0.5/(TiO*TiE)*(Etax*It/TTa+Ir/TRa) # IL = 1 Eq.(64) # -- It from It only with atrue without corrections - for BERTHA (and PollyXTs) #Text1 = " x-axis: IT(LDRtrue) / IT(LDRtrue = 0.004) - 1" #F11corr = It/(TaT*TiT*TiO*TiE) #/(TaT*TiT*TiO*TiE*(GT0+atrue*HT0)) # !!! see below line 1673ff aF11corr[iLDR,iN] = F11corr aA[iLDR,iN] = LDRCorr aX[0,iN] = GR aX[1,iN] = GT aX[2,iN] = HR aX[3,iN] = HT aX[4,iN] = K aLDRCal[iN] = iLDRCal aERaT[iN] = iERaT aERaR[iN] = iERaR aRotaT[iN] = iRotaT aRotaR[iN] = iRotaR aRetT[iN] = iRetT aRetR[iN] = iRetR aRotL[iN] = iRotL aRotE[iN] = iRotE aRetE[iN] = iRetE aRotO[iN] = iRotO aRetO[iN] = iRetO aRotC[iN] = iRotC aRetC[iN] = iRetC aDiO[iN] = iDiO aDiE[iN] = iDiE aDiC[iN] = iDiC aTP[iN] = iTP aTS[iN] = iTS aRP[iN] = iRP aRS[iN] = iRS # --- END loop btime = clock() print("\r done in ", "{0:5.0f}".format(btime-atime), "sec") #, end="\r") # --- Plot ----------------------------------------------------------------- if (sns_loaded): sns.set_style("whitegrid") sns.set_palette("bright", 6) ''' fig2 = plt.figure() plt.plot(aA[2,:],'b.') plt.plot(aA[3,:],'r.') plt.plot(aA[4,:],'g.') #plt.plot(aA[6,:],'c.') plt.show ''' # Plot LDR def PlotSubHist(aVar, aX, X0, daX, iaX, naX): fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) iLDR = -1 for LDRTrue in LDRrange: iLDR = iLDR + 1 LDRmin[iLDR] = np.min(aA[iLDR,:]) LDRmax[iLDR] = np.max(aA[iLDR,:]) Rmin = LDRmin[iLDR] * 0.995 # np.min(aA[iLDR,:]) * 0.995 Rmax = LDRmax[iLDR] * 1.005 # np.max(aA[iLDR,:]) * 1.005 #plt.subplot(5,2,iLDR+1) plt.subplot(1,5,iLDR+1) (n, bins, patches) = plt.hist(aA[iLDR,:], bins=100, log=False, range=[Rmin, Rmax], alpha=0.5, normed=False, color = '0.5', histtype='stepfilled') for iaX in range(-naX,naX+1): plt.hist(aA[iLDR,aX == iaX], range=[Rmin, Rmax], bins=100, log=False, alpha=0.3, normed=False, histtype='stepfilled', label = str(round(X0 + iaX*daX/naX,5))) if (iLDR == 2): plt.legend() plt.tick_params(axis='both', labelsize=9) plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2) #plt.title(LID + ' ' + aVar, fontsize=18) #plt.ylabel('frequency', fontsize=10) #plt.xlabel('LDRcorr', fontsize=10) #fig.tight_layout() fig.suptitle(LID + ' ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar, fontsize=14, y=1.05) #plt.show() #fig.savefig(LID + '_' + aVar + '.png', dpi=150, bbox_inches='tight', pad_inches=0) #plt.close return if (nRotL > 0): PlotSubHist("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) if (nRetE > 0): PlotSubHist("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) if (nRotE > 0): PlotSubHist("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) if (nDiE > 0): PlotSubHist("DiE", aDiE, DiE0, dDiE, iDiE, nDiE) if (nRetO > 0): PlotSubHist("RetO", aRetO, RetO0, dRetO, iRetO, nRetO) if (nRotO > 0): PlotSubHist("RotO", aRotO, RotO0, dRotO, iRotO, nRotO) if (nDiO > 0): PlotSubHist("DiO", aDiO, DiO0, dDiO, iDiO, nDiO) if (nDiC > 0): PlotSubHist("DiC", aDiC, DiC0, dDiC, iDiC, nDiC) if (nRotC > 0): PlotSubHist("RotC", aRotC, RotC0, dRotC, iRotC, nRotC) if (nRetC > 0): PlotSubHist("RetC", aRetC, RetC0, dRetC, iRetC, nRetC) if (nTP > 0): PlotSubHist("TP", aTP, TP0, dTP, iTP, nTP) if (nTS > 0): PlotSubHist("TS", aTS, TS0, dTS, iTS, nTS) if (nRP > 0): PlotSubHist("RP", aRP, RP0, dRP, iRP, nRP) if (nRS > 0): PlotSubHist("RS", aRS, RS0, dRS, iRS, nRS) if (nRetT > 0): PlotSubHist("RetT", aRetT, RetT0, dRetT, iRetT, nRetT) if (nRetR > 0): PlotSubHist("RetR", aRetR, RetR0, dRetR, iRetR, nRetR) if (nERaT > 0): PlotSubHist("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) if (nERaR > 0): PlotSubHist("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) if (nRotaT > 0): PlotSubHist("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) if (nRotaR > 0): PlotSubHist("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) if (nLDRCal > 0): PlotSubHist("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) plt.show() plt.close print() #print("IT(LDRtrue) devided by IT(LDRtrue = 0.004)") print(" ############################################################################## ") print(Text1) print() iLDR = 5 for LDRTrue in LDRrange: iLDR = iLDR - 1 aF11corr[iLDR,:] = aF11corr[iLDR,:] / aF11corr[0,:] - 1.0 # Plot F11 def PlotSubHistF11(aVar, aX, X0, daX, iaX, naX): fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) iLDR = -1 for LDRTrue in LDRrange: iLDR = iLDR + 1 ''' F11min[iLDR] = np.min(aF11corr[iLDR,:]) F11max[iLDR] = np.max(aF11corr[iLDR,:]) Rmin = F11min[iLDR] * 0.995 # np.min(aA[iLDR,:]) * 0.995 Rmax = F11max[iLDR] * 1.005 # np.max(aA[iLDR,:]) * 1.005 ''' #Rmin = 0.8 #Rmax = 1.2 #plt.subplot(5,2,iLDR+1) plt.subplot(1,5,iLDR+1) (n, bins, patches) = plt.hist(aF11corr[iLDR,:], bins=100, log=False, alpha=0.5, normed=False, color = '0.5', histtype='stepfilled') for iaX in range(-naX,naX+1): plt.hist(aF11corr[iLDR,aX == iaX], bins=100, log=False, alpha=0.3, normed=False, histtype='stepfilled', label = str(round(X0 + iaX*daX/naX,5))) if (iLDR == 2): plt.legend() plt.tick_params(axis='both', labelsize=9) #plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2) #plt.title(LID + ' ' + aVar, fontsize=18) #plt.ylabel('frequency', fontsize=10) #plt.xlabel('LDRcorr', fontsize=10) #fig.tight_layout() fig.suptitle(LID + ' ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar, fontsize=14, y=1.05) #plt.show() #fig.savefig(LID + '_' + aVar + '.png', dpi=150, bbox_inches='tight', pad_inches=0) #plt.close return if (nRotL > 0): PlotSubHistF11("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) if (nRetE > 0): PlotSubHistF11("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) if (nRotE > 0): PlotSubHistF11("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) if (nDiE > 0): PlotSubHistF11("DiE", aDiE, DiE0, dDiE, iDiE, nDiE) if (nRetO > 0): PlotSubHistF11("RetO", aRetO, RetO0, dRetO, iRetO, nRetO) if (nRotO > 0): PlotSubHistF11("RotO", aRotO, RotO0, dRotO, iRotO, nRotO) if (nDiO > 0): PlotSubHistF11("DiO", aDiO, DiO0, dDiO, iDiO, nDiO) if (nDiC > 0): PlotSubHistF11("DiC", aDiC, DiC0, dDiC, iDiC, nDiC) if (nRotC > 0): PlotSubHistF11("RotC", aRotC, RotC0, dRotC, iRotC, nRotC) if (nRetC > 0): PlotSubHistF11("RetC", aRetC, RetC0, dRetC, iRetC, nRetC) if (nTP > 0): PlotSubHistF11("TP", aTP, TP0, dTP, iTP, nTP) if (nTS > 0): PlotSubHistF11("TS", aTS, TS0, dTS, iTS, nTS) if (nRP > 0): PlotSubHistF11("RP", aRP, RP0, dRP, iRP, nRP) if (nRS > 0): PlotSubHistF11("RS", aRS, RS0, dRS, iRS, nRS) if (nRetT > 0): PlotSubHistF11("RetT", aRetT, RetT0, dRetT, iRetT, nRetT) if (nRetR > 0): PlotSubHistF11("RetR", aRetR, RetR0, dRetR, iRetR, nRetR) if (nERaT > 0): PlotSubHistF11("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) if (nERaR > 0): PlotSubHistF11("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) if (nRotaT > 0): PlotSubHistF11("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) if (nRotaR > 0): PlotSubHistF11("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) if (nLDRCal > 0): PlotSubHistF11("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) plt.show() plt.close ''' # only histogram #print("******************* " + aVar + " *******************") fig, ax = plt.subplots(nrows=5, ncols=2, sharex=True, sharey=True, figsize=(10, 10)) iLDR = -1 for LDRTrue in LDRrange: iLDR = iLDR + 1 LDRmin[iLDR] = np.min(aA[iLDR,:]) LDRmax[iLDR] = np.max(aA[iLDR,:]) Rmin = np.min(aA[iLDR,:]) * 0.999 Rmax = np.max(aA[iLDR,:]) * 1.001 plt.subplot(5,2,iLDR+1) (n, bins, patches) = plt.hist(aA[iLDR,:], range=[Rmin, Rmax], bins=200, log=False, alpha=0.2, normed=False, color = '0.5', histtype='stepfilled') plt.tick_params(axis='both', labelsize=9) plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2) plt.show() plt.close ''' # --- Plot LDRmin, LDRmax fig2 = plt.figure() plt.plot(LDRrange,LDRmax-LDRrange, linewidth=2.0, color='b') plt.plot(LDRrange,LDRmin-LDRrange, linewidth=2.0, color='g') plt.xlabel('LDRtrue', fontsize=18) plt.ylabel('LDRTrue-LDRmin, LDRTrue-LDRmax', fontsize=14) plt.title(LID + ' ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]), fontsize=18) #plt.ylimit(-0.07, 0.07) plt.show() plt.close # --- Save LDRmin, LDRmax to file # http://stackoverflow.com/questions/4675728/redirect-stdout-to-a-file-in-python with open('LDR_min_max_ver7_' + LID + '.dat', 'w') as f: with redirect_stdout(f): print(LID) print("LDRtrue, LDRmin, LDRmax") for i in range(len(LDRrange)): print("{0:7.4f},{1:7.4f},{2:7.4f}".format(LDRrange[i], LDRmin[i], LDRmax[i])) ''' # --- Plot K over LDRCal fig3 = plt.figure() plt.plot(LDRCal0+aLDRCal*dLDRCal/nLDRCal,aX[4,:], linewidth=2.0, color='b') plt.xlabel('LDRCal', fontsize=18) plt.ylabel('K', fontsize=14) plt.title(LID, fontsize=18) plt.show() plt.close ''' # Additional plot routines ======> ''' #****************************************************************************** # 1. Plot LDRcorrected - LDR(measured Icross/Iparallel) LDRa = np.arange(1.,100.)*0.005 LDRCorra = np.arange(1.,100.) if Y == - 1.: LDRa = 1./LDRa LDRCorra = (1./Eta*LDRa*(GT+HT)-(GR+HR))/((GR-HR)-1./Eta*LDRa*(GT-HT)) if Y == - 1.: LDRa = 1./LDRa # #fig = plt.figure() plt.plot(LDRa,LDRCorra-LDRa) plt.plot([0.,0.5],[0.,0.5]) plt.suptitle('LDRcorrected - LDR(measured Icross/Iparallel)', fontsize=16) plt.xlabel('LDR', fontsize=18) plt.ylabel('LDRCorr - LDR', fontsize=16) #plt.savefig('test.png') # ''' ''' #****************************************************************************** # 2. Plot LDRsim (simulated measurements without corrections = Icross/Iparallel) over LDRtrue LDRa = np.arange(1.,100.)*0.005 LDRsima = np.arange(1.,100.) atruea = (1.-LDRa)/(1+LDRa) Ita = TiT*TiO*IinL*(GT+atruea*HT) Ira = TiR*TiO*IinL*(GR+atruea*HR) LDRsima = Ira/Ita # simulated uncorrected LDR with Y from input file if Y == -1.: LDRsima = 1./LDRsima # #fig = plt.figure() plt.plot(LDRa,LDRsima) plt.plot([0.,0.5],[0.,0.5]) plt.suptitle('LDRsim (simulated measurements without corrections = Icross/Iparallel) over LDRtrue', fontsize=10) plt.xlabel('LDRtrue', fontsize=18) plt.ylabel('LDRsim', fontsize=16) #plt.savefig('test.png') # '''
