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