Mercurial > public > atmospheric_lidar_ghk / file comparison
comparison: lidar_correction_ghk.py
lidar_correction_ghk.py
- changeset 16
- 313ac320b970
- parent 14
- 82dba9904149
- child 19
- 40d55af749b6
equal
deleted
inserted
replaced
| 15:94eac33c6e6e | 16:313ac320b970 |
|---|---|
| 81 | 81 |
| 82 #PrintToOutputFile = True | 82 #PrintToOutputFile = True |
| 83 | 83 |
| 84 sqr05 = 0.5**0.5 | 84 sqr05 = 0.5**0.5 |
| 85 | 85 |
| 86 # Do you want to calculate the errors? If not, just the GHK-parameters are determined. | |
| 87 Error_Calc = True | |
| 86 # ---- Initial definition of variables; the actual values will be read in with exec(open('./optic_input.py').read()) below | 88 # ---- Initial definition of variables; the actual values will be read in with exec(open('./optic_input.py').read()) below |
| 87 LID = "internal" | 89 LID = "internal" |
| 88 EID = "internal" | 90 EID = "internal" |
| 89 # --- IL Laser IL and +-Uncertainty | 91 # --- IL Laser IL and +-Uncertainty |
| 90 bL = 1. #degree of linear polarization; default 1 | 92 bL = 1. #degree of linear polarization; default 1 |
| 170 | 172 |
| 171 # end of initial definition of variables | 173 # end of initial definition of variables |
| 172 # ******************************************************************************************************************************* | 174 # ******************************************************************************************************************************* |
| 173 | 175 |
| 174 # --- Read actual lidar system parameters from ./optic_input.py (must be in the same directory) | 176 # --- Read actual lidar system parameters from ./optic_input.py (must be in the same directory) |
| 175 | |
| 176 InputFile = 'optic_input_example_lidar.py' | 177 InputFile = 'optic_input_example_lidar.py' |
| 177 #InputFile = 'optic_input_ver6e_POLIS_355.py' | 178 #InputFile = 'optic_input_ver6e_POLIS_355.py' |
| 178 #InputFile = 'optic_input_ver6e_POLIS_355_JA.py' | 179 #InputFile = 'optic_input_ver6e_POLIS_355_JA.py' |
| 179 #InputFile = 'optic_input_ver6c_POLIS_532.py' | 180 #InputFile = 'optic_input_ver6c_POLIS_532.py' |
| 180 #InputFile = 'optic_input_ver6e_POLIS_532.py' | 181 #InputFile = 'optic_input_ver6e_POLIS_532.py' |
| 841 print(" Diatt., Tunpol, Retard., Rotation (deg)") | 842 print(" Diatt., Tunpol, Retard., Rotation (deg)") |
| 842 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("Emitter ", DiE0, dDiE, TiE, RetE0, dRetE, RotE0, dRotE, nDiE, nRetE, nRotE)) | 843 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("Emitter ", DiE0, dDiE, TiE, RetE0, dRetE, RotE0, dRotE, nDiE, nRetE, nRotE)) |
| 843 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)) | 844 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)) |
| 844 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)) | 845 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)) |
| 845 print("{0:12}".format(" --- Pol.-filter ---")) | 846 print("{0:12}".format(" --- Pol.-filter ---")) |
| 846 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)) | 847 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)) |
| 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)) | 848 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)) |
| 848 print("{0:12}".format(" --- PBS ---")) | 849 print("{0:12}".format(" --- PBS ---")) |
| 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)) | 850 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)) |
| 851 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)) | |
| 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)) | 852 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)) |
| 851 print("{0:12}".format(" --- Combined PBS + Pol.-filter ---")) | 853 print("{0:12}".format(" --- Combined PBS + Pol.-filter ---")) |
| 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)) | 854 print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}".format("DT,TT,DR,TR :", DTa0, TTa0, DRa0, TRa0)) |
| 853 print() | 855 print() |
| 854 print("Rotation Error Epsilon For Normal Measurements = ", RotationErrorEpsilonForNormalMeasurements) | 856 print("Rotation Error Epsilon For Normal Measurements = ", RotationErrorEpsilonForNormalMeasurements) |
| 864 #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)) | 866 #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)) |
| 865 #print("{0:8} |{1:8}->{2:8}->{3:8}".format(" LDRCal"," LDRtrue", " LDRsimx", " LDRcorr")) | 867 #print("{0:8} |{1:8}->{2:8}->{3:8}".format(" LDRCal"," LDRtrue", " LDRsimx", " LDRcorr")) |
| 866 #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)) | 868 #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)) |
| 867 #print("{0:8} |{1:8}->{2:8}->{3:8}".format(" LDRCal"," LDRtrue", " LDRsimx", " LDRcorr")) | 869 #print("{0:8} |{1:8}->{2:8}->{3:8}".format(" LDRCal"," LDRtrue", " LDRsimx", " LDRcorr")) |
| 868 #print(" --- LDRCal during calibration") | 870 #print(" --- LDRCal during calibration") |
| 869 print("{0:26}: {1:6.3f}±{2:5.3f}/{3:2d}".format("LDRCal during calibration", LDRCal0, dLDRCal, nLDRCal)) | 871 print("{0:26}: {1:6.3f}±{2:5.3f}/{3:2d}".format("LDRCal during calibration in calibration range", LDRCal0, dLDRCal, nLDRCal)) |
| 870 | 872 |
| 871 #print("{0:8}={1:8.5f};{2:8}={3:8.5f}".format(" IinP",IinP," F11sim",F11sim)) | 873 #print("{0:8}={1:8.5f};{2:8}={3:8.5f}".format(" IinP",IinP," F11sim",F11sim)) |
| 872 print() | 874 print() |
| 873 | 875 |
| 874 K0List = np.zeros(3) | 876 K0List = np.zeros(3) |
| 906 if(PrintToOutputFile): | 908 if(PrintToOutputFile): |
| 907 sys.stdout.flush() | 909 sys.stdout.flush() |
| 908 f.close | 910 f.close |
| 909 sys.stdout = old_target | 911 sys.stdout = old_target |
| 910 ''' | 912 ''' |
| 911 # --- CALC again truth with LDRCal0 to reset all 0-values | 913 if(Error_Calc): |
| 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) | 914 # --- CALC again truth with LDRCal0 to reset all 0-values |
| 913 | 915 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) |
| 914 # --- Start Errors calculation with variable parameters ------------------------------------------------------------------ | 916 |
| 915 | 917 # --- Start Errors calculation with variable parameters ------------------------------------------------------------------ |
| 916 iN = -1 | 918 |
| 917 N = ((nRotL*2+1)* | 919 iN = -1 |
| 918 (nRotE*2+1)*(nRetE*2+1)*(nDiE*2+1)* | 920 N = ((nRotL*2+1)* |
| 919 (nRotO*2+1)*(nRetO*2+1)*(nDiO*2+1)* | 921 (nRotE*2+1)*(nRetE*2+1)*(nDiE*2+1)* |
| 920 (nRotC*2+1)*(nRetC*2+1)*(nDiC*2+1)* | 922 (nRotO*2+1)*(nRetO*2+1)*(nDiO*2+1)* |
| 921 (nTP*2+1)*(nTS*2+1)*(nRP*2+1)*(nRS*2+1)*(nERaT*2+1)*(nERaR*2+1)* | 923 (nRotC*2+1)*(nRetC*2+1)*(nDiC*2+1)* |
| 922 (nRotaT*2+1)*(nRotaR*2+1)*(nRetT*2+1)*(nRetR*2+1)*(nLDRCal*2+1)) | 924 (nTP*2+1)*(nTS*2+1)*(nRP*2+1)*(nRS*2+1)*(nERaT*2+1)*(nERaR*2+1)* |
| 923 print("N = ",N ," ", end="") | 925 (nRotaT*2+1)*(nRotaR*2+1)*(nRetT*2+1)*(nRetR*2+1)*(nLDRCal*2+1)) |
| 924 | 926 print("N = ",N ," ", end="") |
| 925 if N > 1e6: | 927 |
| 926 if user_yes_no_query('Warning: processing ' + str(N) + ' samples will take very long. Do you want to proceed?') == 0: sys.exit() | 928 if N > 1e6: |
| 927 if N > 5e6: | 929 if user_yes_no_query('Warning: processing ' + str(N) + ' samples will take very long. Do you want to proceed?') == 0: sys.exit() |
| 928 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() | 930 if N > 5e6: |
| 929 | 931 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() |
| 930 #if user_yes_no_query('Warning: processing' + str(N) + ' samples will take very long. Do you want to proceed?') == 0: sys.exit() | 932 |
| 931 | 933 #if user_yes_no_query('Warning: processing' + str(N) + ' samples will take very long. Do you want to proceed?') == 0: sys.exit() |
| 932 # --- Arrays for plotting ------ | 934 |
| 933 LDRmin = np.zeros(5) | 935 # --- Arrays for plotting ------ |
| 934 LDRmax = np.zeros(5) | 936 LDRmin = np.zeros(5) |
| 935 F11min = np.zeros(5) | 937 LDRmax = np.zeros(5) |
| 936 F11max = np.zeros(5) | 938 F11min = np.zeros(5) |
| 937 | 939 F11max = np.zeros(5) |
| 938 LDRrange = np.zeros(5) | 940 |
| 939 LDRrange = 0.004, 0.02, 0.1, 0.3, 0.45 | 941 LDRrange = np.zeros(5) |
| 940 #aLDRsimx = np.zeros(N) | 942 LDRrange = 0.004, 0.02, 0.1, 0.3, 0.45 |
| 941 #aLDRsimx2 = np.zeros(N) | 943 #aLDRsimx = np.zeros(N) |
| 942 #aLDRcorr = np.zeros(N) | 944 #aLDRsimx2 = np.zeros(N) |
| 943 #aLDRcorr2 = np.zeros(N) | 945 #aLDRcorr = np.zeros(N) |
| 944 aERaT = np.zeros(N) | 946 #aLDRcorr2 = np.zeros(N) |
| 945 aERaR = np.zeros(N) | 947 aERaT = np.zeros(N) |
| 946 aRotaT = np.zeros(N) | 948 aERaR = np.zeros(N) |
| 947 aRotaR = np.zeros(N) | 949 aRotaT = np.zeros(N) |
| 948 aRetT = np.zeros(N) | 950 aRotaR = np.zeros(N) |
| 949 aRetR = np.zeros(N) | 951 aRetT = np.zeros(N) |
| 950 aTP = np.zeros(N) | 952 aRetR = np.zeros(N) |
| 951 aTS = np.zeros(N) | 953 aTP = np.zeros(N) |
| 952 aRP = np.zeros(N) | 954 aTS = np.zeros(N) |
| 953 aRS = np.zeros(N) | 955 aRP = np.zeros(N) |
| 954 aDiE = np.zeros(N) | 956 aRS = np.zeros(N) |
| 955 aDiO = np.zeros(N) | 957 aDiE = np.zeros(N) |
| 956 aDiC = np.zeros(N) | 958 aDiO = np.zeros(N) |
| 957 aRotC = np.zeros(N) | 959 aDiC = np.zeros(N) |
| 958 aRetC = np.zeros(N) | 960 aRotC = np.zeros(N) |
| 959 aRotL = np.zeros(N) | 961 aRetC = np.zeros(N) |
| 960 aRetE = np.zeros(N) | 962 aRotL = np.zeros(N) |
| 961 aRotE = np.zeros(N) | 963 aRetE = np.zeros(N) |
| 962 aRetO = np.zeros(N) | 964 aRotE = np.zeros(N) |
| 963 aRotO = np.zeros(N) | 965 aRetO = np.zeros(N) |
| 964 aLDRCal = np.zeros(N) | 966 aRotO = np.zeros(N) |
| 965 aA = np.zeros((5,N)) | 967 aLDRCal = np.zeros(N) |
| 966 aX = np.zeros((5,N)) | 968 aA = np.zeros((5,N)) |
| 967 aF11corr = np.zeros((5,N)) | 969 aX = np.zeros((5,N)) |
| 968 | 970 aF11corr = np.zeros((5,N)) |
| 969 atime = clock() | 971 |
| 970 dtime = clock() | 972 atime = clock() |
| 971 | 973 dtime = clock() |
| 972 # --- Calc Error signals | 974 |
| 973 #GT, HT, GR, HR, K, Eta, LDRsim = Calc(RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC, TP, TS) | 975 # --- Calc Error signals |
| 974 # ---- Do the calculations of bra-ket vectors | 976 #GT, HT, GR, HR, K, Eta, LDRsim = Calc(RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC, TP, TS) |
| 975 h = -1. if TypeC == 2 else 1 | 977 # ---- Do the calculations of bra-ket vectors |
| 976 | 978 h = -1. if TypeC == 2 else 1 |
| 977 # from input file: measured LDRm and true LDRtrue, LDRtrue2 => | 979 |
| 978 ameas = (1.-LDRmeas)/(1+LDRmeas) | 980 # from input file: measured LDRm and true LDRtrue, LDRtrue2 => |
| 979 atrue = (1.-LDRtrue)/(1+LDRtrue) | 981 ameas = (1.-LDRmeas)/(1+LDRmeas) |
| 980 atrue2 = (1.-LDRtrue2)/(1+LDRtrue2) | 982 atrue = (1.-LDRtrue)/(1+LDRtrue) |
| 981 | 983 atrue2 = (1.-LDRtrue2)/(1+LDRtrue2) |
| 982 for iLDRCal in range(-nLDRCal,nLDRCal+1): | 984 |
| 983 # from input file: assumed LDRCal for calibration measurements | 985 for iLDRCal in range(-nLDRCal,nLDRCal+1): |
| 984 LDRCal = LDRCal0 | 986 # from input file: assumed LDRCal for calibration measurements |
| 985 if nLDRCal > 0: LDRCal = LDRCal0 + iLDRCal*dLDRCal/nLDRCal | 987 LDRCal = LDRCal0 |
| 986 | 988 if nLDRCal > 0: LDRCal = LDRCal0 + iLDRCal*dLDRCal/nLDRCal |
| 987 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) | 989 |
| 988 aCal = (1.-LDRCal)/(1+LDRCal) | 990 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) |
| 989 for iRotL, iRotE, iRetE, iDiE \ | 991 aCal = (1.-LDRCal)/(1+LDRCal) |
| 990 in [(iRotL,iRotE,iRetE,iDiE) | 992 for iRotL, iRotE, iRetE, iDiE \ |
| 991 for iRotL in range(-nRotL,nRotL+1) | 993 in [(iRotL,iRotE,iRetE,iDiE) |
| 992 for iRotE in range(-nRotE,nRotE+1) | 994 for iRotL in range(-nRotL,nRotL+1) |
| 993 for iRetE in range(-nRetE,nRetE+1) | 995 for iRotE in range(-nRotE,nRotE+1) |
| 994 for iDiE in range(-nDiE,nDiE+1)]: | 996 for iRetE in range(-nRetE,nRetE+1) |
| 995 | 997 for iDiE in range(-nDiE,nDiE+1)]: |
| 996 if nRotL > 0: RotL = RotL0 + iRotL*dRotL/nRotL | 998 |
| 997 if nRotE > 0: RotE = RotE0 + iRotE*dRotE/nRotE | 999 if nRotL > 0: RotL = RotL0 + iRotL*dRotL/nRotL |
| 998 if nRetE > 0: RetE = RetE0 + iRetE*dRetE/nRetE | 1000 if nRotE > 0: RotE = RotE0 + iRotE*dRotE/nRotE |
| 999 if nDiE > 0: DiE = DiE0 + iDiE*dDiE/nDiE | 1001 if nRetE > 0: RetE = RetE0 + iRetE*dRetE/nRetE |
| 1000 | 1002 if nDiE > 0: DiE = DiE0 + iDiE*dDiE/nDiE |
| 1001 # angles of emitter and laser and calibrator and receiver optics | 1003 |
| 1002 # RotL = alpha, RotE = beta, RotO = gamma, RotC = epsilon | 1004 # angles of emitter and laser and calibrator and receiver optics |
| 1003 S2a = np.sin(2*np.deg2rad(RotL)) | 1005 # RotL = alpha, RotE = beta, RotO = gamma, RotC = epsilon |
| 1004 C2a = np.cos(2*np.deg2rad(RotL)) | 1006 S2a = np.sin(2*np.deg2rad(RotL)) |
| 1005 S2b = np.sin(2*np.deg2rad(RotE)) | 1007 C2a = np.cos(2*np.deg2rad(RotL)) |
| 1006 C2b = np.cos(2*np.deg2rad(RotE)) | 1008 S2b = np.sin(2*np.deg2rad(RotE)) |
| 1007 S2ab = np.sin(np.deg2rad(2*RotL-2*RotE)) | 1009 C2b = np.cos(2*np.deg2rad(RotE)) |
| 1008 C2ab = np.cos(np.deg2rad(2*RotL-2*RotE)) | 1010 S2ab = np.sin(np.deg2rad(2*RotL-2*RotE)) |
| 1009 | 1011 C2ab = np.cos(np.deg2rad(2*RotL-2*RotE)) |
| 1010 # Laser with Degree of linear polarization DOLP = bL | 1012 |
| 1011 IinL = 1. | 1013 # Laser with Degree of linear polarization DOLP = bL |
| 1012 QinL = bL | 1014 IinL = 1. |
| 1013 UinL = 0. | 1015 QinL = bL |
| 1014 VinL = (1. - bL**2)**0.5 | 1016 UinL = 0. |
| 1015 | 1017 VinL = (1. - bL**2)**0.5 |
| 1016 # Stokes Input Vector rotation Eq. E.4 | 1018 |
| 1017 A = C2a*QinL - S2a*UinL | 1019 # Stokes Input Vector rotation Eq. E.4 |
| 1018 B = S2a*QinL + C2a*UinL | 1020 A = C2a*QinL - S2a*UinL |
| 1019 # Stokes Input Vector rotation Eq. E.9 | 1021 B = S2a*QinL + C2a*UinL |
| 1020 C = C2ab*QinL - S2ab*UinL | 1022 # Stokes Input Vector rotation Eq. E.9 |
| 1021 D = S2ab*QinL + C2ab*UinL | 1023 C = C2ab*QinL - S2ab*UinL |
| 1022 | 1024 D = S2ab*QinL + C2ab*UinL |
| 1023 # emitter optics | 1025 |
| 1024 CosE = np.cos(np.deg2rad(RetE)) | 1026 # emitter optics |
| 1025 SinE = np.sin(np.deg2rad(RetE)) | 1027 CosE = np.cos(np.deg2rad(RetE)) |
| 1026 ZiE = (1. - DiE**2)**0.5 | 1028 SinE = np.sin(np.deg2rad(RetE)) |
| 1027 WiE = (1. - ZiE*CosE) | 1029 ZiE = (1. - DiE**2)**0.5 |
| 1028 | 1030 WiE = (1. - ZiE*CosE) |
| 1029 # Stokes Input Vector after emitter optics equivalent to Eq. E.9 with already rotated input vector from Eq. E.4 | 1031 |
| 1030 # b = beta | 1032 # Stokes Input Vector after emitter optics equivalent to Eq. E.9 with already rotated input vector from Eq. E.4 |
| 1031 IinE = (IinL + DiE*C) | 1033 # b = beta |
| 1032 QinE = (C2b*DiE*IinL + A + S2b*(WiE*D - ZiE*SinE*VinL)) | 1034 IinE = (IinL + DiE*C) |
| 1033 UinE = (S2b*DiE*IinL + B - C2b*(WiE*D - ZiE*SinE*VinL)) | 1035 QinE = (C2b*DiE*IinL + A + S2b*(WiE*D - ZiE*SinE*VinL)) |
| 1034 VinE = (-ZiE*SinE*D + ZiE*CosE*VinL) | 1036 UinE = (S2b*DiE*IinL + B - C2b*(WiE*D - ZiE*SinE*VinL)) |
| 1035 | 1037 VinE = (-ZiE*SinE*D + ZiE*CosE*VinL) |
| 1036 #------------------------- | 1038 |
| 1037 # F11 assuemd to be = 1 => measured: F11m = IinP / IinE with atrue | 1039 #------------------------- |
| 1038 #F11sim = (IinE + DiO*atrue*(C2g*QinE - S2g*UinE))/IinE | 1040 # F11 assuemd to be = 1 => measured: F11m = IinP / IinE with atrue |
| 1039 #------------------------- | 1041 #F11sim = (IinE + DiO*atrue*(C2g*QinE - S2g*UinE))/IinE |
| 1040 | 1042 #------------------------- |
| 1041 for iRotO, iRetO, iDiO, iRotC, iRetC, iDiC, iTP, iTS, iRP, iRS, iERaT, iRotaT, iRetT, iERaR, iRotaR, iRetR \ | 1043 |
| 1042 in [(iRotO,iRetO,iDiO,iRotC,iRetC,iDiC,iTP,iTS,iRP,iRS,iERaT,iRotaT,iRetT,iERaR,iRotaR,iRetR ) | 1044 for iRotO, iRetO, iDiO, iRotC, iRetC, iDiC, iTP, iTS, iRP, iRS, iERaT, iRotaT, iRetT, iERaR, iRotaR, iRetR \ |
| 1043 for iRotO in range(-nRotO,nRotO+1) | 1045 in [(iRotO,iRetO,iDiO,iRotC,iRetC,iDiC,iTP,iTS,iRP,iRS,iERaT,iRotaT,iRetT,iERaR,iRotaR,iRetR ) |
| 1044 for iRetO in range(-nRetO,nRetO+1) | 1046 for iRotO in range(-nRotO,nRotO+1) |
| 1045 for iDiO in range(-nDiO,nDiO+1) | 1047 for iRetO in range(-nRetO,nRetO+1) |
| 1046 for iRotC in range(-nRotC,nRotC+1) | 1048 for iDiO in range(-nDiO,nDiO+1) |
| 1047 for iRetC in range(-nRetC,nRetC+1) | 1049 for iRotC in range(-nRotC,nRotC+1) |
| 1048 for iDiC in range(-nDiC,nDiC+1) | 1050 for iRetC in range(-nRetC,nRetC+1) |
| 1049 for iTP in range(-nTP,nTP+1) | 1051 for iDiC in range(-nDiC,nDiC+1) |
| 1050 for iTS in range(-nTS,nTS+1) | 1052 for iTP in range(-nTP,nTP+1) |
| 1051 for iRP in range(-nRP,nRP+1) | 1053 for iTS in range(-nTS,nTS+1) |
| 1052 for iRS in range(-nRS,nRS+1) | 1054 for iRP in range(-nRP,nRP+1) |
| 1053 for iERaT in range(-nERaT,nERaT+1) | 1055 for iRS in range(-nRS,nRS+1) |
| 1054 for iRotaT in range(-nRotaT,nRotaT+1) | 1056 for iERaT in range(-nERaT,nERaT+1) |
| 1055 for iRetT in range(-nRetT,nRetT+1) | 1057 for iRotaT in range(-nRotaT,nRotaT+1) |
| 1056 for iERaR in range(-nERaR,nERaR+1) | 1058 for iRetT in range(-nRetT,nRetT+1) |
| 1057 for iRotaR in range(-nRotaR,nRotaR+1) | 1059 for iERaR in range(-nERaR,nERaR+1) |
| 1058 for iRetR in range(-nRetR,nRetR+1)]: | 1060 for iRotaR in range(-nRotaR,nRotaR+1) |
| 1059 | 1061 for iRetR in range(-nRetR,nRetR+1)]: |
| 1060 iN = iN + 1 | 1062 |
| 1061 if (iN == 10001): | 1063 iN = iN + 1 |
| 1064 if (iN == 10001): | |
| 1065 ctime = clock() | |
| 1066 print(" estimated time ", "{0:4.2f}".format(N/10000 * (ctime-atime)), "sec ") #, end="") | |
| 1067 print("\r elapsed time ", "{0:5.0f}".format((ctime-atime)), "sec ", end="\r") | |
| 1062 ctime = clock() | 1068 ctime = clock() |
| 1063 print(" estimated time ", "{0:4.2f}".format(N/10000 * (ctime-atime)), "sec ") #, end="") | 1069 if ((ctime - dtime) > 10): |
| 1064 print("\r elapsed time ", "{0:5.0f}".format((ctime-atime)), "sec ", end="\r") | 1070 print("\r elapsed time ", "{0:5.0f}".format((ctime-atime)), "sec ", end="\r") |
| 1065 ctime = clock() | 1071 dtime = ctime |
| 1066 if ((ctime - dtime) > 10): | 1072 |
| 1067 print("\r elapsed time ", "{0:5.0f}".format((ctime-atime)), "sec ", end="\r") | 1073 if nRotO > 0: RotO = RotO0 + iRotO*dRotO/nRotO |
| 1068 dtime = ctime | 1074 if nRetO > 0: RetO = RetO0 + iRetO*dRetO/nRetO |
| 1069 | 1075 if nDiO > 0: DiO = DiO0 + iDiO*dDiO/nDiO |
| 1070 if nRotO > 0: RotO = RotO0 + iRotO*dRotO/nRotO | 1076 if nRotC > 0: RotC = RotC0 + iRotC*dRotC/nRotC |
| 1071 if nRetO > 0: RetO = RetO0 + iRetO*dRetO/nRetO | 1077 if nRetC > 0: RetC = RetC0 + iRetC*dRetC/nRetC |
| 1072 if nDiO > 0: DiO = DiO0 + iDiO*dDiO/nDiO | 1078 if nDiC > 0: DiC = DiC0 + iDiC*dDiC/nDiC |
| 1073 if nRotC > 0: RotC = RotC0 + iRotC*dRotC/nRotC | 1079 if nTP > 0: TP = TP0 + iTP*dTP/nTP |
| 1074 if nRetC > 0: RetC = RetC0 + iRetC*dRetC/nRetC | 1080 if nTS > 0: TS = TS0 + iTS*dTS/nTS |
| 1075 if nDiC > 0: DiC = DiC0 + iDiC*dDiC/nDiC | 1081 if nRP > 0: RP = RP0 + iRP*dRP/nRP |
| 1076 if nTP > 0: TP = TP0 + iTP*dTP/nTP | 1082 if nRS > 0: RS = RS0 + iRS*dRS/nRS |
| 1077 if nTS > 0: TS = TS0 + iTS*dTS/nTS | 1083 if nERaT > 0: ERaT = ERaT0 + iERaT*dERaT/nERaT |
| 1078 if nRP > 0: RP = RP0 + iRP*dRP/nRP | 1084 if nRotaT > 0:RotaT= RotaT0+ iRotaT*dRotaT/nRotaT |
| 1079 if nRS > 0: RS = RS0 + iRS*dRS/nRS | 1085 if nRetT > 0: RetT = RetT0 + iRetT*dRetT/nRetT |
| 1080 if nERaT > 0: ERaT = ERaT0 + iERaT*dERaT/nERaT | 1086 if nERaR > 0: ERaR = ERaR0 + iERaR*dERaR/nERaR |
| 1081 if nRotaT > 0:RotaT= RotaT0+ iRotaT*dRotaT/nRotaT | 1087 if nRotaR > 0:RotaR= RotaR0+ iRotaR*dRotaR/nRotaR |
| 1082 if nRetT > 0: RetT = RetT0 + iRetT*dRetT/nRetT | 1088 if nRetR > 0: RetR = RetR0 + iRetR*dRetR/nRetR |
| 1083 if nERaR > 0: ERaR = ERaR0 + iERaR*dERaR/nERaR | 1089 |
| 1084 if nRotaR > 0:RotaR= RotaR0+ iRotaR*dRotaR/nRotaR | 1090 #print("{0:5.2f}, {1:5.2f}, {2:5.2f}, {3:10d}".format(RotL, RotE, RotO, iN)) |
| 1085 if nRetR > 0: RetR = RetR0 + iRetR*dRetR/nRetR | 1091 |
| 1086 | 1092 # receiver optics |
| 1087 #print("{0:5.2f}, {1:5.2f}, {2:5.2f}, {3:10d}".format(RotL, RotE, RotO, iN)) | 1093 CosO = np.cos(np.deg2rad(RetO)) |
| 1088 | 1094 SinO = np.sin(np.deg2rad(RetO)) |
| 1089 # receiver optics | 1095 ZiO = (1. - DiO**2)**0.5 |
| 1090 CosO = np.cos(np.deg2rad(RetO)) | 1096 WiO = (1. - ZiO*CosO) |
| 1091 SinO = np.sin(np.deg2rad(RetO)) | 1097 S2g = np.sin(np.deg2rad(2*RotO)) |
| 1092 ZiO = (1. - DiO**2)**0.5 | 1098 C2g = np.cos(np.deg2rad(2*RotO)) |
| 1093 WiO = (1. - ZiO*CosO) | 1099 # calibrator |
| 1094 S2g = np.sin(np.deg2rad(2*RotO)) | 1100 CosC = np.cos(np.deg2rad(RetC)) |
| 1095 C2g = np.cos(np.deg2rad(2*RotO)) | 1101 SinC = np.sin(np.deg2rad(RetC)) |
| 1096 # calibrator | 1102 ZiC = (1. - DiC**2)**0.5 |
| 1097 CosC = np.cos(np.deg2rad(RetC)) | 1103 WiC = (1. - ZiC*CosC) |
| 1098 SinC = np.sin(np.deg2rad(RetC)) | 1104 |
| 1099 ZiC = (1. - DiC**2)**0.5 | 1105 # analyser |
| 1100 WiC = (1. - ZiC*CosC) | 1106 #For POLLY_XTs |
| 1101 | 1107 if(RS_RP_depend_on_TS_TP): |
| 1102 # analyser | 1108 RS = 1 - TS |
| 1103 #For POLLY_XTs | 1109 RP = 1 - TP |
| 1104 if(RS_RP_depend_on_TS_TP): | 1110 TiT = 0.5 * (TP + TS) |
| 1105 RS = 1 - TS | 1111 DiT = (TP-TS)/(TP+TS) |
| 1106 RP = 1 - TP | 1112 ZiT = (1. - DiT**2)**0.5 |
| 1107 TiT = 0.5 * (TP + TS) | 1113 TiR = 0.5 * (RP + RS) |
| 1108 DiT = (TP-TS)/(TP+TS) | 1114 DiR = (RP-RS)/(RP+RS) |
| 1109 ZiT = (1. - DiT**2)**0.5 | 1115 ZiR = (1. - DiR**2)**0.5 |
| 1110 TiR = 0.5 * (RP + RS) | 1116 CosT = np.cos(np.deg2rad(RetT)) |
| 1111 DiR = (RP-RS)/(RP+RS) | 1117 SinT = np.sin(np.deg2rad(RetT)) |
| 1112 ZiR = (1. - DiR**2)**0.5 | 1118 CosR = np.cos(np.deg2rad(RetR)) |
| 1113 CosT = np.cos(np.deg2rad(RetT)) | 1119 SinR = np.sin(np.deg2rad(RetR)) |
| 1114 SinT = np.sin(np.deg2rad(RetT)) | 1120 |
| 1115 CosR = np.cos(np.deg2rad(RetR)) | 1121 DaT = (1-ERaT)/(1+ERaT) |
| 1116 SinR = np.sin(np.deg2rad(RetR)) | 1122 DaR = (1-ERaR)/(1+ERaR) |
| 1117 | 1123 TaT = 0.5*(1+ERaT) |
| 1118 DaT = (1-ERaT)/(1+ERaT) | 1124 TaR = 0.5*(1+ERaR) |
| 1119 DaR = (1-ERaR)/(1+ERaR) | 1125 |
| 1120 TaT = 0.5*(1+ERaT) | 1126 S2aT = np.sin(np.deg2rad(h*2*RotaT)) |
| 1121 TaR = 0.5*(1+ERaR) | 1127 C2aT = np.cos(np.deg2rad(2*RotaT)) |
| 1122 | 1128 S2aR = np.sin(np.deg2rad(h*2*RotaR)) |
| 1123 S2aT = np.sin(np.deg2rad(h*2*RotaT)) | 1129 C2aR = np.cos(np.deg2rad(2*RotaR)) |
| 1124 C2aT = np.cos(np.deg2rad(2*RotaT)) | 1130 |
| 1125 S2aR = np.sin(np.deg2rad(h*2*RotaR)) | 1131 # Aanalyzer As before the PBS Eq. D.5 |
| 1126 C2aR = np.cos(np.deg2rad(2*RotaR)) | 1132 ATP1 = (1+C2aT*DaT*DiT) |
| 1127 | 1133 ATP2 = Y*(DiT+C2aT*DaT) |
| 1128 # Aanalyzer As before the PBS Eq. D.5 | 1134 ATP3 = Y*S2aT*DaT*ZiT*CosT |
| 1129 ATP1 = (1+C2aT*DaT*DiT) | 1135 ATP4 = S2aT*DaT*ZiT*SinT |
| 1130 ATP2 = Y*(DiT+C2aT*DaT) | 1136 ATP = np.array([ATP1,ATP2,ATP3,ATP4]) |
| 1131 ATP3 = Y*S2aT*DaT*ZiT*CosT | 1137 |
| 1132 ATP4 = S2aT*DaT*ZiT*SinT | 1138 ARP1 = (1+C2aR*DaR*DiR) |
| 1133 ATP = np.array([ATP1,ATP2,ATP3,ATP4]) | 1139 ARP2 = Y*(DiR+C2aR*DaR) |
| 1134 | 1140 ARP3 = Y*S2aR*DaR*ZiR*CosR |
| 1135 ARP1 = (1+C2aR*DaR*DiR) | 1141 ARP4 = S2aR*DaR*ZiR*SinR |
| 1136 ARP2 = Y*(DiR+C2aR*DaR) | 1142 ARP = np.array([ARP1,ARP2,ARP3,ARP4]) |
| 1137 ARP3 = Y*S2aR*DaR*ZiR*CosR | 1143 |
| 1138 ARP4 = S2aR*DaR*ZiR*SinR | 1144 TTa = TiT*TaT #*ATP1 |
| 1139 ARP = np.array([ARP1,ARP2,ARP3,ARP4]) | 1145 TRa = TiR*TaR #*ARP1 |
| 1140 | 1146 |
| 1141 TTa = TiT*TaT #*ATP1 | 1147 # ---- Calculate signals and correction parameters for diffeent locations and calibrators |
| 1142 TRa = TiR*TaR #*ARP1 | 1148 if LocC == 4: # Calibrator before the PBS |
| 1143 | 1149 #print("Calibrator location not implemented yet") |
| 1144 # ---- Calculate signals and correction parameters for diffeent locations and calibrators | 1150 |
| 1145 if LocC == 4: # Calibrator before the PBS | 1151 #S2ge = np.sin(np.deg2rad(2*RotO + h*2*RotC)) |
| 1146 #print("Calibrator location not implemented yet") | 1152 #C2ge = np.cos(np.deg2rad(2*RotO + h*2*RotC)) |
| 1147 | 1153 S2e = np.sin(np.deg2rad(h*2*RotC)) |
| 1148 #S2ge = np.sin(np.deg2rad(2*RotO + h*2*RotC)) | 1154 C2e = np.cos(np.deg2rad(2*RotC)) |
| 1149 #C2ge = np.cos(np.deg2rad(2*RotO + h*2*RotC)) | 1155 # rotated AinP by epsilon Eq. C.3 |
| 1150 S2e = np.sin(np.deg2rad(h*2*RotC)) | 1156 ATP2e = C2e*ATP2 + S2e*ATP3 |
| 1151 C2e = np.cos(np.deg2rad(2*RotC)) | 1157 ATP3e = C2e*ATP3 - S2e*ATP2 |
| 1152 # rotated AinP by epsilon Eq. C.3 | 1158 ARP2e = C2e*ARP2 + S2e*ARP3 |
| 1153 ATP2e = C2e*ATP2 + S2e*ATP3 | 1159 ARP3e = C2e*ARP3 - S2e*ARP2 |
| 1154 ATP3e = C2e*ATP3 - S2e*ATP2 | 1160 ATPe = np.array([ATP1,ATP2e,ATP3e,ATP4]) |
| 1155 ARP2e = C2e*ARP2 + S2e*ARP3 | 1161 ARPe = np.array([ARP1,ARP2e,ARP3e,ARP4]) |
| 1156 ARP3e = C2e*ARP3 - S2e*ARP2 | 1162 # Stokes Input Vector before the polarising beam splitter Eq. E.31 |
| 1157 ATPe = np.array([ATP1,ATP2e,ATP3e,ATP4]) | 1163 A = C2g*QinE - S2g*UinE |
| 1158 ARPe = np.array([ARP1,ARP2e,ARP3e,ARP4]) | 1164 B = S2g*QinE + C2g*UinE |
| 1159 # Stokes Input Vector before the polarising beam splitter Eq. E.31 | 1165 #C = (WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE) |
| 1160 A = C2g*QinE - S2g*UinE | 1166 Co = ZiO*SinO*VinE |
| 1161 B = S2g*QinE + C2g*UinE | 1167 Ca = (WiO*B - 2*ZiO*SinO*VinE) |
| 1162 #C = (WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE) | 1168 #C = Co + aCal*Ca |
| 1163 Co = ZiO*SinO*VinE | 1169 #IinP = (IinE + DiO*aCal*A) |
| 1164 Ca = (WiO*B - 2*ZiO*SinO*VinE) | 1170 #QinP = (C2g*DiO*IinE + aCal*QinE - S2g*C) |
| 1165 #C = Co + aCal*Ca | 1171 #UinP = (S2g*DiO*IinE - aCal*UinE + C2g*C) |
| 1166 #IinP = (IinE + DiO*aCal*A) | 1172 #VinP = (ZiO*SinO*aCal*B + ZiO*CosO*(1-2*aCal)*VinE) |
| 1167 #QinP = (C2g*DiO*IinE + aCal*QinE - S2g*C) | 1173 IinPo = IinE |
| 1168 #UinP = (S2g*DiO*IinE - aCal*UinE + C2g*C) | 1174 QinPo = (C2g*DiO*IinE - S2g*Co) |
| 1169 #VinP = (ZiO*SinO*aCal*B + ZiO*CosO*(1-2*aCal)*VinE) | 1175 UinPo = (S2g*DiO*IinE + C2g*Co) |
| 1170 IinPo = IinE | 1176 VinPo = ZiO*CosO*VinE |
| 1171 QinPo = (C2g*DiO*IinE - S2g*Co) | 1177 |
| 1172 UinPo = (S2g*DiO*IinE + C2g*Co) | 1178 IinPa = DiO*A |
| 1173 VinPo = ZiO*CosO*VinE | 1179 QinPa = QinE - S2g*Ca |
| 1174 | 1180 UinPa = -UinE + C2g*Ca |
| 1175 IinPa = DiO*A | 1181 VinPa = ZiO*(SinO*B - 2*CosO*VinE) |
| 1176 QinPa = QinE - S2g*Ca | 1182 |
| 1177 UinPa = -UinE + C2g*Ca | 1183 IinP = IinPo + aCal*IinPa |
| 1178 VinPa = ZiO*(SinO*B - 2*CosO*VinE) | 1184 QinP = QinPo + aCal*QinPa |
| 1179 | 1185 UinP = UinPo + aCal*UinPa |
| 1180 IinP = IinPo + aCal*IinPa | 1186 VinP = VinPo + aCal*VinPa |
| 1181 QinP = QinPo + aCal*QinPa | 1187 # Stokes Input Vector before the polarising beam splitter rotated by epsilon Eq. C.3 |
| 1182 UinP = UinPo + aCal*UinPa | 1188 #QinPe = C2e*QinP + S2e*UinP |
| 1183 VinP = VinPo + aCal*VinPa | 1189 #UinPe = C2e*UinP - S2e*QinP |
| 1184 # Stokes Input Vector before the polarising beam splitter rotated by epsilon Eq. C.3 | 1190 QinPoe = C2e*QinPo + S2e*UinPo |
| 1185 #QinPe = C2e*QinP + S2e*UinP | 1191 UinPoe = C2e*UinPo - S2e*QinPo |
| 1186 #UinPe = C2e*UinP - S2e*QinP | 1192 QinPae = C2e*QinPa + S2e*UinPa |
| 1187 QinPoe = C2e*QinPo + S2e*UinPo | 1193 UinPae = C2e*UinPa - S2e*QinPa |
| 1188 UinPoe = C2e*UinPo - S2e*QinPo | 1194 QinPe = C2e*QinP + S2e*UinP |
| 1189 QinPae = C2e*QinPa + S2e*UinPa | 1195 UinPe = C2e*UinP - S2e*QinP |
| 1190 UinPae = C2e*UinPa - S2e*QinPa | 1196 |
| 1191 QinPe = C2e*QinP + S2e*UinP | 1197 # Calibration signals and Calibration correction K from measurements with LDRCal / aCal |
| 1192 UinPe = C2e*UinP - S2e*QinP | 1198 if (TypeC == 2) or (TypeC == 1): # rotator calibration Eq. C.4 |
| 1193 | 1199 # parameters for calibration with aCal |
| 1194 # Calibration signals and Calibration correction K from measurements with LDRCal / aCal | 1200 AT = ATP1*IinP + h*ATP4*VinP |
| 1195 if (TypeC == 2) or (TypeC == 1): # rotator calibration Eq. C.4 | 1201 BT = ATP3e*QinP - h*ATP2e*UinP |
| 1196 # parameters for calibration with aCal | 1202 AR = ARP1*IinP + h*ARP4*VinP |
| 1197 AT = ATP1*IinP + h*ATP4*VinP | 1203 BR = ARP3e*QinP - h*ARP2e*UinP |
| 1198 BT = ATP3e*QinP - h*ATP2e*UinP | 1204 # Correction paremeters for normal measurements; they are independent of LDR |
| 1199 AR = ARP1*IinP + h*ARP4*VinP | 1205 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out |
| 1200 BR = ARP3e*QinP - h*ARP2e*UinP | 1206 IS1 = np.array([IinPo,QinPo,UinPo,VinPo]) |
| 1201 # Correction paremeters for normal measurements; they are independent of LDR | 1207 IS2 = np.array([IinPa,QinPa,UinPa,VinPa]) |
| 1202 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out | 1208 GT = np.dot(ATP,IS1) |
| 1203 IS1 = np.array([IinPo,QinPo,UinPo,VinPo]) | 1209 GR = np.dot(ARP,IS1) |
| 1204 IS2 = np.array([IinPa,QinPa,UinPa,VinPa]) | 1210 HT = np.dot(ATP,IS2) |
| 1205 GT = np.dot(ATP,IS1) | 1211 HR = np.dot(ARP,IS2) |
| 1206 GR = np.dot(ARP,IS1) | 1212 else: |
| 1207 HT = np.dot(ATP,IS2) | 1213 IS1 = np.array([IinPo,QinPo,UinPo,VinPo]) |
| 1208 HR = np.dot(ARP,IS2) | 1214 IS2 = np.array([IinPa,QinPa,UinPa,VinPa]) |
| 1215 GT = np.dot(ATPe,IS1) | |
| 1216 GR = np.dot(ARPe,IS1) | |
| 1217 HT = np.dot(ATPe,IS2) | |
| 1218 HR = np.dot(ARPe,IS2) | |
| 1219 elif (TypeC == 3) or (TypeC == 4): # linear polariser calibration Eq. C.5 | |
| 1220 # parameters for calibration with aCal | |
| 1221 AT = ATP1*IinP + ATP3e*UinPe + ZiC*CosC*(ATP2e*QinPe + ATP4*VinP) | |
| 1222 BT = DiC*(ATP1*UinPe + ATP3e*IinP) - ZiC*SinC*(ATP2e*VinP - ATP4*QinPe) | |
| 1223 AR = ARP1*IinP + ARP3e*UinPe + ZiC*CosC*(ARP2e*QinPe + ARP4*VinP) | |
| 1224 BR = DiC*(ARP1*UinPe + ARP3e*IinP) - ZiC*SinC*(ARP2e*VinP - ARP4*QinPe) | |
| 1225 # Correction paremeters for normal measurements; they are independent of LDR | |
| 1226 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out | |
| 1227 IS1 = np.array([IinPo,QinPo,UinPo,VinPo]) | |
| 1228 IS2 = np.array([IinPa,QinPa,UinPa,VinPa]) | |
| 1229 GT = np.dot(ATP,IS1) | |
| 1230 GR = np.dot(ARP,IS1) | |
| 1231 HT = np.dot(ATP,IS2) | |
| 1232 HR = np.dot(ARP,IS2) | |
| 1233 else: | |
| 1234 IS1e = np.array([IinPo+DiC*QinPoe,DiC*IinPo+QinPoe,ZiC*(CosC*UinPoe+SinC*VinPo),-ZiC*(SinC*UinPoe-CosC*VinPo)]) | |
| 1235 IS2e = np.array([IinPa+DiC*QinPae,DiC*IinPa+QinPae,ZiC*(CosC*UinPae+SinC*VinPa),-ZiC*(SinC*UinPae-CosC*VinPa)]) | |
| 1236 GT = np.dot(ATPe,IS1e) | |
| 1237 GR = np.dot(ARPe,IS1e) | |
| 1238 HT = np.dot(ATPe,IS2e) | |
| 1239 HR = np.dot(ARPe,IS2e) | |
| 1240 elif (TypeC == 6): # diattenuator calibration +-22.5° rotated_diattenuator_X22x5deg.odt | |
| 1241 # parameters for calibration with aCal | |
| 1242 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) | |
| 1243 BT = sqr05*DiC*(ATP1*UinPe + ATP3e*IinP) + 0.5*WiC*(ATP2e*UinPe + ATP3e*QinPe) - sqr05*ZiC*SinC*(ATP2e*VinP - ATP4*QinPe) | |
| 1244 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) | |
| 1245 BR = sqr05*DiC*(ARP1*UinPe + ARP3e*IinP) + 0.5*WiC*(ARP2e*UinPe + ARP3e*QinPe) - sqr05*ZiC*SinC*(ARP2e*VinP - ARP4*QinPe) | |
| 1246 # Correction paremeters for normal measurements; they are independent of LDR | |
| 1247 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out | |
| 1248 IS1 = np.array([IinPo,QinPo,UinPo,VinPo]) | |
| 1249 IS2 = np.array([IinPa,QinPa,UinPa,VinPa]) | |
| 1250 GT = np.dot(ATP,IS1) | |
| 1251 GR = np.dot(ARP,IS1) | |
| 1252 HT = np.dot(ATP,IS2) | |
| 1253 HR = np.dot(ARP,IS2) | |
| 1254 else: | |
| 1255 IS1e = np.array([IinPo+DiC*QinPoe,DiC*IinPo+QinPoe,ZiC*(CosC*UinPoe+SinC*VinPo),-ZiC*(SinC*UinPoe-CosC*VinPo)]) | |
| 1256 IS2e = np.array([IinPa+DiC*QinPae,DiC*IinPa+QinPae,ZiC*(CosC*UinPae+SinC*VinPa),-ZiC*(SinC*UinPae-CosC*VinPa)]) | |
| 1257 GT = np.dot(ATPe,IS1e) | |
| 1258 GR = np.dot(ARPe,IS1e) | |
| 1259 HT = np.dot(ATPe,IS2e) | |
| 1260 HR = np.dot(ARPe,IS2e) | |
| 1209 else: | 1261 else: |
| 1210 IS1 = np.array([IinPo,QinPo,UinPo,VinPo]) | 1262 print("Calibrator not implemented yet") |
| 1211 IS2 = np.array([IinPa,QinPa,UinPa,VinPa]) | 1263 sys.exit() |
| 1212 GT = np.dot(ATPe,IS1) | 1264 |
| 1213 GR = np.dot(ARPe,IS1) | 1265 elif LocC == 3: # C before receiver optics Eq.57 |
| 1214 HT = np.dot(ATPe,IS2) | 1266 |
| 1215 HR = np.dot(ARPe,IS2) | 1267 #S2ge = np.sin(np.deg2rad(2*RotO - 2*RotC)) |
| 1216 elif (TypeC == 3) or (TypeC == 4): # linear polariser calibration Eq. C.5 | 1268 #C2ge = np.cos(np.deg2rad(2*RotO - 2*RotC)) |
| 1217 # parameters for calibration with aCal | 1269 S2e = np.sin(np.deg2rad(2*RotC)) |
| 1218 AT = ATP1*IinP + ATP3e*UinPe + ZiC*CosC*(ATP2e*QinPe + ATP4*VinP) | 1270 C2e = np.cos(np.deg2rad(2*RotC)) |
| 1219 BT = DiC*(ATP1*UinPe + ATP3e*IinP) - ZiC*SinC*(ATP2e*VinP - ATP4*QinPe) | 1271 |
| 1220 AR = ARP1*IinP + ARP3e*UinPe + ZiC*CosC*(ARP2e*QinPe + ARP4*VinP) | 1272 # AS with C before the receiver optics (see document rotated_diattenuator_X22x5deg.odt) |
| 1221 BR = DiC*(ARP1*UinPe + ARP3e*IinP) - ZiC*SinC*(ARP2e*VinP - ARP4*QinPe) | 1273 AF1 = np.array([1,C2g*DiO,S2g*DiO,0]) |
| 1222 # Correction paremeters for normal measurements; they are independent of LDR | 1274 AF2 = np.array([C2g*DiO,1-S2g**2*WiO,S2g*C2g*WiO,-S2g*ZiO*SinO]) |
| 1223 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out | 1275 AF3 = np.array([S2g*DiO, S2g*C2g*WiO, 1-C2g**2*WiO, C2g*ZiO*SinO]) |
| 1224 IS1 = np.array([IinPo,QinPo,UinPo,VinPo]) | 1276 AF4 = np.array([0, S2g*SinO, -C2g*SinO, CosO]) |
| 1225 IS2 = np.array([IinPa,QinPa,UinPa,VinPa]) | 1277 |
| 1226 GT = np.dot(ATP,IS1) | 1278 ATF = (ATP1*AF1+ATP2*AF2+ATP3*AF3+ATP4*AF4) |
| 1227 GR = np.dot(ARP,IS1) | 1279 ARF = (ARP1*AF1+ARP2*AF2+ARP3*AF3+ARP4*AF4) |
| 1228 HT = np.dot(ATP,IS2) | 1280 ATF1 = ATF[0] |
| 1229 HR = np.dot(ARP,IS2) | 1281 ATF2 = ATF[1] |
| 1282 ATF3 = ATF[2] | |
| 1283 ATF4 = ATF[3] | |
| 1284 ARF1 = ARF[0] | |
| 1285 ARF2 = ARF[1] | |
| 1286 ARF3 = ARF[2] | |
| 1287 ARF4 = ARF[3] | |
| 1288 | |
| 1289 # rotated AinF by epsilon | |
| 1290 ATF2e = C2e*ATF[1] + S2e*ATF[2] | |
| 1291 ATF3e = C2e*ATF[2] - S2e*ATF[1] | |
| 1292 ARF2e = C2e*ARF[1] + S2e*ARF[2] | |
| 1293 ARF3e = C2e*ARF[2] - S2e*ARF[1] | |
| 1294 | |
| 1295 ATFe = np.array([ATF1,ATF2e,ATF3e,ATF4]) | |
| 1296 ARFe = np.array([ARF1,ARF2e,ARF3e,ARF4]) | |
| 1297 | |
| 1298 QinEe = C2e*QinE + S2e*UinE | |
| 1299 UinEe = C2e*UinE - S2e*QinE | |
| 1300 | |
| 1301 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) | |
| 1302 IinF = IinE | |
| 1303 QinF = aCal*QinE | |
| 1304 UinF = -aCal*UinE | |
| 1305 VinF = (1.-2.*aCal)*VinE | |
| 1306 | |
| 1307 IinFo = IinE | |
| 1308 QinFo = 0. | |
| 1309 UinFo = 0. | |
| 1310 VinFo = VinE | |
| 1311 | |
| 1312 IinFa = 0. | |
| 1313 QinFa = QinE | |
| 1314 UinFa = -UinE | |
| 1315 VinFa = -2.*VinE | |
| 1316 | |
| 1317 # Stokes Input Vector before receiver optics rotated by epsilon Eq. C.3 | |
| 1318 QinFe = C2e*QinF + S2e*UinF | |
| 1319 UinFe = C2e*UinF - S2e*QinF | |
| 1320 QinFoe = C2e*QinFo + S2e*UinFo | |
| 1321 UinFoe = C2e*UinFo - S2e*QinFo | |
| 1322 QinFae = C2e*QinFa + S2e*UinFa | |
| 1323 UinFae = C2e*UinFa - S2e*QinFa | |
| 1324 | |
| 1325 # Calibration signals and Calibration correction K from measurements with LDRCal / aCal | |
| 1326 if (TypeC == 2) or (TypeC == 1): # rotator calibration Eq. C.4 | |
| 1327 AT = ATF1*IinF + ATF4*h*VinF | |
| 1328 BT = ATF3e*QinF - ATF2e*h*UinF | |
| 1329 AR = ARF1*IinF + ARF4*h*VinF | |
| 1330 BR = ARF3e*QinF - ARF2e*h*UinF | |
| 1331 | |
| 1332 # Correction paremeters for normal measurements; they are independent of LDR | |
| 1333 if (not RotationErrorEpsilonForNormalMeasurements): | |
| 1334 GT = ATF1*IinE + ATF4*VinE | |
| 1335 GR = ARF1*IinE + ARF4*VinE | |
| 1336 HT = ATF2*QinE - ATF3*UinE - ATF4*2*VinE | |
| 1337 HR = ARF2*QinE - ARF3*UinE - ARF4*2*VinE | |
| 1338 else: | |
| 1339 GT = ATF1*IinE + ATF4*h*VinE | |
| 1340 GR = ARF1*IinE + ARF4*h*VinE | |
| 1341 HT = ATF2e*QinE - ATF3e*h*UinE - ATF4*h*2*VinE | |
| 1342 HR = ARF2e*QinE - ARF3e*h*UinE - ARF4*h*2*VinE | |
| 1343 | |
| 1344 elif (TypeC == 3) or (TypeC == 4): # linear polariser calibration Eq. C.5 | |
| 1345 # p = +45°, m = -45° | |
| 1346 IF1e = np.array([IinF, ZiC*CosC*QinFe, UinFe, ZiC*CosC*VinF]) | |
| 1347 IF2e = np.array([DiC*UinFe, -ZiC*SinC*VinF, DiC*IinF, ZiC*SinC*QinFe]) | |
| 1348 | |
| 1349 AT = np.dot(ATFe,IF1e) | |
| 1350 AR = np.dot(ARFe,IF1e) | |
| 1351 BT = np.dot(ATFe,IF2e) | |
| 1352 BR = np.dot(ARFe,IF2e) | |
| 1353 | |
| 1354 # Correction paremeters for normal measurements; they are independent of LDR --- the same as for TypeC = 6 | |
| 1355 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out | |
| 1356 IS1 = np.array([IinE,0,0,VinE]) | |
| 1357 IS2 = np.array([0,QinE,-UinE,-2*VinE]) | |
| 1358 | |
| 1359 GT = np.dot(ATF,IS1) | |
| 1360 GR = np.dot(ARF,IS1) | |
| 1361 HT = np.dot(ATF,IS2) | |
| 1362 HR = np.dot(ARF,IS2) | |
| 1363 else: | |
| 1364 IS1e = np.array([IinFo+DiC*QinFoe,DiC*IinFo+QinFoe,ZiC*(CosC*UinFoe+SinC*VinFo),-ZiC*(SinC*UinFoe-CosC*VinFo)]) | |
| 1365 IS2e = np.array([IinFa+DiC*QinFae,DiC*IinFa+QinFae,ZiC*(CosC*UinFae+SinC*VinFa),-ZiC*(SinC*UinFae-CosC*VinFa)]) | |
| 1366 GT = np.dot(ATFe,IS1e) | |
| 1367 GR = np.dot(ARFe,IS1e) | |
| 1368 HT = np.dot(ATFe,IS2e) | |
| 1369 HR = np.dot(ARFe,IS2e) | |
| 1370 | |
| 1371 elif (TypeC == 6): # diattenuator calibration +-22.5° rotated_diattenuator_X22x5deg.odt | |
| 1372 # p = +22.5°, m = -22.5° | |
| 1373 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]) | |
| 1374 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]) | |
| 1375 | |
| 1376 AT = np.dot(ATFe,IF1e) | |
| 1377 AR = np.dot(ARFe,IF1e) | |
| 1378 BT = np.dot(ATFe,IF2e) | |
| 1379 BR = np.dot(ARFe,IF2e) | |
| 1380 | |
| 1381 # Correction paremeters for normal measurements; they are independent of LDR | |
| 1382 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out | |
| 1383 #IS1 = np.array([IinE,0,0,VinE]) | |
| 1384 #IS2 = np.array([0,QinE,-UinE,-2*VinE]) | |
| 1385 IS1 = np.array([IinFo,0,0,VinFo]) | |
| 1386 IS2 = np.array([0,QinFa,UinFa,VinFa]) | |
| 1387 GT = np.dot(ATF,IS1) | |
| 1388 GR = np.dot(ARF,IS1) | |
| 1389 HT = np.dot(ATF,IS2) | |
| 1390 HR = np.dot(ARF,IS2) | |
| 1391 else: | |
| 1392 #IS1e = np.array([IinE,DiC*IinE,ZiC*SinC*VinE,ZiC*CosC*VinE]) | |
| 1393 #IS2e = np.array([DiC*QinEe,QinEe,-ZiC*(CosC*UinEe+2*SinC*VinE),ZiC*(SinC*UinEe-2*CosC*VinE)]) | |
| 1394 IS1e = np.array([IinFo+DiC*QinFoe,DiC*IinFo+QinFoe,ZiC*(CosC*UinFoe+SinC*VinFo),-ZiC*(SinC*UinFoe-CosC*VinFo)]) | |
| 1395 IS2e = np.array([IinFa+DiC*QinFae,DiC*IinFa+QinFae,ZiC*(CosC*UinFae+SinC*VinFa),-ZiC*(SinC*UinFae-CosC*VinFa)]) | |
| 1396 GT = np.dot(ATFe,IS1e) | |
| 1397 GR = np.dot(ARFe,IS1e) | |
| 1398 HT = np.dot(ATFe,IS2e) | |
| 1399 HR = np.dot(ARFe,IS2e) | |
| 1400 | |
| 1401 | |
| 1230 else: | 1402 else: |
| 1231 IS1e = np.array([IinPo+DiC*QinPoe,DiC*IinPo+QinPoe,ZiC*(CosC*UinPoe+SinC*VinPo),-ZiC*(SinC*UinPoe-CosC*VinPo)]) | 1403 print('Calibrator not implemented yet') |
| 1232 IS2e = np.array([IinPa+DiC*QinPae,DiC*IinPa+QinPae,ZiC*(CosC*UinPae+SinC*VinPa),-ZiC*(SinC*UinPae-CosC*VinPa)]) | 1404 sys.exit() |
| 1233 GT = np.dot(ATPe,IS1e) | 1405 |
| 1234 GR = np.dot(ARPe,IS1e) | 1406 elif LocC == 2: # C behind emitter optics Eq.57 |
| 1235 HT = np.dot(ATPe,IS2e) | 1407 #print("Calibrator location not implemented yet") |
| 1236 HR = np.dot(ARPe,IS2e) | 1408 S2e = np.sin(np.deg2rad(2*RotC)) |
| 1237 elif (TypeC == 6): # diattenuator calibration +-22.5° rotated_diattenuator_X22x5deg.odt | 1409 C2e = np.cos(np.deg2rad(2*RotC)) |
| 1238 # parameters for calibration with aCal | 1410 |
| 1239 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) | 1411 # AS with C before the receiver optics (see document rotated_diattenuator_X22x5deg.odt) |
| 1240 BT = sqr05*DiC*(ATP1*UinPe + ATP3e*IinP) + 0.5*WiC*(ATP2e*UinPe + ATP3e*QinPe) - sqr05*ZiC*SinC*(ATP2e*VinP - ATP4*QinPe) | 1412 AF1 = np.array([1,C2g*DiO,S2g*DiO,0]) |
| 1241 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) | 1413 AF2 = np.array([C2g*DiO,1-S2g**2*WiO,S2g*C2g*WiO,-S2g*ZiO*SinO]) |
| 1242 BR = sqr05*DiC*(ARP1*UinPe + ARP3e*IinP) + 0.5*WiC*(ARP2e*UinPe + ARP3e*QinPe) - sqr05*ZiC*SinC*(ARP2e*VinP - ARP4*QinPe) | 1414 AF3 = np.array([S2g*DiO, S2g*C2g*WiO, 1-C2g**2*WiO, C2g*ZiO*SinO]) |
| 1243 # Correction paremeters for normal measurements; they are independent of LDR | 1415 AF4 = np.array([0, S2g*SinO, -C2g*SinO, CosO]) |
| 1244 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out | 1416 |
| 1245 IS1 = np.array([IinPo,QinPo,UinPo,VinPo]) | 1417 ATF = (ATP1*AF1+ATP2*AF2+ATP3*AF3+ATP4*AF4) |
| 1246 IS2 = np.array([IinPa,QinPa,UinPa,VinPa]) | 1418 ARF = (ARP1*AF1+ARP2*AF2+ARP3*AF3+ARP4*AF4) |
| 1247 GT = np.dot(ATP,IS1) | 1419 ATF1 = ATF[0] |
| 1248 GR = np.dot(ARP,IS1) | 1420 ATF2 = ATF[1] |
| 1249 HT = np.dot(ATP,IS2) | 1421 ATF3 = ATF[2] |
| 1250 HR = np.dot(ARP,IS2) | 1422 ATF4 = ATF[3] |
| 1423 ARF1 = ARF[0] | |
| 1424 ARF2 = ARF[1] | |
| 1425 ARF3 = ARF[2] | |
| 1426 ARF4 = ARF[3] | |
| 1427 | |
| 1428 # AS with C behind the emitter -------------------------------------------- | |
| 1429 # terms without aCal | |
| 1430 ATE1o, ARE1o = ATF1, ARF1 | |
| 1431 ATE2o, ARE2o = 0., 0. | |
| 1432 ATE3o, ARE3o = 0., 0. | |
| 1433 ATE4o, ARE4o = ATF4, ARF4 | |
| 1434 # terms with aCal | |
| 1435 ATE1a, ARE1a = 0. , 0. | |
| 1436 ATE2a, ARE2a = ATF2, ARF2 | |
| 1437 ATE3a, ARE3a = -ATF3, -ARF3 | |
| 1438 ATE4a, ARE4a = -2*ATF4, -2*ARF4 | |
| 1439 # rotated AinEa by epsilon | |
| 1440 ATE2ae = C2e*ATF2 + S2e*ATF3 | |
| 1441 ATE3ae = -S2e*ATF2 - C2e*ATF3 | |
| 1442 ARE2ae = C2e*ARF2 + S2e*ARF3 | |
| 1443 ARE3ae = -S2e*ARF2 - C2e*ARF3 | |
| 1444 | |
| 1445 ATE1 = ATE1o | |
| 1446 ATE2e = aCal*ATE2ae | |
| 1447 ATE3e = aCal*ATE3ae | |
| 1448 ATE4 = (1-2*aCal)*ATF4 | |
| 1449 ARE1 = ARE1o | |
| 1450 ARE2e = aCal*ARE2ae | |
| 1451 ARE3e = aCal*ARE3ae | |
| 1452 ARE4 = (1-2*aCal)*ARF4 | |
| 1453 | |
| 1454 # rotated IinE | |
| 1455 QinEe = C2e*QinE + S2e*UinE | |
| 1456 UinEe = C2e*UinE - S2e*QinE | |
| 1457 | |
| 1458 # --- Calibration signals and Calibration correction K from measurements with LDRCal / aCal | |
| 1459 if (TypeC == 2) or (TypeC == 1): # +++++++++ rotator calibration Eq. C.4 | |
| 1460 AT = ATE1o*IinE + (ATE4o+aCal*ATE4a)*h*VinE | |
| 1461 BT = aCal * (ATE3ae*QinEe - ATE2ae*h*UinEe) | |
| 1462 AR = ARE1o*IinE + (ARE4o+aCal*ARE4a)*h*VinE | |
| 1463 BR = aCal * (ARE3ae*QinEe - ARE2ae*h*UinEe) | |
| 1464 | |
| 1465 # Correction paremeters for normal measurements; they are independent of LDR | |
| 1466 if (not RotationErrorEpsilonForNormalMeasurements): | |
| 1467 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) | |
| 1468 GT = ATE1o*IinE + ATE4o*h*VinE | |
| 1469 GR = ARE1o*IinE + ARE4o*h*VinE | |
| 1470 HT = ATE2a*QinE + ATE3a*h*UinEe + ATE4a*h*VinE | |
| 1471 HR = ARE2a*QinE + ARE3a*h*UinEe + ARE4a*h*VinE | |
| 1472 else: | |
| 1473 GT = ATE1o*IinE + ATE4o*h*VinE | |
| 1474 GR = ARE1o*IinE + ARE4o*h*VinE | |
| 1475 HT = ATE2ae*QinE + ATE3ae*h*UinEe + ATE4a*h*VinE | |
| 1476 HR = ARE2ae*QinE + ARE3ae*h*UinEe + ARE4a*h*VinE | |
| 1477 | |
| 1478 elif (TypeC == 3) or (TypeC == 4): # +++++++++ linear polariser calibration Eq. C.5 | |
| 1479 # p = +45°, m = -45° | |
| 1480 AT = ATE1*IinE + ZiC*CosC*(ATE2e*QinEe + ATE4*VinE) + ATE3e*UinEe | |
| 1481 BT = DiC*(ATE1*UinEe + ATE3e*IinE) + ZiC*SinC*(ATE4*QinEe - ATE2e*VinE) | |
| 1482 AR = ARE1*IinE + ZiC*CosC*(ARE2e*QinEe + ARE4*VinE) + ARE3e*UinEe | |
| 1483 BR = DiC*(ARE1*UinEe + ARE3e*IinE) + ZiC*SinC*(ARE4*QinEe - ARE2e*VinE) | |
| 1484 | |
| 1485 # Correction paremeters for normal measurements; they are independent of LDR | |
| 1486 if (not RotationErrorEpsilonForNormalMeasurements): | |
| 1487 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) | |
| 1488 GT = ATE1o*IinE + ATE4o*VinE | |
| 1489 GR = ARE1o*IinE + ARE4o*VinE | |
| 1490 HT = ATE2a*QinE + ATE3a*UinE + ATE4a*VinE | |
| 1491 HR = ARE2a*QinE + ARE3a*UinE + ARE4a*VinE | |
| 1492 else: | |
| 1493 D = IinE + DiC*QinEe | |
| 1494 A = DiC*IinE + QinEe | |
| 1495 B = ZiC*(CosC*UinEe + SinC*VinE) | |
| 1496 C = -ZiC*(SinC*UinEe - CosC*VinE) | |
| 1497 GT = ATE1o*D + ATE4o*C | |
| 1498 GR = ARE1o*D + ARE4o*C | |
| 1499 HT = ATE2a*A + ATE3a*B + ATE4a*C | |
| 1500 HR = ARE2a*A + ARE3a*B + ARE4a*C | |
| 1501 | |
| 1502 elif (TypeC == 6): # real HWP calibration +-22.5° rotated_diattenuator_X22x5deg.odt | |
| 1503 # p = +22.5°, m = -22.5° | |
| 1504 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]) | |
| 1505 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]) | |
| 1506 ATEe = np.array([ATE1,ATE2e,ATE3e,ATE4]) | |
| 1507 AREe = np.array([ARE1,ARE2e,ARE3e,ARE4]) | |
| 1508 AT = np.dot(ATEe,IE1e) | |
| 1509 AR = np.dot(AREe,IE1e) | |
| 1510 BT = np.dot(ATEe,IE2e) | |
| 1511 BR = np.dot(AREe,IE2e) | |
| 1512 | |
| 1513 # Correction paremeters for normal measurements; they are independent of LDR | |
| 1514 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out | |
| 1515 GT = ATE1o*IinE + ATE4o*VinE | |
| 1516 GR = ARE1o*IinE + ARE4o*VinE | |
| 1517 HT = ATE2a*QinE + ATE3a*UinE + ATE4a*VinE | |
| 1518 HR = ARE2a*QinE + ARE3a*UinE + ARE4a*VinE | |
| 1519 else: | |
| 1520 D = IinE + DiC*QinEe | |
| 1521 A = DiC*IinE + QinEe | |
| 1522 B = ZiC*(CosC*UinEe + SinC*VinE) | |
| 1523 C = -ZiC*(SinC*UinEe - CosC*VinE) | |
| 1524 GT = ATE1o*D + ATE4o*C | |
| 1525 GR = ARE1o*D + ARE4o*C | |
| 1526 HT = ATE2a*A + ATE3a*B + ATE4a*C | |
| 1527 HR = ARE2a*A + ARE3a*B + ARE4a*C | |
| 1528 | |
| 1251 else: | 1529 else: |
| 1252 IS1e = np.array([IinPo+DiC*QinPoe,DiC*IinPo+QinPoe,ZiC*(CosC*UinPoe+SinC*VinPo),-ZiC*(SinC*UinPoe-CosC*VinPo)]) | 1530 print('Calibrator not implemented yet') |
| 1253 IS2e = np.array([IinPa+DiC*QinPae,DiC*IinPa+QinPae,ZiC*(CosC*UinPae+SinC*VinPa),-ZiC*(SinC*UinPae-CosC*VinPa)]) | 1531 sys.exit() |
| 1254 GT = np.dot(ATPe,IS1e) | 1532 |
| 1255 GR = np.dot(ARPe,IS1e) | 1533 # Calibration signals with aCal => Determination of the correction K of the real calibration factor |
| 1256 HT = np.dot(ATPe,IS2e) | 1534 IoutTp = TaT*TiT*TiO*TiE*(AT + BT) |
| 1257 HR = np.dot(ARPe,IS2e) | 1535 IoutTm = TaT*TiT*TiO*TiE*(AT - BT) |
| 1258 else: | 1536 IoutRp = TaR*TiR*TiO*TiE*(AR + BR) |
| 1259 print("Calibrator not implemented yet") | 1537 IoutRm = TaR*TiR*TiO*TiE*(AR - BR) |
| 1260 sys.exit() | 1538 # --- Results and Corrections; electronic etaR and etaT are assumed to be 1 |
| 1261 | 1539 #Eta = TiR/TiT # Eta = Eta*/K Eq. 84 |
| 1262 elif LocC == 3: # C before receiver optics Eq.57 | 1540 Etapx = IoutRp/IoutTp |
| 1263 | 1541 Etamx = IoutRm/IoutTm |
| 1264 #S2ge = np.sin(np.deg2rad(2*RotO - 2*RotC)) | 1542 Etax = (Etapx*Etamx)**0.5 |
| 1265 #C2ge = np.cos(np.deg2rad(2*RotO - 2*RotC)) | 1543 K = Etax / Eta0 |
| 1266 S2e = np.sin(np.deg2rad(2*RotC)) | 1544 #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)) |
| 1267 C2e = np.cos(np.deg2rad(2*RotC)) | 1545 #print("{0:6.3f},{1:6.3f},{2:6.3f},{3:6.3f}".format(DiC, ZiC, Kp, Km)) |
| 1268 | 1546 |
| 1269 # AS with C before the receiver optics (see document rotated_diattenuator_X22x5deg.odt) | 1547 # For comparison with Volkers Libreoffice Müller Matrix spreadsheet |
| 1270 AF1 = np.array([1,C2g*DiO,S2g*DiO,0]) | 1548 #Eta_test_p = (IoutRp/IoutTp) |
| 1271 AF2 = np.array([C2g*DiO,1-S2g**2*WiO,S2g*C2g*WiO,-S2g*ZiO*SinO]) | 1549 #Eta_test_m = (IoutRm/IoutTm) |
| 1272 AF3 = np.array([S2g*DiO, S2g*C2g*WiO, 1-C2g**2*WiO, C2g*ZiO*SinO]) | 1550 #Eta_test = (Eta_test_p*Eta_test_m)**0.5 |
| 1273 AF4 = np.array([0, S2g*SinO, -C2g*SinO, CosO]) | 1551 |
| 1274 | 1552 # ************************************************************************* |
| 1275 ATF = (ATP1*AF1+ATP2*AF2+ATP3*AF3+ATP4*AF4) | 1553 iLDR = -1 |
| 1276 ARF = (ARP1*AF1+ARP2*AF2+ARP3*AF3+ARP4*AF4) | 1554 for LDRTrue in LDRrange: |
| 1277 ATF1 = ATF[0] | 1555 iLDR = iLDR + 1 |
| 1278 ATF2 = ATF[1] | 1556 atrue = (1-LDRTrue)/(1+LDRTrue) |
| 1279 ATF3 = ATF[2] | 1557 # ----- Forward simulated signals and LDRsim with atrue; from input file |
| 1280 ATF4 = ATF[3] | 1558 It = TaT*TiT*TiO*TiE*(GT+atrue*HT) # TaT*TiT*TiC*TiO*IinL*(GT+atrue*HT) |
| 1281 ARF1 = ARF[0] | 1559 Ir = TaR*TiR*TiO*TiE*(GR+atrue*HR) # TaR*TiR*TiC*TiO*IinL*(GR+atrue*HR) |
| 1282 ARF2 = ARF[1] | 1560 |
| 1283 ARF3 = ARF[2] | 1561 # LDRsim = 1/Eta*Ir/It # simulated LDR* with Y from input file |
| 1284 ARF4 = ARF[3] | 1562 LDRsim = Ir/It # simulated uncorrected LDR with Y from input file |
| 1285 | 1563 ''' |
| 1286 # rotated AinF by epsilon | 1564 if Y == 1.: |
| 1287 ATF2e = C2e*ATF[1] + S2e*ATF[2] | 1565 LDRsimx = LDRsim |
| 1288 ATF3e = C2e*ATF[2] - S2e*ATF[1] | 1566 LDRsimx2 = LDRsim2 |
| 1289 ARF2e = C2e*ARF[1] + S2e*ARF[2] | |
| 1290 ARF3e = C2e*ARF[2] - S2e*ARF[1] | |
| 1291 | |
| 1292 ATFe = np.array([ATF1,ATF2e,ATF3e,ATF4]) | |
| 1293 ARFe = np.array([ARF1,ARF2e,ARF3e,ARF4]) | |
| 1294 | |
| 1295 QinEe = C2e*QinE + S2e*UinE | |
| 1296 UinEe = C2e*UinE - S2e*QinE | |
| 1297 | |
| 1298 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) | |
| 1299 IinF = IinE | |
| 1300 QinF = aCal*QinE | |
| 1301 UinF = -aCal*UinE | |
| 1302 VinF = (1.-2.*aCal)*VinE | |
| 1303 | |
| 1304 IinFo = IinE | |
| 1305 QinFo = 0. | |
| 1306 UinFo = 0. | |
| 1307 VinFo = VinE | |
| 1308 | |
| 1309 IinFa = 0. | |
| 1310 QinFa = QinE | |
| 1311 UinFa = -UinE | |
| 1312 VinFa = -2.*VinE | |
| 1313 | |
| 1314 # Stokes Input Vector before receiver optics rotated by epsilon Eq. C.3 | |
| 1315 QinFe = C2e*QinF + S2e*UinF | |
| 1316 UinFe = C2e*UinF - S2e*QinF | |
| 1317 QinFoe = C2e*QinFo + S2e*UinFo | |
| 1318 UinFoe = C2e*UinFo - S2e*QinFo | |
| 1319 QinFae = C2e*QinFa + S2e*UinFa | |
| 1320 UinFae = C2e*UinFa - S2e*QinFa | |
| 1321 | |
| 1322 # Calibration signals and Calibration correction K from measurements with LDRCal / aCal | |
| 1323 if (TypeC == 2) or (TypeC == 1): # rotator calibration Eq. C.4 | |
| 1324 AT = ATF1*IinF + ATF4*h*VinF | |
| 1325 BT = ATF3e*QinF - ATF2e*h*UinF | |
| 1326 AR = ARF1*IinF + ARF4*h*VinF | |
| 1327 BR = ARF3e*QinF - ARF2e*h*UinF | |
| 1328 | |
| 1329 # Correction paremeters for normal measurements; they are independent of LDR | |
| 1330 if (not RotationErrorEpsilonForNormalMeasurements): | |
| 1331 GT = ATF1*IinE + ATF4*VinE | |
| 1332 GR = ARF1*IinE + ARF4*VinE | |
| 1333 HT = ATF2*QinE - ATF3*UinE - ATF4*2*VinE | |
| 1334 HR = ARF2*QinE - ARF3*UinE - ARF4*2*VinE | |
| 1335 else: | 1567 else: |
| 1336 GT = ATF1*IinE + ATF4*h*VinE | 1568 LDRsimx = 1./LDRsim |
| 1337 GR = ARF1*IinE + ARF4*h*VinE | 1569 LDRsimx2 = 1./LDRsim2 |
| 1338 HT = ATF2e*QinE - ATF3e*h*UinE - ATF4*h*2*VinE | 1570 ''' |
| 1339 HR = ARF2e*QinE - ARF3e*h*UinE - ARF4*h*2*VinE | 1571 # ----- Backward correction |
| 1340 | 1572 # Corrected LDRCorr from forward simulated LDRsim (atrue) with assumed true G0,H0,K0 |
| 1341 elif (TypeC == 3) or (TypeC == 4): # linear polariser calibration Eq. C.5 | 1573 LDRCorr = (LDRsim*K0/Etax*(GT0+HT0)-(GR0+HR0))/((GR0-HR0)-LDRsim*K0/Etax*(GT0-HT0)) |
| 1342 # p = +45°, m = -45° | 1574 |
| 1343 IF1e = np.array([IinF, ZiC*CosC*QinFe, UinFe, ZiC*CosC*VinF]) | 1575 # -- F11corr from It and Ir and calibration EtaX |
| 1344 IF2e = np.array([DiC*UinFe, -ZiC*SinC*VinF, DiC*IinF, ZiC*SinC*QinFe]) | 1576 Text1 = "!!! EXPERIMENTAL !!! F11corr from It and Ir with calibration EtaX: x-axis: F11corr(LDRtrue) / F11corr(LDRtrue = 0.004) - 1" |
| 1345 | 1577 F11corr = 1/(TiO*TiE)*((HR0*Etax/K0*It/TTa-HT0*Ir/TRa)/(HR0*GT0-HT0*GR0)) # IL = 1 Eq.(64) |
| 1346 AT = np.dot(ATFe,IF1e) | 1578 |
| 1347 AR = np.dot(ARFe,IF1e) | 1579 #Text1 = "F11corr from It and Ir without corrections but with calibration EtaX: x-axis: F11corr(LDRtrue) devided by F11corr(LDRtrue = 0.004)" |
| 1348 BT = np.dot(ATFe,IF2e) | 1580 #F11corr = 0.5/(TiO*TiE)*(Etax*It/TTa+Ir/TRa) # IL = 1 Eq.(64) |
| 1349 BR = np.dot(ARFe,IF2e) | 1581 |
| 1350 | 1582 # -- It from It only with atrue without corrections - for BERTHA (and PollyXTs) |
| 1351 # Correction paremeters for normal measurements; they are independent of LDR --- the same as for TypeC = 6 | 1583 #Text1 = " x-axis: IT(LDRtrue) / IT(LDRtrue = 0.004) - 1" |
| 1352 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out | 1584 #F11corr = It/(TaT*TiT*TiO*TiE) #/(TaT*TiT*TiO*TiE*(GT0+atrue*HT0)) |
| 1353 IS1 = np.array([IinE,0,0,VinE]) | 1585 # !!! see below line 1673ff |
| 1354 IS2 = np.array([0,QinE,-UinE,-2*VinE]) | 1586 |
| 1355 | 1587 aF11corr[iLDR,iN] = F11corr |
| 1356 GT = np.dot(ATF,IS1) | 1588 aA[iLDR,iN] = LDRCorr |
| 1357 GR = np.dot(ARF,IS1) | 1589 |
| 1358 HT = np.dot(ATF,IS2) | 1590 aX[0,iN] = GR |
| 1359 HR = np.dot(ARF,IS2) | 1591 aX[1,iN] = GT |
| 1360 else: | 1592 aX[2,iN] = HR |
| 1361 IS1e = np.array([IinFo+DiC*QinFoe,DiC*IinFo+QinFoe,ZiC*(CosC*UinFoe+SinC*VinFo),-ZiC*(SinC*UinFoe-CosC*VinFo)]) | 1593 aX[3,iN] = HT |
| 1362 IS2e = np.array([IinFa+DiC*QinFae,DiC*IinFa+QinFae,ZiC*(CosC*UinFae+SinC*VinFa),-ZiC*(SinC*UinFae-CosC*VinFa)]) | 1594 aX[4,iN] = K |
| 1363 GT = np.dot(ATFe,IS1e) | 1595 |
| 1364 GR = np.dot(ARFe,IS1e) | 1596 aLDRCal[iN] = iLDRCal |
| 1365 HT = np.dot(ATFe,IS2e) | 1597 aERaT[iN] = iERaT |
| 1366 HR = np.dot(ARFe,IS2e) | 1598 aERaR[iN] = iERaR |
| 1367 | 1599 aRotaT[iN] = iRotaT |
| 1368 elif (TypeC == 6): # diattenuator calibration +-22.5° rotated_diattenuator_X22x5deg.odt | 1600 aRotaR[iN] = iRotaR |
| 1369 # p = +22.5°, m = -22.5° | 1601 aRetT[iN] = iRetT |
| 1370 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]) | 1602 aRetR[iN] = iRetR |
| 1371 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]) | 1603 |
| 1372 | 1604 aRotL[iN] = iRotL |
| 1373 AT = np.dot(ATFe,IF1e) | 1605 aRotE[iN] = iRotE |
| 1374 AR = np.dot(ARFe,IF1e) | 1606 aRetE[iN] = iRetE |
| 1375 BT = np.dot(ATFe,IF2e) | 1607 aRotO[iN] = iRotO |
| 1376 BR = np.dot(ARFe,IF2e) | 1608 aRetO[iN] = iRetO |
| 1377 | 1609 aRotC[iN] = iRotC |
| 1378 # Correction paremeters for normal measurements; they are independent of LDR | 1610 aRetC[iN] = iRetC |
| 1379 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out | 1611 aDiO[iN] = iDiO |
| 1380 #IS1 = np.array([IinE,0,0,VinE]) | 1612 aDiE[iN] = iDiE |
| 1381 #IS2 = np.array([0,QinE,-UinE,-2*VinE]) | 1613 aDiC[iN] = iDiC |
| 1382 IS1 = np.array([IinFo,0,0,VinFo]) | 1614 aTP[iN] = iTP |
| 1383 IS2 = np.array([0,QinFa,UinFa,VinFa]) | 1615 aTS[iN] = iTS |
| 1384 GT = np.dot(ATF,IS1) | 1616 aRP[iN] = iRP |
| 1385 GR = np.dot(ARF,IS1) | 1617 aRS[iN] = iRS |
| 1386 HT = np.dot(ATF,IS2) | 1618 |
| 1387 HR = np.dot(ARF,IS2) | 1619 # --- END loop |
| 1388 else: | 1620 btime = clock() |
| 1389 #IS1e = np.array([IinE,DiC*IinE,ZiC*SinC*VinE,ZiC*CosC*VinE]) | 1621 print("\r done in ", "{0:5.0f}".format(btime-atime), "sec") #, end="\r") |
| 1390 #IS2e = np.array([DiC*QinEe,QinEe,-ZiC*(CosC*UinEe+2*SinC*VinE),ZiC*(SinC*UinEe-2*CosC*VinE)]) | 1622 |
| 1391 IS1e = np.array([IinFo+DiC*QinFoe,DiC*IinFo+QinFoe,ZiC*(CosC*UinFoe+SinC*VinFo),-ZiC*(SinC*UinFoe-CosC*VinFo)]) | 1623 # --- Plot ----------------------------------------------------------------- |
| 1392 IS2e = np.array([IinFa+DiC*QinFae,DiC*IinFa+QinFae,ZiC*(CosC*UinFae+SinC*VinFa),-ZiC*(SinC*UinFae-CosC*VinFa)]) | 1624 if (sns_loaded): |
| 1393 GT = np.dot(ATFe,IS1e) | 1625 sns.set_style("whitegrid") |
| 1394 GR = np.dot(ARFe,IS1e) | 1626 sns.set_palette("bright", 6) |
| 1395 HT = np.dot(ATFe,IS2e) | 1627 |
| 1396 HR = np.dot(ARFe,IS2e) | 1628 ''' |
| 1397 | 1629 fig2 = plt.figure() |
| 1398 | 1630 plt.plot(aA[2,:],'b.') |
| 1399 else: | 1631 plt.plot(aA[3,:],'r.') |
| 1400 print('Calibrator not implemented yet') | 1632 plt.plot(aA[4,:],'g.') |
| 1401 sys.exit() | 1633 #plt.plot(aA[6,:],'c.') |
| 1402 | 1634 plt.show |
| 1403 elif LocC == 2: # C behind emitter optics Eq.57 | 1635 ''' |
| 1404 #print("Calibrator location not implemented yet") | 1636 # Plot LDR |
| 1405 S2e = np.sin(np.deg2rad(2*RotC)) | 1637 def PlotSubHist(aVar, aX, X0, daX, iaX, naX): |
| 1406 C2e = np.cos(np.deg2rad(2*RotC)) | 1638 fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) |
| 1407 | 1639 iLDR = -1 |
| 1408 # AS with C before the receiver optics (see document rotated_diattenuator_X22x5deg.odt) | 1640 for LDRTrue in LDRrange: |
| 1409 AF1 = np.array([1,C2g*DiO,S2g*DiO,0]) | 1641 iLDR = iLDR + 1 |
| 1410 AF2 = np.array([C2g*DiO,1-S2g**2*WiO,S2g*C2g*WiO,-S2g*ZiO*SinO]) | 1642 |
| 1411 AF3 = np.array([S2g*DiO, S2g*C2g*WiO, 1-C2g**2*WiO, C2g*ZiO*SinO]) | 1643 LDRmin[iLDR] = np.min(aA[iLDR,:]) |
| 1412 AF4 = np.array([0, S2g*SinO, -C2g*SinO, CosO]) | 1644 LDRmax[iLDR] = np.max(aA[iLDR,:]) |
| 1413 | 1645 Rmin = LDRmin[iLDR] * 0.995 # np.min(aA[iLDR,:]) * 0.995 |
| 1414 ATF = (ATP1*AF1+ATP2*AF2+ATP3*AF3+ATP4*AF4) | 1646 Rmax = LDRmax[iLDR] * 1.005 # np.max(aA[iLDR,:]) * 1.005 |
| 1415 ARF = (ARP1*AF1+ARP2*AF2+ARP3*AF3+ARP4*AF4) | 1647 |
| 1416 ATF1 = ATF[0] | 1648 #plt.subplot(5,2,iLDR+1) |
| 1417 ATF2 = ATF[1] | 1649 plt.subplot(1,5,iLDR+1) |
| 1418 ATF3 = ATF[2] | 1650 (n, bins, patches) = plt.hist(aA[iLDR,:], |
| 1419 ATF4 = ATF[3] | 1651 bins=100, log=False, |
| 1420 ARF1 = ARF[0] | 1652 range=[Rmin, Rmax], |
| 1421 ARF2 = ARF[1] | 1653 alpha=0.5, normed=False, color = '0.5', histtype='stepfilled') |
| 1422 ARF3 = ARF[2] | 1654 |
| 1423 ARF4 = ARF[3] | 1655 for iaX in range(-naX,naX+1): |
| 1424 | 1656 plt.hist(aA[iLDR,aX == iaX], |
| 1425 # AS with C behind the emitter -------------------------------------------- | 1657 range=[Rmin, Rmax], |
| 1426 # terms without aCal | 1658 bins=100, log=False, alpha=0.3, normed=False, histtype='stepfilled', label = str(round(X0 + iaX*daX/naX,5))) |
| 1427 ATE1o, ARE1o = ATF1, ARF1 | 1659 |
| 1428 ATE2o, ARE2o = 0., 0. | 1660 if (iLDR == 2): plt.legend() |
| 1429 ATE3o, ARE3o = 0., 0. | 1661 |
| 1430 ATE4o, ARE4o = ATF4, ARF4 | 1662 plt.tick_params(axis='both', labelsize=9) |
| 1431 # terms with aCal | 1663 plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2) |
| 1432 ATE1a, ARE1a = 0. , 0. | 1664 |
| 1433 ATE2a, ARE2a = ATF2, ARF2 | 1665 #plt.title(LID + ' ' + aVar, fontsize=18) |
| 1434 ATE3a, ARE3a = -ATF3, -ARF3 | 1666 #plt.ylabel('frequency', fontsize=10) |
| 1435 ATE4a, ARE4a = -2*ATF4, -2*ARF4 | 1667 #plt.xlabel('LDRcorr', fontsize=10) |
| 1436 # rotated AinEa by epsilon | 1668 #fig.tight_layout() |
| 1437 ATE2ae = C2e*ATF2 + S2e*ATF3 | 1669 fig.suptitle(LID + ' with ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar, fontsize=14, y=1.05) |
| 1438 ATE3ae = -S2e*ATF2 - C2e*ATF3 | 1670 #plt.show() |
| 1439 ARE2ae = C2e*ARF2 + S2e*ARF3 | 1671 #fig.savefig(LID + '_' + aVar + '.png', dpi=150, bbox_inches='tight', pad_inches=0) |
| 1440 ARE3ae = -S2e*ARF2 - C2e*ARF3 | 1672 #plt.close |
| 1441 | 1673 return |
| 1442 ATE1 = ATE1o | 1674 |
| 1443 ATE2e = aCal*ATE2ae | 1675 if (nRotL > 0): PlotSubHist("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) |
| 1444 ATE3e = aCal*ATE3ae | 1676 if (nRetE > 0): PlotSubHist("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) |
| 1445 ATE4 = (1-2*aCal)*ATF4 | 1677 if (nRotE > 0): PlotSubHist("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) |
| 1446 ARE1 = ARE1o | 1678 if (nDiE > 0): PlotSubHist("DiE", aDiE, DiE0, dDiE, iDiE, nDiE) |
| 1447 ARE2e = aCal*ARE2ae | 1679 if (nRetO > 0): PlotSubHist("RetO", aRetO, RetO0, dRetO, iRetO, nRetO) |
| 1448 ARE3e = aCal*ARE3ae | 1680 if (nRotO > 0): PlotSubHist("RotO", aRotO, RotO0, dRotO, iRotO, nRotO) |
| 1449 ARE4 = (1-2*aCal)*ARF4 | 1681 if (nDiO > 0): PlotSubHist("DiO", aDiO, DiO0, dDiO, iDiO, nDiO) |
| 1450 | 1682 if (nDiC > 0): PlotSubHist("DiC", aDiC, DiC0, dDiC, iDiC, nDiC) |
| 1451 # rotated IinE | 1683 if (nRotC > 0): PlotSubHist("RotC", aRotC, RotC0, dRotC, iRotC, nRotC) |
| 1452 QinEe = C2e*QinE + S2e*UinE | 1684 if (nRetC > 0): PlotSubHist("RetC", aRetC, RetC0, dRetC, iRetC, nRetC) |
| 1453 UinEe = C2e*UinE - S2e*QinE | 1685 if (nTP > 0): PlotSubHist("TP", aTP, TP0, dTP, iTP, nTP) |
| 1454 | 1686 if (nTS > 0): PlotSubHist("TS", aTS, TS0, dTS, iTS, nTS) |
| 1455 # --- Calibration signals and Calibration correction K from measurements with LDRCal / aCal | 1687 if (nRP > 0): PlotSubHist("RP", aRP, RP0, dRP, iRP, nRP) |
| 1456 if (TypeC == 2) or (TypeC == 1): # +++++++++ rotator calibration Eq. C.4 | 1688 if (nRS > 0): PlotSubHist("RS", aRS, RS0, dRS, iRS, nRS) |
| 1457 AT = ATE1o*IinE + (ATE4o+aCal*ATE4a)*h*VinE | 1689 if (nRetT > 0): PlotSubHist("RetT", aRetT, RetT0, dRetT, iRetT, nRetT) |
| 1458 BT = aCal * (ATE3ae*QinEe - ATE2ae*h*UinEe) | 1690 if (nRetR > 0): PlotSubHist("RetR", aRetR, RetR0, dRetR, iRetR, nRetR) |
| 1459 AR = ARE1o*IinE + (ARE4o+aCal*ARE4a)*h*VinE | 1691 if (nERaT > 0): PlotSubHist("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) |
| 1460 BR = aCal * (ARE3ae*QinEe - ARE2ae*h*UinEe) | 1692 if (nERaR > 0): PlotSubHist("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) |
| 1461 | 1693 if (nRotaT > 0): PlotSubHist("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) |
| 1462 # Correction paremeters for normal measurements; they are independent of LDR | 1694 if (nRotaR > 0): PlotSubHist("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) |
| 1463 if (not RotationErrorEpsilonForNormalMeasurements): | 1695 if (nLDRCal > 0): PlotSubHist("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) |
| 1464 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) | 1696 |
| 1465 GT = ATE1o*IinE + ATE4o*h*VinE | 1697 plt.show() |
| 1466 GR = ARE1o*IinE + ARE4o*h*VinE | 1698 plt.close |
| 1467 HT = ATE2a*QinE + ATE3a*h*UinEe + ATE4a*h*VinE | 1699 ''' |
| 1468 HR = ARE2a*QinE + ARE3a*h*UinEe + ARE4a*h*VinE | 1700 print() |
| 1469 else: | 1701 #print("IT(LDRtrue) devided by IT(LDRtrue = 0.004)") |
| 1470 GT = ATE1o*IinE + ATE4o*h*VinE | 1702 print(" ############################################################################## ") |
| 1471 GR = ARE1o*IinE + ARE4o*h*VinE | 1703 print(Text1) |
| 1472 HT = ATE2ae*QinE + ATE3ae*h*UinEe + ATE4a*h*VinE | 1704 print() |
| 1473 HR = ARE2ae*QinE + ARE3ae*h*UinEe + ARE4a*h*VinE | 1705 |
| 1474 | 1706 iLDR = 5 |
| 1475 elif (TypeC == 3) or (TypeC == 4): # +++++++++ linear polariser calibration Eq. C.5 | 1707 for LDRTrue in LDRrange: |
| 1476 # p = +45°, m = -45° | 1708 iLDR = iLDR - 1 |
| 1477 AT = ATE1*IinE + ZiC*CosC*(ATE2e*QinEe + ATE4*VinE) + ATE3e*UinEe | 1709 aF11corr[iLDR,:] = aF11corr[iLDR,:] / aF11corr[0,:] - 1.0 |
| 1478 BT = DiC*(ATE1*UinEe + ATE3e*IinE) + ZiC*SinC*(ATE4*QinEe - ATE2e*VinE) | 1710 |
| 1479 AR = ARE1*IinE + ZiC*CosC*(ARE2e*QinEe + ARE4*VinE) + ARE3e*UinEe | 1711 # Plot F11 |
| 1480 BR = DiC*(ARE1*UinEe + ARE3e*IinE) + ZiC*SinC*(ARE4*QinEe - ARE2e*VinE) | 1712 def PlotSubHistF11(aVar, aX, X0, daX, iaX, naX): |
| 1481 | 1713 fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) |
| 1482 # Correction paremeters for normal measurements; they are independent of LDR | 1714 iLDR = -1 |
| 1483 if (not RotationErrorEpsilonForNormalMeasurements): | 1715 for LDRTrue in LDRrange: |
| 1484 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F) | 1716 iLDR = iLDR + 1 |
| 1485 GT = ATE1o*IinE + ATE4o*VinE | 1717 |
| 1486 GR = ARE1o*IinE + ARE4o*VinE | 1718 #F11min[iLDR] = np.min(aF11corr[iLDR,:]) |
| 1487 HT = ATE2a*QinE + ATE3a*UinE + ATE4a*VinE | 1719 #F11max[iLDR] = np.max(aF11corr[iLDR,:]) |
| 1488 HR = ARE2a*QinE + ARE3a*UinE + ARE4a*VinE | 1720 #Rmin = F11min[iLDR] * 0.995 # np.min(aA[iLDR,:]) * 0.995 |
| 1489 else: | 1721 #Rmax = F11max[iLDR] * 1.005 # np.max(aA[iLDR,:]) * 1.005 |
| 1490 D = IinE + DiC*QinEe | 1722 |
| 1491 A = DiC*IinE + QinEe | 1723 #Rmin = 0.8 |
| 1492 B = ZiC*(CosC*UinEe + SinC*VinE) | 1724 #Rmax = 1.2 |
| 1493 C = -ZiC*(SinC*UinEe - CosC*VinE) | 1725 |
| 1494 GT = ATE1o*D + ATE4o*C | 1726 #plt.subplot(5,2,iLDR+1) |
| 1495 GR = ARE1o*D + ARE4o*C | 1727 plt.subplot(1,5,iLDR+1) |
| 1496 HT = ATE2a*A + ATE3a*B + ATE4a*C | 1728 (n, bins, patches) = plt.hist(aF11corr[iLDR,:], |
| 1497 HR = ARE2a*A + ARE3a*B + ARE4a*C | 1729 bins=100, log=False, |
| 1498 | 1730 alpha=0.5, normed=False, color = '0.5', histtype='stepfilled') |
| 1499 elif (TypeC == 6): # real HWP calibration +-22.5° rotated_diattenuator_X22x5deg.odt | 1731 |
| 1500 # p = +22.5°, m = -22.5° | 1732 for iaX in range(-naX,naX+1): |
| 1501 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]) | 1733 plt.hist(aF11corr[iLDR,aX == iaX], |
| 1502 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]) | 1734 bins=100, log=False, alpha=0.3, normed=False, histtype='stepfilled', label = str(round(X0 + iaX*daX/naX,5))) |
| 1503 ATEe = np.array([ATE1,ATE2e,ATE3e,ATE4]) | 1735 |
| 1504 AREe = np.array([ARE1,ARE2e,ARE3e,ARE4]) | 1736 if (iLDR == 2): plt.legend() |
| 1505 AT = np.dot(ATEe,IE1e) | 1737 |
| 1506 AR = np.dot(AREe,IE1e) | 1738 plt.tick_params(axis='both', labelsize=9) |
| 1507 BT = np.dot(ATEe,IE2e) | 1739 #plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2) |
| 1508 BR = np.dot(AREe,IE2e) | 1740 |
| 1509 | 1741 #plt.title(LID + ' ' + aVar, fontsize=18) |
| 1510 # Correction paremeters for normal measurements; they are independent of LDR | 1742 #plt.ylabel('frequency', fontsize=10) |
| 1511 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out | 1743 #plt.xlabel('LDRcorr', fontsize=10) |
| 1512 GT = ATE1o*IinE + ATE4o*VinE | 1744 #fig.tight_layout() |
| 1513 GR = ARE1o*IinE + ARE4o*VinE | 1745 fig.suptitle(LID + ' ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar, fontsize=14, y=1.05) |
| 1514 HT = ATE2a*QinE + ATE3a*UinE + ATE4a*VinE | 1746 #plt.show() |
| 1515 HR = ARE2a*QinE + ARE3a*UinE + ARE4a*VinE | 1747 #fig.savefig(LID + '_' + aVar + '.png', dpi=150, bbox_inches='tight', pad_inches=0) |
| 1516 else: | 1748 #plt.close |
| 1517 D = IinE + DiC*QinEe | 1749 return |
| 1518 A = DiC*IinE + QinEe | 1750 |
| 1519 B = ZiC*(CosC*UinEe + SinC*VinE) | 1751 if (nRotL > 0): PlotSubHistF11("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) |
| 1520 C = -ZiC*(SinC*UinEe - CosC*VinE) | 1752 if (nRetE > 0): PlotSubHistF11("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) |
| 1521 GT = ATE1o*D + ATE4o*C | 1753 if (nRotE > 0): PlotSubHistF11("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) |
| 1522 GR = ARE1o*D + ARE4o*C | 1754 if (nDiE > 0): PlotSubHistF11("DiE", aDiE, DiE0, dDiE, iDiE, nDiE) |
| 1523 HT = ATE2a*A + ATE3a*B + ATE4a*C | 1755 if (nRetO > 0): PlotSubHistF11("RetO", aRetO, RetO0, dRetO, iRetO, nRetO) |
| 1524 HR = ARE2a*A + ARE3a*B + ARE4a*C | 1756 if (nRotO > 0): PlotSubHistF11("RotO", aRotO, RotO0, dRotO, iRotO, nRotO) |
| 1525 | 1757 if (nDiO > 0): PlotSubHistF11("DiO", aDiO, DiO0, dDiO, iDiO, nDiO) |
| 1526 else: | 1758 if (nDiC > 0): PlotSubHistF11("DiC", aDiC, DiC0, dDiC, iDiC, nDiC) |
| 1527 print('Calibrator not implemented yet') | 1759 if (nRotC > 0): PlotSubHistF11("RotC", aRotC, RotC0, dRotC, iRotC, nRotC) |
| 1528 sys.exit() | 1760 if (nRetC > 0): PlotSubHistF11("RetC", aRetC, RetC0, dRetC, iRetC, nRetC) |
| 1529 | 1761 if (nTP > 0): PlotSubHistF11("TP", aTP, TP0, dTP, iTP, nTP) |
| 1530 # Calibration signals with aCal => Determination of the correction K of the real calibration factor | 1762 if (nTS > 0): PlotSubHistF11("TS", aTS, TS0, dTS, iTS, nTS) |
| 1531 IoutTp = TaT*TiT*TiO*TiE*(AT + BT) | 1763 if (nRP > 0): PlotSubHistF11("RP", aRP, RP0, dRP, iRP, nRP) |
| 1532 IoutTm = TaT*TiT*TiO*TiE*(AT - BT) | 1764 if (nRS > 0): PlotSubHistF11("RS", aRS, RS0, dRS, iRS, nRS) |
| 1533 IoutRp = TaR*TiR*TiO*TiE*(AR + BR) | 1765 if (nRetT > 0): PlotSubHistF11("RetT", aRetT, RetT0, dRetT, iRetT, nRetT) |
| 1534 IoutRm = TaR*TiR*TiO*TiE*(AR - BR) | 1766 if (nRetR > 0): PlotSubHistF11("RetR", aRetR, RetR0, dRetR, iRetR, nRetR) |
| 1535 # --- Results and Corrections; electronic etaR and etaT are assumed to be 1 | 1767 if (nERaT > 0): PlotSubHistF11("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) |
| 1536 #Eta = TiR/TiT # Eta = Eta*/K Eq. 84 | 1768 if (nERaR > 0): PlotSubHistF11("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) |
| 1537 Etapx = IoutRp/IoutTp | 1769 if (nRotaT > 0): PlotSubHistF11("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) |
| 1538 Etamx = IoutRm/IoutTm | 1770 if (nRotaR > 0): PlotSubHistF11("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) |
| 1539 Etax = (Etapx*Etamx)**0.5 | 1771 if (nLDRCal > 0): PlotSubHistF11("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) |
| 1540 K = Etax / Eta0 | 1772 |
| 1541 #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)) | 1773 plt.show() |
| 1542 #print("{0:6.3f},{1:6.3f},{2:6.3f},{3:6.3f}".format(DiC, ZiC, Kp, Km)) | 1774 plt.close |
| 1543 | 1775 ''' |
| 1544 # For comparison with Volkers Libreoffice Müller Matrix spreadsheet | 1776 |
| 1545 #Eta_test_p = (IoutRp/IoutTp) | 1777 ''' |
| 1546 #Eta_test_m = (IoutRm/IoutTm) | 1778 # only histogram |
| 1547 #Eta_test = (Eta_test_p*Eta_test_m)**0.5 | 1779 #print("******************* " + aVar + " *******************") |
| 1548 | 1780 fig, ax = plt.subplots(nrows=5, ncols=2, sharex=True, sharey=True, figsize=(10, 10)) |
| 1549 # ************************************************************************* | |
| 1550 iLDR = -1 | |
| 1551 for LDRTrue in LDRrange: | |
| 1552 iLDR = iLDR + 1 | |
| 1553 atrue = (1-LDRTrue)/(1+LDRTrue) | |
| 1554 # ----- Forward simulated signals and LDRsim with atrue; from input file | |
| 1555 It = TaT*TiT*TiO*TiE*(GT+atrue*HT) # TaT*TiT*TiC*TiO*IinL*(GT+atrue*HT) | |
| 1556 Ir = TaR*TiR*TiO*TiE*(GR+atrue*HR) # TaR*TiR*TiC*TiO*IinL*(GR+atrue*HR) | |
| 1557 | |
| 1558 # LDRsim = 1/Eta*Ir/It # simulated LDR* with Y from input file | |
| 1559 LDRsim = Ir/It # simulated uncorrected LDR with Y from input file | |
| 1560 ''' | |
| 1561 if Y == 1.: | |
| 1562 LDRsimx = LDRsim | |
| 1563 LDRsimx2 = LDRsim2 | |
| 1564 else: | |
| 1565 LDRsimx = 1./LDRsim | |
| 1566 LDRsimx2 = 1./LDRsim2 | |
| 1567 ''' | |
| 1568 # ----- Backward correction | |
| 1569 # Corrected LDRCorr from forward simulated LDRsim (atrue) with assumed true G0,H0,K0 | |
| 1570 LDRCorr = (LDRsim*K0/Etax*(GT0+HT0)-(GR0+HR0))/((GR0-HR0)-LDRsim*K0/Etax*(GT0-HT0)) | |
| 1571 | |
| 1572 # -- F11corr from It and Ir and calibration EtaX | |
| 1573 Text1 = "!!! EXPERIMENTAL !!! F11corr from It and Ir with calibration EtaX: x-axis: F11corr(LDRtrue) / F11corr(LDRtrue = 0.004) - 1" | |
| 1574 F11corr = 1/(TiO*TiE)*((HR0*Etax/K0*It/TTa-HT0*Ir/TRa)/(HR0*GT0-HT0*GR0)) # IL = 1 Eq.(64) | |
| 1575 | |
| 1576 #Text1 = "F11corr from It and Ir without corrections but with calibration EtaX: x-axis: F11corr(LDRtrue) devided by F11corr(LDRtrue = 0.004)" | |
| 1577 #F11corr = 0.5/(TiO*TiE)*(Etax*It/TTa+Ir/TRa) # IL = 1 Eq.(64) | |
| 1578 | |
| 1579 # -- It from It only with atrue without corrections - for BERTHA (and PollyXTs) | |
| 1580 #Text1 = " x-axis: IT(LDRtrue) / IT(LDRtrue = 0.004) - 1" | |
| 1581 #F11corr = It/(TaT*TiT*TiO*TiE) #/(TaT*TiT*TiO*TiE*(GT0+atrue*HT0)) | |
| 1582 # !!! see below line 1673ff | |
| 1583 | |
| 1584 aF11corr[iLDR,iN] = F11corr | |
| 1585 aA[iLDR,iN] = LDRCorr | |
| 1586 | |
| 1587 aX[0,iN] = GR | |
| 1588 aX[1,iN] = GT | |
| 1589 aX[2,iN] = HR | |
| 1590 aX[3,iN] = HT | |
| 1591 aX[4,iN] = K | |
| 1592 | |
| 1593 aLDRCal[iN] = iLDRCal | |
| 1594 aERaT[iN] = iERaT | |
| 1595 aERaR[iN] = iERaR | |
| 1596 aRotaT[iN] = iRotaT | |
| 1597 aRotaR[iN] = iRotaR | |
| 1598 aRetT[iN] = iRetT | |
| 1599 aRetR[iN] = iRetR | |
| 1600 | |
| 1601 aRotL[iN] = iRotL | |
| 1602 aRotE[iN] = iRotE | |
| 1603 aRetE[iN] = iRetE | |
| 1604 aRotO[iN] = iRotO | |
| 1605 aRetO[iN] = iRetO | |
| 1606 aRotC[iN] = iRotC | |
| 1607 aRetC[iN] = iRetC | |
| 1608 aDiO[iN] = iDiO | |
| 1609 aDiE[iN] = iDiE | |
| 1610 aDiC[iN] = iDiC | |
| 1611 aTP[iN] = iTP | |
| 1612 aTS[iN] = iTS | |
| 1613 aRP[iN] = iRP | |
| 1614 aRS[iN] = iRS | |
| 1615 | |
| 1616 # --- END loop | |
| 1617 btime = clock() | |
| 1618 print("\r done in ", "{0:5.0f}".format(btime-atime), "sec") #, end="\r") | |
| 1619 | |
| 1620 # --- Plot ----------------------------------------------------------------- | |
| 1621 if (sns_loaded): | |
| 1622 sns.set_style("whitegrid") | |
| 1623 sns.set_palette("bright", 6) | |
| 1624 | |
| 1625 ''' | |
| 1626 fig2 = plt.figure() | |
| 1627 plt.plot(aA[2,:],'b.') | |
| 1628 plt.plot(aA[3,:],'r.') | |
| 1629 plt.plot(aA[4,:],'g.') | |
| 1630 #plt.plot(aA[6,:],'c.') | |
| 1631 plt.show | |
| 1632 ''' | |
| 1633 # Plot LDR | |
| 1634 def PlotSubHist(aVar, aX, X0, daX, iaX, naX): | |
| 1635 fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) | |
| 1636 iLDR = -1 | 1781 iLDR = -1 |
| 1637 for LDRTrue in LDRrange: | 1782 for LDRTrue in LDRrange: |
| 1638 iLDR = iLDR + 1 | 1783 iLDR = iLDR + 1 |
| 1639 | |
| 1640 LDRmin[iLDR] = np.min(aA[iLDR,:]) | 1784 LDRmin[iLDR] = np.min(aA[iLDR,:]) |
| 1641 LDRmax[iLDR] = np.max(aA[iLDR,:]) | 1785 LDRmax[iLDR] = np.max(aA[iLDR,:]) |
| 1642 Rmin = LDRmin[iLDR] * 0.995 # np.min(aA[iLDR,:]) * 0.995 | 1786 Rmin = np.min(aA[iLDR,:]) * 0.999 |
| 1643 Rmax = LDRmax[iLDR] * 1.005 # np.max(aA[iLDR,:]) * 1.005 | 1787 Rmax = np.max(aA[iLDR,:]) * 1.001 |
| 1644 | 1788 plt.subplot(5,2,iLDR+1) |
| 1645 #plt.subplot(5,2,iLDR+1) | |
| 1646 plt.subplot(1,5,iLDR+1) | |
| 1647 (n, bins, patches) = plt.hist(aA[iLDR,:], | 1789 (n, bins, patches) = plt.hist(aA[iLDR,:], |
| 1648 bins=100, log=False, | |
| 1649 range=[Rmin, Rmax], | 1790 range=[Rmin, Rmax], |
| 1650 alpha=0.5, normed=False, color = '0.5', histtype='stepfilled') | 1791 bins=200, log=False, alpha=0.2, normed=False, color = '0.5', histtype='stepfilled') |
| 1651 | |
| 1652 for iaX in range(-naX,naX+1): | |
| 1653 plt.hist(aA[iLDR,aX == iaX], | |
| 1654 range=[Rmin, Rmax], | |
| 1655 bins=100, log=False, alpha=0.3, normed=False, histtype='stepfilled', label = str(round(X0 + iaX*daX/naX,5))) | |
| 1656 | |
| 1657 if (iLDR == 2): plt.legend() | |
| 1658 | |
| 1659 plt.tick_params(axis='both', labelsize=9) | 1792 plt.tick_params(axis='both', labelsize=9) |
| 1660 plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2) | 1793 plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2) |
| 1661 | 1794 plt.show() |
| 1662 #plt.title(LID + ' ' + aVar, fontsize=18) | 1795 plt.close |
| 1663 #plt.ylabel('frequency', fontsize=10) | 1796 ''' |
| 1664 #plt.xlabel('LDRcorr', fontsize=10) | 1797 |
| 1665 #fig.tight_layout() | 1798 # --- Plot LDRmin, LDRmax |
| 1666 fig.suptitle(LID + ' ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar, fontsize=14, y=1.05) | 1799 fig2 = plt.figure() |
| 1667 #plt.show() | 1800 plt.plot(LDRrange,LDRmax-LDRrange, linewidth=2.0, color='b') |
| 1668 #fig.savefig(LID + '_' + aVar + '.png', dpi=150, bbox_inches='tight', pad_inches=0) | 1801 plt.plot(LDRrange,LDRmin-LDRrange, linewidth=2.0, color='g') |
| 1669 #plt.close | 1802 |
| 1670 return | 1803 plt.xlabel('LDRtrue', fontsize=18) |
| 1671 | 1804 plt.ylabel('LDRTrue-LDRmin, LDRTrue-LDRmax', fontsize=14) |
| 1672 if (nRotL > 0): PlotSubHist("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) | 1805 plt.title(LID + ' ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]), fontsize=18) |
| 1673 if (nRetE > 0): PlotSubHist("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) | 1806 #plt.ylimit(-0.07, 0.07) |
| 1674 if (nRotE > 0): PlotSubHist("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) | 1807 plt.show() |
| 1675 if (nDiE > 0): PlotSubHist("DiE", aDiE, DiE0, dDiE, iDiE, nDiE) | 1808 plt.close |
| 1676 if (nRetO > 0): PlotSubHist("RetO", aRetO, RetO0, dRetO, iRetO, nRetO) | 1809 |
| 1677 if (nRotO > 0): PlotSubHist("RotO", aRotO, RotO0, dRotO, iRotO, nRotO) | 1810 # --- Save LDRmin, LDRmax to file |
| 1678 if (nDiO > 0): PlotSubHist("DiO", aDiO, DiO0, dDiO, iDiO, nDiO) | 1811 # http://stackoverflow.com/questions/4675728/redirect-stdout-to-a-file-in-python |
| 1679 if (nDiC > 0): PlotSubHist("DiC", aDiC, DiC0, dDiC, iDiC, nDiC) | 1812 with open('output_files\LDR_min_max_' + LID + '.dat', 'w') as f: |
| 1680 if (nRotC > 0): PlotSubHist("RotC", aRotC, RotC0, dRotC, iRotC, nRotC) | 1813 with redirect_stdout(f): |
| 1681 if (nRetC > 0): PlotSubHist("RetC", aRetC, RetC0, dRetC, iRetC, nRetC) | 1814 print(LID) |
| 1682 if (nTP > 0): PlotSubHist("TP", aTP, TP0, dTP, iTP, nTP) | 1815 print("LDRtrue, LDRmin, LDRmax") |
| 1683 if (nTS > 0): PlotSubHist("TS", aTS, TS0, dTS, iTS, nTS) | 1816 for i in range(len(LDRrange)): |
| 1684 if (nRP > 0): PlotSubHist("RP", aRP, RP0, dRP, iRP, nRP) | 1817 print("{0:7.4f},{1:7.4f},{2:7.4f}".format(LDRrange[i], LDRmin[i], LDRmax[i])) |
| 1685 if (nRS > 0): PlotSubHist("RS", aRS, RS0, dRS, iRS, nRS) | 1818 |
| 1686 if (nRetT > 0): PlotSubHist("RetT", aRetT, RetT0, dRetT, iRetT, nRetT) | 1819 ''' |
| 1687 if (nRetR > 0): PlotSubHist("RetR", aRetR, RetR0, dRetR, iRetR, nRetR) | 1820 # --- Plot K over LDRCal |
| 1688 if (nERaT > 0): PlotSubHist("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) | 1821 fig3 = plt.figure() |
| 1689 if (nERaR > 0): PlotSubHist("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) | 1822 plt.plot(LDRCal0+aLDRCal*dLDRCal/nLDRCal,aX[4,:], linewidth=2.0, color='b') |
| 1690 if (nRotaT > 0): PlotSubHist("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) | 1823 |
| 1691 if (nRotaR > 0): PlotSubHist("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) | 1824 plt.xlabel('LDRCal', fontsize=18) |
| 1692 if (nLDRCal > 0): PlotSubHist("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) | 1825 plt.ylabel('K', fontsize=14) |
| 1693 | 1826 plt.title(LID, fontsize=18) |
| 1694 plt.show() | 1827 plt.show() |
| 1695 plt.close | 1828 plt.close |
| 1696 | 1829 ''' |
| 1697 print() | |
| 1698 #print("IT(LDRtrue) devided by IT(LDRtrue = 0.004)") | |
| 1699 print(" ############################################################################## ") | |
| 1700 print(Text1) | |
| 1701 print() | |
| 1702 | |
| 1703 iLDR = 5 | |
| 1704 for LDRTrue in LDRrange: | |
| 1705 iLDR = iLDR - 1 | |
| 1706 aF11corr[iLDR,:] = aF11corr[iLDR,:] / aF11corr[0,:] - 1.0 | |
| 1707 | |
| 1708 # Plot F11 | |
| 1709 def PlotSubHistF11(aVar, aX, X0, daX, iaX, naX): | |
| 1710 fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2)) | |
| 1711 iLDR = -1 | |
| 1712 for LDRTrue in LDRrange: | |
| 1713 iLDR = iLDR + 1 | |
| 1714 | |
| 1715 ''' | |
| 1716 F11min[iLDR] = np.min(aF11corr[iLDR,:]) | |
| 1717 F11max[iLDR] = np.max(aF11corr[iLDR,:]) | |
| 1718 Rmin = F11min[iLDR] * 0.995 # np.min(aA[iLDR,:]) * 0.995 | |
| 1719 Rmax = F11max[iLDR] * 1.005 # np.max(aA[iLDR,:]) * 1.005 | |
| 1720 ''' | |
| 1721 #Rmin = 0.8 | |
| 1722 #Rmax = 1.2 | |
| 1723 | |
| 1724 #plt.subplot(5,2,iLDR+1) | |
| 1725 plt.subplot(1,5,iLDR+1) | |
| 1726 (n, bins, patches) = plt.hist(aF11corr[iLDR,:], | |
| 1727 bins=100, log=False, | |
| 1728 alpha=0.5, normed=False, color = '0.5', histtype='stepfilled') | |
| 1729 | |
| 1730 for iaX in range(-naX,naX+1): | |
| 1731 plt.hist(aF11corr[iLDR,aX == iaX], | |
| 1732 bins=100, log=False, alpha=0.3, normed=False, histtype='stepfilled', label = str(round(X0 + iaX*daX/naX,5))) | |
| 1733 | |
| 1734 if (iLDR == 2): plt.legend() | |
| 1735 | |
| 1736 plt.tick_params(axis='both', labelsize=9) | |
| 1737 #plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2) | |
| 1738 | |
| 1739 #plt.title(LID + ' ' + aVar, fontsize=18) | |
| 1740 #plt.ylabel('frequency', fontsize=10) | |
| 1741 #plt.xlabel('LDRcorr', fontsize=10) | |
| 1742 #fig.tight_layout() | |
| 1743 fig.suptitle(LID + ' ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar, fontsize=14, y=1.05) | |
| 1744 #plt.show() | |
| 1745 #fig.savefig(LID + '_' + aVar + '.png', dpi=150, bbox_inches='tight', pad_inches=0) | |
| 1746 #plt.close | |
| 1747 return | |
| 1748 | |
| 1749 if (nRotL > 0): PlotSubHistF11("RotL", aRotL, RotL0, dRotL, iRotL, nRotL) | |
| 1750 if (nRetE > 0): PlotSubHistF11("RetE", aRetE, RetE0, dRetE, iRetE, nRetE) | |
| 1751 if (nRotE > 0): PlotSubHistF11("RotE", aRotE, RotE0, dRotE, iRotE, nRotE) | |
| 1752 if (nDiE > 0): PlotSubHistF11("DiE", aDiE, DiE0, dDiE, iDiE, nDiE) | |
| 1753 if (nRetO > 0): PlotSubHistF11("RetO", aRetO, RetO0, dRetO, iRetO, nRetO) | |
| 1754 if (nRotO > 0): PlotSubHistF11("RotO", aRotO, RotO0, dRotO, iRotO, nRotO) | |
| 1755 if (nDiO > 0): PlotSubHistF11("DiO", aDiO, DiO0, dDiO, iDiO, nDiO) | |
| 1756 if (nDiC > 0): PlotSubHistF11("DiC", aDiC, DiC0, dDiC, iDiC, nDiC) | |
| 1757 if (nRotC > 0): PlotSubHistF11("RotC", aRotC, RotC0, dRotC, iRotC, nRotC) | |
| 1758 if (nRetC > 0): PlotSubHistF11("RetC", aRetC, RetC0, dRetC, iRetC, nRetC) | |
| 1759 if (nTP > 0): PlotSubHistF11("TP", aTP, TP0, dTP, iTP, nTP) | |
| 1760 if (nTS > 0): PlotSubHistF11("TS", aTS, TS0, dTS, iTS, nTS) | |
| 1761 if (nRP > 0): PlotSubHistF11("RP", aRP, RP0, dRP, iRP, nRP) | |
| 1762 if (nRS > 0): PlotSubHistF11("RS", aRS, RS0, dRS, iRS, nRS) | |
| 1763 if (nRetT > 0): PlotSubHistF11("RetT", aRetT, RetT0, dRetT, iRetT, nRetT) | |
| 1764 if (nRetR > 0): PlotSubHistF11("RetR", aRetR, RetR0, dRetR, iRetR, nRetR) | |
| 1765 if (nERaT > 0): PlotSubHistF11("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT) | |
| 1766 if (nERaR > 0): PlotSubHistF11("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR) | |
| 1767 if (nRotaT > 0): PlotSubHistF11("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT) | |
| 1768 if (nRotaR > 0): PlotSubHistF11("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR) | |
| 1769 if (nLDRCal > 0): PlotSubHistF11("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal) | |
| 1770 | |
| 1771 plt.show() | |
| 1772 plt.close | |
| 1773 ''' | |
| 1774 # only histogram | |
| 1775 #print("******************* " + aVar + " *******************") | |
| 1776 fig, ax = plt.subplots(nrows=5, ncols=2, sharex=True, sharey=True, figsize=(10, 10)) | |
| 1777 iLDR = -1 | |
| 1778 for LDRTrue in LDRrange: | |
| 1779 iLDR = iLDR + 1 | |
| 1780 LDRmin[iLDR] = np.min(aA[iLDR,:]) | |
| 1781 LDRmax[iLDR] = np.max(aA[iLDR,:]) | |
| 1782 Rmin = np.min(aA[iLDR,:]) * 0.999 | |
| 1783 Rmax = np.max(aA[iLDR,:]) * 1.001 | |
| 1784 plt.subplot(5,2,iLDR+1) | |
| 1785 (n, bins, patches) = plt.hist(aA[iLDR,:], | |
| 1786 range=[Rmin, Rmax], | |
| 1787 bins=200, log=False, alpha=0.2, normed=False, color = '0.5', histtype='stepfilled') | |
| 1788 plt.tick_params(axis='both', labelsize=9) | |
| 1789 plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2) | |
| 1790 plt.show() | |
| 1791 plt.close | |
| 1792 ''' | |
| 1793 | |
| 1794 # --- Plot LDRmin, LDRmax | |
| 1795 fig2 = plt.figure() | |
| 1796 plt.plot(LDRrange,LDRmax-LDRrange, linewidth=2.0, color='b') | |
| 1797 plt.plot(LDRrange,LDRmin-LDRrange, linewidth=2.0, color='g') | |
| 1798 | |
| 1799 plt.xlabel('LDRtrue', fontsize=18) | |
| 1800 plt.ylabel('LDRTrue-LDRmin, LDRTrue-LDRmax', fontsize=14) | |
| 1801 plt.title(LID + ' ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]), fontsize=18) | |
| 1802 #plt.ylimit(-0.07, 0.07) | |
| 1803 plt.show() | |
| 1804 plt.close | |
| 1805 | |
| 1806 # --- Save LDRmin, LDRmax to file | |
| 1807 # http://stackoverflow.com/questions/4675728/redirect-stdout-to-a-file-in-python | |
| 1808 with open('output_files\LDR_min_max_' + LID + '.dat', 'w') as f: | |
| 1809 with redirect_stdout(f): | |
| 1810 print(LID) | |
| 1811 print("LDRtrue, LDRmin, LDRmax") | |
| 1812 for i in range(len(LDRrange)): | |
| 1813 print("{0:7.4f},{1:7.4f},{2:7.4f}".format(LDRrange[i], LDRmin[i], LDRmax[i])) | |
| 1814 | |
| 1815 ''' | |
| 1816 # --- Plot K over LDRCal | |
| 1817 fig3 = plt.figure() | |
| 1818 plt.plot(LDRCal0+aLDRCal*dLDRCal/nLDRCal,aX[4,:], linewidth=2.0, color='b') | |
| 1819 | |
| 1820 plt.xlabel('LDRCal', fontsize=18) | |
| 1821 plt.ylabel('K', fontsize=14) | |
| 1822 plt.title(LID, fontsize=18) | |
| 1823 plt.show() | |
| 1824 plt.close | |
| 1825 ''' | |
| 1826 | 1830 |
| 1827 # Additional plot routines ======> | 1831 # Additional plot routines ======> |
| 1828 ''' | 1832 ''' |
| 1829 #****************************************************************************** | 1833 #****************************************************************************** |
| 1830 # 1. Plot LDRcorrected - LDR(measured Icross/Iparallel) | 1834 # 1. Plot LDRcorrected - LDR(measured Icross/Iparallel) |
