main.py 6.17 KB
 Orkun Şensebat committed May 06, 2021 1 import numpy as np Orkun Şensebat committed May 09, 2021 2 import matplotlib.pyplot as plt Tobias Jörg Theodor Kempe committed May 11, 2021 3 4 5 import time showPlots = True Orkun Şensebat committed May 06, 2021 6 Tobias Jörg Theodor Kempe committed May 11, 2021 7 # initialization Orkun Şensebat committed May 07, 2021 8 9 wxyz = 1234 np.random.seed(wxyz) Tobias Jörg Theodor Kempe committed May 10, 2021 10 Tobias Jörg Theodor Kempe committed May 11, 2021 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 # global settings n_configs = 32 # number of detector configs n_samples = 100000 # number of simulations per config T0 = 1000 # (ns), maximum time delay W = 1 # (ns), time coincidence window # orientation of half-wave plate in observation station 2 #HWP2 = 0 phi2 = 0 ## cosine and sine components of HWP2 #cHWP2 = np.cos(HWP2*np.pi/180) #sHWP2 = np.sin(HWP2*np.pi/180) # simulation station module # ta def simulation_station(theta, phi): # c, s, cHWP, sHWP): Orkun Şensebat committed May 07, 2021 27 # EOM : plane rotation Tobias Jörg Theodor Kempe committed May 10, 2021 28 # this is basically the application of a rotation matrix! Tobias Jörg Theodor Kempe committed May 11, 2021 29 30 31 32 33 34 35 36 # #c2 = cHWP*c + sHWP*s #s2 = -sHWP*c + cHWP*s #x = c2*c2 - s2*s2 # cos(2(x-a)) #y = 2*c2*s2 # sin(2(x-a)) c = np.cos(2*(phi-theta)) s = np.sin(2*(phi-theta)) Tobias Jörg Theodor Kempe committed May 10, 2021 37 38 # when the x component is larger than a certain value, we count it as +1 and else as -1 event Orkun Şensebat committed May 07, 2021 39 # Malus law Tobias Jörg Theodor Kempe committed May 11, 2021 40 41 42 43 44 45 46 47 48 49 50 #if(x>2*np.random.rand() - 1): # x = 0 # +1 event #else: # x = 1 #-1 event x = int(c > 2*np.random.rand() - 1) # computes time delay based on #t = y**4*T0*np.random.rand() # delay time: T_0 sin(2*(theta1 - x))**4 t = s**4*T0*np.random.rand() # delay time: T_0 sin(2*(theta1 - x))**4 return x, t Orkun Şensebat committed May 07, 2021 51 52 Tobias Jörg Theodor Kempe committed May 11, 2021 53 start = time.time() Tobias Jörg Theodor Kempe committed May 11, 2021 54 55 Tobias Jörg Theodor Kempe committed May 10, 2021 56 57 58 59 60 61 # this is our result array, in which we store, how detection combinations we detected during the entire experiment # 1st index: arriving in station1 at bottom (0) or at top (1) # 2nd index: arriving in station2 at bottom (0) or at top (1) # 3rd index: non-coincidences included (0) or not (1) # 4th index: configuration of the polarizors/detectors Tobias Jörg Theodor Kempe committed May 11, 2021 62 counts = np.zeros(shape=(2, 2, 2, n_configs), dtype=np.int32) # set all counts to zero Tobias Jörg Theodor Kempe committed May 10, 2021 63 64 65 66 67 # we split a whole rotation (2PI) into nsteps and make experiments for each step # TODO What do these angles represent? # it seems to be that the angles is the difference between the polarizations angles of both analyzer stations Tobias Jörg Theodor Kempe committed May 11, 2021 68 for delta_psi in range(n_configs): # loop over different settings of EOM1 Tobias Jörg Theodor Kempe committed May 10, 2021 69 Tobias Jörg Theodor Kempe committed May 11, 2021 70 71 72 #cHWP1 = np.cos(delta_psi*2*np.pi/n_configs) #sHWP1 = np.sin(delta_psi*2*np.pi/n_configs) phi1 = delta_psi * 2 * np.pi / n_configs Tobias Jörg Theodor Kempe committed May 10, 2021 73 74 # we want probabilities, basically, so we iterate over a number of random samples Tobias Jörg Theodor Kempe committed May 11, 2021 75 for i in range(n_samples): Tobias Jörg Theodor Kempe committed May 10, 2021 76 Orkun Şensebat committed May 06, 2021 77 # source Tobias Jörg Theodor Kempe committed May 10, 2021 78 79 # we randomize the polarization of the photons, # however they're always entangle, so they will have a phase of pi/2 between them Tobias Jörg Theodor Kempe committed May 11, 2021 80 81 82 83 #r0 = np.random.rand() theta1 = np.random.rand() * 2 * np.pi theta2 = theta1 + np.pi / 2 Tobias Jörg Theodor Kempe committed May 10, 2021 84 85 86 # but these aren't angles, right? # why don't we input the angles into the analzer? # would be more declarative... Tobias Jörg Theodor Kempe committed May 11, 2021 87 88 89 90 #c1 = np.cos(r0*2*np.pi) # polarization angle x of particle going to station 1 #s1 = np.sin(r0*2*np.pi) # polarization angle x + pi/2 of particle going to station 2 #c2 = -s1 #s2 = c1 Tobias Jörg Theodor Kempe committed May 10, 2021 91 92 # produces binary x values along with a timestamp for the current event Tobias Jörg Theodor Kempe committed May 11, 2021 93 94 95 96 97 # for both observation stations #x1, t1 = simulation_station(c1, s1, cHWP1, sHWP1) #x2, t2 = simulation_station(c2, s2, cHWP2, sHWP2) x1, t1 = simulation_station(theta1, phi1) x2, t2 = simulation_station(theta2, phi2) Tobias Jörg Theodor Kempe committed May 10, 2021 98 99 Tobias Jörg Theodor Kempe committed May 11, 2021 100 101 102 103 104 # increases the count of the associated (x1, x2) counter element # (with and without time-window) counts[x1, x2, 0, delta_psi] = counts[x1, x2, 0, delta_psi] + 1 # Malus law model if(abs(t1 - t2)0): E_12[i, j] = E_12[i, j]/E_tot[i, j] E_1[i, j] = E_1[i, j]/E_tot[i, j] E_2[i, j] = E_2[i, j]/E_tot[i, j] # plots the graphs # 1. as expected # 2. from simulation results Tobias Jörg Theodor Kempe committed May 12, 2021 135 136 plot_i = 1 Tobias Jörg Theodor Kempe committed May 11, 2021 137 # correlation Tobias Jörg Theodor Kempe committed May 12, 2021 138 139 phi = np.linspace(0, 360/n_configs*(n_configs-1), n_configs) phi_theory = np.linspace(0, 360/n_configs*(n_configs-1), n_configs*100) Tobias Jörg Theodor Kempe committed May 11, 2021 140 141 theory = -np.cos(2 * phi_theory * (2 * np.pi / 360)) plt.figure() Tobias Jörg Theodor Kempe committed May 11, 2021 142 plt.xlabel(r'$\Delta \varphi$') Tobias Jörg Theodor Kempe committed May 11, 2021 143 144 plt.ylabel(r'$$') plt.plot(phi_theory, theory, '.', markersize=1, color='orange') Tobias Jörg Theodor Kempe committed May 12, 2021 145 plt.plot(phi, E_12[plot_i, :], 'o', color='blue') Tobias Jörg Theodor Kempe committed May 11, 2021 146 plt.savefig('correlation.pdf') Tobias Jörg Theodor Kempe committed May 11, 2021 147 148 if showPlots: plt.show() Tobias Jörg Theodor Kempe committed May 10, 2021 149 150 plt.close() Tobias Jörg Theodor Kempe committed May 11, 2021 151 152 # product of expectation values line = np.zeros(phi_theory.shape) Tobias Jörg Theodor Kempe committed May 10, 2021 153 plt.figure() Tobias Jörg Theodor Kempe committed May 11, 2021 154 plt.xlabel(r'\Delta \varphi') Tobias Jörg Theodor Kempe committed May 11, 2021 155 156 157 158 plt.ylabel(r'$$') plt.ylim([-1, 1]) plt.xticks([0, 45, 90, 135, 180, 225, 270, 305, 360]) plt.plot(phi_theory, line, '.', markersize=1, color='orange') Tobias Jörg Theodor Kempe committed May 12, 2021 159 plt.plot(phi, E_1[plot_i, :]*E_2[plot_i, :], 'o', color='blue') Tobias Jörg Theodor Kempe committed May 11, 2021 160 plt.savefig('expectation_value.pdf') Tobias Jörg Theodor Kempe committed May 11, 2021 161 162 if showPlots: plt.show() Tobias Jörg Theodor Kempe committed May 11, 2021 163 164 165 plt.close() # displays stopwatch Tobias Jörg Theodor Kempe committed May 11, 2021 166 167 duration = time.time() - start print(f'simulation took {round(duration, 2)} seconds')