16 permissions and limitations under the Licence. |
16 permissions and limitations under the Licence. |
17 |
17 |
18 Equation reference: http://www.atmos-meas-tech-discuss.net/amt-2015-338/amt-2015-338.pdf |
18 Equation reference: http://www.atmos-meas-tech-discuss.net/amt-2015-338/amt-2015-338.pdf |
19 With equations code from Appendix C |
19 With equations code from Appendix C |
20 Python 3.4.2 |
20 Python 3.4.2 |
|
21 |
|
22 Code description: |
|
23 |
|
24 From measured lidar signals we cannot directly determine the desired backscatter coefficient (F11) and the linear depolarization ratio (LDR) |
|
25 because of the cross talk between the channles and systematic errors of a lidar system. |
|
26 http://www.atmos-meas-tech-discuss.net/amt-2015-338/amt-2015-338.pdf provides an analytical model for the description of these errors, |
|
27 with which the measured signals can be corrected. |
|
28 This code simulates the lidar measurements with "assumed true" model parameters from an input file, and calculates the correction parameters (G,H, and K). |
|
29 The "assumed true" system parameters are the ones we think are the right ones, but in reality these parameters probably deviate from the assumed truth due to |
|
30 uncertainties. The uncertainties of the "assumed true" parameters can be described in the input file. Then this code calculates the lidar signals and the |
|
31 gain ratio eta* with all possible combinations of "errors", which represents the distribution of "possibly real" signals, and "corrects" them with the "assumed true" |
|
32 GHK parameters (GT0, GR0, HT0, HR0, and K0) to derive finally the distributions of "possibly real" linear depolarization ratios (LDRcorr), |
|
33 which are plotted for five different input linear depolarization ratios (LDRtrue). The red bars in the plots represent the input values of LDRtrue. |
|
34 A complication arises from the fact that the correction parameter K = eta*/eta (Eq. 83) can depend on the LDR during the calibration measurement, i.e. LDRcal or aCal |
|
35 in the code (see e.g. Eqs. (103), (115), and (141); mind the mistake in Eq. (116)). Therefor values of K for LDRcal = 0.004, 0.2, and 0.45 are calculated for |
|
36 "assumed true" system parameters and printed in the output file behind the GH parameters. The full impact of the LDRcal dependent K can be considered in the error |
|
37 calculation by specifying a range of possible LDRcal values in the input file. For the real calibration measurements a calibration range with low or no aerosol |
|
38 content should be chosen, and the default in the input file is a range of LDRcal between 0.004 and 0.014 (i.e. 0.009 +-0.005). |
|
39 |
|
40 Tip: In case you run the code with Spyder, all output text and plots can be displayed together in an IPython console, which can be saved as an html file. |
21 """ |
41 """ |
22 # !/usr/bin/env python3 |
|
23 |
42 |
24 # Comment: The code works with Python 2.7 with the help of following line, which enables Python2 to correctly interpret the Python 3 print statements. |
43 # Comment: The code works with Python 2.7 with the help of following line, which enables Python2 to correctly interpret the Python 3 print statements. |
25 from __future__ import print_function |
44 from __future__ import print_function |
|
45 # !/usr/bin/env python3 |
26 |
46 |
27 import os |
47 import os |
28 import sys |
48 import sys |
29 |
49 |
30 import numpy as np |
50 import numpy as np |
46 ## pp.savefig can be called multiple times to save to multiple pages |
66 ## pp.savefig can be called multiple times to save to multiple pages |
47 # pp.savefig() |
67 # pp.savefig() |
48 # pp.close() |
68 # pp.close() |
49 |
69 |
50 from contextlib import contextmanager |
70 from contextlib import contextmanager |
51 |
|
52 |
71 |
53 @contextmanager |
72 @contextmanager |
54 def redirect_stdout(new_target): |
73 def redirect_stdout(new_target): |
55 old_target, sys.stdout = sys.stdout, new_target # replace sys.stdout |
74 old_target, sys.stdout = sys.stdout, new_target # replace sys.stdout |
56 try: |
75 try: |
57 yield new_target # run some code with the replaced stdout |
76 yield new_target # run some code with the replaced stdout |
58 finally: |
77 finally: |
59 sys.stdout.flush() |
78 sys.stdout.flush() |
60 sys.stdout = old_target # restore to the previous value |
79 sys.stdout = old_target # restore to the previous value |
61 |
80 |
62 |
|
63 ''' |
81 ''' |
64 real_raw_input = vars(__builtins__).get('raw_input',input) |
82 real_raw_input = vars(__builtins__).get('raw_input',input) |
65 ''' |
83 ''' |
66 try: |
84 try: |
67 import __builtin__ |
85 import __builtin__ |
91 |
109 |
92 # PrintToOutputFile = True |
110 # PrintToOutputFile = True |
93 |
111 |
94 sqr05 = 0.5 ** 0.5 |
112 sqr05 = 0.5 ** 0.5 |
95 |
113 |
|
114 # ---- Initial definition of variables; the actual values will be read in with exec(open('./optic_input.py').read()) below |
96 # Do you want to calculate the errors? If not, just the GHK-parameters are determined. |
115 # Do you want to calculate the errors? If not, just the GHK-parameters are determined. |
97 Error_Calc = True |
116 Error_Calc = True |
98 # ---- Initial definition of variables; the actual values will be read in with exec(open('./optic_input.py').read()) below |
|
99 LID = "internal" |
117 LID = "internal" |
100 EID = "internal" |
118 EID = "internal" |
101 # --- IL Laser IL and +-Uncertainty |
119 # --- IL Laser IL and +-Uncertainty |
102 bL = 1. # degree of linear polarization; default 1 |
120 bL = 1. # degree of linear polarization; default 1 |
103 RotL, dRotL, nRotL = 0.0, 0.0, 1 # alpha; rotation of laser polarization in degrees; default 0 |
121 RotL, dRotL, nRotL = 0.0, 0.0, 1 # alpha; rotation of laser polarization in degrees; default 0 |
182 dY = ['reflected channel', '', 'transmitted channel'] |
200 dY = ['reflected channel', '', 'transmitted channel'] |
183 |
201 |
184 # end of initial definition of variables |
202 # end of initial definition of variables |
185 # ******************************************************************************************************************************* |
203 # ******************************************************************************************************************************* |
186 |
204 |
187 # --- Read actual lidar system parameters from ./optic_input.py (must be in the same directory) |
205 # --- Read actual lidar system parameters from optic_input.py (should be in sub-directory 'system_settings') |
188 InputFile = 'optic_input_example_lidar.py' |
206 InputFile = 'optic_input_example_lidar.py' |
189 # InputFile = 'optic_input_ver6e_POLIS_355.py' |
|
190 # InputFile = 'optic_input_ver6e_POLIS_355_JA.py' |
|
191 # InputFile = 'optic_input_ver6c_POLIS_532.py' |
|
192 # InputFile = 'optic_input_ver6e_POLIS_532.py' |
|
193 # InputFile = 'optic_input_ver8c_POLIS_532.py' |
|
194 # InputFile = 'optic_input_ver6e_MUSA.py' |
|
195 # InputFile = 'optic_input_ver6e_MUSA_JA.py' |
|
196 # InputFile = 'optic_input_ver6e_PollyXTSea.py' |
|
197 # InputFile = 'optic_input_ver6e_PollyXTSea_JA.py' |
|
198 # InputFile = 'optic_input_ver6e_PollyXT_RALPH.py' |
|
199 # InputFile = 'optic_input_ver8c_PollyXT_RALPH.py' |
|
200 # InputFile = 'optic_input_ver8c_PollyXT_RALPH_2.py' |
|
201 # InputFile = 'optic_input_ver8c_PollyXT_RALPH_3.py' |
|
202 # InputFile = 'optic_input_ver8c_PollyXT_RALPH_4.py' |
|
203 # InputFile = 'optic_input_ver8c_PollyXT_RALPH_5.py' |
|
204 # InputFile = 'optic_input_ver8c_PollyXT_RALPH_6.py' |
|
205 # InputFile = 'optic_input_ver8c_PollyXT_RALPH_7.py' |
|
206 # InputFile = 'optic_input_ver8a_MOHP_DPL_355.py' |
|
207 # InputFile = 'optic_input_ver9_MOHP_DPL_355.py' |
|
208 # InputFile = 'optic_input_ver6e_RALI.py' |
|
209 # InputFile = 'optic_input_ver6e_RALI_JA.py' |
|
210 # InputFile = 'optic_input_ver6e_RALI_new.py' |
|
211 # InputFile = 'optic_input_ver6e_RALI_act.py' |
|
212 # InputFile = 'optic_input_ver6e_MULHACEN.py' |
|
213 # InputFile = 'optic_input_ver6e_MULHACEN_JA.py' |
|
214 # InputFile = 'optic_input_ver6e-IPRAL.py' |
|
215 # InputFile = 'optic_input_ver6e-IPRAL_JA.py' |
|
216 # InputFile = 'optic_input_ver6e-LB21.py' |
|
217 # InputFile = 'optic_input_ver6e-LB21_JA.py' |
|
218 # InputFile = 'optic_input_ver6e_Bertha_b_355.py' |
|
219 # InputFile = 'optic_input_ver6e_Bertha_b_532.py' |
|
220 # InputFile = 'optic_input_ver6e_Bertha_b_1064.py' |
|
221 |
207 |
222 ''' |
208 ''' |
223 print("From ", dname) |
209 print("From ", dname) |
224 print("Running ", fname) |
210 print("Running ", fname) |
225 print("Reading input file ", InputFile, " for") |
211 print("Reading input file ", InputFile, " for") |
226 ''' |
212 ''' |
227 input_path = os.path.join('.', 'system_settings', InputFile) |
213 input_path = os.path.join('.', 'system_settings', InputFile) |
228 # this works with Python 2 - and 3? |
214 # this works with Python 2 and 3! |
229 exec (open(input_path).read(), globals()) |
215 exec (open(input_path).read(), globals()) |
230 # end of read actual system parameters |
216 # end of read actual system parameters |
|
217 |
231 |
218 |
232 # --- Manual Parameter Change --- |
219 # --- Manual Parameter Change --- |
233 # (use for quick parameter changes without changing the input file ) |
220 # (use for quick parameter changes without changing the input file ) |
234 # DiO = 0. |
221 # DiO = 0. |
235 # LDRtrue = 0.45 |
222 # LDRtrue = 0.45 |
366 # F11 assuemd to be = 1 => measured: F11m = IinP / IinE with atrue |
353 # F11 assuemd to be = 1 => measured: F11m = IinP / IinE with atrue |
367 # F11sim = TiO*(IinE + DiO*atrue*A)/IinE |
354 # F11sim = TiO*(IinE + DiO*atrue*A)/IinE |
368 # ------------------------- |
355 # ------------------------- |
369 |
356 |
370 # analyser |
357 # analyser |
371 # For POLLY_XTs |
|
372 if (RS_RP_depend_on_TS_TP): |
358 if (RS_RP_depend_on_TS_TP): |
373 RS = 1 - TS |
359 RS = 1 - TS |
374 RP = 1 - TP |
360 RP = 1 - TP |
|
361 |
375 TiT = 0.5 * (TP + TS) |
362 TiT = 0.5 * (TP + TS) |
376 DiT = (TP - TS) / (TP + TS) |
363 DiT = (TP - TS) / (TP + TS) |
377 ZiT = (1. - DiT ** 2) ** 0.5 |
364 ZiT = (1. - DiT ** 2) ** 0.5 |
378 TiR = 0.5 * (RP + RS) |
365 TiR = 0.5 * (RP + RS) |
379 DiR = (RP - RS) / (RP + RS) |
366 DiR = (RP - RS) / (RP + RS) |
851 |
838 |
852 TTa = TiT * TaT # *ATP1 |
839 TTa = TiT * TaT # *ATP1 |
853 TRa = TiR * TaR # *ARP1 |
840 TRa = TiR * TaR # *ARP1 |
854 |
841 |
855 F11sim = 1 / (TiO * TiE) * ( |
842 F11sim = 1 / (TiO * TiE) * ( |
856 (HR * Etax / K * It / TTa - HT * Ir / TRa) / (HR * GT - HT * GR)) # IL = 1, Etat = Etar = 1 |
843 (HR * Etax / K * It / TTa - HT * Ir / TRa) / (HR * GT - HT * GR)) # IL = 1, Etat = Etar = 1 ; AMT Eq.64; what is Etax/K? => see about 20 lines above: = Eta |
857 |
844 |
858 return (GT, HT, GR, HR, K, Eta, LDRsimx, LDRCorr, DTa, DRa, TTa, TRa, F11sim) |
845 return (GT, HT, GR, HR, K, Eta, LDRsimx, LDRCorr, DTa, DRa, TTa, TRa, F11sim) |
859 |
846 |
860 |
847 |
861 # ******************************************************************************************************************************* |
848 # ******************************************************************************************************************************* |
862 |
849 |
863 # --- CALC truth |
850 # --- CALC truth with assumed true parameters from the input file |
864 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, |
851 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, |
865 RetO0, DiO0, RotC0, RetC0, DiC0, |
852 RetO0, DiO0, RotC0, RetC0, DiC0, |
866 TP0, TS0, RP0, RS0, ERaT0, |
853 TP0, TS0, RP0, RS0, ERaT0, |
867 RotaT0, RetT0, ERaR0, RotaR0, |
854 RotaT0, RetT0, ERaR0, RotaR0, |
868 RetR0, LDRCal0) |
855 RetR0, LDRCal0) |
899 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, |
886 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, |
900 dRS, nRS)) |
887 dRS, nRS)) |
901 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)) |
888 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)) |
902 print("{0:12}".format(" --- Combined PBS + Pol.-filter ---")) |
889 print("{0:12}".format(" --- Combined PBS + Pol.-filter ---")) |
903 print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}".format("DT,TT,DR,TR :", DTa0, TTa0, DRa0, TRa0)) |
890 print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}".format("DT,TT,DR,TR :", DTa0, TTa0, DRa0, TRa0)) |
|
891 print("{0:26}: {1:6.3f}± {2:5.3f}/{3:2d}".format("LDRCal during calibration in calibration range", LDRCal0, |
|
892 dLDRCal, nLDRCal)) |
904 print() |
893 print() |
905 print("Rotation Error Epsilon For Normal Measurements = ", RotationErrorEpsilonForNormalMeasurements) |
894 print("Rotation Error Epsilon For Normal Measurements = ", RotationErrorEpsilonForNormalMeasurements) |
906 print(Type[TypeC], Loc[LocC]) |
895 print(Type[TypeC], Loc[LocC]) |
907 print("Parallel signal detected in", dY[int(Y + 1)]) |
896 print("Parallel signal detected in", dY[int(Y + 1)]) |
908 print("RS_RP_depend_on_TS_TP = ", RS_RP_depend_on_TS_TP) |
897 print("RS_RP_depend_on_TS_TP = ", RS_RP_depend_on_TS_TP) |
915 # print("{0:8.5f} |{1:8.5f}->{2:8.5f},{3:9.5f}->{4:9.5f} |{5:8.5f}->{6:8.5f}".format(LDRCal, LDRtrue, LDRsim, LDRtrue2, LDRsim2, LDRmeas, LDRCorr)) |
904 # print("{0:8.5f} |{1:8.5f}->{2:8.5f},{3:9.5f}->{4:9.5f} |{5:8.5f}->{6:8.5f}".format(LDRCal, LDRtrue, LDRsim, LDRtrue2, LDRsim2, LDRmeas, LDRCorr)) |
916 # print("{0:8} |{1:8}->{2:8}->{3:8}".format(" LDRCal"," LDRtrue", " LDRsimx", " LDRcorr")) |
905 # print("{0:8} |{1:8}->{2:8}->{3:8}".format(" LDRCal"," LDRtrue", " LDRsimx", " LDRcorr")) |
917 # 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)) |
906 # 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)) |
918 # print("{0:8} |{1:8}->{2:8}->{3:8}".format(" LDRCal"," LDRtrue", " LDRsimx", " LDRcorr")) |
907 # print("{0:8} |{1:8}->{2:8}->{3:8}".format(" LDRCal"," LDRtrue", " LDRsimx", " LDRcorr")) |
919 # print(" --- LDRCal during calibration") |
908 # print(" --- LDRCal during calibration") |
920 print("{0:26}: {1:6.3f}±{2:5.3f}/{3:2d}".format("LDRCal during calibration in calibration range", LDRCal0, |
|
921 dLDRCal, nLDRCal)) |
|
922 |
909 |
923 # print("{0:8}={1:8.5f};{2:8}={3:8.5f}".format(" IinP",IinP," F11sim",F11sim)) |
910 # print("{0:8}={1:8.5f};{2:8}={3:8.5f}".format(" IinP",IinP," F11sim",F11sim)) |
924 print() |
911 print() |
925 |
912 |
926 K0List = np.zeros(3) |
913 K0List = np.zeros(3) |
927 LDRsimxList = np.zeros(3) |
914 LDRsimxList = np.zeros(3) |
928 LDRCalList = 0.004, 0.2, 0.45 |
915 LDRCalList = 0.004, 0.2, 0.45 |
|
916 # 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. |
|
917 # Still with assumed true parameters in input file |
929 for i, LDRCal in enumerate(LDRCalList): |
918 for i, LDRCal in enumerate(LDRCalList): |
930 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, |
919 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, |
931 DiE0, RotO0, RetO0, |
920 DiE0, RotO0, RetO0, |
932 DiO0, RotC0, RetC0, |
921 DiO0, RotC0, RetC0, |
933 DiC0, TP0, TS0, RP0, |
922 DiC0, TP0, TS0, RP0, |
941 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)", |
930 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)", |
942 " K(0.45)")) |
931 " K(0.45)")) |
943 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], |
932 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], |
944 K0List[1], K0List[2])) |
933 K0List[1], K0List[2])) |
945 print('========================================================================') |
934 print('========================================================================') |
946 |
|
947 print("{0:9},{1:9},{2:9}".format(" LDRtrue", " LDRsimx", " LDRCorr")) |
935 print("{0:9},{1:9},{2:9}".format(" LDRtrue", " LDRsimx", " LDRCorr")) |
948 LDRtrueList = 0.004, 0.02, 0.2, 0.45 |
936 |
949 for i, LDRtrue in enumerate(LDRtrueList): |
937 #LDRtrueList = 0.004, 0.02, 0.2, 0.45 |
|
938 aF11sim0 = np.zeros(5) |
|
939 LDRrange = np.zeros(5) |
|
940 LDRrange = 0.004, 0.02, 0.1, 0.3, 0.45 |
|
941 |
|
942 # The loop over LDRtrueList is ony for checking how much the uncorrected LDRsimx deviates from LDRtrue ... and whether the corrections work. |
|
943 # LDRsimx = LDRsim = Ir / It or 1/LDRsim |
|
944 # Still with assumed true parameters in input file |
|
945 for i, LDRtrue in enumerate(LDRrange): |
|
946 #for LDRtrue in LDRrange: |
950 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, |
947 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, |
951 DiE0, RotO0, RetO0, |
948 DiE0, RotO0, RetO0, |
952 DiO0, RotC0, RetC0, |
949 DiO0, RotC0, RetC0, |
953 DiC0, TP0, TS0, RP0, |
950 DiC0, TP0, TS0, RP0, |
954 RS0, ERaT0, RotaT0, |
951 RS0, ERaT0, RotaT0, |
955 RetT0, ERaR0, RotaR0, |
952 RetT0, ERaR0, RotaR0, |
956 RetR0, LDRCal0) |
953 RetR0, LDRCal0) |
957 print("{0:9.5f},{1:9.5f},{2:9.5f}".format(LDRtrue, LDRsimx, LDRCorr)) |
954 print("{0:9.5f},{1:9.5f},{2:9.5f}".format(LDRtrue, LDRsimx, LDRCorr)) |
|
955 aF11sim0[i] = F11sim0 |
|
956 # the assumed true aF11sim0 results will be used below to calc the deviation from the real signals |
958 |
957 |
959 file = open('output_files\output_' + LID + '.dat', 'r') |
958 file = open('output_files\output_' + LID + '.dat', 'r') |
960 print(file.read()) |
959 print(file.read()) |
961 file.close() |
960 file.close() |
962 |
961 |
972 sys.stdout.flush() |
971 sys.stdout.flush() |
973 f.close |
972 f.close |
974 sys.stdout = old_target |
973 sys.stdout = old_target |
975 ''' |
974 ''' |
976 if (Error_Calc): |
975 if (Error_Calc): |
977 # --- CALC again truth with LDRCal0 to reset all 0-values |
976 # --- CALC again assumed truth with LDRCal0 and with assumed true parameters in input file to reset all 0-values |
978 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, |
977 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, |
979 RotO0, RetO0, DiO0, RotC0, |
978 RotO0, RetO0, DiO0, RotC0, |
980 RetC0, DiC0, TP0, TS0, RP0, |
979 RetC0, DiC0, TP0, TS0, RP0, |
981 RS0, ERaT0, RotaT0, RetT0, |
980 RS0, ERaT0, RotaT0, RetT0, |
982 ERaR0, RotaR0, RetR0, |
981 ERaR0, RotaR0, RetR0, |
983 LDRCal0) |
982 LDRCal0) |
984 |
|
985 # --- Start Errors calculation with variable parameters ------------------------------------------------------------------ |
983 # --- Start Errors calculation with variable parameters ------------------------------------------------------------------ |
986 |
984 |
987 iN = -1 |
985 iN = -1 |
988 N = ((nRotL * 2 + 1) * |
986 N = ((nRotL * 2 + 1) * |
989 (nRotE * 2 + 1) * (nRetE * 2 + 1) * (nDiE * 2 + 1) * |
987 (nRotE * 2 + 1) * (nRetE * 2 + 1) * (nDiE * 2 + 1) * |
1006 LDRmin = np.zeros(5) |
1004 LDRmin = np.zeros(5) |
1007 LDRmax = np.zeros(5) |
1005 LDRmax = np.zeros(5) |
1008 F11min = np.zeros(5) |
1006 F11min = np.zeros(5) |
1009 F11max = np.zeros(5) |
1007 F11max = np.zeros(5) |
1010 |
1008 |
1011 LDRrange = np.zeros(5) |
1009 #LDRrange = np.zeros(5) |
1012 LDRrange = 0.004, 0.02, 0.1, 0.3, 0.45 |
1010 #LDRrange = 0.004, 0.02, 0.1, 0.3, 0.45 |
1013 # aLDRsimx = np.zeros(N) |
1011 # aLDRsimx = np.zeros(N) |
1014 # aLDRsimx2 = np.zeros(N) |
1012 # aLDRsimx2 = np.zeros(N) |
1015 # aLDRcorr = np.zeros(N) |
1013 # aLDRcorr = np.zeros(N) |
1016 # aLDRcorr2 = np.zeros(N) |
1014 # aLDRcorr2 = np.zeros(N) |
1017 aERaT = np.zeros(N) |
1015 aERaT = np.zeros(N) |
1041 |
1039 |
1042 atime = clock() |
1040 atime = clock() |
1043 dtime = clock() |
1041 dtime = clock() |
1044 |
1042 |
1045 # --- Calc Error signals |
1043 # --- Calc Error signals |
1046 # GT, HT, GR, HR, K, Eta, LDRsim = Calc(RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC, TP, TS) |
|
1047 # ---- Do the calculations of bra-ket vectors |
1044 # ---- Do the calculations of bra-ket vectors |
1048 h = -1. if TypeC == 2 else 1 |
1045 h = -1. if TypeC == 2 else 1 |
1049 |
1046 |
|
1047 ''' |
1050 # from input file: measured LDRm and true LDRtrue, LDRtrue2 => |
1048 # from input file: measured LDRm and true LDRtrue, LDRtrue2 => |
1051 ameas = (1. - LDRmeas) / (1 + LDRmeas) |
1049 ameas = (1. - LDRmeas) / (1 + LDRmeas) |
1052 atrue = (1. - LDRtrue) / (1 + LDRtrue) |
1050 atrue = (1. - LDRtrue) / (1 + LDRtrue) |
1053 atrue2 = (1. - LDRtrue2) / (1 + LDRtrue2) |
1051 atrue2 = (1. - LDRtrue2) / (1 + LDRtrue2) |
|
1052 ''' |
1054 |
1053 |
1055 for iLDRCal in range(-nLDRCal, nLDRCal + 1): |
1054 for iLDRCal in range(-nLDRCal, nLDRCal + 1): |
1056 # from input file: assumed LDRCal for calibration measurements |
1055 # from input file: assumed LDRCal for calibration measurements |
1057 LDRCal = LDRCal0 |
1056 LDRCal = LDRCal0 |
1058 if nLDRCal > 0: LDRCal = LDRCal0 + iLDRCal * dLDRCal / nLDRCal |
1057 if nLDRCal > 0: LDRCal = LDRCal0 + iLDRCal * dLDRCal / nLDRCal |
1059 |
1058 |
1060 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, |
1059 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim = Calc(RotL0, RotE0, RetE0, |
1061 DiE0, RotO0, RetO0, DiO0, |
1060 DiE0, RotO0, RetO0, DiO0, |
1062 RotC0, RetC0, DiC0, TP0, |
1061 RotC0, RetC0, DiC0, TP0, |
1063 TS0, RP0, RS0, ERaT0, |
1062 TS0, RP0, RS0, ERaT0, |
1064 RotaT0, RetT0, ERaR0, |
1063 RotaT0, RetT0, ERaR0, |
1065 RotaR0, RetR0, LDRCal) |
1064 RotaR0, RetR0, LDRCal) |
1671 else: |
1670 else: |
1672 LDRsimx = 1./LDRsim |
1671 LDRsimx = 1./LDRsim |
1673 LDRsimx2 = 1./LDRsim2 |
1672 LDRsimx2 = 1./LDRsim2 |
1674 ''' |
1673 ''' |
1675 # ----- Backward correction |
1674 # ----- Backward correction |
1676 # Corrected LDRCorr from forward simulated LDRsim (atrue) with assumed true G0,H0,K0 |
1675 # Corrected LDRCorr with assumed true G0,H0,K0 from forward simulated (real) LDRsim (atrue) |
1677 LDRCorr = (LDRsim * K0 / Etax * (GT0 + HT0) - (GR0 + HR0)) / ( |
1676 LDRCorr = (LDRsim * K0 / Etax * (GT0 + HT0) - (GR0 + HR0)) / ( |
1678 (GR0 - HR0) - LDRsim * K0 / Etax * (GT0 - HT0)) |
1677 (GR0 - HR0) - LDRsim * K0 / Etax * (GT0 - HT0)) |
1679 |
1678 ''' |
1680 # -- F11corr from It and Ir and calibration EtaX |
1679 # -- F11corr from It and Ir and calibration EtaX |
1681 Text1 = "!!! EXPERIMENTAL !!! F11corr from It and Ir with calibration EtaX: x-axis: F11corr(LDRtrue) / F11corr(LDRtrue = 0.004) - 1" |
1680 Text1 = "!!! EXPERIMENTAL !!! F11corr from It and Ir with calibration EtaX: x-axis: F11corr(LDRtrue) / F11corr(LDRtrue = 0.004) - 1" |
1682 F11corr = 1 / (TiO * TiE) * ( |
1681 F11corr = 1 / (TiO * TiE) * ( |
1683 (HR0 * Etax / K0 * It / TTa - HT0 * Ir / TRa) / (HR0 * GT0 - HT0 * GR0)) # IL = 1 Eq.(64) |
1682 (HR0 * Etax / K0 * It / TTa - HT0 * Ir / TRa) / (HR0 * GT0 - HT0 * GR0)) # IL = 1 Eq.(64); Etax/K0 = Eta0. |
|
1683 ''' |
|
1684 # Corrected F11corr with assumed true G0,H0,K0 from forward simulated (real) It and Ir (atrue) |
|
1685 Text1 = "!!! EXPERIMENTAL !!! F11corr from real It and Ir with real calibration EtaX: x-axis: F11corr(LDRtrue) / aF11sim0(LDRtrue) - 1" |
|
1686 F11corr = 1 / (TiO * TiE) * ( |
|
1687 (HR0 * Etax / K0 * It / TTa - HT0 * Ir / TRa) / (HR0 * GT0 - HT0 * GR0)) # IL = 1 Eq.(64); Etax/K0 = Eta0. |
1684 |
1688 |
1685 # Text1 = "F11corr from It and Ir without corrections but with calibration EtaX: x-axis: F11corr(LDRtrue) devided by F11corr(LDRtrue = 0.004)" |
1689 # Text1 = "F11corr from It and Ir without corrections but with calibration EtaX: x-axis: F11corr(LDRtrue) devided by F11corr(LDRtrue = 0.004)" |
1686 # F11corr = 0.5/(TiO*TiE)*(Etax*It/TTa+Ir/TRa) # IL = 1 Eq.(64) |
1690 # F11corr = 0.5/(TiO*TiE)*(Etax*It/TTa+Ir/TRa) # IL = 1 Eq.(64) |
1687 |
1691 |
1688 # -- It from It only with atrue without corrections - for BERTHA (and PollyXTs) |
1692 # -- It from It only with atrue without corrections - for BERTHA (and PollyXTs) |
1804 if (nRotaR > 0): PlotSubHist("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) |
1808 if (nRotaR > 0): PlotSubHist("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) |
1805 if (nLDRCal > 0): PlotSubHist("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) |
1809 if (nLDRCal > 0): PlotSubHist("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) |
1806 |
1810 |
1807 plt.show() |
1811 plt.show() |
1808 plt.close |
1812 plt.close |
|
1813 |
1809 ''' |
1814 ''' |
|
1815 # --- Plot F11 histograms |
1810 print() |
1816 print() |
1811 #print("IT(LDRtrue) devided by IT(LDRtrue = 0.004)") |
|
1812 print(" ############################################################################## ") |
1817 print(" ############################################################################## ") |
1813 print(Text1) |
1818 print(Text1) |
1814 print() |
1819 print() |
1815 |
1820 |
1816 iLDR = 5 |
1821 iLDR = 5 |
1817 for LDRTrue in LDRrange: |
1822 for LDRTrue in LDRrange: |
1818 iLDR = iLDR - 1 |
1823 iLDR = iLDR - 1 |
1819 aF11corr[iLDR,:] = aF11corr[iLDR,:] / aF11corr[0,:] - 1.0 |
1824 #aF11corr[iLDR,:] = aF11corr[iLDR,:] / aF11corr[0,:] - 1.0 |
1820 |
1825 aF11corr[iLDR,:] = aF11corr[iLDR,:] / aF11sim0[iLDR] - 1.0 |
1821 # Plot F11 |
1826 # Plot F11 |
1822 def PlotSubHistF11(aVar, aX, X0, daX, iaX, naX): |
1827 def PlotSubHistF11(aVar, aX, X0, daX, iaX, naX): |
1823 fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) |
1828 fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) |
1824 iLDR = -1 |
1829 iLDR = -1 |
1825 for LDRTrue in LDRrange: |
1830 for LDRTrue in LDRrange: |
1880 if (nRotaR > 0): PlotSubHistF11("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) |
1885 if (nRotaR > 0): PlotSubHistF11("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) |
1881 if (nLDRCal > 0): PlotSubHistF11("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) |
1886 if (nLDRCal > 0): PlotSubHistF11("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) |
1882 |
1887 |
1883 plt.show() |
1888 plt.show() |
1884 plt.close |
1889 plt.close |
|
1890 |
1885 ''' |
1891 ''' |
1886 |
|
1887 ''' |
1892 ''' |
1888 # only histogram |
1893 # only histogram |
1889 #print("******************* " + aVar + " *******************") |
1894 #print("******************* " + aVar + " *******************") |
1890 fig, ax = plt.subplots(nrows=5, ncols=2, sharex=True, sharey=True, figsize=(10, 10)) |
1895 fig, ax = plt.subplots(nrows=5, ncols=2, sharex=True, sharey=True, figsize=(10, 10)) |
1891 iLDR = -1 |
1896 iLDR = -1 |
1901 bins=200, log=False, alpha=0.2, normed=False, color = '0.5', histtype='stepfilled') |
1906 bins=200, log=False, alpha=0.2, normed=False, color = '0.5', histtype='stepfilled') |
1902 plt.tick_params(axis='both', labelsize=9) |
1907 plt.tick_params(axis='both', labelsize=9) |
1903 plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2) |
1908 plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2) |
1904 plt.show() |
1909 plt.show() |
1905 plt.close |
1910 plt.close |
|
1911 # --- End of Plot F11 histograms |
1906 ''' |
1912 ''' |
1907 |
1913 |
1908 # --- Plot LDRmin, LDRmax |
1914 # --- Plot LDRmin, LDRmax |
1909 fig2 = plt.figure() |
1915 fig2 = plt.figure() |
1910 plt.plot(LDRrange, LDRmax - LDRrange, linewidth=2.0, color='b') |
1916 plt.plot(LDRrange, LDRmax - LDRrange, linewidth=2.0, color='b') |