# HG changeset patch # User Volker Freudenthaler # Date 1625665350 -7200 # Node ID 7c14c9a085baa5d1c23ac3690994b9b3bf180676 # Parent a882c113284fd564241af14a8cfe23384b037d81 2nd order polynomial fits to - K(LDRCal) and - LDR-uncertainty max. and min depending on LDRtrue diff -r a882c113284f -r 7c14c9a085ba GHK_0.9.8h_Py3.7.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/GHK_0.9.8h_Py3.7.py Wed Jul 07 15:42:30 2021 +0200 @@ -0,0 +1,2963 @@ +# -*- coding: utf-8 -*- +""" +Copyright 2016...2021 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.7, seaborn 0.9.0 + +Code description: + +From measured lidar signals we cannot directly determine the desired backscatter coefficient (F11) and the linear depolarization ratio (LDR) +because of the cross talk between the channles and systematic errors of a lidar system. +http://www.atmos-meas-tech-discuss.net/amt-2015-338/amt-2015-338.pdf provides an analytical model for the description of these errors, +with which the measured signals can be corrected. +This code simulates the lidar measurements with "assumed true" model parameters from an input file, and calculates the correction parameters (G,H, and K). +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 +uncertainties. The uncertainties of the "assumed true" parameters can be described in the input file. Then this code calculates the lidar signals and the +gain ratio eta* with all possible combinations of "errors", which represents the distribution of "possibly real" signals, and "corrects" them with the "assumed true" +GHK parameters (GT0, GR0, HT0, HR0, and K0) to derive finally the distributions of "possibly real" linear depolarization ratios (LDRCorr), +which are plotted for five different input linear depolarization ratios (LDRtrue). The red bars in the plots represent the input values of LDRtrue. +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 +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 +"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 +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 +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). + +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. + +Ver. 0.9.7: includes the random error (signal noise) of the calibration and standard measurements +Changes: + Line 1687 Eta = (TaR * TiR) / (TaT * TiT) + Line 1691 K = Etax / Eta # K of the real system; but correction in Line 1721 with K0 / Etax + should work with nTCalT = nTCalR = 0 +Ver. 0.9.7b: + ToDo: include error due to TCalT und TCalR => determination of NCalT and NCalR etc. in error calculation line 1741ff + combined error loops iNI and INCal for signals +Ver. 0.9.7c: individual error loops for each of the six signals +Ver. 0.9.7c2: different calculation of the signal noise errors +Ver. 0.9.7c3: n.a.different calculation of the signal noise errors +Ver. 0.9.7c4: test to speed up the loops for error calculation by moving them just before the actual calculation: still some code errors +Ver. 0.9.8: + - correct calculation of Eta for cleaned anaylsers considering the combined transmission Eta = (TaT* TiT)(1 + cos2RotaT * DaT * DiT) and (TaR * TiR)(1 + cos2RotaR * DaR * DiR) according to the papers supplement Eqs. (S.10.10.1) ff + - calculation of the PLDR from LDR and BSR, BSR, and LDRm + - ND-filters can be added for the calibration measurements in the transmitted (TCalT) and the reflected path (TCalR) in order to include their uncertainties in the error calculation. + The ND-filters effect is not included in the K-parameter. To do that, you have to multiply the eta, which is determined from the "calibration measurements with ND-filter", with the attenuation of the ND-filter "by hand", in order to get the calibration factor for the "standard measurements without ND-filter". Furthermore, the additional error of the calibration factor eta due to the uncertainty of the attenuation of the ND-filter has to be considered "by hand". +Ver. 0.9.8b: change from "TTa = TiT * TaT" to "TTa = TiT * TaT * ATPT" etc. (compare ver 0.9.8 with 0.9.8b) removes + - the strong Tp dependence of the errors + - the factor 2 in the GH parameters + - see c:\technik\Optik\Polarizers\DepCal\ApplOpt\GH-parameters-190114.odt +Ver. 0.9.8c: includes error of Etax +Ver. 0.9.8d: Eta0, K0 etc in error loop replaced by Eta0y, K0y etc. Changes in signal noise calculations +Ver. 0.9.8e: ambiguous laser spec. DOLP (no discrimination between left and right circular polarisation) replaced by Stokes parameters Qin, Uin +Ver. 0.9.8e2: Added plot of LDRsim, Etax, Etapx, Etamx; LDRCorr and aLDRcorr consistently named +Ver. 0.9.8e3: Change of OutputFile name; Change of Ir and It noise if (CalcFrom0deg) = False; (Different calculation of error contributions tested but not implemented) +Ver. 0.9.8e4: text changed for y=+-1 (see line 274 ff and line 1044 ff +Ver. 0.9.8e5: changed: LDRunCorr = LDRsim / Etax +Ver. 0.9.8e6: K(0.05) instead K(0.02) +Ver. 0.9.8f: Ioannis Binietoglou => for LINUX compatibility: Using os.path.join should work in all operating systems. + # After line 1007: + output_path = os.path.join('output_files', OutputFile) + with open(output_path, 'w') as f: + # Line 1147 + file = open(output_path, 'r') +Ver. 0.9.8g: The LDRCal in the optic_input file is taken as first value for K(LDRCal) +Ver. 0.9.8h: 1. K(LDRCal) second order polynomial fit to values calculated for the listed LDRCal = 0.004 to 0.45 (without error calculation). + 2. Second order polynomial fits to upper and lower systematic error limits of the retrieved LDR depending on the true LDR. + + ======================================================== +simulation: LDRsim = Ir / It with variable parameters (possible truths) + G,H,Eta,Etax,K + It = TaT * TiT * ATP1 * TiO * TiE * (GT + atrue * HT) + LDRsim = Ir / It +consistency test: is forward simulation and correction consistent? + LDRCorr = (LDRsim / Eta * (GT + HT) - (GR + HR)) / ((GR - HR) - LDRsim / Eta * (GT - HT)) => atrue? +assumed true: G0,H0,Eta0,Etax0,K0 => actual retrievals of LDRCorr + => correct possible truths with assumed true G0,H0,Eta0 + measure: It, Ir, EtaX + LDRunCorr = LDRsim / Etax + correct it with G0,H0,K0: + LDRCorr = (LDRsim / (Etax / K0) * (GT0 + HT0) - (GR0 + HR0)) / ((GR0 - HR0) - LDRsim0 / (Etax / K0) * (GT0 - HT0)) +""" +# Comment: The code might 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 +# !/usr/bin/env python3 + +import os +import sys + +from scipy.stats import kurtosis +from scipy.stats import skew +# use: kurtosis(data, fisher=True,bias=False) => 0; skew(data,bias=False) => 0 +# Comment: the seaborn library makes nicer plots, but the code works also without it. +import numpy as np +import matplotlib.pyplot as plt + +try: + import seaborn as sns + + sns_loaded = True +except ImportError: + sns_loaded = False + +# from time import clock # python 2 +from timeit import default_timer as 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 +# Do you want to calculate the errors? If not, just the GHK-parameters are determined. +Error_Calc = True +LID = "internal" +EID = "internal" +# --- IL Laser IL and +-Uncertainty +Qin, dQin, nQin = 1., 0.0, 0 # second Stokes vector parameter; default 1 => linear polarization +Vin, dVin, nVin = 0., 0.0, 0 # fourth Stokes vector parameter +RotL, dRotL, nRotL = 0.0, 0.0, 1 # alpha; rotation of laser polarization in degrees; default 0 +# IL = 1e5 #photons in the laser beam, including detection efficiency of the telescope, atmodspheric and r^2 attenuation +# --- 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_RP_depend_on_TS_TP = False +if (RS_RP_depend_on_TS_TP): + RP, dRP, nRP = 1 - TP, 0.0, 0 + RS, dRS, nRS = 1 - TS, 0.0, 0 +else: + RP, dRP, nRP = 0.05, 0.01, 1 + RS, dRS, nRS = 0.98, 0.01, 1 +TiR = 0.5 * (RP + RS) +DiR = (RP - RS) / (RP + RS) +# 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) + +# +++ Orientation of the PBS with respect to the reference plane (see Polarisation-orientation.png and Polarisation-orientation-2.png in /system_settings) +# Y = +1: PBS incidence plane is parallel to the reference plane and polarisation in reference plane is finally transmitted. +# Y = -1: PBS incidence plane is perpendicular to the reference plane and polarisation in reference plane is finally reflected. +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 + +# --- Additional attenuation (transmission of the ND-filter) during the calibration +TCalT, dTCalT, nTCalT = 1, 0., 0 # transmitting path; error calc not working yet +TCalR, dTCalR, nTCalR = 1, 0., 0 # reflecting path; error calc not working yet + +# *** signal noise error calculation +# --- number of photon counts in the signal summed up in the calibration range during the calibration measurements +NCalT = 1e6 # default 1e6, assumed the same in +45° and -45° signals +NCalR = 1e6 # default 1e6, assumed the same in +45° and -45° signals +NILfac = 1.0 # duration of standard (0°) measurement relative to calibration measurements +nNCal = 0 # error nNCal: one-sigma in steps to left and right for calibration signals +nNI = 0 # error nNI: one-sigma in steps to left and right for 0° signals +NI = 50000 #number of photon counts in the parallel 0°-signal +eFacT = 1.0 # rel. amplification of transmitted channel, approximate values are sufficient; def. = 1 +eFacR = 10.0 +IoutTp0, IoutTp, dIoutTp0 = 0.5, 0.5, 0.0 +IoutTm0, IoutTm, dIoutTm0 = 0.5, 0.5, 0.0 +IoutRp0, IoutRp, dIoutRp0 = 0.5, 0.5, 0.0 +IoutRm0, IoutRm, dIoutRm0 = 0.5, 0.5, 0.0 +It0, It, dIt0 = 1 , 1, 0 +Ir0, Ir, dTr0 = 1 , 1, 0 +CalcFrom0deg = True + +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.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 +# BSR backscatter ratio +# BSR, dBSR, nBSR = 10, 0.05, 1 +BSR = np.zeros(5) +BSR = [1.1, 2, 5, 10., 50.] +# theoretical molecular LDR LDRm +LDRm, dLDRm, nLDRm = 0.004, 0.001, 1 +# 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.004 +LDRtrue2 = 0.004 +LDRunCorr = 1. +# 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 the laser', 'behind the emitter', 'before the receiver', 'before the PBS'] +Type = ['', 'Mechanical rotator', 'HWP rotator', 'Linear polarizer', 'QWP rotator', 'Circular polarizer', + 'real HWP +-22.5°'] + +bPlotEtax = False + +# end of initial definition of variables +# ******************************************************************************************************************************* +# ---- Read actual lidar system parameters from optic_input.py (must be in the programs sub-directory 'system_settings') +# ******************************************************************************************************************************* +# InputFile = 'optic_input_0.9.8e4-PollyXT_Lacros.py' +InputFile = 'optic_input_example_lidar_ver0.9.8e.py' +InputFile = 'calibrator-test-1-ver0.9.8e.py' +InputFile = 'optic_input_0.9.8e4-BRC-532.py' +InputFile = '' +InputFile = '' +InputFile = '' +InputFile = '' +InputFile = '' +InputFile = 'CIEMAT-200929-test-ver0.9.8f.py' +InputFile = 'calibrator-HWP-retard-error-ver0.9.8e.py' +InputFile = 'MUSA C1A-ver0.98e.py' +InputFile = 'EVE-210411-ver0.9.8f.py' +InputFile = 'Cilbolton-210609-ver0.9.8f.py' +InputFile = 'optic_input_0.9.8e4-PollyXT_Cyprus_355_210115_vf_4.py' +InputFile = 'MUSA C1B-ver0.98e.py' +InputFile = 'RALI-may2020-X.py' +InputFile = 'optic_input_0.9.8h-PollyXT_Cyprus_532_210429_vf_4.py' + +# InputFile = '' + +# ******************************************************************************************************************************* + +''' +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 +# # --- Errors +Qin0, dQin, nQin = Qin, dQin, nQin +Vin0, dVin, nVin = Vin, dVin, nVin +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 + +# BSR0, dBSR, nBSR = BSR, dBSR, nBSR +LDRm0, dLDRm, nLDRm = LDRm, dLDRm, nLDRm +# ---------- 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. +TCalT0, TCalR0 = TCalT, TCalR + +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 + +C2aT = np.cos(np.deg2rad(2. * RotaT)) +C2aR = np.cos(np.deg2rad(2. * RotaR)) +ATPT = float(1. + C2aT * DaT * DiT) +ARPT = float(1. + C2aR * DaR * DiR) +TTa = TiT * TaT * ATPT # unpolarized transmission +TRa = TiR * TaR * ARPT # unpolarized transmission +Eta0 = TRa / TTa + +# --- alternative texts for output +dY = ['perpendicular', '', 'parallel'] +dY2 = ['reflected', '', 'transmitted'] +if ((abs(RotL) < 45 and Y == 1) or (abs(RotL) >= 45 and Y == -1)): + dY3 = "the parallel laser polarisation is detected in the transmitted channel." +else: + dY3 = "the parallel laser polarisation is detected in the reflected channel." + +# --- check input errors +if ((Qin ** 2 + Vin ** 2) ** 0.5) > 1: + print("Error: degree of polarisation of laser > 1. Check Qin and Vin! ") + sys.exit() + +# --- this subroutine is for the calculation of the PLDR from LDR, BSR, and LDRm ------------------- +def CalcPLDR(LDR, BSR, LDRm): + PLDR = (BSR * (1. + LDRm) * LDR - LDRm * (1. + LDR)) / (BSR * (1. + LDRm) - (1. + LDR)) + return (PLDR) +# --- this subroutine is for the calculation with certain fixed parameters ------------------------ +def Calc(TCalT, TCalR, NCalT, NCalR, Qin, Vin, 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) + atrue = (1. - LDRtrue) / (1. + LDRtrue) + + # 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 + IinL = 1. + QinL = Qin + UinL = 0. + VinL = Vin + # VinL = (1. - DOLP ** 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 + # ------------------------- + + # analyser + if (RS_RP_depend_on_TS_TP): + RS = 1. - TS + RP = 1. - TP + + TiT = 0.5 * (TP + TS) + DiT = (TP - TS) / (TP + TS) + ZiT = (1. - DiT ** 2) ** 0.5 + 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)) + + # Analyzer As before the PBS Eq. D.5; combined PBS and cleaning pol-filter + ATPT = (1. + C2aT * DaT * DiT) # unpolarized transmission correction + TTa = TiT * TaT * ATPT # unpolarized transmission + ATP1 = 1. + ATP2 = Y * (DiT + C2aT * DaT) / ATPT + ATP3 = Y * S2aT * DaT * ZiT * CosT / ATPT + ATP4 = S2aT * DaT * ZiT * SinT / ATPT + ATP = np.array([ATP1, ATP2, ATP3, ATP4]) + DTa = ATP2 * Y + + ARPT = (1 + C2aR * DaR * DiR) # unpolarized transmission correction + TRa = TiR * TaR * ARPT # unpolarized transmission + ARP1 = 1 + ARP2 = Y * (DiR + C2aR * DaR) / ARPT + ARP3 = Y * S2aR * DaR * ZiR * CosR / ARPT + ARP4 = S2aR * DaR * ZiR * SinR / ARPT + ARP = np.array([ARP1, ARP2, ARP3, ARP4]) + DRa = ARP2 * Y + + + # ---- 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 parameters 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 parameters 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 parameters 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 parameters 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 parameters 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 parameters 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 parameters 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 parameters 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 parameters 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 = TTa * TiC * TiO * TiE * (AT + BT) + IoutTm = TTa * TiC * TiO * TiE * (AT - BT) + IoutRp = TRa * TiC * TiO * TiE * (AR + BR) + IoutRm = TRa * TiC * 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 = (TRa / TTa) # = TRa / TTa; Eta = Eta*/K Eq. 84 => K = Eta* / Eta; equation corrected according to the papers supplement Eqs. (S.10.10.1) ff + 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 + + # ----- random error calculation ---------- + # noise must be calculated with the photon counts of measured signals; + # relative standard deviation of calibration signals with LDRcal; assumed to be statisitcally independent + # normalised noise errors + if (CalcFrom0deg): + dIoutTp = (NCalT * IoutTp) ** -0.5 + dIoutTm = (NCalT * IoutTm) ** -0.5 + dIoutRp = (NCalR * IoutRp) ** -0.5 + dIoutRm = (NCalR * IoutRm) ** -0.5 + else: + dIoutTp = (NCalT ** -0.5) + dIoutTm = (NCalT ** -0.5) + dIoutRp = (NCalR ** -0.5) + dIoutRm = (NCalR ** -0.5) + # Forward simulated 0°-signals with LDRCal with atrue; from input file + + It = TTa * TiO * TiE * (GT + atrue * HT) + Ir = TRa * TiO * TiE * (GR + atrue * HR) + # relative standard deviation of standard signals with LDRmeas; assumed to be statisitcally independent + if (CalcFrom0deg): # this works! + dIt = ((It * NI * eFacT) ** -0.5) + dIr = ((Ir * NI * eFacR) ** -0.5) + ''' + dIt = ((NCalT * It / IoutTp * NILfac / TCalT) ** -0.5) + dIr = ((NCalR * Ir / IoutRp * NILfac / TCalR) ** -0.5) + ''' + else: # does this work? Why not as above? + dIt = ((NCalT * 2 * NILfac / TCalT ) ** -0.5) + dIr = ((NCalR * 2 * NILfac / TCalR) ** -0.5) + + # ----- Forward simulated 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.) and (abs(RotL0) < 45)) or ((Y == +1.) and (abs(RotL0) > 45)): + LDRsimx = 1. / LDRsim / Etax + else: + LDRsimx = LDRsim / Etax + ''' + LDRsimx = LDRsim + + # The following is correct without doubt + # LDRCorr = (LDRsim/(Etax/K)*(GT+HT)-(GR+HR))/((GR-HR)-LDRsim/(Etax/K)*(GT-HT)) + + # The following is a test whether the equations for calibration Etax and normal signal (GHK, LDRsim) are consistent + LDRCorr = (LDRsim / (Etax / K) * (GT + HT) - (GR + HR)) / ((GR - HR) - LDRsim / (Etax / K) * (GT - HT)) + # here we could also use Eta instead of Etax / K => how to test whether Etax is correct? => comparison with MüllerMatrix simulation! + # Without any correction: only measured It, Ir, EtaX are used + LDRunCorr = LDRsim / Etax + # LDRunCorr = (LDRsim / Etax * (GT / abs(GT) + HT / abs(HT)) - (GR / abs(GR) + HR / abs(HR))) / ((GR / abs(GR) - HR / abs(HR)) - LDRsim / Etax * (GT / abs(GT) - HT / abs(HT))) + + #LDRCorr = LDRsimx # for test only + + F11sim = 1 / (TiO * TiE) * ((HR * Eta * It - HT * Ir) / (HR * GT - HT * GR)) # IL = 1, Etat = Etar = 1 ; AMT Eq.64; what is Etax/K? => see about 20 lines above: = Eta + + return (IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr, + GT, HT, GR, HR, K, Eta, LDRsimx, LDRCorr, DTa, DRa, TTa, TRa, F11sim, LDRunCorr) + + + +# ******************************************************************************************************************************* + +# --- CALC with assumed true parameters from the input file +LDRtrue = LDRtrue2 +IoutTp0, IoutTm0, IoutRp0, IoutRm0, It0, Ir0, dIoutTp0, dIoutTm0, dIoutRp0, dIoutRm0, dIt0, dIr0, \ +GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \ +Calc(TCalT, TCalR, NCalT, NCalR, Qin0, Vin0, RotL0, RotE0, RetE0, DiE0, + RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, + ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0) +Etax0 = K0 * Eta0 +Etapx0 = IoutRp0 / IoutTp0 +Etamx0 = IoutRm0 / IoutTm0 +# --- Print parameters to console and output file +OutputFile = 'output_' + InputFile[0:-3] + '_' + fname[0:-3] +'.dat' +output_path = os.path.join('output_files', OutputFile) +# with open('output_files\\' + OutputFile, 'w') as f: +with open(output_path, '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:7}{1:17} {2:6.4f}±{3:7.4f}/{4:2d}".format("Laser: ", "Qin =", Qin0, dQin, nQin)) + print("{0:7}{1:17} {2:6.4f}±{3:7.4f}/{4:2d}".format("", "Vin =", Vin0, dVin, nVin)) + print("{0:7}{1:17} {2:6.4f}±{3:7.4f}/{4:2d}".format("", "Rotation alpha = ", RotL0, dRotL, nRotL)) + print("{0:7}{1:15} {2:8.4f} {3:17}".format("", "=> DOP", ((Qin ** 2 + Vin ** 2) ** 0.5), " (degree of polarisation)")) + + print("Optic: 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:9.6f}±{2:9.6f}/{8:2d}, {3:7.4f}, {4:3.0f}±{5:3.0f}/{9:2d}, {6:7.4f}±{7:7.4f}/{10:2d}".format( + "Calibrator ", DiC0, dDiC, TiC, RetC0, dRetC, RotC0, dRotC, nDiC, nRetC, nRotC)) + print("{0:12}".format(" Pol.-filter ------ ")) + print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format( + "ERT, RotT :", ERaT0, dERaT, nERaT, RotaT0, dRotaT, nRotaT)) + print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format( + "ERR, RotR :", ERaR0, dERaR, nERaR, RotaR0, dRotaR, nRotaR)) + print("{0:12}".format(" PBS ------ ")) + print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format( + "TP,TS :", TP0, dTP, nTP, TS0, dTS, nTS)) + print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format( + "RP,RS :", RP0, dRP, nRP, RS0, dRS, nRS)) + print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}, {5:1.0f}".format( + "DT,TT,DR,TR,Y :", DiT, TiT, DiR, TiR, Y)) + print("{0:12}".format(" Combined PBS + Pol.-filter ------ ")) + print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}".format( + "DT,TT,DR,TR :", DTa0, TTa0, DRa0, TRa0)) + print("{0:26}: {1:6.3f}± {2:5.3f}/{3:2d}".format( + "LDRCal during calibration in calibration range", LDRCal0, dLDRCal, nLDRCal)) + print("{0:12}".format(" --- Additional ND filter attenuation (transmission) during the calibration ---")) + print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format( + "TCalT,TCalR :", TCalT0, dTCalT, nTCalT, TCalR0, dTCalR, nTCalR)) + print() + print(Type[TypeC],"calibrator is located", Loc[LocC]) + print("Rotation error epsilon is considered also for normal measurements = ", RotationErrorEpsilonForNormalMeasurements) + print("The PBS incidence plane is ", dY[int(Y + 1)], "to the reference plane" ) + print("The laser polarisation in the reference plane is finally", dY2[int(Y + 1)], "=>", dY3) + print("RS_RP_depend_on_TS_TP = ", RS_RP_depend_on_TS_TP) + # end of print actual system parameters + # ****************************************************************************** + + + print() + + K0List = np.zeros(7) + LDRsimxList = np.zeros(7) + LDRCalList = LDRCal0, 0.004, 0.05, 0.1, 0.2, 0.3, 0.45 + # 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. + # Still with assumed true parameters in input file + + ''' + facIt = NCalT / TCalT0 * NILfac + facIr = NCalR / TCalR0 * NILfac + ''' + facIt = NI * eFacT + facIr = NI * eFacR + if (bPlotEtax): + # check error signals + # dIs are relative stdevs + print("LDRCal, IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp,dIoutTm,dIoutRp,dIoutRm,dIt, dIr") + + for i, LDRCal in enumerate(LDRCalList): + IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr, \ + GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \ + Calc(TCalT0, TCalR0, NCalT, NCalR, Qin0, Vin0, 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 + + if (bPlotEtax): + # check error signals + print( "{:0.2f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}".format(LDRCal, IoutTp * NCalT, IoutTm * NCalT, IoutRp * NCalR, IoutRm * NCalR, It * facIt, Ir * facIr, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr)) + #print( "{:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}".format(IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr)) + # end check error signals + print('===========================================================================================================') + print("{0:8},{1:8},{2:8},{3:8},{4:4}{5:5.3f}{6:1},{7:8},{8:9},{9:9},{10:9},{11:9},{12:9}".format( + " GR", " GT", " HR", " HT", " K(", LDRCalList[0], ")", " K(0.004)", " K(0.05)", " K(0.1)", " K(0.2)", " K(0.3)", " 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},{7:10.5f},{8:9.5f},{9:9.5f},{10:9.5f}".format( + GR0, GT0, HR0, HT0, K0List[0], K0List[1], K0List[2], K0List[3], K0List[4], K0List[5], K0List[6])) + print('===========================================================================================================') + print() + # --- Calculate 2nd order polynom-fit to K(LDRCalList) + K_poly = np.polyfit(LDRCalList, K0List, 2) + # p = np.poly1d(K_poly) + # f-string formatting + print("2nd order polyfit: K(LDRcal) = a + (b * LDRcal) + (c * LDRcal^2)") + print("(LDRcal is the volume linear depolarisation ratio in the range used for the polarisation calibration)") + print(" a = ", f"{K_poly[2]:.5f}") + print(" b = ", f"{K_poly[1]:.5f}") + print(" c = ", f"{K_poly[0]:.5f}") + print() + + print("Errors from neglecting GHK corrections and/or calibration:") + print("{0:>10},{1:>10},{2:>10},{3:>10},{4:>10},{5:>10}".format( + "LDRtrue", "LDRunCorr", "1/LDRunCorr", "LDRsimx", "1/LDRsimx", "LDRCorr")) + + aF11sim0 = np.zeros(5) + LDRrange = np.zeros(5) + LDRsim0 = np.zeros(5) + LDRrange = [0.004, 0.02, 0.1, 0.3, 0.45] # list + LDRrange[0] = LDRtrue2 # value in the input file; default 0.004 + + # The loop over LDRtrueList is only for checking how much the uncorrected LDRsimx deviates from LDRtrue ... and whether the corrections work. + # LDRsimx = LDRsim = Ir / It or 1/LDRsim + # Still with assumed true parameters in input file + for i, LDRtrue in enumerate(LDRrange): + #for LDRtrue in LDRrange: + IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr, \ + GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \ + Calc(TCalT0, TCalR0, NCalT, NCalR, Qin0, Vin0, RotL0, RotE0, RetE0, DiE0, + RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, + ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0) + print("{0:10.5f},{1:10.5f},{2:10.5f},{3:10.5f},{4:10.5f},{5:10.5f}".format(LDRtrue, LDRunCorr, 1/LDRunCorr, LDRsimx, 1/LDRsimx, LDRCorr)) + aF11sim0[i] = F11sim0 + LDRsim0[i] = Ir / It + # the assumed true aF11sim0 results will be used below to calc the deviation from the real signals + print("LDRsimx = LDR of the nominal system directly from measured signals without calibration and GHK-corrections") + print("LDRunCorr = LDR of the nominal system directly from measured signals with calibration but without GHK-corrections; electronic amplifications = 1 assumed") + print("LDRCorr = LDR calibrated and GHK-corrected") + print() + print("Errors from signal noise:") + print("Signal counts: NI, NCalT, NCalR, NILfac, nNCal, nNI, stdev(NI)/NI = {0:10.0f},{1:10.0f},{2:10.0f},{3:3.0f},{4:2.0f},{5:2.0f},{6:8.5f}".format( + NI, NCalT, NCalR, NILfac, nNCal, nNI, 1.0 / NI ** 0.5)) + print() + print() + '''# das muß wieder weg + print("IoutTp, IoutTm, IoutRp, IoutRm, It , Ir , dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr") + LDRCal = 0.01 + for i, LDRtrue in enumerate(LDRrange): + IoutTp, IoutTm, IoutRp, IoutRm, It, Ir, dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr, \ + GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \ + Calc(TCalT0, TCalR0, NCalT, NCalR, DOLP0, RotL0, RotE0, RetE0, DiE0, + RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, + ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0) + print( "{:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}, {:0.4f}".format( + IoutTp * NCalT, IoutTm * NCalT, IoutRp * NCalR, IoutRm * NCalR, It * facIt, Ir * facIr, + dIoutTp, dIoutTm, dIoutRp, dIoutRm, dIt, dIr)) + aF11sim0[i] = F11sim0 + # the assumed true aF11sim0 results will be used below to calc the deviation from the real signals + # bis hierher weg + ''' + +# file = open('output_files\\' + OutputFile, 'r') +file = open(output_path, 'r') +print(file.read()) +file.close() + +# --- CALC again assumed truth with LDRCal0 and with assumed true parameters in input file to reset all 0-values +LDRtrue = LDRtrue2 +IoutTp0, IoutTm0, IoutRp0, IoutRm0, It0, Ir0, dIoutTp0, dIoutTm0, dIoutRp0, dIoutRm0, dIt0, dIr0, \ +GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \ +Calc(TCalT0, TCalR0, NCalT, NCalR, Qin0, Vin0, RotL0, RotE0, RetE0, DiE0, + RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, + ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0) +Etax0 = K0 * Eta0 +Etapx0 = IoutRp0 / IoutTp0 +Etamx0 = IoutRm0 / IoutTm0 +''' +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 +''' +if (Error_Calc): + print("-----------------------------------------") + print("Errors from GHK correction uncertainties:") + print("-----------------------------------------") + # --- CALC again assumed truth with LDRCal0 and with assumed true parameters in input file to reset all 0-values + LDRtrue = LDRtrue2 + IoutTp0, IoutTm0, IoutRp0, IoutRm0, It0, Ir0, dIoutTp0, dIoutTm0, dIoutRp0, dIoutRm0, dIt0, dIr0, \ + GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0, LDRunCorr = \ + Calc(TCalT0, TCalR0, NCalT, NCalR, Qin0, Vin0, RotL0, RotE0, RetE0, DiE0, + RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, + ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0) + Etax0 = K0 * Eta0 + Etapx0 = IoutRp0 / IoutTp0 + Etamx0 = IoutRm0 / IoutTm0 + + # --- Start Error calculation with variable parameters ------------------------------------------------------------------ + # error nNCal: one-sigma in steps to left and right for calibration signals + # error nNI: one-sigma in steps to left and right for 0° signals + + iN = -1 + N = ((nTCalT * 2 + 1) * (nTCalR * 2 + 1) * + (nNCal * 2 + 1) ** 4 * (nNI * 2 + 1) ** 2 * + (nQin * 2 + 1) * (nVin * 2 + 1) * (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("number of system variations 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) + LDRstd = np.zeros(5) + LDRmean = np.zeros(5) + LDRmedian = np.zeros(5) + LDRskew = np.zeros(5) + LDRkurt = np.zeros(5) + LDRsimmin = np.zeros(5) + LDRsimmax = np.zeros(5) + LDRsimmean = np.zeros(5) + + F11min = np.zeros(5) + F11max = np.zeros(5) + Etaxmin = np.zeros(5) + Etaxmax = np.zeros(5) + + aQin = np.zeros(N) + aVin = 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) + aNCalTp = np.zeros(N) + aNCalTm = np.zeros(N) + aNCalRp = np.zeros(N) + aNCalRm = np.zeros(N) + aNIt = np.zeros(N) + aNIr = np.zeros(N) + aTCalT = np.zeros(N) + aTCalR = np.zeros(N) + + # each np.zeros((LDRrange, N)) array has the same N-dependency + aLDRcorr = np.zeros((5, N)) + aLDRsim = np.zeros((5, N)) + aF11corr = np.zeros((5, N)) + aPLDR = np.zeros((5, N)) + aEtax = np.zeros((5, N)) + aEtapx = np.zeros((5, N)) + aEtamx = np.zeros((5, N)) + + # np.zeros((GHKs, N)) + aGHK = np.zeros((5, N)) + + atime = clock() + dtime = clock() + + # --- Calc Error signals + # ---- Do the calculations of bra-ket vectors + h = -1. if TypeC == 2 else 1 + + for iLDRCal in range(-nLDRCal, nLDRCal + 1): + # from input file: LDRCal for calibration measurements + LDRCal = LDRCal0 + if nLDRCal > 0: + LDRCal = LDRCal0 + iLDRCal * dLDRCal / nLDRCal + # provides the intensities of the calibration measurements at various LDRCal for signal noise errors + # IoutTp, IoutTm, IoutRp, IoutRm, dIoutTp, dIoutTm, dIoutRp, dIoutRm + + aCal = (1. - LDRCal) / (1. + LDRCal) + for iQin, iVin, iRotL, iRotE, iRetE, iDiE \ + in [(iQin, iVin, iRotL, iRotE, iRetE, iDiE) + for iQin in range(-nQin, nQin + 1) + for iVin in range(-nVin, nVin + 1) + 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 nQin > 0: Qin = Qin0 + iQin * dQin / nQin + if nVin > 0: Vin = Vin0 + iVin * dVin / nVin + if nRotL > 0: RotL = RotL0 + iRotL * dRotL / nRotL + if nRotE > 0: RotE = RotE0 + iRotE * dRotE / nRotE + if nRetE > 0: RetE = RetE0 + iRetE * dRetE / nRetE + if nDiE > 0: DiE = DiE0 + iDiE * dDiE / nDiE + + if ((Qin ** 2 + Vin ** 2) ** 0.5) > 1.0: + print("Error: degree of polarisation of laser > 1. Check Qin and Vin! ") + sys.exit() + # 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 + IinL = 1. + QinL = Qin + UinL = 0. + VinL = Vin + # VinL = (1. - DOLP ** 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)]: + + 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) + + # analyser + # For POLLY_XTs + if (RS_RP_depend_on_TS_TP): + RS = 1.0 - TS + RP = 1.0 - 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)) + + # cleaning pol-filter + DaT = (1.0 - ERaT) / (1.0 + ERaT) + DaR = (1.0 - ERaR) / (1.0 + ERaR) + TaT = 0.5 * (1.0 + ERaT) + TaR = 0.5 * (1.0 + ERaR) + + S2aT = np.sin(np.deg2rad(h * 2.0 * RotaT)) + C2aT = np.cos(np.deg2rad(2.0 * RotaT)) + S2aR = np.sin(np.deg2rad(h * 2.0 * RotaR)) + C2aR = np.cos(np.deg2rad(2.0 * RotaR)) + + # Analyzer As before the PBS Eq. D.5; combined PBS and cleaning pol-filter + ATPT = (1 + C2aT * DaT * DiT) # unpolarized transmission correction + TTa = TiT * TaT * ATPT # unpolarized transmission + ATP1 = 1.0 + ATP2 = Y * (DiT + C2aT * DaT) / ATPT + ATP3 = Y * S2aT * DaT * ZiT * CosT / ATPT + ATP4 = S2aT * DaT * ZiT * SinT / ATPT + ATP = np.array([ATP1, ATP2, ATP3, ATP4]) + DTa = ATP2 * Y + + ARPT = (1 + C2aR * DaR * DiR) # unpolarized transmission correction + TRa = TiR * TaR * ARPT # unpolarized transmission + ARP1 = 1 + ARP2 = Y * (DiR + C2aR * DaR) / ARPT + ARP3 = Y * S2aR * DaR * ZiR * CosR / ARPT + ARP4 = S2aR * DaR * ZiR * SinR / ARPT + ARP = np.array([ARP1, ARP2, ARP3, ARP4]) + DRa = ARP2 * Y + + # ---- 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 parameters 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 parameters 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 parameters 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 parameters 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 parameters 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 parameters 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 parameters 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 parameters 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 parameters 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() + + for iTCalT, iTCalR, iNCalTp, iNCalTm, iNCalRp, iNCalRm, iNIt, iNIr \ + in [ + (iTCalT, iTCalR, iNCalTp, iNCalTm, iNCalRp, iNCalRm, iNIt, iNIr) + for iTCalT in range(-nTCalT, nTCalT + 1) # Etax + for iTCalR in range(-nTCalR, nTCalR + 1) # Etax + for iNCalTp in range(-nNCal, nNCal + 1) # noise error of calibration signals => Etax + for iNCalTm in range(-nNCal, nNCal + 1) # noise error of calibration signals => Etax + for iNCalRp in range(-nNCal, nNCal + 1) # noise error of calibration signals => Etax + for iNCalRm in range(-nNCal, nNCal + 1) # noise error of calibration signals => Etax + for iNIt in range(-nNI, nNI + 1) + for iNIr in range(-nNI, nNI + 1)]: + + # Calibration signals with aCal => Determination of the correction K of the real calibration factor + IoutTp = TTa * TiC * TiO * TiE * (AT + BT) + IoutTm = TTa * TiC * TiO * TiE * (AT - BT) + IoutRp = TRa * TiC * TiO * TiE * (AR + BR) + IoutRm = TRa * TiC * TiO * TiE * (AR - BR) + + if nTCalT > 0: TCalT = TCalT0 + iTCalT * dTCalT / nTCalT + if nTCalR > 0: TCalR = TCalR0 + iTCalR * dTCalR / nTCalR + # signal noise errors + # ----- random error calculation ---------- + # noise must be calculated from/with the actually measured signals; influence of TCalT, TCalR errors on noise are not considered ? + # actually measured signal counts are in input file and don't change + # relative standard deviation of calibration signals with LDRcal; assumed to be statisitcally independent + # error nNCal: one-sigma in steps to left and right for calibration signals + if nNCal > 0: + if (CalcFrom0deg): + dIoutTp = (NCalT * IoutTp) ** -0.5 + dIoutTm = (NCalT * IoutTm) ** -0.5 + dIoutRp = (NCalR * IoutRp) ** -0.5 + dIoutRm = (NCalR * IoutRm) ** -0.5 + else: + dIoutTp = dIoutTp0 * (IoutTp / IoutTp0) + dIoutTm = dIoutTm0 * (IoutTm / IoutTm0) + dIoutRp = dIoutRp0 * (IoutRp / IoutRp0) + dIoutRm = dIoutRm0 * (IoutRm / IoutRm0) + # print(iTCalT, iTCalR, iNCalTp, iNCalTm, iNCalRp, iNCalRm, iNIt, iNIr, IoutTp, dIoutTp) + IoutTp = IoutTp * (1. + iNCalTp * dIoutTp / nNCal) + IoutTm = IoutTm * (1. + iNCalTm * dIoutTm / nNCal) + IoutRp = IoutRp * (1. + iNCalRp * dIoutRp / nNCal) + IoutRm = IoutRm * (1. + iNCalRm * dIoutRm / nNCal) + + IoutTp = IoutTp * TCalT / TCalT0 + IoutTm = IoutTm * TCalT / TCalT0 + IoutRp = IoutRp * TCalR / TCalR0 + IoutRm = IoutRm * TCalR / TCalR0 + # --- Results and Corrections; electronic etaR and etaT are assumed to be 1 for true and assumed true systems + # calibration factor + Eta = (TRa / TTa) # = TRa / TTa; Eta = Eta*/K Eq. 84; corrected according to the papers supplement Eqs. (S.10.10.1) ff + # possibly real calibration factor + Etapx = IoutRp / IoutTp + Etamx = IoutRm / IoutTm + Etax = (Etapx * Etamx) ** 0.5 + K = Etax / Eta + # 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 + ''' + for iIt, iIr \ + in [(iIt, iIr) + for iIt in range(-nNI, nNI + 1) + for iIr in range(-nNI, nNI + 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 + + # *** loop for different real LDRs ********************************************************************** + iLDR = -1 + for LDRTrue in LDRrange: + iLDR = iLDR + 1 + atrue = (1. - LDRTrue) / (1. + LDRTrue) + # ----- Forward simulated signals and LDRsim with atrue; from input file; not considering TiC. + It = TTa * TiO * TiE * (GT + atrue * HT) # TaT*TiT*TiC*TiO*IinL*(GT+atrue*HT) + Ir = TRa * TiO * TiE * (GR + atrue * HR) # TaR*TiR*TiC*TiO*IinL*(GR+atrue*HR) + # # signal noise errors; standard deviation of signals; assumed to be statisitcally independent + # because the signals depend on LDRtrue, the errors dIt and dIr must be calculated for each LDRtrue + if (CalcFrom0deg): + ''' + dIt = ((NCalT * It / IoutTp * NILfac / TCalT) ** -0.5) + dIr = ((NCalR * Ir / IoutRp * NILfac / TCalR) ** -0.5) + ''' + dIt = ((It * NI * eFacT) ** -0.5) + dIr = ((Ir * NI * eFacR) ** -0.5) + else: + dIt = ((It * NI * eFacT) ** -0.5) + dIr = ((Ir * NI * eFacR) ** -0.5) + ''' + # does this work? Why not as above? + dIt = ((NCalT * 2. * NILfac / TCalT ) ** -0.5) + dIr = ((NCalR * 2. * NILfac / TCalR) ** -0.5) + ''' + # error nNI: one-sigma in steps to left and right for 0° signals + if nNI > 0: + It = It * (1. + iNIt * dIt / nNI) + Ir = Ir * (1. + iNIr * dIr / nNI) + + # LDRsim = 1/Eta*Ir/It # simulated LDR* with Y from input file + LDRsim = Ir / It # simulated uncorrected LDR with Y from input file + + # ----- Backward correction + # Corrected LDRCorr with assumed true G0,H0,K0,Eta0 from forward simulated (real) LDRsim(atrue) + LDRCorr = (LDRsim / (Etax / K0) * (GT0 + HT0) - (GR0 + HR0)) / ((GR0 - HR0) - LDRsim / (Etax / K0) * (GT0 - HT0)) + + # 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 / Eta * (GT - HT)) + # Without any correction + LDRunCorr = LDRsim / Etax + # LDRunCorr = (LDRsim / Etax * (GT / abs(GT) + HT / abs(HT)) - (GR / abs(GR) + HR / abs(HR))) / ((GR / abs(GR) - HR / abs(HR)) - LDRsim / Etax * (GT / abs(GT) - HT / abs(HT))) + + + ''' + # -- 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); Etax/K0 = Eta0. + ''' + # Corrected F11corr with assumed true G0,H0,K0 from forward simulated (real) It and Ir (atrue) + Text1 = "!!! EXPERIMENTAL !!! F11corr from real It and Ir with real calibration EtaX: x-axis: F11corr(LDRtrue) / aF11sim0(LDRtrue) - 1" + F11corr = 1 / (TiO * TiE) * ( + (HR0 * Etax / K0 * It / TTa - HT0 * Ir / TRa) / (HR0 * GT0 - HT0 * GR0)) # IL = 1 Eq.(64); Etax/K0 = Eta0. + + # 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 + aLDRcorr[iLDR, iN] = LDRCorr # LDRCorr # LDRsim # for test only + aLDRsim[iLDR, iN] = LDRsim # LDRCorr # LDRsim # for test only + # aPLDR[iLDR, iN] = CalcPLDR(LDRCorr, BSR[iLDR], LDRm0) + aEtax[iLDR, iN] = Etax + aEtapx[iLDR, iN] = Etapx + aEtamx[iLDR, iN] = Etamx + + aGHK[0, iN] = GR + aGHK[1, iN] = GT + aGHK[2, iN] = HR + aGHK[3, iN] = HT + aGHK[4, iN] = K + + aLDRCal[iN] = iLDRCal + aQin[iN] = iQin + aVin[iN] = iVin + 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 + aTCalT[iN] = iTCalT + aTCalR[iN] = iTCalR + + aNCalTp[iN] = iNCalTp # IoutTp, IoutTm, IoutRp, IoutRm => Etax + aNCalTm[iN] = iNCalTm # IoutTp, IoutTm, IoutRp, IoutRm => Etax + aNCalRp[iN] = iNCalRp # IoutTp, IoutTm, IoutRp, IoutRm => Etax + aNCalRm[iN] = iNCalRm # IoutTp, IoutTm, IoutRp, IoutRm => Etax + aNIt[iN] = iNIt # It, Tr + aNIr[iN] = iNIr # It, Tr + + # --- END loop + btime = clock() + # print("\r done in ", "{0:5.0f}".format(btime - atime), "sec. => producing plots now .... some more seconds ..."), # , end="\r"); + print(" done in ", "{0:5.0f}".format(btime - atime), "sec. => producing plots now .... some more seconds ...") + # --- Plot ----------------------------------------------------------------- + # print("Errors from GHK correction uncertainties:") + if (sns_loaded): + sns.set_style("whitegrid") + sns.set_palette("bright6", 6) + # for older seaborn versions use: + # sns.set_palette("bright", 6) + + ''' + fig2 = plt.figure() + plt.plot(aLDRcorr[2,:],'b.') + plt.plot(aLDRcorr[3,:],'r.') + plt.plot(aLDRcorr[4,:],'g.') + #plt.plot(aLDRcorr[6,:],'c.') + plt.show + ''' + + # Plot LDR + def PlotSubHist(aVar, aX, X0, daX, iaX, naX): + # aVar is the name of the parameter and aX is the subset of aLDRcorr which is coloured in the plot + # example: PlotSubHist("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP) + fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) + iLDR = -1 + for LDRTrue in LDRrange: + aXmean = np.zeros(2 * naX + 1) + iLDR = iLDR + 1 + LDRmin[iLDR] = np.amin(aLDRcorr[iLDR, :]) + LDRmax[iLDR] = np.amax(aLDRcorr[iLDR, :]) + if (LDRmax[iLDR] > 10): LDRmax[iLDR] = 10 + if (LDRmin[iLDR] < -10): LDRmin[iLDR] = -10 + Rmin = LDRmin[iLDR] * 0.995 # np.min(aLDRcorr[iLDR,:]) * 0.995 + Rmax = LDRmax[iLDR] * 1.005 # np.max(aLDRcorr[iLDR,:]) * 1.005 + + # Determine mean distance of all aXmean from each other for each iLDR + meanDist = 0.0 + for iaX in range(-naX, naX + 1): + # mean LDRCorr value for certain error (iaX) of parameter aVar + aXmean[iaX + naX] = np.mean(aLDRcorr[iLDR, aX == iaX]) + # relative to absolute spread of LDRCorrs + meanDist = (np.max(aXmean) - np.min(aXmean)) / (LDRmax[iLDR] - LDRmin[iLDR]) * 100 + + plt.subplot(1, 5, iLDR + 1) + (n, bins, patches) = plt.hist(aLDRcorr[iLDR, :], + bins=100, log=False, + range=[Rmin, Rmax], + alpha=0.5, density=False, color='0.5', histtype='stepfilled') + + for iaX in range(-naX, naX + 1): + # mean LDRCorr value for certain error (iaX) of parameter aVar + plt.hist(aLDRcorr[iLDR, aX == iaX], + range=[Rmin, Rmax], + bins=100, log=False, alpha=0.3, density=False, histtype='stepfilled', + label=str(round(X0 + iaX * daX / naX, 5))) + + if (iLDR == 2): + leg = plt.legend() + leg.get_frame().set_alpha(0.1) + + plt.tick_params(axis='both', labelsize=10) + plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2) + plt.gca().set_title("{0:3.0f}%".format(meanDist)) + plt.gca().set_xlabel('LDRtrue', color="red") + + # plt.ylabel('frequency', fontsize=10) + # plt.xlabel('LDRCorr', fontsize=10) + # fig.tight_layout() + fig.suptitle(LID + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar + ' error contribution', fontsize=14, y=1.10) + # plt.show() + # fig.savefig(LID + '_' + aVar + '.png', dpi=150, bbox_inches='tight', pad_inches=0) + # plt.close + return + + def PlotLDRsim(aVar, aX, X0, daX, iaX, naX): + # aVar is the name of the parameter and aX is the subset of aLDRsim which is coloured in the plot + # example: PlotSubHist("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP) + fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) + iLDR = -1 + for LDRTrue in LDRrange: + aXmean = np.zeros(2 * naX + 1) + iLDR = iLDR + 1 + LDRsimmin[iLDR] = np.amin(aLDRsim[iLDR, :]) + LDRsimmax[iLDR] = np.amax(aLDRsim[iLDR, :]) + # print("LDRsimmin[iLDR], LDRsimmax[iLDR] = ", LDRsimmin[iLDR], LDRsimmax[iLDR]) + # if (LDRsimmax[iLDR] > 10): LDRsimmax[iLDR] = 10 + # if (LDRsimmin[iLDR] < -10): LDRsimmin[iLDR] = -10 + Rmin = LDRsimmin[iLDR] * 0.995 # np.min(aLDRsim[iLDR,:]) * 0.995 + Rmax = LDRsimmax[iLDR] * 1.005 # np.max(aLDRsim[iLDR,:]) * 1.005 + # print("Rmin, Rmax = ", Rmin, Rmax) + + # Determine mean distance of all aXmean from each other for each iLDR + meanDist = 0.0 + for iaX in range(-naX, naX + 1): + # mean LDRCorr value for certain error (iaX) of parameter aVar + aXmean[iaX + naX] = np.mean(aLDRsim[iLDR, aX == iaX]) + # relative to absolute spread of LDRCorrs + meanDist = (np.max(aXmean) - np.min(aXmean)) / (LDRsimmax[iLDR] - LDRsimmin[iLDR]) * 100 + + plt.subplot(1, 5, iLDR + 1) + (n, bins, patches) = plt.hist(aLDRsim[iLDR, :], + bins=100, log=False, + range=[Rmin, Rmax], + alpha=0.5, density=False, color='0.5', histtype='stepfilled') + + for iaX in range(-naX, naX + 1): + # mean LDRCorr value for certain error (iaX) of parameter aVar + plt.hist(aLDRsim[iLDR, aX == iaX], + range=[Rmin, Rmax], + bins=100, log=False, alpha=0.3, density=False, histtype='stepfilled', + label=str(round(X0 + iaX * daX / naX, 5))) + + if (iLDR == 2): + leg = plt.legend() + leg.get_frame().set_alpha(0.1) + + plt.tick_params(axis='both', labelsize=10) + plt.plot([LDRsim0[iLDR], LDRsim0[iLDR]], [0, np.max(n)], 'r-', lw=2) + plt.gca().set_title("{0:3.0f}%".format(meanDist)) + plt.gca().set_xlabel('LDRsim0', color="red") + + fig.suptitle('LDRsim - ' +LID + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar + ' error contribution', fontsize=14, y=1.10) + return + + + # Plot Etax + def PlotEtax(aVar, aX, X0, daX, iaX, naX): + # aVar is the name of the parameter and aX is the subset of aLDRcorr which is coloured in the plot + # example: PlotSubHist("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP) + fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) + iLDR = -1 + for LDRTrue in LDRrange: + aXmean = np.zeros(2 * naX + 1) + iLDR = iLDR + 1 + Etaxmin = np.amin(aEtax[iLDR, :]) + Etaxmax = np.amax(aEtax[iLDR, :]) + Rmin = Etaxmin * 0.995 # np.min(aLDRcorr[iLDR,:]) * 0.995 + Rmax = Etaxmax * 1.005 # np.max(aLDRcorr[iLDR,:]) * 1.005 + + # Determine mean distance of all aXmean from each other for each iLDR + meanDist = 0.0 + for iaX in range(-naX, naX + 1): + # mean Etax value for certain error (iaX) of parameter aVar + aXmean[iaX + naX] = np.mean(aEtax[iLDR, aX == iaX]) + # relative to absolute spread of Etax + meanDist = (np.max(aXmean) - np.min(aXmean)) / (Etaxmax - Etaxmin) * 100 + + plt.subplot(1, 5, iLDR + 1) + (n, bins, patches) = plt.hist(aEtax[iLDR, :], + bins=50, log=False, + range=[Rmin, Rmax], + alpha=0.5, density=False, color='0.5', histtype='stepfilled') + for iaX in range(-naX, naX + 1): + plt.hist(aEtax[iLDR, aX == iaX], + range=[Rmin, Rmax], + bins=50, log=False, alpha=0.3, density=False, histtype='stepfilled', + label=str(round(X0 + iaX * daX / naX, 5))) + if (iLDR == 2): + leg = plt.legend() + leg.get_frame().set_alpha(0.1) + plt.tick_params(axis='both', labelsize=10) + plt.plot([Etax0, Etax0], [0, np.max(n)], 'r-', lw=2) + plt.gca().set_title("{0:3.0f}%".format(meanDist)) + plt.gca().set_xlabel('Etax0', color="red") + fig.suptitle('Etax - ' + LID + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar + ' error contribution', fontsize=14, y=1.10) + return + + def PlotEtapx(aVar, aX, X0, daX, iaX, naX): + # aVar is the name of the parameter and aX is the subset of aLDRcorr which is coloured in the plot + # example: PlotSubHist("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP) + fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) + iLDR = -1 + for LDRTrue in LDRrange: + aXmean = np.zeros(2 * naX + 1) + iLDR = iLDR + 1 + Etapxmin = np.amin(aEtapx[iLDR, :]) + Etapxmax = np.amax(aEtapx[iLDR, :]) + Rmin = Etapxmin * 0.995 # np.min(aLDRcorr[iLDR,:]) * 0.995 + Rmax = Etapxmax * 1.005 # np.max(aLDRcorr[iLDR,:]) * 1.005 + + # Determine mean distance of all aXmean from each other for each iLDR + meanDist = 0.0 + for iaX in range(-naX, naX + 1): + # mean Etapx value for certain error (iaX) of parameter aVar + aXmean[iaX + naX] = np.mean(aEtapx[iLDR, aX == iaX]) + # relative to absolute spread of Etapx + meanDist = (np.max(aXmean) - np.min(aXmean)) / (Etapxmax - Etapxmin) * 100 + + plt.subplot(1, 5, iLDR + 1) + (n, bins, patches) = plt.hist(aEtapx[iLDR, :], + bins=50, log=False, + range=[Rmin, Rmax], + alpha=0.5, density=False, color='0.5', histtype='stepfilled') + for iaX in range(-naX, naX + 1): + plt.hist(aEtapx[iLDR, aX == iaX], + range=[Rmin, Rmax], + bins=50, log=False, alpha=0.3, density=False, histtype='stepfilled', + label=str(round(X0 + iaX * daX / naX, 5))) + if (iLDR == 2): + leg = plt.legend() + leg.get_frame().set_alpha(0.1) + plt.tick_params(axis='both', labelsize=10) + plt.plot([Etapx0, Etapx0], [0, np.max(n)], 'r-', lw=2) + plt.gca().set_title("{0:3.0f}%".format(meanDist)) + plt.gca().set_xlabel('Etapx0', color="red") + fig.suptitle('Etapx - ' + LID + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar + ' error contribution', fontsize=14, y=1.10) + return + + def PlotEtamx(aVar, aX, X0, daX, iaX, naX): + # aVar is the name of the parameter and aX is the subset of aLDRcorr which is coloured in the plot + # example: PlotSubHist("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP) + fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) + iLDR = -1 + for LDRTrue in LDRrange: + aXmean = np.zeros(2 * naX + 1) + iLDR = iLDR + 1 + Etamxmin = np.amin(aEtamx[iLDR, :]) + Etamxmax = np.amax(aEtamx[iLDR, :]) + Rmin = Etamxmin * 0.995 # np.min(aLDRcorr[iLDR,:]) * 0.995 + Rmax = Etamxmax * 1.005 # np.max(aLDRcorr[iLDR,:]) * 1.005 + + # Determine mean distance of all aXmean from each other for each iLDR + meanDist = 0.0 + for iaX in range(-naX, naX + 1): + # mean Etamx value for certain error (iaX) of parameter aVar + aXmean[iaX + naX] = np.mean(aEtamx[iLDR, aX == iaX]) + # relative to absolute spread of Etamx + meanDist = (np.max(aXmean) - np.min(aXmean)) / (Etamxmax - Etamxmin) * 100 + + plt.subplot(1, 5, iLDR + 1) + (n, bins, patches) = plt.hist(aEtamx[iLDR, :], + bins=50, log=False, + range=[Rmin, Rmax], + alpha=0.5, density=False, color='0.5', histtype='stepfilled') + for iaX in range(-naX, naX + 1): + plt.hist(aEtamx[iLDR, aX == iaX], + range=[Rmin, Rmax], + bins=50, log=False, alpha=0.3, density=False, histtype='stepfilled', + label=str(round(X0 + iaX * daX / naX, 5))) + if (iLDR == 2): + leg = plt.legend() + leg.get_frame().set_alpha(0.1) + plt.tick_params(axis='both', labelsize=10) + plt.plot([Etamx0, Etamx0], [0, np.max(n)], 'r-', lw=2) + plt.gca().set_title("{0:3.0f}%".format(meanDist)) + plt.gca().set_xlabel('Etamx0', color="red") + fig.suptitle('Etamx - ' + LID + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar + ' error contribution', fontsize=14, y=1.10) + return + + # calc contribution of the error of aVar = aX to aY for each LDRtrue + def Contribution(aVar, aX, X0, daX, iaX, naX, aY, Ysum, widthSum): + # aVar is the name of the parameter and aX is the subset of aY which is coloured in the plot + # example: Contribution("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP, aLDRcorr, DOLPcontr) + iLDR = -1 + # Ysum, widthSum = np.zeros(5) + meanDist = np.zeros(5) # iLDR + widthDist = np.zeros(5) # iLDR + for LDRTrue in LDRrange: + aXmean = np.zeros(2 * naX + 1) + aXwidth = np.zeros(2 * naX + 1) + iLDR = iLDR + 1 + # total width of distribution + aYmin = np.amin(aY[iLDR, :]) + aYmax = np.amax(aY[iLDR, :]) + aYwidth = aYmax - aYmin + # Determine mean distance of all aXmean from each other for each iLDR + for iaX in range(-naX, naX + 1): + # mean LDRCorr value for all errors iaX of parameter aVar + aXmean[iaX + naX] = np.mean(aY[iLDR, aX == iaX]) + aXwidth[iaX + naX] = np.max(aY[iLDR, aX == iaX]) - np.min(aY[iLDR, aX == iaX]) + # relative to absolute spread of LDRCorrs + meanDist[iLDR] = (np.max(aXmean) - np.min(aXmean)) / aYwidth * 1000 + # meanDist[iLDR] = (aYwidth - aXwidth[naX]) / aYwidth * 1000 + widthDist[iLDR] = (np.max(aXwidth) - aXwidth[naX]) / aYwidth * 1000 + + print("{:12}{:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f}"\ + .format(aVar,meanDist[0],meanDist[1],meanDist[2],meanDist[3],meanDist[4],widthDist[0],widthDist[1],widthDist[2],widthDist[3],widthDist[4])) + Ysum = Ysum + meanDist + widthSum = widthSum + widthDist + return(Ysum, widthSum) + + # print(.format(LDRrangeA[iLDR],)) + + # error contributions to a certain output aY; loop over all variables + def Contribution_aY(aYvar, aY): + Ysum = np.zeros(5) + widthSum = np.zeros(5) + # meanDist = np.zeros(5) # iLDR + LDRrangeA = np.array(LDRrange) + print() + print(aYvar + ": contribution to the total error (per mill)") + print(" of individual parameter errors of combined parameter errors") + print(" at LDRtrue {:5.3f} {:5.3f} {:5.3f} {:5.3f} {:5.3f} {:5.3f} {:5.3f} {:5.3f} {:5.3f} {:5.3f}"\ + .format(LDRrangeA[0],LDRrangeA[1],LDRrangeA[2],LDRrangeA[3],LDRrangeA[4],LDRrangeA[0],LDRrangeA[1],LDRrangeA[2],LDRrangeA[3],LDRrangeA[4])) + print() + if (nQin > 0): Ysum, widthSum = Contribution("Qin", aQin, Qin0, dQin, iQin, nQin, aY, Ysum, widthSum) + if (nVin > 0): Ysum, widthSum = Contribution("Vin", aVin, Vin0, dVin, iVin, nVin, aY, Ysum, widthSum) + if (nRotL > 0): Ysum, widthSum = Contribution("RotL", aRotL, RotL0, dRotL, iRotL, nRotL, aY, Ysum, widthSum) + if (nRetE > 0): Ysum, widthSum = Contribution("RetE", aRetE, RetE0, dRetE, iRetE, nRetE, aY, Ysum, widthSum) + if (nRotE > 0): Ysum, widthSum = Contribution("RotE", aRotE, RotE0, dRotE, iRotE, nRotE, aY, Ysum, widthSum) + if (nDiE > 0): Ysum, widthSum = Contribution("DiE", aDiE, DiE0, dDiE, iDiE, nDiE, aY, Ysum, widthSum) + if (nRetO > 0): Ysum, widthSum = Contribution("RetO", aRetO, RetO0, dRetO, iRetO, nRetO, aY, Ysum, widthSum) + if (nRotO > 0): Ysum, widthSum = Contribution("RotO", aRotO, RotO0, dRotO, iRotO, nRotO, aY, Ysum, widthSum) + if (nDiO > 0): Ysum, widthSum = Contribution("DiO", aDiO, DiO0, dDiO, iDiO, nDiO, aY, Ysum, widthSum) + if (nDiC > 0): Ysum, widthSum = Contribution("DiC", aDiC, DiC0, dDiC, iDiC, nDiC, aY, Ysum, widthSum) + if (nRotC > 0): Ysum, widthSum = Contribution("RotC", aRotC, RotC0, dRotC, iRotC, nRotC, aY, Ysum, widthSum) + if (nRetC > 0): Ysum, widthSum = Contribution("RetC", aRetC, RetC0, dRetC, iRetC, nRetC, aY, Ysum, widthSum) + if (nTP > 0): Ysum, widthSum = Contribution("TP", aTP, TP0, dTP, iTP, nTP, aY, Ysum, widthSum) + if (nTS > 0): Ysum, widthSum = Contribution("TS", aTS, TS0, dTS, iTS, nTS, aY, Ysum, widthSum) + if (nRP > 0): Ysum, widthSum = Contribution("RP", aRP, RP0, dRP, iRP, nRP, aY, Ysum, widthSum) + if (nRS > 0): Ysum, widthSum = Contribution("RS", aRS, RS0, dRS, iRS, nRS, aY, Ysum, widthSum) + if (nRetT > 0): Ysum, widthSum = Contribution("RetT", aRetT, RetT0, dRetT, iRetT, nRetT, aY, Ysum, widthSum) + if (nRetR > 0): Ysum, widthSum = Contribution("RetR", aRetR, RetR0, dRetR, iRetR, nRetR, aY, Ysum, widthSum) + if (nERaT > 0): Ysum, widthSum = Contribution("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT, aY, Ysum, widthSum) + if (nERaR > 0): Ysum, widthSum = Contribution("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR, aY, Ysum, widthSum) + if (nRotaT > 0): Ysum, widthSum = Contribution("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT, aY, Ysum, widthSum) + if (nRotaR > 0): Ysum, widthSum = Contribution("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR, aY, Ysum, widthSum) + if (nLDRCal > 0): Ysum, widthSum = Contribution("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal, aY, Ysum, widthSum) + if (nTCalT > 0): Ysum, widthSum = Contribution("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT, aY, Ysum, widthSum) + if (nTCalR > 0): Ysum, widthSum = Contribution("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR, aY, Ysum, widthSum) + if (nNCal > 0): Ysum, widthSum = Contribution("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal, aY, Ysum, widthSum) + if (nNCal > 0): Ysum, widthSum = Contribution("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal, aY, Ysum, widthSum) + if (nNCal > 0): Ysum, widthSum = Contribution("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal, aY, Ysum, widthSum) + if (nNCal > 0): Ysum, widthSum = Contribution("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal, aY, Ysum, widthSum) + if (nNI > 0): Ysum, widthSum = Contribution("SigNoiseIt", aNIt, 0, 1, iNIt, nNI, aY, Ysum, widthSum) + if (nNI > 0): Ysum, widthSum = Contribution("SigNoiseIr", aNIr, 0, 1, iNIr, nNI, aY, Ysum, widthSum) + print("{:12}{:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f} {:5.0f}"\ + .format("Sum ",Ysum[0],Ysum[1],Ysum[2],Ysum[3],Ysum[4],widthSum[0],widthSum[1],widthSum[2],widthSum[3],widthSum[4])) + + + # Plot LDR histograms + if (nQin > 0): PlotSubHist("Qin", aQin, Qin0, dQin, iQin, nQin) + if (nVin > 0): PlotSubHist("Vin", aVin, Vin0, dVin, iVin, nVin) + 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) + if (nTCalT > 0): PlotSubHist("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT) + if (nTCalR > 0): PlotSubHist("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR) + if (nNCal > 0): PlotSubHist("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal) + if (nNCal > 0): PlotSubHist("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal) + if (nNCal > 0): PlotSubHist("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal) + if (nNCal > 0): PlotSubHist("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal) + if (nNI > 0): PlotSubHist("SigNoiseIt", aNIt, 0, 1, iNIt, nNI) + if (nNI > 0): PlotSubHist("SigNoiseIr", aNIr, 0, 1, iNIr, nNI) + plt.show() + plt.close + + + + # --- Plot LDRmin, LDRmax + iLDR = -1 + for LDRTrue in LDRrange: + iLDR = iLDR + 1 + LDRmin[iLDR] = np.amin(aLDRcorr[iLDR, :]) + LDRmax[iLDR] = np.amax(aLDRcorr[iLDR, :]) + LDRstd[iLDR] = np.std(aLDRcorr[iLDR, :]) + LDRmean[iLDR] = np.mean(aLDRcorr[iLDR, :]) + LDRmedian[iLDR] = np.median(aLDRcorr[iLDR, :]) + LDRskew[iLDR] = skew(aLDRcorr[iLDR, :],bias=False) + LDRkurt[iLDR] = kurtosis(aLDRcorr[iLDR, :],fisher=True,bias=False) + + fig2 = plt.figure() + LDRrangeA = np.array(LDRrange) + if((np.amax(LDRmax - LDRrangeA)-np.amin(LDRmin - LDRrangeA)) < 0.001): + plt.ylim(-0.001,0.001) + plt.plot(LDRrangeA, LDRmax - LDRrangeA, linewidth=2.0, color='b') + plt.plot(LDRrangeA, LDRmin - LDRrangeA, linewidth=2.0, color='g') + + plt.xlabel('LDRtrue', fontsize=18) + plt.ylabel('LDR: min - true, max - true', fontsize=14) + plt.title(LID + ' ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]), fontsize=18) + # plt.ylimit(-0.07, 0.07) + plt.show() + plt.close + + # Print LDR statistics + print("LDRtrue , mean , median, max-true, min-true, std, excess_kurtosis, skewness") + iLDR = -1 + LDRrangeA = np.array(LDRrange) + for LDRTrue in LDRrange: + iLDR = iLDR + 1 + print("{0:8.5f},{1:8.5f},{2:8.5f}, {3:8.5f},{4:8.5f},{5:8.5f}, {6:8.5f},{7:8.5f}".format(LDRrangeA[iLDR], LDRmean[iLDR], LDRmedian[iLDR], LDRmax[iLDR]-LDRrangeA[iLDR], LDRmin[iLDR]-LDRrangeA[iLDR], LDRstd[iLDR],LDRkurt[iLDR],LDRskew[iLDR])) + print() + # --- Calculate 2nd order polynom fit to lower and upper LDR errors between LDR = 0.004 and 0.45 and also 2nd order fit to K(LDRcal) + LDRmin_poly = np.polyfit(LDRrangeA, (LDRmin-LDRrangeA), 2) + LDRmax_poly = np.polyfit(LDRrangeA, (LDRmax-LDRrangeA), 2) + print("LDRCal is the volume linear depolarisation ratio in the range used for the polarisation calibration. K(LDRCal) is the variability of K depending on LDRCal.") + print("If LDRCal is not known, its uncertainty has to be included in the input file for the error calcualtion below.") + print("LDR_min-true(LDR_true) and LDR_max-true(LDR_true) are the difference between the lowest and highest possible retrieved LDR, respectively, and the true LDR depending on the true LDR. They represent the negative and positive uncertainty of the retrieved LDR resulting from all known uncertainties of the optical input parameters.") + print("2nd order polyfit = a + (b * LDR) + (c * LDR^2)") + print() + print("K(LDRCal), LDR_min-true(LDRtrue), LDR_max-true(LDRtrue)") + print(" a, ", f"{K_poly[2]:+.5f},", f"{LDRmin_poly[2]:+.5f},", f"{LDRmax_poly[2]:+.5f}") + print(" b, ", f"{K_poly[1]:+.5f},", f"{LDRmin_poly[1]:+.5f},", f"{LDRmax_poly[1]:+.5f}") + print(" c, ", f"{K_poly[0]:+.5f},", f"{LDRmin_poly[0]:+.5f},", f"{LDRmax_poly[0]:+.5f}") + print() + + + # --- Save LDRmin, LDRmax to file + # http://stackoverflow.com/questions/4675728/redirect-stdout-to-a-file-in-python + with open('output_files\\' + OutputFile, 'a') as f: + # with open('output_files\\' + LID + '-' + InputFile[0:-3] + '-LDR_min_max.dat', 'w') as f: + with redirect_stdout(f): + print("Lidar ID: " + LID) + print() + print("Correction K(LDRCal in calibration range)") + print("LDRCal, K(LDRCal)") + for i in range(1, 7): + # print("{0:5.3f},{1:9.5f},{2:9.5f}".format(LDRCalList[i], K0List[i]), p(LDRCalList[i])) + print("{0:5.3f},{1:9.5f}".format(LDRCalList[i], K0List[i])) + print() + print("minimum and maximum values of the distributions of possibly measured LDR for different LDRtrue") + print("LDRtrue , LDRmin, LDRmax") + for i in range(len(LDRrangeA)): + print("{0:7.4f},{1:7.4f},{2:7.4f}".format(LDRrangeA[i], LDRmin[i], LDRmax[i])) + print() + print("LDRCal is the volume linear depolarisation ratio in the range used for the polarisation calibration. K(LDRCal) is the variability of K depending on LDRCal.") + print("If LDRCal is not known, its uncertainty has to be included in the input file for the error calcualtion below.") + print("LDR_min-true(LDR_true) and LDR_max-true(LDR_true) are the difference between the lowest and highest possible retrieved LDR, respectively, and the true LDR depending on the true LDR. They represent the negative and positive uncertainty of the retrieved LDR resulting from all known uncertainties of the optical input parameters.") + print("-----------------------------------------------") + print("2nd order polyfit = a + (b * LDR) + (c * LDR^2)") + print("fit-parameter, K(LDRCal), LDR_min-true(LDRtrue), LDR_max-true(LDRtrue)") + print(" a, ", f"{K_poly[2]:+.5f},", f"{LDRmin_poly[2]:+.5f},", f"{LDRmax_poly[2]:+.5f}") + print(" b, ", f"{K_poly[1]:+.5f},", f"{LDRmin_poly[1]:+.5f},", f"{LDRmax_poly[1]:+.5f}") + print(" c, ", f"{K_poly[0]:+.5f},", f"{LDRmin_poly[0]:+.5f},", f"{LDRmax_poly[0]:+.5f}") + print() + # Print LDR statistics + print("##################################################################################") + print("Further statistical analyses of the uncertainties, their distribution and sources:") + print() + print("LDRtrue , mean , median, max-mean, min-mean, std, excess_kurtosis, skewness") + iLDR = -1 + LDRrangeA = np.array(LDRrange) + for LDRTrue in LDRrange: + iLDR = iLDR + 1 + print("{0:8.5f},{1:8.5f},{2:8.5f}, {3:8.5f},{4:8.5f},{5:8.5f}, {6:8.5f},{7:8.5f}"\ + .format(LDRrangeA[iLDR], LDRmean[iLDR], LDRmedian[iLDR], LDRmax[iLDR]-LDRrangeA[iLDR], \ + LDRmin[iLDR]-LDRrangeA[iLDR], LDRstd[iLDR], LDRkurt[iLDR], LDRskew[iLDR])) + print() + # Calculate and print statistics for calibration factors + print("minimum and maximum values of the distributions of signal ratios and calibration factors for different LDRtrue") + iLDR = -1 + LDRrangeA = np.array(LDRrange) + print("LDRtrue , LDRsim, (max-min)/2, relerr") + for LDRTrue in LDRrange: + iLDR = iLDR + 1 + LDRsimmin[iLDR] = np.amin(aLDRsim[iLDR, :]) + LDRsimmax[iLDR] = np.amax(aLDRsim[iLDR, :]) + # LDRsimstd = np.std(aLDRsim[iLDR, :]) + LDRsimmean[iLDR] = np.mean(aLDRsim[iLDR, :]) + # LDRsimmedian = np.median(aLDRsim[iLDR, :]) + print("{0:8.5f}, {1:8.5f}, {2:8.5f}, {3:8.5f}".format(LDRrangeA[iLDR],LDRsimmean[iLDR],(LDRsimmax[iLDR]-LDRsimmin[iLDR])/2,(LDRsimmax[iLDR]-LDRsimmin[iLDR])/2/LDRsimmean[iLDR])) + iLDR = -1 + print("LDRtrue , Etax , (max-min)/2, relerr") + for LDRTrue in LDRrange: + iLDR = iLDR + 1 + Etaxmin = np.amin(aEtax[iLDR, :]) + Etaxmax = np.amax(aEtax[iLDR, :]) + # Etaxstd = np.std(aEtax[iLDR, :]) + Etaxmean = np.mean(aEtax[iLDR, :]) + # Etaxmedian = np.median(aEtax[iLDR, :]) + print("{0:8.5f}, {1:8.5f}, {2:8.5f}, {3:8.5f}".format(LDRrangeA[iLDR], Etaxmean, (Etaxmax-Etaxmin)/2, (Etaxmax-Etaxmin)/2/Etaxmean)) + iLDR = -1 + print("LDRtrue , Etapx , (max-min)/2, relerr") + for LDRTrue in LDRrange: + iLDR = iLDR + 1 + Etapxmin = np.amin(aEtapx[iLDR, :]) + Etapxmax = np.amax(aEtapx[iLDR, :]) + # Etapxstd = np.std(aEtapx[iLDR, :]) + Etapxmean = np.mean(aEtapx[iLDR, :]) + # Etapxmedian = np.median(aEtapx[iLDR, :]) + print("{0:8.5f}, {1:8.5f}, {2:8.5f}, {3:8.5f}".format(LDRrangeA[iLDR], Etapxmean, (Etapxmax-Etapxmin)/2, (Etapxmax-Etapxmin)/2/Etapxmean)) + iLDR = -1 + print("LDRtrue , Etamx , (max-min)/2, relerr") + for LDRTrue in LDRrange: + iLDR = iLDR + 1 + Etamxmin = np.amin(aEtamx[iLDR, :]) + Etamxmax = np.amax(aEtamx[iLDR, :]) + # Etamxstd = np.std(aEtamx[iLDR, :]) + Etamxmean = np.mean(aEtamx[iLDR, :]) + # Etamxmedian = np.median(aEtamx[iLDR, :]) + print("{0:8.5f}, {1:8.5f}, {2:8.5f}, {3:8.5f}".format(LDRrangeA[iLDR], Etamxmean, (Etamxmax-Etamxmin)/2, (Etamxmax-Etamxmin)/2/Etamxmean)) + + with open('output_files\\' + OutputFile, 'a') as f: + # with open('output_files\\' + LID + '-' + InputFile[0:-3] + '-LDR_min_max.dat', 'a') as f: + with redirect_stdout(f): + Contribution_aY("LDRCorr", aLDRcorr) + Contribution_aY("LDRsim", aLDRsim) + Contribution_aY("EtaX, D90", aEtax) + Contribution_aY("Etapx, +45°", aEtapx) + Contribution_aY("Etamx -45°", aEtamx) + + + # Plot other histograms + if (bPlotEtax): + + if (nQin > 0): PlotLDRsim("Qin", aQin, Qin0, dQin, iQin, nQin) + if (nVin > 0): PlotLDRsim("Vin", aVin, Vin0, dVin, iVin, nVin) + if (nRotL > 0): PlotLDRsim("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) + if (nRetE > 0): PlotLDRsim("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) + if (nRotE > 0): PlotLDRsim("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) + if (nDiE > 0): PlotLDRsim("DiE", aDiE, DiE0, dDiE, iDiE, nDiE) + if (nRetO > 0): PlotLDRsim("RetO", aRetO, RetO0, dRetO, iRetO, nRetO) + if (nRotO > 0): PlotLDRsim("RotO", aRotO, RotO0, dRotO, iRotO, nRotO) + if (nDiO > 0): PlotLDRsim("DiO", aDiO, DiO0, dDiO, iDiO, nDiO) + if (nDiC > 0): PlotLDRsim("DiC", aDiC, DiC0, dDiC, iDiC, nDiC) + if (nRotC > 0): PlotLDRsim("RotC", aRotC, RotC0, dRotC, iRotC, nRotC) + if (nRetC > 0): PlotLDRsim("RetC", aRetC, RetC0, dRetC, iRetC, nRetC) + if (nTP > 0): PlotLDRsim("TP", aTP, TP0, dTP, iTP, nTP) + if (nTS > 0): PlotLDRsim("TS", aTS, TS0, dTS, iTS, nTS) + if (nRP > 0): PlotLDRsim("RP", aRP, RP0, dRP, iRP, nRP) + if (nRS > 0): PlotLDRsim("RS", aRS, RS0, dRS, iRS, nRS) + if (nRetT > 0): PlotLDRsim("RetT", aRetT, RetT0, dRetT, iRetT, nRetT) + if (nRetR > 0): PlotLDRsim("RetR", aRetR, RetR0, dRetR, iRetR, nRetR) + if (nERaT > 0): PlotLDRsim("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) + if (nERaR > 0): PlotLDRsim("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) + if (nRotaT > 0): PlotLDRsim("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) + if (nRotaR > 0): PlotLDRsim("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) + if (nLDRCal > 0): PlotLDRsim("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) + if (nTCalT > 0): PlotLDRsim("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT) + if (nTCalR > 0): PlotLDRsim("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR) + if (nNCal > 0): PlotLDRsim("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal) + if (nNCal > 0): PlotLDRsim("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal) + if (nNCal > 0): PlotLDRsim("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal) + if (nNCal > 0): PlotLDRsim("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal) + if (nNI > 0): PlotLDRsim("SigNoiseIt", aNIt, 0, 1, iNIt, nNI) + if (nNI > 0): PlotLDRsim("SigNoiseIr", aNIr, 0, 1, iNIr, nNI) + plt.show() + plt.close + print("---------------------------------------...producing more plots...------------------------------------------------------------------") + + if (nQin > 0): PlotEtax("Qin", aQin, Qin0, dQin, iQin, nQin) + if (nVin > 0): PlotEtax("Vin", aVin, Vin0, dVin, iVin, nVin) + if (nRotL > 0): PlotEtax("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) + if (nRetE > 0): PlotEtax("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) + if (nRotE > 0): PlotEtax("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) + if (nDiE > 0): PlotEtax("DiE", aDiE, DiE0, dDiE, iDiE, nDiE) + if (nRetO > 0): PlotEtax("RetO", aRetO, RetO0, dRetO, iRetO, nRetO) + if (nRotO > 0): PlotEtax("RotO", aRotO, RotO0, dRotO, iRotO, nRotO) + if (nDiO > 0): PlotEtax("DiO", aDiO, DiO0, dDiO, iDiO, nDiO) + if (nDiC > 0): PlotEtax("DiC", aDiC, DiC0, dDiC, iDiC, nDiC) + if (nRotC > 0): PlotEtax("RotC", aRotC, RotC0, dRotC, iRotC, nRotC) + if (nRetC > 0): PlotEtax("RetC", aRetC, RetC0, dRetC, iRetC, nRetC) + if (nTP > 0): PlotEtax("TP", aTP, TP0, dTP, iTP, nTP) + if (nTS > 0): PlotEtax("TS", aTS, TS0, dTS, iTS, nTS) + if (nRP > 0): PlotEtax("RP", aRP, RP0, dRP, iRP, nRP) + if (nRS > 0): PlotEtax("RS", aRS, RS0, dRS, iRS, nRS) + if (nRetT > 0): PlotEtax("RetT", aRetT, RetT0, dRetT, iRetT, nRetT) + if (nRetR > 0): PlotEtax("RetR", aRetR, RetR0, dRetR, iRetR, nRetR) + if (nERaT > 0): PlotEtax("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) + if (nERaR > 0): PlotEtax("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) + if (nRotaT > 0): PlotEtax("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) + if (nRotaR > 0): PlotEtax("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) + if (nLDRCal > 0): PlotEtax("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) + if (nTCalT > 0): PlotEtax("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT) + if (nTCalR > 0): PlotEtax("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR) + if (nNCal > 0): PlotEtax("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal) + if (nNCal > 0): PlotEtax("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal) + if (nNCal > 0): PlotEtax("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal) + if (nNCal > 0): PlotEtax("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal) + if (nNI > 0): PlotEtax("SigNoiseIt", aNIt, 0, 1, iNIt, nNI) + if (nNI > 0): PlotEtax("SigNoiseIr", aNIr, 0, 1, iNIr, nNI) + plt.show() + plt.close + print("---------------------------------------...producing more plots...------------------------------------------------------------------") + + if (nQin > 0): PlotEtapx("Qin", aQin, Qin0, dQin, iQin, nQin) + if (nVin > 0): PlotEtapx("Vin", aVin, Vin0, dVin, iVin, nVin) + if (nRotL > 0): PlotEtapx("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) + if (nRetE > 0): PlotEtapx("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) + if (nRotE > 0): PlotEtapx("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) + if (nDiE > 0): PlotEtapx("DiE", aDiE, DiE0, dDiE, iDiE, nDiE) + if (nRetO > 0): PlotEtapx("RetO", aRetO, RetO0, dRetO, iRetO, nRetO) + if (nRotO > 0): PlotEtapx("RotO", aRotO, RotO0, dRotO, iRotO, nRotO) + if (nDiO > 0): PlotEtapx("DiO", aDiO, DiO0, dDiO, iDiO, nDiO) + if (nDiC > 0): PlotEtapx("DiC", aDiC, DiC0, dDiC, iDiC, nDiC) + if (nRotC > 0): PlotEtapx("RotC", aRotC, RotC0, dRotC, iRotC, nRotC) + if (nRetC > 0): PlotEtapx("RetC", aRetC, RetC0, dRetC, iRetC, nRetC) + if (nTP > 0): PlotEtapx("TP", aTP, TP0, dTP, iTP, nTP) + if (nTS > 0): PlotEtapx("TS", aTS, TS0, dTS, iTS, nTS) + if (nRP > 0): PlotEtapx("RP", aRP, RP0, dRP, iRP, nRP) + if (nRS > 0): PlotEtapx("RS", aRS, RS0, dRS, iRS, nRS) + if (nRetT > 0): PlotEtapx("RetT", aRetT, RetT0, dRetT, iRetT, nRetT) + if (nRetR > 0): PlotEtapx("RetR", aRetR, RetR0, dRetR, iRetR, nRetR) + if (nERaT > 0): PlotEtapx("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) + if (nERaR > 0): PlotEtapx("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) + if (nRotaT > 0): PlotEtapx("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) + if (nRotaR > 0): PlotEtapx("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) + if (nLDRCal > 0): PlotEtapx("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) + if (nTCalT > 0): PlotEtapx("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT) + if (nTCalR > 0): PlotEtapx("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR) + if (nNCal > 0): PlotEtapx("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal) + if (nNCal > 0): PlotEtapx("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal) + if (nNCal > 0): PlotEtapx("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal) + if (nNCal > 0): PlotEtapx("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal) + if (nNI > 0): PlotEtapx("SigNoiseIt", aNIt, 0, 1, iNIt, nNI) + if (nNI > 0): PlotEtapx("SigNoiseIr", aNIr, 0, 1, iNIr, nNI) + plt.show() + plt.close + print("---------------------------------------...producing more plots...------------------------------------------------------------------") + + if (nQin > 0): PlotEtamx("Qin", aQin, Qin0, dQin, iQin, nQin) + if (nVin > 0): PlotEtamx("Vin", aVin, Vin0, dVin, iVin, nVin) + if (nRotL > 0): PlotEtamx("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) + if (nRetE > 0): PlotEtamx("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) + if (nRotE > 0): PlotEtamx("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) + if (nDiE > 0): PlotEtamx("DiE", aDiE, DiE0, dDiE, iDiE, nDiE) + if (nRetO > 0): PlotEtamx("RetO", aRetO, RetO0, dRetO, iRetO, nRetO) + if (nRotO > 0): PlotEtamx("RotO", aRotO, RotO0, dRotO, iRotO, nRotO) + if (nDiO > 0): PlotEtamx("DiO", aDiO, DiO0, dDiO, iDiO, nDiO) + if (nDiC > 0): PlotEtamx("DiC", aDiC, DiC0, dDiC, iDiC, nDiC) + if (nRotC > 0): PlotEtamx("RotC", aRotC, RotC0, dRotC, iRotC, nRotC) + if (nRetC > 0): PlotEtamx("RetC", aRetC, RetC0, dRetC, iRetC, nRetC) + if (nTP > 0): PlotEtamx("TP", aTP, TP0, dTP, iTP, nTP) + if (nTS > 0): PlotEtamx("TS", aTS, TS0, dTS, iTS, nTS) + if (nRP > 0): PlotEtamx("RP", aRP, RP0, dRP, iRP, nRP) + if (nRS > 0): PlotEtamx("RS", aRS, RS0, dRS, iRS, nRS) + if (nRetT > 0): PlotEtamx("RetT", aRetT, RetT0, dRetT, iRetT, nRetT) + if (nRetR > 0): PlotEtamx("RetR", aRetR, RetR0, dRetR, iRetR, nRetR) + if (nERaT > 0): PlotEtamx("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) + if (nERaR > 0): PlotEtamx("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) + if (nRotaT > 0): PlotEtamx("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) + if (nRotaR > 0): PlotEtamx("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) + if (nLDRCal > 0): PlotEtamx("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) + if (nTCalT > 0): PlotEtamx("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT) + if (nTCalR > 0): PlotEtamx("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR) + if (nNCal > 0): PlotEtamx("CalNoiseTp", aNCalTp, 0, 1, iNCalTp, nNCal) + if (nNCal > 0): PlotEtamx("CalNoiseTm", aNCalTm, 0, 1, iNCalTm, nNCal) + if (nNCal > 0): PlotEtamx("CalNoiseRp", aNCalRp, 0, 1, iNCalRp, nNCal) + if (nNCal > 0): PlotEtamx("CalNoiseRm", aNCalRm, 0, 1, iNCalRm, nNCal) + if (nNI > 0): PlotEtamx("SigNoiseIt", aNIt, 0, 1, iNIt, nNI) + if (nNI > 0): PlotEtamx("SigNoiseIr", aNIr, 0, 1, iNIr, nNI) + plt.show() + plt.close + + # Print Etax statistics + Etaxmin = np.amin(aEtax[1, :]) + Etaxmax = np.amax(aEtax[1, :]) + Etaxstd = np.std(aEtax[1, :]) + Etaxmean = np.mean(aEtax[1, :]) + Etaxmedian = np.median(aEtax[1, :]) + print("Etax , max-mean, min-mean, median, mean ± std, eta") + print("{0:8.5f} ±({1:8.5f},{2:8.5f}),{3:8.5f},{4:8.5f}±{5:8.5f},{6:8.5f}".format(Etax0, Etaxmax-Etax0, Etaxmin-Etax0, Etaxmedian, Etaxmean, Etaxstd, Etax0 / K0)) + print() + + # Calculate and print statistics for calibration factors + iLDR = -1 + LDRrangeA = np.array(LDRrange) + print("LDR...., LDRsim, (max-min)/2, relerr") + for LDRTrue in LDRrange: + iLDR = iLDR + 1 + LDRsimmin[iLDR] = np.amin(aLDRsim[iLDR, :]) + LDRsimmax[iLDR] = np.amax(aLDRsim[iLDR, :]) + # LDRsimstd = np.std(aLDRsim[iLDR, :]) + LDRsimmean[iLDR] = np.mean(aLDRsim[iLDR, :]) + # LDRsimmedian = np.median(aLDRsim[iLDR, :]) + print("{0:8.5f}, {1:8.5f}, {2:8.5f}, {3:8.5f}".format(LDRrangeA[iLDR], LDRsimmean[iLDR], (LDRsimmax[iLDR]-LDRsimmin[iLDR])/2, (LDRsimmax[iLDR]-LDRsimmin[iLDR])/2/LDRsimmean[iLDR])) + iLDR = -1 + print("LDR...., Etax , (max-min)/2, relerr") + for LDRTrue in LDRrange: + iLDR = iLDR + 1 + Etaxmin = np.amin(aEtax[iLDR, :]) + Etaxmax = np.amax(aEtax[iLDR, :]) + # Etaxstd = np.std(aEtax[iLDR, :]) + Etaxmean = np.mean(aEtax[iLDR, :]) + # Etaxmedian = np.median(aEtax[iLDR, :]) + print("{0:8.5f}, {1:8.5f}, {2:8.5f}, {3:8.5f}".format(LDRrangeA[iLDR], Etaxmean, (Etaxmax-Etaxmin)/2, (Etaxmax-Etaxmin)/2/Etaxmean)) + iLDR = -1 + print("LDR...., Etapx , (max-min)/2, relerr") + for LDRTrue in LDRrange: + iLDR = iLDR + 1 + Etapxmin = np.amin(aEtapx[iLDR, :]) + Etapxmax = np.amax(aEtapx[iLDR, :]) + # Etapxstd = np.std(aEtapx[iLDR, :]) + Etapxmean = np.mean(aEtapx[iLDR, :]) + # Etapxmedian = np.median(aEtapx[iLDR, :]) + print("{0:8.5f}, {1:8.5f}, {2:8.5f}, {3:8.5f}".format(LDRrangeA[iLDR], Etapxmean, (Etapxmax-Etapxmin)/2, (Etapxmax-Etapxmin)/2/Etapxmean)) + iLDR = -1 + print("LDR...., Etamx , (max-min)/2, relerr") + for LDRTrue in LDRrange: + iLDR = iLDR + 1 + Etamxmin = np.amin(aEtamx[iLDR, :]) + Etamxmax = np.amax(aEtamx[iLDR, :]) + # Etamxstd = np.std(aEtamx[iLDR, :]) + Etamxmean = np.mean(aEtamx[iLDR, :]) + # Etamxmedian = np.median(aEtamx[iLDR, :]) + print("{0:8.5f}, {1:8.5f}, {2:8.5f}, {3:8.5f}".format(LDRrangeA[iLDR], Etamxmean, (Etamxmax-Etamxmin)/2, (Etamxmax-Etamxmin)/2/Etamxmean)) + + f.close() + + +''' + # --- Plot F11 histograms + print() + print(" ############################################################################## ") + print(Text1) + print() + + iLDR = 5 + for LDRTrue in LDRrange: + iLDR = iLDR - 1 + #aF11corr[iLDR,:] = aF11corr[iLDR,:] / aF11corr[0,:] - 1.0 + aF11corr[iLDR,:] = aF11corr[iLDR,:] / aF11sim0[iLDR] - 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(aLDRcorr[iLDR,:]) * 0.995 + #Rmax = F11max[iLDR] * 1.005 # np.max(aLDRcorr[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, density=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, density=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 (nQin > 0): PlotSubHistF11("Qin", aQin, Qin0, dQin, iQin, nQin) + if (nVin > 0): PlotSubHistF11("Vin", aVin, Vin0, dVin, iVin, nVin) + 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) + if (nTCalT > 0): PlotSubHistF11("TCalT", aTCalT, TCalT0, dTCalT, iTCalT, nTCalT) + if (nTCalR > 0): PlotSubHistF11("TCalR", aTCalR, TCalR0, dTCalR, iTCalR, nTCalR) + if (nNCal > 0): PlotSubHistF11("CalNoise", aNCal, 0, 1/nNCal, iNCal, nNCal) + if (nNI > 0): PlotSubHistF11("SigNoise", aNI, 0, 1/nNI, iNI, nNI) + + + 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(aLDRcorr[iLDR,:]) + LDRmax[iLDR] = np.max(aLDRcorr[iLDR,:]) + Rmin = np.min(aLDRcorr[iLDR,:]) * 0.999 + Rmax = np.max(aLDRcorr[iLDR,:]) * 1.001 + plt.subplot(5,2,iLDR+1) + (n, bins, patches) = plt.hist(aLDRcorr[iLDR,:], + range=[Rmin, Rmax], + bins=200, log=False, alpha=0.2, density=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 + # --- End of Plot F11 histograms + ''' + + +''' + # --- Plot K over LDRCal + fig3 = plt.figure() + plt.plot(LDRCal0+aLDRCal*dLDRCal/nLDRCal,aGHK[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') +# +''' \ No newline at end of file diff -r a882c113284f -r 7c14c9a085ba output_files/output_optic_input_0.9.8h-PollyXT_Cyprus_532_210429_vf_4_GHK_0.9.8h_Py3.7.dat --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/output_files/output_optic_input_0.9.8h-PollyXT_Cyprus_532_210429_vf_4_GHK_0.9.8h_Py3.7.dat Wed Jul 07 15:42:30 2021 +0200 @@ -0,0 +1,197 @@ +From C:\Projects\Earlinet-Asos\Lidarsystems\calc_lidar_correction_parameters +Running GHK_0.9.8h_Py3.7.py +Reading input file optic_input_0.9.8h-PollyXT_Cyprus_532_210429_vf_4.py +for Lidar system: li , PollyXT Cyprus 210429 532 approx 4 + --- Input parameters: value ±error / ±steps ---------------------- +Laser: Qin = 0.9672± 0.0100/ 1 + Vin = 0.0000± 0.0000/ 0 + Rotation alpha = 91.6500± 0.2400/ 1 + => DOP 0.9672 (degree of polarisation) +Optic: Diatt., Tunpol, Retard., Rotation (deg) +Emitter 0.0000 ± 0.0200 / 0, 1.0000, 0±180/ 0, 0.0000± 1.0000/ 0 +Receiver 0.0000 ± 0.0200 / 0, 1.0000, 0±180/ 0, 0.0000± 0.5000/ 0 +Calibrator 0.999280± 0.000210/ 1, 0.4000, 0± 0/ 0, 0.0000± 0.1000/ 0 + Pol.-filter ------ +ERT, RotT : 1.0000± 0.0000/ 0, 0.0000± 1.0000/ 1 +ERR, RotR : 1.0000± 0.0000/ 0, 0.0000± 1.0000/ 0 + PBS ------ +TP,TS : 0.5000± 0.0100/ 1, 0.5000± 0.0100/ 1 +RP,RS : 0.0020± 0.0010/ 1, 0.9980± 0.0010/ 1 +DT,TT,DR,TR,Y : 0.0000, 0.5000, -0.9960, 0.5000, -1 + Combined PBS + Pol.-filter ------ +DT,TT,DR,TR : 0.0000, 0.5000, -0.9960, 0.5000 +LDRCal during calibration in calibration range: 0.110± 0.100/ 1 + --- Additional ND filter attenuation (transmission) during the calibration --- +TCalT,TCalR : 1.0000± 0.0100/ 0, 1.0000± 0.0001/ 0 + +Linear polarizer calibrator is located before the receiver +Rotation error epsilon is considered also for normal measurements = False +The PBS incidence plane is perpendicular to the reference plane +The laser polarisation in the reference plane is finally reflected => the parallel laser polarisation is detected in the transmitted channel. +RS_RP_depend_on_TS_TP = False + +=========================================================================================================== + GR , GT , HR , HT , K(0.110), K(0.004), K(0.05) , K(0.1) , K(0.2) , K(0.3) , K(0.45) + 1.00000, 1.00000,-0.96173, 0.00000, 0.97068 , 0.96369, 0.96690, 0.97008, 0.97564, 0.98033, 0.98615 +=========================================================================================================== + +2nd order polyfit: K(LDRcal) = a + (b * LDRcal) + (c * LDRcal^2) +(LDRcal is the volume linear depolarisation ratio in the range used for the polarisation calibration) + a = 0.96352 + b = 0.06904 + c = -0.04188 + +Errors from neglecting GHK corrections and/or calibration: + LDRtrue, LDRunCorr,1/LDRunCorr, LDRsimx, 1/LDRsimx, LDRCorr + 0.00400, 0.04732, 21.13427, 0.04593, 21.77255, 0.00400 + 0.02000, 0.07828, 12.77531, 0.07598, 13.16114, 0.02000 + 0.10000, 0.21956, 4.55449, 0.21313, 4.69204, 0.10000 + 0.30000, 0.49670, 2.01327, 0.48214, 2.07407, 0.30000 + 0.45000, 0.65439, 1.52814, 0.63520, 1.57430, 0.45000 +LDRsimx = LDR of the nominal system directly from measured signals without calibration and GHK-corrections +LDRunCorr = LDR of the nominal system directly from measured signals with calibration but without GHK-corrections; electronic amplifications = 1 assumed +LDRCorr = LDR calibrated and GHK-corrected + +Errors from signal noise: +Signal counts: NI, NCalT, NCalR, NILfac, nNCal, nNI, stdev(NI)/NI = 100000000, 100000000, 100000000, 1, 0, 0, 0.00010 + + +Lidar ID: PollyXT Cyprus 210429 532 approx 4 + +Correction K(LDRCal in calibration range) +LDRCal, K(LDRCal) +0.004, 0.96369 +0.050, 0.96690 +0.100, 0.97008 +0.200, 0.97564 +0.300, 0.98033 +0.450, 0.98615 + +minimum and maximum values of the distributions of possibly measured LDR for different LDRtrue +LDRtrue , LDRmin, LDRmax + 0.0040,-0.0029, 0.0114 + 0.0200, 0.0127, 0.0279 + 0.1000, 0.0906, 0.1102 + 0.3000, 0.2863, 0.3151 + 0.4500, 0.4340, 0.4679 + +LDRCal is the volume linear depolarisation ratio in the range used for the polarisation calibration. K(LDRCal) is the variability of K depending on LDRCal. +If LDRCal is not known, its uncertainty has to be included in the input file for the error calcualtion below. +LDR_min-true(LDR_true) and LDR_max-true(LDR_true) are the difference between the lowest and highest possible retrieved LDR, respectively, and the true LDR depending on the true LDR. They represent the negative and positive uncertainty of the retrieved LDR resulting from all known uncertainties of the optical input parameters. +----------------------------------------------- +2nd order polyfit = a + (b * LDR) + (c * LDR^2) +fit-parameter, K(LDRCal), LDR_min-true(LDRtrue), LDR_max-true(LDRtrue) + a, +0.96352, -0.00674, +0.00727 + b, +0.06904, -0.02805, +0.03080 + c, -0.04188, +0.01659, -0.01587 + +################################################################################## +Further statistical analyses of the uncertainties, their distribution and sources: + +LDRtrue , mean , median, max-mean, min-mean, std, excess_kurtosis, skewness + 0.00400, 0.00404, 0.00402, 0.00739,-0.00685, 0.00431, -1.37505, 0.01060 + 0.02000, 0.02004, 0.02003, 0.00788,-0.00729, 0.00432, -1.35168, 0.01361 + 0.10000, 0.10006, 0.10001, 0.01019,-0.00938, 0.00447, -1.11574, 0.02519 + 0.30000, 0.30013, 0.30010, 0.01508,-0.01366, 0.00523, -0.47128, 0.06082 + 0.45000, 0.45018, 0.45000, 0.01792,-0.01600, 0.00607, -0.34137, 0.12718 + +minimum and maximum values of the distributions of signal ratios and calibration factors for different LDRtrue +LDRtrue , LDRsim, (max-min)/2, relerr + 0.00400, 0.04596, 0.01325, 0.28826 + 0.02000, 0.07602, 0.01353, 0.17791 + 0.10000, 0.21319, 0.01478, 0.06935 + 0.30000, 0.48224, 0.01725, 0.03578 + 0.45000, 0.63531, 0.01866, 0.02937 +LDRtrue , Etax , (max-min)/2, relerr + 0.00400, 0.97066, 0.03207, 0.03304 + 0.02000, 0.97066, 0.03207, 0.03304 + 0.10000, 0.97066, 0.03207, 0.03304 + 0.30000, 0.97066, 0.03207, 0.03304 + 0.45000, 0.97066, 0.03207, 0.03304 +LDRtrue , Etapx , (max-min)/2, relerr + 0.00400, 0.97203, 0.03152, 0.03243 + 0.02000, 0.97203, 0.03152, 0.03243 + 0.10000, 0.97203, 0.03152, 0.03243 + 0.30000, 0.97203, 0.03152, 0.03243 + 0.45000, 0.97203, 0.03152, 0.03243 +LDRtrue , Etamx , (max-min)/2, relerr + 0.00400, 0.96930, 0.03303, 0.03408 + 0.02000, 0.96930, 0.03303, 0.03408 + 0.10000, 0.96930, 0.03303, 0.03408 + 0.30000, 0.96930, 0.03303, 0.03408 + 0.45000, 0.96930, 0.03303, 0.03408 + +LDRCorr: contribution to the total error (per mill) + of individual parameter errors of combined parameter errors + at LDRtrue 0.004 0.020 0.100 0.300 0.450 0.004 0.020 0.100 0.300 0.450 + +Qin 725 680 519 318 230 21 19 13 4 0 +RotL 34 32 25 15 12 1 1 1 1 1 +DiC 15 24 62 133 185 7 8 15 28 37 +TP 31 49 104 145 138 9 9 7 3 0 +TS 31 49 104 145 138 9 9 7 3 0 +RP 141 132 101 62 45 4 4 3 1 0 +RS 0 0 0 0 0 0 0 0 0 0 +RotaT 0 0 0 0 0 0 0 0 0 0 +LDRCal 21 33 84 182 252 10 12 19 34 44 +Sum 999 999 999 999 1000 61 62 64 74 83 + +LDRsim: contribution to the total error (per mill) + of individual parameter errors of combined parameter errors + at LDRtrue 0.004 0.020 0.100 0.300 0.450 0.004 0.020 0.100 0.300 0.450 + +Qin 745 706 550 310 202 15 14 11 6 4 +RotL 35 33 26 14 9 0 0 0 0 0 +DiC 0 0 0 0 0 0 0 0 0 0 +TP 68 108 258 425 465 18 17 13 9 9 +TS 2 4 30 134 216 1 1 4 8 9 +RP 148 143 121 88 73 2 2 2 1 1 +RS 3 5 14 28 34 1 1 1 1 1 +RotaT 0 0 0 0 0 0 0 0 0 0 +LDRCal 0 0 0 0 0 0 0 0 0 0 +Sum 1000 1000 1000 1000 1000 37 35 30 26 25 + +EtaX, D90: contribution to the total error (per mill) + of individual parameter errors of combined parameter errors + at LDRtrue 0.004 0.020 0.100 0.300 0.450 0.004 0.020 0.100 0.300 0.450 + +Qin 10 10 10 10 10 3 3 3 3 3 +RotL 0 0 0 0 0 0 0 0 0 0 +DiC 136 136 136 136 136 24 24 24 24 24 +TP 312 312 312 312 312 7 7 7 7 7 +TS 294 294 294 294 294 13 13 13 13 13 +RP 32 32 32 32 32 0 0 0 0 0 +RS 30 30 30 30 30 1 1 1 1 1 +RotaT 0 0 0 0 0 0 0 0 0 0 +LDRCal 186 186 186 186 186 29 29 29 29 29 +Sum 1000 1000 1000 1000 1000 78 78 78 78 78 + +Etapx, +45°: contribution to the total error (per mill) + of individual parameter errors of combined parameter errors + at LDRtrue 0.004 0.020 0.100 0.300 0.450 0.004 0.020 0.100 0.300 0.450 + +Qin 9 9 9 9 9 3 3 3 3 3 +RotL 6 6 6 6 6 3 3 3 3 3 +DiC 132 132 132 132 132 23 23 23 23 23 +TP 317 317 317 317 317 7 7 7 7 7 +TS 300 300 300 300 300 13 13 13 13 13 +RP 33 33 33 33 33 0 0 0 0 0 +RS 31 31 31 31 31 1 1 1 1 1 +RotaT 0 0 0 0 0 0 0 0 0 0 +LDRCal 172 172 172 172 172 28 28 28 28 28 +Sum 1000 1000 1000 1000 1000 79 79 79 79 79 + +Etamx -45°: contribution to the total error (per mill) + of individual parameter errors of combined parameter errors + at LDRtrue 0.004 0.020 0.100 0.300 0.450 0.004 0.020 0.100 0.300 0.450 + +Qin 10 10 10 10 10 4 4 4 4 4 +RotL 6 6 6 6 6 4 4 4 4 4 +DiC 139 139 139 139 139 27 27 27 27 27 +TP 303 303 303 303 303 7 7 7 7 7 +TS 284 284 284 284 284 13 13 13 13 13 +RP 31 31 31 31 31 0 0 0 0 0 +RS 29 29 29 29 29 1 1 1 1 1 +RotaT 0 0 0 0 0 0 0 0 0 0 +LDRCal 197 197 197 197 197 33 33 33 33 33 +Sum 999 999 999 999 999 89 89 89 89 89 diff -r a882c113284f -r 7c14c9a085ba system_settings/output_optic_input_0.9.8h-PollyXT_Cyprus_532_210429_vf_4_GHK_0.9.8h_Py3.7.dat --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/system_settings/output_optic_input_0.9.8h-PollyXT_Cyprus_532_210429_vf_4_GHK_0.9.8h_Py3.7.dat Wed Jul 07 15:42:30 2021 +0200 @@ -0,0 +1,197 @@ +From C:\Projects\Earlinet-Asos\Lidarsystems\calc_lidar_correction_parameters +Running GHK_0.9.8h_Py3.7.py +Reading input file optic_input_0.9.8h-PollyXT_Cyprus_532_210429_vf_4.py +for Lidar system: li , PollyXT Cyprus 210429 532 approx 4 + --- Input parameters: value ±error / ±steps ---------------------- +Laser: Qin = 0.9672± 0.0100/ 1 + Vin = 0.0000± 0.0000/ 0 + Rotation alpha = 91.6500± 0.2400/ 1 + => DOP 0.9672 (degree of polarisation) +Optic: Diatt., Tunpol, Retard., Rotation (deg) +Emitter 0.0000 ± 0.0200 / 0, 1.0000, 0±180/ 0, 0.0000± 1.0000/ 0 +Receiver 0.0000 ± 0.0200 / 0, 1.0000, 0±180/ 0, 0.0000± 0.5000/ 0 +Calibrator 0.999280± 0.000210/ 1, 0.4000, 0± 0/ 0, 0.0000± 0.1000/ 0 + Pol.-filter ------ +ERT, RotT : 1.0000± 0.0000/ 0, 0.0000± 1.0000/ 1 +ERR, RotR : 1.0000± 0.0000/ 0, 0.0000± 1.0000/ 0 + PBS ------ +TP,TS : 0.5000± 0.0100/ 1, 0.5000± 0.0100/ 1 +RP,RS : 0.0020± 0.0010/ 1, 0.9980± 0.0010/ 1 +DT,TT,DR,TR,Y : 0.0000, 0.5000, -0.9960, 0.5000, -1 + Combined PBS + Pol.-filter ------ +DT,TT,DR,TR : 0.0000, 0.5000, -0.9960, 0.5000 +LDRCal during calibration in calibration range: 0.110± 0.100/ 1 + --- Additional ND filter attenuation (transmission) during the calibration --- +TCalT,TCalR : 1.0000± 0.0100/ 0, 1.0000± 0.0001/ 0 + +Linear polarizer calibrator is located before the receiver +Rotation error epsilon is considered also for normal measurements = False +The PBS incidence plane is perpendicular to the reference plane +The laser polarisation in the reference plane is finally reflected => the parallel laser polarisation is detected in the transmitted channel. +RS_RP_depend_on_TS_TP = False + +=========================================================================================================== + GR , GT , HR , HT , K(0.110), K(0.004), K(0.05) , K(0.1) , K(0.2) , K(0.3) , K(0.45) + 1.00000, 1.00000,-0.96173, 0.00000, 0.97068 , 0.96369, 0.96690, 0.97008, 0.97564, 0.98033, 0.98615 +=========================================================================================================== + +2nd order polyfit: K(LDRcal) = a + (b * LDRcal) + (c * LDRcal^2) +(LDRcal is the volume linear depolarisation ratio in the range used for the polarisation calibration) + a = 0.96352 + b = 0.06904 + c = -0.04188 + +Errors from neglecting GHK corrections and/or calibration: + LDRtrue, LDRunCorr,1/LDRunCorr, LDRsimx, 1/LDRsimx, LDRCorr + 0.00400, 0.04732, 21.13427, 0.04593, 21.77255, 0.00400 + 0.02000, 0.07828, 12.77531, 0.07598, 13.16114, 0.02000 + 0.10000, 0.21956, 4.55449, 0.21313, 4.69204, 0.10000 + 0.30000, 0.49670, 2.01327, 0.48214, 2.07407, 0.30000 + 0.45000, 0.65439, 1.52814, 0.63520, 1.57430, 0.45000 +LDRsimx = LDR of the nominal system directly from measured signals without calibration and GHK-corrections +LDRunCorr = LDR of the nominal system directly from measured signals with calibration but without GHK-corrections; electronic amplifications = 1 assumed +LDRCorr = LDR calibrated and GHK-corrected + +Errors from signal noise: +Signal counts: NI, NCalT, NCalR, NILfac, nNCal, nNI, stdev(NI)/NI = 100000000, 100000000, 100000000, 1, 0, 0, 0.00010 + + +Lidar ID: PollyXT Cyprus 210429 532 approx 4 + +Correction K(LDRCal in calibration range) +LDRCal, K(LDRCal) +0.004, 0.96369 +0.050, 0.96690 +0.100, 0.97008 +0.200, 0.97564 +0.300, 0.98033 +0.450, 0.98615 + +minimum and maximum values of the distributions of possibly measured LDR for different LDRtrue +LDRtrue , LDRmin, LDRmax + 0.0040,-0.0029, 0.0114 + 0.0200, 0.0127, 0.0279 + 0.1000, 0.0906, 0.1102 + 0.3000, 0.2863, 0.3151 + 0.4500, 0.4340, 0.4679 + +LDRCal is the volume linear depolarisation ratio in the range used for the polarisation calibration. K(LDRCal) is the variability of K depending on LDRCal. +If LDRCal is not known, its uncertainty has to be included in the input file for the error calcualtion below. +LDR_min-true(LDR_true) and LDR_max-true(LDR_true) are the difference between the lowest and highest possible retrieved LDR, respectively, and the true LDR depending on the true LDR. They represent the negative and positive uncertainty of the retrieved LDR resulting from all known uncertainties of the optical input parameters. +----------------------------------------------- +2nd order polyfit = a + (b * LDR) + (c * LDR^2) +fit-parameter, K(LDRCal), LDR_min-true(LDRtrue), LDR_max-true(LDRtrue) + a, +0.96352, -0.00674, +0.00727 + b, +0.06904, -0.02805, +0.03080 + c, -0.04188, +0.01659, -0.01587 + +################################################################################## +Further statistical analyses of the uncertainties, their distribution and sources: + +LDRtrue , mean , median, max-mean, min-mean, std, excess_kurtosis, skewness + 0.00400, 0.00404, 0.00402, 0.00739,-0.00685, 0.00431, -1.37505, 0.01060 + 0.02000, 0.02004, 0.02003, 0.00788,-0.00729, 0.00432, -1.35168, 0.01361 + 0.10000, 0.10006, 0.10001, 0.01019,-0.00938, 0.00447, -1.11574, 0.02519 + 0.30000, 0.30013, 0.30010, 0.01508,-0.01366, 0.00523, -0.47128, 0.06082 + 0.45000, 0.45018, 0.45000, 0.01792,-0.01600, 0.00607, -0.34137, 0.12718 + +minimum and maximum values of the distributions of signal ratios and calibration factors for different LDRtrue +LDRtrue , LDRsim, (max-min)/2, relerr + 0.00400, 0.04596, 0.01325, 0.28826 + 0.02000, 0.07602, 0.01353, 0.17791 + 0.10000, 0.21319, 0.01478, 0.06935 + 0.30000, 0.48224, 0.01725, 0.03578 + 0.45000, 0.63531, 0.01866, 0.02937 +LDRtrue , Etax , (max-min)/2, relerr + 0.00400, 0.97066, 0.03207, 0.03304 + 0.02000, 0.97066, 0.03207, 0.03304 + 0.10000, 0.97066, 0.03207, 0.03304 + 0.30000, 0.97066, 0.03207, 0.03304 + 0.45000, 0.97066, 0.03207, 0.03304 +LDRtrue , Etapx , (max-min)/2, relerr + 0.00400, 0.97203, 0.03152, 0.03243 + 0.02000, 0.97203, 0.03152, 0.03243 + 0.10000, 0.97203, 0.03152, 0.03243 + 0.30000, 0.97203, 0.03152, 0.03243 + 0.45000, 0.97203, 0.03152, 0.03243 +LDRtrue , Etamx , (max-min)/2, relerr + 0.00400, 0.96930, 0.03303, 0.03408 + 0.02000, 0.96930, 0.03303, 0.03408 + 0.10000, 0.96930, 0.03303, 0.03408 + 0.30000, 0.96930, 0.03303, 0.03408 + 0.45000, 0.96930, 0.03303, 0.03408 + +LDRCorr: contribution to the total error (per mill) + of individual parameter errors of combined parameter errors + at LDRtrue 0.004 0.020 0.100 0.300 0.450 0.004 0.020 0.100 0.300 0.450 + +Qin 725 680 519 318 230 21 19 13 4 0 +RotL 34 32 25 15 12 1 1 1 1 1 +DiC 15 24 62 133 185 7 8 15 28 37 +TP 31 49 104 145 138 9 9 7 3 0 +TS 31 49 104 145 138 9 9 7 3 0 +RP 141 132 101 62 45 4 4 3 1 0 +RS 0 0 0 0 0 0 0 0 0 0 +RotaT 0 0 0 0 0 0 0 0 0 0 +LDRCal 21 33 84 182 252 10 12 19 34 44 +Sum 999 999 999 999 1000 61 62 64 74 83 + +LDRsim: contribution to the total error (per mill) + of individual parameter errors of combined parameter errors + at LDRtrue 0.004 0.020 0.100 0.300 0.450 0.004 0.020 0.100 0.300 0.450 + +Qin 745 706 550 310 202 15 14 11 6 4 +RotL 35 33 26 14 9 0 0 0 0 0 +DiC 0 0 0 0 0 0 0 0 0 0 +TP 68 108 258 425 465 18 17 13 9 9 +TS 2 4 30 134 216 1 1 4 8 9 +RP 148 143 121 88 73 2 2 2 1 1 +RS 3 5 14 28 34 1 1 1 1 1 +RotaT 0 0 0 0 0 0 0 0 0 0 +LDRCal 0 0 0 0 0 0 0 0 0 0 +Sum 1000 1000 1000 1000 1000 37 35 30 26 25 + +EtaX, D90: contribution to the total error (per mill) + of individual parameter errors of combined parameter errors + at LDRtrue 0.004 0.020 0.100 0.300 0.450 0.004 0.020 0.100 0.300 0.450 + +Qin 10 10 10 10 10 3 3 3 3 3 +RotL 0 0 0 0 0 0 0 0 0 0 +DiC 136 136 136 136 136 24 24 24 24 24 +TP 312 312 312 312 312 7 7 7 7 7 +TS 294 294 294 294 294 13 13 13 13 13 +RP 32 32 32 32 32 0 0 0 0 0 +RS 30 30 30 30 30 1 1 1 1 1 +RotaT 0 0 0 0 0 0 0 0 0 0 +LDRCal 186 186 186 186 186 29 29 29 29 29 +Sum 1000 1000 1000 1000 1000 78 78 78 78 78 + +Etapx, +45°: contribution to the total error (per mill) + of individual parameter errors of combined parameter errors + at LDRtrue 0.004 0.020 0.100 0.300 0.450 0.004 0.020 0.100 0.300 0.450 + +Qin 9 9 9 9 9 3 3 3 3 3 +RotL 6 6 6 6 6 3 3 3 3 3 +DiC 132 132 132 132 132 23 23 23 23 23 +TP 317 317 317 317 317 7 7 7 7 7 +TS 300 300 300 300 300 13 13 13 13 13 +RP 33 33 33 33 33 0 0 0 0 0 +RS 31 31 31 31 31 1 1 1 1 1 +RotaT 0 0 0 0 0 0 0 0 0 0 +LDRCal 172 172 172 172 172 28 28 28 28 28 +Sum 1000 1000 1000 1000 1000 79 79 79 79 79 + +Etamx -45°: contribution to the total error (per mill) + of individual parameter errors of combined parameter errors + at LDRtrue 0.004 0.020 0.100 0.300 0.450 0.004 0.020 0.100 0.300 0.450 + +Qin 10 10 10 10 10 4 4 4 4 4 +RotL 6 6 6 6 6 4 4 4 4 4 +DiC 139 139 139 139 139 27 27 27 27 27 +TP 303 303 303 303 303 7 7 7 7 7 +TS 284 284 284 284 284 13 13 13 13 13 +RP 31 31 31 31 31 0 0 0 0 0 +RS 29 29 29 29 29 1 1 1 1 1 +RotaT 0 0 0 0 0 0 0 0 0 0 +LDRCal 197 197 197 197 197 33 33 33 33 33 +Sum 999 999 999 999 999 89 89 89 89 89