-
Notifications
You must be signed in to change notification settings - Fork 0
/
ExcelonCamera.py
58 lines (49 loc) · 1.94 KB
/
ExcelonCamera.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
# we're going to model the data coming from an
# Excelon Camera at 1340x400 pixels. This is to get familiar
# with data from the quantum array detection experiments.
# coding=utf-8
from pylab import *
from scipy import *
import BeamOptics as bopt
from numpy import random, real, imag
import Qfunction as qf
NX = 1024 #select an FFT-friendly subset of the array
NY = 400
pitch = 20e-6 # 20 micron CCD pitch in meters
px, py = ogrid[0:NX,0:NY] # pixel index
center = [NX/2, NY/2] # center of the CCD sensor
#physical coordinates, referenced from 0,0 (center of CCD), in meters
xmin = -(NX/2)*pitch
xmax = (NX/2)*pitch
ymin = -(NY/2)*pitch
ymax = (NY/2)*pitch
x,y = ogrid[xmin:xmax:pitch, ymin:ymax:pitch]
k = 2*pi/780e-9 # in radians per meter
theta = 0.005 #radian measure of beam angle
k1 = [0,0,k]
k2 = [k*sin(theta),0,k*cos(theta)]
pos = array([x,y])
amp = 380
darkcts = 0.0 # no idea what is reasonable here, just tinkering
# the ccd specifies 2-3 electrons per pixel per hour
# seems best way to model (classical) detector is with partial loss (2%)
# and added noise (dark counts)
values = []
for i in range(5000): #this is the loop to parallelize
print i
total = bopt.plane_wave_beam(x,y,0,amp,k1) +
exp(1j*1)*bopt.plane_wave_beam(x,y,0,0.001*amp,k2)
intensity = total * total.conjugate() # +
#darkcts*(random.random([max(shape(x)),max(shape(y))]) +
#1j*random.random([max(shape(x)),max(shape(y))]))
# add dark noise and QE
K = fftshift(fft(intensity[:,200])) # complex intensity after FFT2
values.append(K[643]/1e5)
# pixel of interest in FFT is 643 (this is 131 pixels away from the
# center @ 512) 131 pixels corresponds to the off-axis component of
# the wave:
# ∆K is 2π / (20e-6 * 1024) = 306.796 rad/m
# 306.796 * 131 is 40190 rad/m which is k * sin(0.005)
# where k = 2π/780e-9 = 8055365 rad/m TADA!
qfuncoutput = qf.qfuncimage(real(values),imag(values),30)
show()