Thu, 16 Feb 2017 21:34:49 +0100
Changed parameter name for laser degree of linear polarization from bL to DOLP.
Added loop for DOLP uncertainty.
Input file changed accordingly.
lidar_correction_ghk.py | file | annotate | diff | comparison | revisions | |
system_settings/optic_input_example_lidar.py | file | annotate | diff | comparison | revisions |
--- 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)
--- a/system_settings/optic_input_example_lidar.py Tue Jan 24 02:39:39 2017 +0100 +++ b/system_settings/optic_input_example_lidar.py Thu Feb 16 21:34:49 2017 +0100 @@ -3,6 +3,7 @@ # which might improve the portability of the code within an executable. # Due to problems I had with some two letter variables, most variables are now with at least # three letters mixed small and capital. +# To be used with lidar_correction_ghk.py ver. 0.9.5 and larger # Do you want to calculate the errors? If not, just the GHK-parameters are determined. Error_Calc = True @@ -14,7 +15,7 @@ print(" Lidar system :", EID, ", ", LID) # +++ IL Laser 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., 2., 1 #alpha; rotation of laser polarization in degrees; default 0 # +++ ME Emitter optics and +-Uncertainty; default = no emitter optics @@ -122,6 +123,7 @@ LDRmeas = 0.3 # --- this is just for correct transfer of the variables to the main file +DOLP0, dDOLP, nDOLP = DOLP, dDOLP, nDOLP RotL0, dRotL, nRotL = RotL, dRotL, nRotL # Emitter DiE0, dDiE, nDiE = DiE, dDiE, nDiE