main.py 6.17 KB
Newer Older
1
import numpy as np
Orkun Şensebat's avatar
Orkun Şensebat committed
2
import matplotlib.pyplot as plt
Tobias Jörg Theodor Kempe's avatar
2nd try    
Tobias Jörg Theodor Kempe committed
3
4
5
import time

showPlots = True
6

7
# initialization
Orkun Şensebat's avatar
Orkun Şensebat committed
8
9
wxyz = 1234
np.random.seed(wxyz)
Tobias Jörg Theodor Kempe's avatar
Tobias Jörg Theodor Kempe committed
10

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's avatar
Orkun Şensebat committed
27
    # EOM : plane rotation
Tobias Jörg Theodor Kempe's avatar
Tobias Jörg Theodor Kempe committed
28
    # this is basically the application of a rotation matrix!
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's avatar
Tobias Jörg Theodor Kempe committed
37
38
    
    # when the x component is larger than a certain value, we count it as +1 and else as -1 event
Orkun Şensebat's avatar
Orkun Şensebat committed
39
    # Malus law
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's avatar
Orkun Şensebat committed
51
52


Tobias Jörg Theodor Kempe's avatar
2nd try    
Tobias Jörg Theodor Kempe committed
53
start = time.time()
54
55


Tobias Jörg Theodor Kempe's avatar
Tobias Jörg Theodor Kempe committed
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

62
counts = np.zeros(shape=(2, 2, 2, n_configs), dtype=np.int32) # set all counts to zero
Tobias Jörg Theodor Kempe's avatar
Tobias Jörg Theodor Kempe committed
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

68
for delta_psi in range(n_configs): # loop over different settings of EOM1
Tobias Jörg Theodor Kempe's avatar
Tobias Jörg Theodor Kempe committed
69

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's avatar
Tobias Jörg Theodor Kempe committed
73
74
    
    # we want probabilities, basically, so we iterate over a number of random samples
75
    for i in range(n_samples):
Tobias Jörg Theodor Kempe's avatar
Tobias Jörg Theodor Kempe committed
76
        
77
        # source
Tobias Jörg Theodor Kempe's avatar
Tobias Jörg Theodor Kempe committed
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
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's avatar
Tobias Jörg Theodor Kempe committed
84
85
86
        # but these aren't angles, right?
        # why don't we input the angles into the analzer?
        # would be more declarative...
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's avatar
Tobias Jörg Theodor Kempe committed
91
92
        
        # produces binary x values along with a timestamp for the current event
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's avatar
Tobias Jörg Theodor Kempe committed
98
99
        
        
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)<W):
            counts[x1, x2, 1, delta_psi] = counts[x1, x2, 1, delta_psi] + 1 # Malus law model + time window
Orkun Şensebat's avatar
Orkun Şensebat committed
105

106
# data analysis
107
108
109
110
111
112
113
114
115
116
# 1st index: no time window (0), use time coincidences (1)
E_tot = np.zeros(shape=(2, n_configs), dtype=np.int32)
E_12 = np.zeros(shape=(2, n_configs), dtype=np.float64)
E_1 = np.zeros(shape=(2, n_configs), dtype=np.float64)
E_2 = np.zeros(shape=(2, n_configs), dtype=np.float64)

# builds the expectation values for each station as well as the correlation between them;
# does this for all configurations and both with and without time-window
# uses formulas from lecture 5, page 25
for j in range(n_configs):
Orkun Şensebat's avatar
Orkun Şensebat committed
117
    for i in [0, 1]:
118
119
120
121
122
123

        # calculates numerators
        E_tot[i, j] = counts[0, 0, i, j] + counts[1, 1, i, j] + counts[1, 0, i, j] + counts[0, 1, i, j]
        E_12[i, j] = counts[0, 0, i, j] + counts[1, 1, i, j] - counts[1, 0, i, j] - counts[0, 1, i, j]
        E_1[i, j] = counts[0, 0, i, j] + counts[0, 1, i, j] - counts[1, 1, i, j] - counts[1, 0, i, j]
        E_2[i, j] = counts[0, 0, i, j] + counts[1, 0, i, j] - counts[1, 1, i, j] - counts[0, 1, i, j]
Tobias Jörg Theodor Kempe's avatar
Tobias Jörg Theodor Kempe committed
124
        
125
126
127
128
129
130
131
132
133
134
        # normalizes (while preventing 0-division)
        if(E_tot[i, j]>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

135
136
plot_i = 1

137
# correlation
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)
140
141
theory = -np.cos(2 * phi_theory * (2 * np.pi / 360))
plt.figure()
Tobias Jörg Theodor Kempe's avatar
minor    
Tobias Jörg Theodor Kempe committed
142
plt.xlabel(r'$\Delta \varphi$')
143
144
plt.ylabel(r'$<x_1 x_2>$')
plt.plot(phi_theory, theory, '.', markersize=1, color='orange')
145
plt.plot(phi, E_12[plot_i, :], 'o', color='blue')
146
plt.savefig('correlation.pdf')
Tobias Jörg Theodor Kempe's avatar
2nd try    
Tobias Jörg Theodor Kempe committed
147
148
if showPlots:
    plt.show()
Tobias Jörg Theodor Kempe's avatar
Tobias Jörg Theodor Kempe committed
149
150
plt.close()

151
152
# product of expectation values
line = np.zeros(phi_theory.shape)
Tobias Jörg Theodor Kempe's avatar
Tobias Jörg Theodor Kempe committed
153
plt.figure()
Tobias Jörg Theodor Kempe's avatar
minor    
Tobias Jörg Theodor Kempe committed
154
plt.xlabel(r'$\Delta \varphi$')
155
156
157
158
plt.ylabel(r'$<x_1><x_2>$')
plt.ylim([-1, 1])
plt.xticks([0, 45, 90, 135, 180, 225, 270, 305, 360])
plt.plot(phi_theory, line, '.', markersize=1, color='orange')
159
plt.plot(phi, E_1[plot_i, :]*E_2[plot_i, :], 'o', color='blue')
160
plt.savefig('expectation_value.pdf')
Tobias Jörg Theodor Kempe's avatar
2nd try    
Tobias Jörg Theodor Kempe committed
161
162
if showPlots:
    plt.show()
163
164
165
plt.close()

# displays stopwatch
Tobias Jörg Theodor Kempe's avatar
2nd try    
Tobias Jörg Theodor Kempe committed
166
167
duration = time.time() - start
print(f'simulation took {round(duration, 2)} seconds')