Mercurial > public > atmospheric_lidar_ghk / file comparison

comparison: lidar_correction_ghk.py

lidar_correction_ghk.py

changeset 13
f08818615e3a
parent 11
453b23dd7f94
child 14
82dba9904149
equal deleted inserted replaced
12:8badc005e347 13:f08818615e3a
109 ERaT, dERaT, nERaT = 0.001, 0.001, 1 109 ERaT, dERaT, nERaT = 0.001, 0.001, 1
110 RotaT, dRotaT, nRotaT = 0., 3., 1 110 RotaT, dRotaT, nRotaT = 0., 3., 1
111 DaT = (1-ERaT)/(1+ERaT) 111 DaT = (1-ERaT)/(1+ERaT)
112 TaT = 0.5*(1+ERaT) 112 TaT = 0.5*(1+ERaT)
113 # --- PBS MR reflecting path defined with (RS,RP); and +-Uncertainty 113 # --- PBS MR reflecting path defined with (RS,RP); and +-Uncertainty
114 RS, dRS, nRS = 1 - TS, 0., 0 114 RS_RP_depend_on_TS_TP = False
115 RP, dRP, nRP = 1 - TP, 0., 0 115 if(RS_RP_depend_on_TS_TP):
116 RP, dRP, nRP = 1-TP, 0.00, 0
117 RS, dRS, nRS = 1-TS, 0.00, 0
118 else:
119 RP, dRP, nRP = 0.05, 0.01, 1
120 RS, dRS, nRS = 0.98, 0.01, 1
116 TiR = 0.5 * (RP + RS) 121 TiR = 0.5 * (RP + RS)
117 DiR = (RP-RS)/(RP+RS) 122 DiR = (RP-RS)/(RP+RS)
118 # PolFilter 123 # PolFilter
119 RetR, dRetR, nRetR = 0., 180., 0 124 RetR, dRetR, nRetR = 0., 180., 0
120 ERaR, dERaR, nERaR = 0.001, 0.001, 1 125 ERaR, dERaR, nERaR = 0.001, 0.001, 1
266 ZiT = (1. - DiT**2)**0.5 271 ZiT = (1. - DiT**2)**0.5
267 TiR = 0.5 * (RP + RS) 272 TiR = 0.5 * (RP + RS)
268 DiR = (RP-RS)/(RP+RS) 273 DiR = (RP-RS)/(RP+RS)
269 ZiR = (1. - DiR**2)**0.5 274 ZiR = (1. - DiR**2)**0.5
270 275
271 # -------------------------------------------------------- 276 # --- this subroutine is for the calculation with certain fixed parameters -----------------------------------------------------
272 def Calc(RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC, TP, TS, RP, RS, ERaT, RotaT, RetT, ERaR, RotaR, RetR, LDRCal): 277 def Calc(RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC, TP, TS, RP, RS, ERaT, RotaT, RetT, ERaR, RotaR, RetR, LDRCal):
273 # ---- Do the calculations of bra-ket vectors 278 # ---- Do the calculations of bra-ket vectors
274 h = -1. if TypeC == 2 else 1 279 h = -1. if TypeC == 2 else 1
275 # from input file: assumed LDRCal for calibration measurements 280 # from input file: assumed LDRCal for calibration measurements
276 aCal = (1.-LDRCal)/(1+LDRCal) 281 aCal = (1.-LDRCal)/(1+LDRCal)
346 #------------------------- 351 #-------------------------
347 # F11 assuemd to be = 1 => measured: F11m = IinP / IinE with atrue 352 # F11 assuemd to be = 1 => measured: F11m = IinP / IinE with atrue
348 #F11sim = TiO*(IinE + DiO*atrue*A)/IinE 353 #F11sim = TiO*(IinE + DiO*atrue*A)/IinE
349 #------------------------- 354 #-------------------------
350 355
351 # For PollyXT
352 # analyser 356 # analyser
353 #RS = 1 - TS 357 #For POLLY_XTs
354 #RP = 1 - TP 358 if(RS_RP_depend_on_TS_TP):
355 359 RS = 1 - TS
360 RP = 1 - TP
356 TiT = 0.5 * (TP + TS) 361 TiT = 0.5 * (TP + TS)
357 DiT = (TP-TS)/(TP+TS) 362 DiT = (TP-TS)/(TP+TS)
358 ZiT = (1. - DiT**2)**0.5 363 ZiT = (1. - DiT**2)**0.5
359 TiR = 0.5 * (RP + RS) 364 TiR = 0.5 * (RP + RS)
360 DiR = (RP-RS)/(RP+RS) 365 DiR = (RP-RS)/(RP+RS)
821 # ******************************************************************************************************************************* 826 # *******************************************************************************************************************************
822 827
823 # --- CALC truth 828 # --- CALC truth
824 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) 829 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)
825 830
826 # -------------------------------------------------------- 831 # --- Print parameters to console and output file
827 with open('output_' + LID + '.dat', 'w') as f: 832 with open('output_files\output_' + LID + '.dat', 'w') as f:
828 with redirect_stdout(f): 833 with redirect_stdout(f):
829 print("From ", dname) 834 print("From ", dname)
830 print("Running ", fname) 835 print("Running ", fname)
831 print("Reading input file ", InputFile) #, " for Lidar system :", EID, ", ", LID) 836 print("Reading input file ", InputFile) #, " for Lidar system :", EID, ", ", LID)
832 print("for Lidar system: ", EID, ", ", LID) 837 print("for Lidar system: ", EID, ", ", LID)
842 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)) 847 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))
843 print("{0:12}".format(" --- PBS ---")) 848 print("{0:12}".format(" --- PBS ---"))
844 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)) 849 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))
845 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)) 850 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))
846 print("{0:12}".format(" --- Combined PBS + Pol.-filter ---")) 851 print("{0:12}".format(" --- Combined PBS + Pol.-filter ---"))
847 print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}".format("DTa,TTa,DRa,TRa: ", DTa0, TTa0, DRa0, TRa0)) 852 print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}".format("DT,TT,DR,TR :", DTa0, TTa0, DRa0, TRa0))
848 print() 853 print()
849 print("Rotation Error Epsilon For Normal Measurements = ", RotationErrorEpsilonForNormalMeasurements) 854 print("Rotation Error Epsilon For Normal Measurements = ", RotationErrorEpsilonForNormalMeasurements)
850 #print ('LocC = ', LocC, Loc[LocC], '; TypeC = ',TypeC, Type[TypeC]) 855 print(Type[TypeC], Loc[LocC])
851 print(Type[TypeC], Loc[LocC], "; Parallel signal detected in", dY[int(Y+1)]) 856 print("Parallel signal detected in", dY[int(Y+1)])
857 print("RS_RP_depend_on_TS_TP = ", RS_RP_depend_on_TS_TP)
852 # end of print actual system parameters 858 # end of print actual system parameters
853 # ****************************************************************************** 859 # ******************************************************************************
854 860
855 #print() 861 #print()
856 #print(" --- LDRCal during calibration | simulated and corrected LDRs -------------") 862 #print(" --- LDRCal during calibration | simulated and corrected LDRs -------------")
871 for i,LDRCal in enumerate(LDRCalList): 877 for i,LDRCal in enumerate(LDRCalList):
872 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) 878 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)
873 K0List[i] = K0 879 K0List[i] = K0
874 LDRsimxList[i] = LDRsimx 880 LDRsimxList[i] = LDRsimx
875 881
876 print("{0:8},{1:8},{2:8},{3:8},{4:9},{5:9},{6:9}".format(" GR", " GT", " HR", " HT", " K(0.004)", " K(0.2)", " K(0.45)")) 882 print('========================================================================')
883 print("{0:8},{1:8},{2:8},{3:8},{4:9},{5:8},{6:9}".format(" GR", " GT", " HR", " HT", " K(0.004)", " K(0.2)", " K(0.45)"))
877 print("{0:8.5f},{1:8.5f},{2:8.5f},{3:8.5f},{4:9.5f},{5:9.5f},{6:9.5f}".format(GR0, GT0, HR0, HT0, K0List[0], K0List[1], K0List[2])) 884 print("{0:8.5f},{1:8.5f},{2:8.5f},{3:8.5f},{4:9.5f},{5:9.5f},{6:9.5f}".format(GR0, GT0, HR0, HT0, K0List[0], K0List[1], K0List[2]))
878 print('========================================================================') 885 print('========================================================================')
879 886
880 print("{0:9},{1:9},{2:9}".format(" LDRtrue", " LDRsimx", " LDRCorr")) 887 print("{0:9},{1:9},{2:9}".format(" LDRtrue", " LDRsimx", " LDRCorr"))
881 LDRtrueList = 0.004, 0.02, 0.2, 0.45 888 LDRtrueList = 0.004, 0.02, 0.2, 0.45
882 for i,LDRtrue in enumerate(LDRtrueList): 889 for i,LDRtrue in enumerate(LDRtrueList):
883 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) 890 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)
884 print("{0:9.5f},{1:9.5f},{2:9.5f}".format(LDRtrue, LDRsimx, LDRCorr)) 891 print("{0:9.5f},{1:9.5f},{2:9.5f}".format(LDRtrue, LDRsimx, LDRCorr))
885 892
886 893
887 file = open('output_' + LID + '.dat', 'r') 894 file = open('output_files\output_' + LID + '.dat', 'r')
888 print (file.read()) 895 print (file.read())
889 file.close() 896 file.close()
890 897
891 ''' 898 '''
892 if(PrintToOutputFile): 899 if(PrintToOutputFile):
902 sys.stdout = old_target 909 sys.stdout = old_target
903 ''' 910 '''
904 # --- CALC again truth with LDRCal0 to reset all 0-values 911 # --- CALC again truth with LDRCal0 to reset all 0-values
905 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) 912 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)
906 913
907 # --- Start Errors calculation 914 # --- Start Errors calculation with variable parameters ------------------------------------------------------------------
908 915
909 iN = -1 916 iN = -1
910 N = ((nRotL*2+1)* 917 N = ((nRotL*2+1)*
911 (nRotE*2+1)*(nRetE*2+1)*(nDiE*2+1)* 918 (nRotE*2+1)*(nRetE*2+1)*(nDiE*2+1)*
912 (nRotO*2+1)*(nRetO*2+1)*(nDiO*2+1)* 919 (nRotO*2+1)*(nRetO*2+1)*(nDiO*2+1)*
1090 CosC = np.cos(np.deg2rad(RetC)) 1097 CosC = np.cos(np.deg2rad(RetC))
1091 SinC = np.sin(np.deg2rad(RetC)) 1098 SinC = np.sin(np.deg2rad(RetC))
1092 ZiC = (1. - DiC**2)**0.5 1099 ZiC = (1. - DiC**2)**0.5
1093 WiC = (1. - ZiC*CosC) 1100 WiC = (1. - ZiC*CosC)
1094 1101
1095 #For POLLY_XT
1096 # analyser 1102 # analyser
1097 #RS = 1 - TS 1103 #For POLLY_XTs
1098 #RP = 1 - TP 1104 if(RS_RP_depend_on_TS_TP):
1105 RS = 1 - TS
1106 RP = 1 - TP
1099 TiT = 0.5 * (TP + TS) 1107 TiT = 0.5 * (TP + TS)
1100 DiT = (TP-TS)/(TP+TS) 1108 DiT = (TP-TS)/(TP+TS)
1101 ZiT = (1. - DiT**2)**0.5 1109 ZiT = (1. - DiT**2)**0.5
1102 TiR = 0.5 * (RP + RS) 1110 TiR = 0.5 * (RP + RS)
1103 DiR = (RP-RS)/(RP+RS) 1111 DiR = (RP-RS)/(RP+RS)
1795 plt.show() 1803 plt.show()
1796 plt.close 1804 plt.close
1797 1805
1798 # --- Save LDRmin, LDRmax to file 1806 # --- Save LDRmin, LDRmax to file
1799 # http://stackoverflow.com/questions/4675728/redirect-stdout-to-a-file-in-python 1807 # http://stackoverflow.com/questions/4675728/redirect-stdout-to-a-file-in-python
1800 with open('LDR_min_max_ver7_' + LID + '.dat', 'w') as f: 1808 with open('LDR_min_max_' + LID + '.dat', 'w') as f:
1801 with redirect_stdout(f): 1809 with redirect_stdout(f):
1802 print(LID) 1810 print(LID)
1803 print("LDRtrue, LDRmin, LDRmax") 1811 print("LDRtrue, LDRmin, LDRmax")
1804 for i in range(len(LDRrange)): 1812 for i in range(len(LDRrange)):
1805 print("{0:7.4f},{1:7.4f},{2:7.4f}".format(LDRrange[i], LDRmin[i], LDRmax[i])) 1813 print("{0:7.4f},{1:7.4f},{2:7.4f}".format(LDRrange[i], LDRmin[i], LDRmax[i]))

EARLINET Repository: Lidar polarisation GHK-parameter and error calculation

mercurial