--- a/lidar_correction_ghk.py Tue Jan 24 02:39:39 2017 +0100
+++ b/lidar_correction_ghk.py Thu Feb 16 21:34:49 2017 +0100
@@ -1,6 +1,6 @@
# -*- coding: utf-8 -*-
"""
-Copyright 2016 Volker Freudenthaler
+Copyright 2016, 2017 Volker Freudenthaler
Licensed under the EUPL, Version 1.1 only (the "Licence").
@@ -117,7 +117,7 @@
LID = "internal"
EID = "internal"
# --- IL Laser IL and +-Uncertainty
-bL = 1. # degree of linear polarization; default 1
+DOLP, dDOLP, nDOLP = 0.995, 0.005, 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
@@ -204,6 +204,8 @@
# --- Read actual lidar system parameters from optic_input.py (should be in sub-directory 'system_settings')
InputFile = 'optic_input_example_lidar.py'
+InputFile = 'optic_input_ver8c_POLIS_532_Sep2016.py'
+InputFile = 'optic_input_Bertha_Saltrace_532.py'
'''
print("From ", dname)
@@ -226,8 +228,8 @@
##TypeC = 6 Don't change the TypeC here
# RotationErrorEpsilonForNormalMeasurements = True
# LDRCal = 0.25
-# bL = 0.8
## --- Errors
+DOLP0, dDOLP, nDOLP = DOLP, dDOLP, nDOLP
RotL0, dRotL, nRotL = RotL, dRotL, nRotL
DiE0, dDiE, nDiE = DiE, dDiE, nDiE
@@ -274,7 +276,7 @@
# --- 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,
+def Calc(DOLP, 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
@@ -296,11 +298,11 @@
S2g = np.sin(np.deg2rad(2 * RotO))
C2g = np.cos(np.deg2rad(2 * RotO))
- # Laser with Degree of linear polarization DOLP = bL
+ # Laser with Degree of linear polarization DOLP
IinL = 1.
- QinL = bL
+ QinL = DOLP
UinL = 0.
- VinL = (1. - bL ** 2) ** 0.5
+ VinL = (1. - DOLP ** 2) ** 0.5
# Stokes Input Vector rotation Eq. E.4
A = C2a * QinL - S2a * UinL
@@ -380,7 +382,7 @@
S2aR = np.sin(np.deg2rad(h * 2 * RotaR))
C2aR = np.cos(np.deg2rad(2 * RotaR))
- # Aanalyzer As before the PBS Eq. D.5
+ # Analyzer As before the PBS Eq. D.5
ATP1 = (1 + C2aT * DaT * DiT)
ATP2 = Y * (DiT + C2aT * DaT)
ATP3 = Y * S2aT * DaT * ZiT * CosT
@@ -453,7 +455,7 @@
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
+ # 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])
@@ -474,7 +476,7 @@
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
+ # 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])
@@ -501,7 +503,7 @@
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
+ # 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])
@@ -589,7 +591,7 @@
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
+ # Correction parameters for normal measurements; they are independent of LDR
if (not RotationErrorEpsilonForNormalMeasurements):
GT = ATF1 * IinE + ATF4 * VinE
GR = ARF1 * IinE + ARF4 * VinE
@@ -609,7 +611,7 @@
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
+ # 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])
@@ -639,7 +641,7 @@
BT = np.dot(ATFe, IF2e)
BR = np.dot(ARFe, IF2e)
- # Correction paremeters for normal measurements; they are independent of LDR
+ # 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])
@@ -724,7 +726,7 @@
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
+ # 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
@@ -744,7 +746,7 @@
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
+ # 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
@@ -775,7 +777,7 @@
BT = np.dot(ATEe, IE2e)
BR = np.dot(AREe, IE2e)
- # Correction paremeters for normal measurements; they are independent of LDR
+ # 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
@@ -848,7 +850,7 @@
# *******************************************************************************************************************************
# --- CALC truth with assumed true parameters from the input file
-GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0,
+GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(DOLP0, RotL0, RotE0, RetE0, DiE0, RotO0,
RetO0, DiO0, RotC0, RetC0, DiC0,
TP0, TS0, RP0, RS0, ERaT0,
RotaT0, RetT0, ERaR0, RotaR0,
@@ -863,8 +865,8 @@
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,
+ print("{0:5}{1:5} {2:6.4f}±{3:7.4f}/{4:2d}; {5:8} {6:8.4f}±{7:7.4f}/{8:2d}".format("Laser: ", "DOLP =", DOLP0, dDOLP, nDOLP,
+ " 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(
@@ -916,7 +918,7 @@
# 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
for i, LDRCal in enumerate(LDRCalList):
- GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0,
+ GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(DOLP0, RotL0, RotE0, RetE0,
DiE0, RotO0, RetO0,
DiO0, RotC0, RetC0,
DiC0, TP0, TS0, RP0,
@@ -944,7 +946,7 @@
# Still with assumed true parameters in input file
for i, LDRtrue in enumerate(LDRrange):
#for LDRtrue in LDRrange:
- GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0,
+ GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(DOLP0, RotL0, RotE0, RetE0,
DiE0, RotO0, RetO0,
DiO0, RotC0, RetC0,
DiC0, TP0, TS0, RP0,
@@ -974,7 +976,7 @@
'''
if (Error_Calc):
# --- CALC again assumed truth with LDRCal0 and with assumed true parameters in input file to reset all 0-values
- GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0,
+ GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(DOLP0, RotL0, RotE0, RetE0, DiE0,
RotO0, RetO0, DiO0, RotC0,
RetC0, DiC0, TP0, TS0, RP0,
RS0, ERaT0, RotaT0, RetT0,
@@ -983,7 +985,7 @@
# --- Start Errors calculation with variable parameters ------------------------------------------------------------------
iN = -1
- N = ((nRotL * 2 + 1) *
+ N = ((nDOLP * 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) *
@@ -1012,6 +1014,7 @@
# aLDRsimx2 = np.zeros(N)
# aLDRcorr = np.zeros(N)
# aLDRcorr2 = np.zeros(N)
+ aDOLP = np.zeros(N)
aERaT = np.zeros(N)
aERaR = np.zeros(N)
aRotaT = np.zeros(N)
@@ -1056,20 +1059,22 @@
LDRCal = LDRCal0
if nLDRCal > 0: LDRCal = LDRCal0 + iLDRCal * dLDRCal / nLDRCal
- GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim = Calc(RotL0, RotE0, RetE0,
+ GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim = Calc(DOLP0, 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 iDOLP, iRotL, iRotE, iRetE, iDiE \
+ in [(iDOLP, iRotL, iRotE, iRetE, iDiE)
+ for iDOLP in range(-nDOLP, nDOLP + 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 nDOLP > 0: DOLP = DOLP0 + iDOLP * dDOLP / nDOLP
if nRotL > 0: RotL = RotL0 + iRotL * dRotL / nRotL
if nRotE > 0: RotE = RotE0 + iRotE * dRotE / nRotE
if nRetE > 0: RetE = RetE0 + iRetE * dRetE / nRetE
@@ -1084,11 +1089,11 @@
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
+ # Laser with Degree of linear polarization DOLP
IinL = 1.
- QinL = bL
+ QinL = DOLP
UinL = 0.
- VinL = (1. - bL ** 2) ** 0.5
+ VinL = (1. - DOLP ** 2) ** 0.5
# Stokes Input Vector rotation Eq. E.4
A = C2a * QinL - S2a * UinL
@@ -1276,7 +1281,7 @@
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
+ # 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])
@@ -1297,7 +1302,7 @@
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
+ # 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])
@@ -1328,7 +1333,7 @@
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
+ # 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])
@@ -1418,7 +1423,7 @@
AR = ARF1 * IinF + ARF4 * h * VinF
BR = ARF3e * QinF - ARF2e * h * UinF
- # Correction paremeters for normal measurements; they are independent of LDR
+ # Correction parameters for normal measurements; they are independent of LDR
if (not RotationErrorEpsilonForNormalMeasurements):
GT = ATF1 * IinE + ATF4 * VinE
GR = ARF1 * IinE + ARF4 * VinE
@@ -1440,7 +1445,7 @@
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
+ # 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])
@@ -1474,7 +1479,7 @@
BT = np.dot(ATFe, IF2e)
BR = np.dot(ARFe, IF2e)
- # Correction paremeters for normal measurements; they are independent of LDR
+ # 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])
@@ -1562,7 +1567,7 @@
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
+ # 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
@@ -1582,7 +1587,7 @@
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
+ # 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
@@ -1613,7 +1618,7 @@
BT = np.dot(ATEe, IE2e)
BR = np.dot(AREe, IE2e)
- # Correction paremeters for normal measurements; they are independent of LDR
+ # 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
@@ -1704,6 +1709,7 @@
aX[4, iN] = K
aLDRCal[iN] = iLDRCal
+ aDOLP[iN] = iDOLP
aERaT[iN] = iERaT
aERaR[iN] = iERaR
aRotaT[iN] = iRotaT
@@ -1786,6 +1792,7 @@
return
+ if (nDOLP > 0): PlotSubHist("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP)
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)
@@ -1863,6 +1870,7 @@
#plt.close
return
+ if (nDOLP > 0): PlotSubHistF11("DOLP", aDOLP, DOLP0, dDOLP, iDOLP, nDOLP)
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)
