-
Notifications
You must be signed in to change notification settings - Fork 0
/
exposure.py
93 lines (76 loc) · 3.46 KB
/
exposure.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
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
# coding: utf-8
# 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.
# from IPython.parallel import require
# import BeamOptics as bopt
# @require(bopt)
def exposure(NX=600, NY=400, pitch=20e-6,
wavelength=780e-9, theta=0.005, amp=100, phase=0, trim=1):
from scipy import pi, sin, cos, exp, conjugate, ogrid
from numpy import random, real, imag, array, absolute, log, log10, arange
from numpy.fft import fft, fftshift
from BeamOptics import plane_wave_beam, gaussian_beam
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.backends.backend_pdf
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)
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/wavelength
k1 = [0, 0, k]
k2 = [k*sin(theta), 0, k*cos(theta)]
# pos = array([x, y])
# darkcts = 0.0 # no idea what is reasonable here, just tinkering
# seems best way to model detector is with partial loss (2%) and
# added noise (dark counts)
#total = plane_wave_beam(x, y, 0, amp, k1) + \
#exp(1j*phase)*plane_wave_beam(x, y, 0, 1e-5*amp, k2)
#total = plane_wave_beam(x, y, 0, amp, k1) + \
#exp(1j*phase)*gaussian_beam(x, y, 0, 1e-2*amp, -1.0, xmin/20.0, k2)
total = gaussian_beam(x, y, 0.99, amp, 0.117e-3, 5.397e-6, k1) + \
exp(1j*phase)*gaussian_beam(x, y, 1.07, 1e-1*amp, 0.209e-3, 7.197e-6, k2)
noise = random.random(total.shape) # add some detector noise?
# really not sure how best to add noise, but this gives the FFT
# a realistic look.
intensity = total * total.conjugate() + 0*noise # 0-2 counts per pixel
#+
#darkcts*(random.random([max(shape(x)),max(shape(y))]) +
#1j*random.random([max(shape(x)),max(shape(y))]))
#add dark noise and
# Calculate the relevant pixels:
deltaK = 2*pi/(20e-6 * (1340-trim))
Kx = k*sin(theta)
pixel = Kx / deltaK
print pixel + (1340-trim)/2.0
# print intensity[0,200]
# trim the last few elements before taking the FFT:
# print "sum: ", log10(intensity[:-trim, 200].sum())
# print "xN: ", log10(intensity[:-trim, 200] * (1340-trim))
# complex intensity after FFT2, trim last few elements:
K = fftshift(fft(intensity[:-trim, 200]))
j = arange(len(K))-len(K)/2.0 # frequency index for plot
# print K[512] used to verify that the sum is equal to DC
# normlzd power spectrum
plt.plot(j, log10(absolute(K)/absolute(K).max()),".-")
plt.ylabel("$\log_{10}(K)$")
plt.xlabel("$j$")
#pp = matplotlib.backends.backend_pdf.PdfPages("figure_1.pdf")
#pp.savefig()
#pp.close()
plt.show()
#return K
# Can find relative amplitude of sideband by dividing fourier amplitude
# F(side) / F(DC) is equal to the amplitude ratio SB/DC.
# pixel of interest in FFT is 816.667 (this is 166.667 pixels away from
# the center @ 650) 166.667 pixels corresponds to the off-axis component
# of the wave: ∆K is 2π / (20e-6 * 1300) = 241.660973353061 rad/m
# 241.660973353061 * 166.667 is 40276.90 rad/m which is k * sin(0.005)
# where k = 2π/780e-9 = 8055365 rad/m TADA!
if __name__ == '__main__':
print exposure()