Mercurial > public > atmospheric_lidar_ghk / changeset
changeset
Error calculation is now optional with parameter "Error_Calc = True/False" in input file.
Tue, 15 Nov 2016 14:17:34 +0100
- author
- Volker Freudenthaler <volker.freudenthaler@lmu.de>
- date
- Tue, 15 Nov 2016 14:17:34 +0100
- changeset 16
- 313ac320b970
- parent 15
- 94eac33c6e6e
- child 17
- 43fe065e63b6
Error calculation is now optional with parameter "Error_Calc = True/False" in input file.
Error plots of experimental F11corr suppressed.
Example lidar input values changed.
Output slightly reformated.
--- a/lidar_correction_ghk.py Tue Nov 15 03:47:25 2016 +0100 +++ b/lidar_correction_ghk.py Tue Nov 15 14:17:34 2016 +0100 @@ -83,6 +83,8 @@ sqr05 = 0.5**0.5 +# Do you want to calculate the errors? If not, just the GHK-parameters are determined. +Error_Calc = True # ---- Initial definition of variables; the actual values will be read in with exec(open('./optic_input.py').read()) below LID = "internal" EID = "internal" @@ -172,7 +174,6 @@ # ******************************************************************************************************************************* # --- Read actual lidar system parameters from ./optic_input.py (must be in the same directory) - InputFile = 'optic_input_example_lidar.py' #InputFile = 'optic_input_ver6e_POLIS_355.py' #InputFile = 'optic_input_ver6e_POLIS_355_JA.py' @@ -843,10 +844,11 @@ 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}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format("ERT, RotT :", ERaT0, dERaT, nERaT, RotaT0, dRotaT, nRotaT)) + print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format("ERR, RotR :", ERaR0, dERaR, nERaR, RotaR0, dRotaR, nRotaR)) print("{0:12}".format(" --- PBS ---")) - print("{0:12}{1:7.4f}±{2:7.4f}/{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:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format("TP,TS :", TP0, dTP, nTP, TS0, dTS, nTS)) + print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format("RP,RS :", RP0, dRP, nRP, RS0, dRS, nRS)) print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}, {5:1.0f}".format("DT,TT,DR,TR,Y :", DiT, TiT, DiR, TiR, Y)) print("{0:12}".format(" --- Combined PBS + Pol.-filter ---")) print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}".format("DT,TT,DR,TR :", DTa0, TTa0, DRa0, TRa0)) @@ -866,7 +868,7 @@ #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:26}: {1:6.3f}±{2:5.3f}/{3:2d}".format("LDRCal during calibration in calibration range", LDRCal0, dLDRCal, nLDRCal)) #print("{0:8}={1:8.5f};{2:8}={3:8.5f}".format(" IinP",IinP," F11sim",F11sim)) print() @@ -908,921 +910,923 @@ 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) +if(Error_Calc): + # --- 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 with variable parameters ------------------------------------------------------------------ + # --- Start Errors calculation with variable parameters ------------------------------------------------------------------ -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="") + 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 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() + #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) + # --- 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)) + 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() + 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 + # --- 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) + # 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 + 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)]: + 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 + 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)) + # 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 + # 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 + # 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) + # 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 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 - #------------------------- + #------------------------- + # 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)]: + 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): + 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() - 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 ((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 + 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)) + #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) + # receiver optics + CosO = np.cos(np.deg2rad(RetO)) + SinO = np.sin(np.deg2rad(RetO)) + ZiO = (1. - DiO**2)**0.5 + WiO = (1. - ZiO*CosO) + S2g = np.sin(np.deg2rad(2*RotO)) + C2g = np.cos(np.deg2rad(2*RotO)) + # calibrator + CosC = np.cos(np.deg2rad(RetC)) + SinC = np.sin(np.deg2rad(RetC)) + ZiC = (1. - DiC**2)**0.5 + WiC = (1. - ZiC*CosC) - # analyser - #For POLLY_XTs - if(RS_RP_depend_on_TS_TP): - RS = 1 - 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)) + # analyser + #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 + 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) + 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)) + 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]) + # 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]) + 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 + 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") + # ---- 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 + #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) + 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 + 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) + # 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: - 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() + print("Calibrator not implemented yet") + sys.exit() + + elif LocC == 3: # C before receiver optics Eq.57 - 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)) + #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]) + # 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] + 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] + # 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]) + ATFe = np.array([ATF1,ATF2e,ATF3e,ATF4]) + ARFe = np.array([ARF1,ARF2e,ARF3e,ARF4]) - QinEe = C2e*QinE + S2e*UinE - UinEe = C2e*UinE - S2e*QinE + 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 + # 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 + IinFo = IinE + QinFo = 0. + UinFo = 0. + VinFo = VinE - IinFa = 0. - QinFa = QinE - UinFa = -UinE - VinFa = -2.*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 + # 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 + # 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 + # 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]) + 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) + 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]) + # 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) + 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]) + 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) + 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) + # 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() + 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)) + 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]) + # 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] + 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 + # 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 - 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) - # rotated IinE - QinEe = C2e*QinE + S2e*UinE - UinEe = C2e*UinE - S2e*QinE + # 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) - # --- 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*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 - # 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) + 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): - # 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 + # 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: - 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 + print('Calibrator not implemented yet') + sys.exit() - 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) + # 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)) - # 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 + # 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: - 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 + LDRsimx = 1./LDRsim + LDRsimx2 = 1./LDRsim2 + ''' + # ----- Backward correction + # Corrected LDRCorr from forward simulated LDRsim (atrue) with assumed true G0,H0,K0 + LDRCorr = (LDRsim*K0/Etax*(GT0+HT0)-(GR0+HR0))/((GR0-HR0)-LDRsim*K0/Etax*(GT0-HT0)) + + # -- F11corr from It and Ir and calibration EtaX + Text1 = "!!! EXPERIMENTAL !!! F11corr from It and Ir with calibration EtaX: x-axis: F11corr(LDRtrue) / F11corr(LDRtrue = 0.004) - 1" + F11corr = 1/(TiO*TiE)*((HR0*Etax/K0*It/TTa-HT0*Ir/TRa)/(HR0*GT0-HT0*GR0)) # IL = 1 Eq.(64) - else: - print('Calibrator not implemented yet') - sys.exit() + #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 - # 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)) + 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 - # 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 + # --- END loop + btime = clock() + print("\r done in ", "{0:5.0f}".format(btime-atime), "sec") #, end="\r") + + # --- Plot ----------------------------------------------------------------- + if (sns_loaded): + sns.set_style("whitegrid") + sns.set_palette("bright", 6) - # ************************************************************************* - 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) + ''' + 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 - # 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)) + 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') - # -- F11corr from It and Ir and calibration EtaX - Text1 = "!!! EXPERIMENTAL !!! F11corr from It and Ir with calibration EtaX: x-axis: F11corr(LDRtrue) / F11corr(LDRtrue = 0.004) - 1" - F11corr = 1/(TiO*TiE)*((HR0*Etax/K0*It/TTa-HT0*Ir/TRa)/(HR0*GT0-HT0*GR0)) # IL = 1 Eq.(64) + 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 + ' with ' + 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 - #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 + 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) - aX[0,iN] = GR - aX[1,iN] = GT - aX[2,iN] = HR - aX[3,iN] = HT - aX[4,iN] = K + plt.show() + plt.close + ''' + print() + #print("IT(LDRtrue) devided by IT(LDRtrue = 0.004)") + print(" ############################################################################## ") + print(Text1) + print() + + iLDR = 5 + for LDRTrue in LDRrange: + iLDR = iLDR - 1 + aF11corr[iLDR,:] = aF11corr[iLDR,:] / aF11corr[0,:] - 1.0 - aLDRCal[iN] = iLDRCal - aERaT[iN] = iERaT - aERaR[iN] = iERaR - aRotaT[iN] = iRotaT - aRotaR[iN] = iRotaR - aRetT[iN] = iRetT - aRetR[iN] = iRetR + # 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 - 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 + #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 -# --- END loop -btime = clock() -print("\r done in ", "{0:5.0f}".format(btime-atime), "sec") #, end="\r") - -# --- Plot ----------------------------------------------------------------- -if (sns_loaded): - sns.set_style("whitegrid") - sns.set_palette("bright", 6) + 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) -''' -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)) + 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 = 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) + 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,:], - 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() - + 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.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 + plt.show() + plt.close + ''' -print() -#print("IT(LDRtrue) devided by IT(LDRtrue = 0.004)") -print(" ############################################################################## ") -print(Text1) -print() - -iLDR = 5 -for LDRTrue in LDRrange: - iLDR = iLDR - 1 - aF11corr[iLDR,:] = aF11corr[iLDR,:] / aF11corr[0,:] - 1.0 - -# Plot F11 -def PlotSubHistF11(aVar, aX, X0, daX, iaX, naX): - fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) - iLDR = -1 - for LDRTrue in LDRrange: - iLDR = iLDR + 1 + # --- 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') - ''' - 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.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 - #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) + # --- Save LDRmin, LDRmax to file + # http://stackoverflow.com/questions/4675728/redirect-stdout-to-a-file-in-python + with open('output_files\LDR_min_max_' + 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])) -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') + ''' + # --- Plot K over LDRCal + fig3 = plt.figure() + plt.plot(LDRCal0+aLDRCal*dLDRCal/nLDRCal,aX[4,:], linewidth=2.0, color='b') -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('output_files\LDR_min_max_' + 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 -''' + plt.xlabel('LDRCal', fontsize=18) + plt.ylabel('K', fontsize=14) + plt.title(LID, fontsize=18) + plt.show() + plt.close + ''' # Additional plot routines ======> '''
--- a/output_files/LDR_min_max_example lidar.dat Tue Nov 15 03:47:25 2016 +0100 +++ b/output_files/LDR_min_max_example lidar.dat Tue Nov 15 14:17:34 2016 +0100 @@ -1,7 +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 + 0.0040, 0.0039, 0.0057 + 0.0200, 0.0197, 0.0219 + 0.1000, 0.0988, 0.1026 + 0.3000, 0.2965, 0.3042 + 0.4500, 0.4448, 0.4554
--- a/output_files/output_example lidar.dat Tue Nov 15 03:47:25 2016 +0100 +++ b/output_files/output_example lidar.dat Tue Nov 15 14:17:34 2016 +0100 @@ -3,32 +3,33 @@ 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 +Laser: DOLP = 1.00000; rotation alpha = 0.0000± 2.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 +Receiver 0.0000± 0.1000/ 1, 1.0000, 0±180/ 0, 0.0000± 0.5000/ 0 +Calibrator 0.9998± 0.0001/ 1, 0.4000, 0± 0/ 0, 0.0000± 0.1000/ 0 --- 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 +ERT, RotT : 0.0010± 0.0010/ 1, 0.0000± 1.0000/ 1 +ERR, RotR : 0.0010± 0.0010/ 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 +TP,TS : 0.9500± 0.0100/ 1, 0.0200± 0.0100/ 1 +RP,RS : 0.0500± 0.0100/ 1, 0.9800± 0.0100/ 1 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 +DT,TT,DR,TR : 1.0000, 0.2427, -0.9999, 0.2578 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 +RS_RP_depend_on_TS_TP = False +LDRCal during calibration in calibration range: 0.009±0.005/ 1 ======================================================================== 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 + 1.90111, 1.95685,-1.90091, 1.95676, 0.93372, 0.94595, 0.95689 ======================================================================== 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 + 0.00400, 0.00418, 0.00400 + 0.02000, 0.02068, 0.02000 + 0.20000, 0.20637, 0.20000 + 0.45000, 0.46426, 0.45000
--- a/system_settings/optic_input_example_lidar.py Tue Nov 15 03:47:25 2016 +0100 +++ b/system_settings/optic_input_example_lidar.py Tue Nov 15 14:17:34 2016 +0100 @@ -4,17 +4,18 @@ # Due to problems I had with some two letter variables, most variables are now with at least # three letters mixed small and capital. +# Do you want to calculate the errors? If not, just the GHK-parameters are determined. +Error_Calc = True + # Header to identify the lidar system # Values of DO, DT, and DR etc. from fit to lamp calibrations in Leipzig (LampCalib_2_invers_c_D0=0.opj) EID = "xx" # Earlinet station ID LID = "example lidar" # Additional lidar ID (short descriptive text) -# firet fit intern (FITLN1) => DO = 0, DT fixed -0.9998, eta and DR fitted, -# => internal calib with LinPol before the receiver print(" Lidar system :", EID, ", ", LID) # +++ IL Laser and +-Uncertainty bL = 1. #degree of linear polarization; default 1 -RotL, dRotL, nRotL = 0., 1., 1 #alpha; rotation of laser polarization in degrees; default 0 +RotL, dRotL, nRotL = 0., 2., 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 @@ -23,9 +24,9 @@ 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.01, 1 +DiO, dDiO, nDiO = 0.0, 0.1, 1 TiO = 1.0 -RetO, dRetO, nRetO = 0., 180., 2 +RetO, dRetO, nRetO = 0., 180., 0 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 @@ -34,8 +35,8 @@ 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° +ERaT, dERaT, nERaT = 0.001, 0.001, 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) @@ -44,7 +45,7 @@ # +++++ 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 +RS_RP_depend_on_TS_TP = False # --- Polarizing beam splitter reflecting path if(RS_RP_depend_on_TS_TP): RP, dRP, nRP = 1-TP, 0.00, 0 # do not change this @@ -52,10 +53,10 @@ 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 +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° +ERaR, dERaR, nERaR = 0.001, 0.001, 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) @@ -84,11 +85,11 @@ #NOTE: use here twice the HWP-rotation-angle RotC, dRotC, nRotC = 0.0, 0.1, 1 #constant calibrator offset epsilon RotationErrorEpsilonForNormalMeasurements = True # is in general True for TypeC == 2 calibrator -elif TypeC == 3: # linear polarizer calibrator +elif TypeC == 3: # linear polarizer calibrator. Diattenuation DiC = (1-ERC)/(1+ERC); ERC = extinction ratio of calibrator DiC, dDiC, nDiC = 0.9998, 0.0001, 1 # ideal 1.0 - TiC = 0.505 # ideal 0.5 + TiC = 0.4 # ideal 0.5 RetC, dRetC, nRetC = 0., 0., 0 - RotC, dRotC, nRotC = 0.0, 0.1, 1 #constant calibrator offset epsilon + RotC, dRotC, nRotC = 0.0, 0.1, 0 #constant calibrator offset epsilon RotationErrorEpsilonForNormalMeasurements = False # is in general False for TypeC == 3 calibrator elif TypeC == 4: # QWP calibrator DiC, dDiC, nDiC = 0.0, 0., 0 # ideal 1.0 @@ -101,14 +102,14 @@ TiC = 1. RetC, dRetC, nRetC = 180., 0., 0 #Note: use real HWP angles here - RotC, dRotC, nRotC = 0.0, 0.1, 1 #constant calibrator offset epsilon -1.15 + RotC, dRotC, nRotC = 0.0, 0.1, 1 #constant calibrator offset epsilon RotationErrorEpsilonForNormalMeasurements = True # is in general True for TypeC == 6 calibrator else: print ('calibrator not implemented yet') sys.exit() -# --- LDRCal assumed atmospheric linear depolarization ratio during the calibration measurements (first guess) -LDRCal,dLDRCal,nLDRCal= 0.008, 0.003, 0 +# --- LDRCal assumed atmospheric linear depolarization ratio during the calibration measurements in calibration range with almost clean air (first guess) +LDRCal,dLDRCal,nLDRCal= 0.009, 0.005, 1 # spans the interference filter influence # ==================================================== # NOTE: there is no need to change anything below.
