### Factor out example implementation into separate module

parent 20efa2cd
Pipeline #466066 passed with stage
in 1 minute and 8 seconds
python/example.py 0 → 100644
 import numpy as np import matplotlib.pyplot as plt # some initialization wxyz = 1234 np.random.seed(wxyz) # input data nsteps = 32 nsamples = 100000 T0 = 1000 # (ns), maximum time delay W = 1 # (ns), time coincidence window HWP2 = 0 # degrees, orientation of wave plate 2 (EOM2) cHWP2 = np.cos(HWP2*np.pi/180) sHWP2 = np.sin(HWP2*np.pi/180) def analyzer(c, s, cHWP, sHWP, T0): # EOM : plane rotation 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)) # Malus law r0 = np.random.rand() if(x>2*r0 - 1): j = 0 # +1 event else: j = 1 #-1 event # time delay r1 = np.random.rand() l = y**4*T0*r1 # delay time: T_0 sin(2*(theta1 - x))**4 return j, l count = np.zeros(shape=(2, 2, 2, nsteps), dtype=np.int32) # set all counts to zero for ipsi0 in range(nsteps): # loop over different settings of EOM1 cHWP1 = np.cos(ipsi0*2*np.pi/nsteps) sHWP1 = np.sin(ipsi0*2*np.pi/nsteps) for i in range(nsamples): # source r0 = np.random.rand() 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 # station 1 j1, l1 = analyzer(c1, s1, cHWP1, sHWP1, T0) # station 2 j2, l2 = analyzer(c2, s2, cHWP2, sHWP2, T0) # count count[j1, j2, 0, ipsi0] = count[j1, j2, 0, ipsi0] + 1 # Malus law model if(abs(l1 - l2) no time window, i = 1 <=> use time coincidences #tot[i, j] = np.sum(count[:, :, i, j]) tot[i, j] = count[0, 0, i, j] + count[1, 1, i, j] + count[1, 0, i, j] + count[0, 1, i, j] E12[i, j] = count[0, 0, i, j] + count[1, 1, i, j] - count[1, 0, i, j] - count[0, 1, i, j] E1[i, j] = count[0, 0, i, j] + count[0, 1, i, j] - count[1, 1, i, j] - count[1, 0, i, j] E2[i, j] = count[0, 0, i, j] + count[1, 0, i, j] - count[1, 1, i, j] - count[0, 1, i, j] if(tot[i, j]>0): E12[i, j] = E12[i, j]/tot[i, j] E1[i, j] = E1[i, j]/tot[i, j] E2[i, j] = E2[i, j]/tot[i, j] theta = np.linspace(0, 360, nsteps) theta_theory = np.linspace(0, 360, nsteps*100) theory = [] for j in range(nsteps*100): theory.append(-np.cos(2 * j * 2 * np.pi / nsteps / 100)) plt.plot(theta, E12[0, :], 'o') plt.plot(theta_theory, theory, '.', markersize=1) plt.savefig('e12.pdf') plt.close() \ No newline at end of file
 import numpy as np import matplotlib.pyplot as plt # some initialization wxyz = 1234 np.random.seed(wxyz) # input data nsteps = 32 nsamples = 100000 T0 = 1000 # (ns), maximum time delay W = 1 # (ns), time coincidence window HWP2 = 0 # degrees, orientation of wave plate 2 (EOM2) cHWP2 = np.cos(HWP2*np.pi/180) sHWP2 = np.sin(HWP2*np.pi/180) def analyzer(c, s, cHWP, sHWP, T0): # EOM : plane rotation 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)) # Malus law r0 = np.random.rand() if(x>2*r0 - 1): j = 0 # +1 event else: j = 1 #-1 event # time delay r1 = np.random.rand() l = y**4*T0*r1 # delay time: T_0 sin(2*(theta1 - x))**4 return j, l count = np.zeros(shape=(2, 2, 2, nsteps), dtype=np.int32) # set all counts to zero for ipsi0 in range(nsteps): # loop over different settings of EOM1 cHWP1 = np.cos(ipsi0*2*np.pi/nsteps) sHWP1 = np.sin(ipsi0*2*np.pi/nsteps) for i in range(nsamples): # source r0 = np.random.rand() 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 # station 1 j1, l1 = analyzer(c1, s1, cHWP1, sHWP1, T0) # station 2 j2, l2 = analyzer(c2, s2, cHWP2, sHWP2, T0) # count count[j1, j2, 0, ipsi0] = count[j1, j2, 0, ipsi0] + 1 # Malus law model if(abs(l1 - l2) no time window, i = 1 <=> use time coincidences #tot[i, j] = np.sum(count[:, :, i, j]) tot[i, j] = count[0, 0, i, j] + count[1, 1, i, j] + count[1, 0, i, j] + count[0, 1, i, j] E12[i, j] = count[0, 0, i, j] + count[1, 1, i, j] - count[1, 0, i, j] - count[0, 1, i, j] E1[i, j] = count[0, 0, i, j] + count[0, 1, i, j] - count[1, 1, i, j] - count[1, 0, i, j] E2[i, j] = count[0, 0, i, j] + count[1, 0, i, j] - count[1, 1, i, j] - count[0, 1, i, j] if(tot[i, j]>0): E12[i, j] = E12[i, j]/tot[i, j] E1[i, j] = E1[i, j]/tot[i, j] E2[i, j] = E2[i, j]/tot[i, j] theta = np.linspace(0, 360, nsteps) theta_theory = np.linspace(0, 360, nsteps*100) theory = [] for j in range(nsteps*100): theory.append(-np.cos(2 * j * 2 * np.pi / nsteps / 100)) plt.plot(theta, E12[0, :], 'o') plt.plot(theta_theory, theory, '.', markersize=1) plt.savefig('e12.pdf') plt.close() \ No newline at end of file import example
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!