I am trying to make a circular diffraction pattern, which has a central spot surrounded by a series of rings. It involves a Bessel integral to do it which is defined in the code.
My problems is that it takes too long like I waited 10 min for the code to run but didn't get anything to display. I understand it is because my Bessel integral has 1000 iterations per point can any one help with this ?
Am I on the right track?
I am trying to self teach myself python and computational physics via Mark Newmans book Computational Physics the exercise is 5.4 of Computational Physics.Here is a link to the chapter. It is on page 9. http://www-personal.umich.edu/~mejn/cp/chapters/int.pdf
Here is the Image I am trying to make.
.
My Code:
import numpy as np
import pylab as plt
import math as mathy
#N = number of slicices
#h = b-a/N
def J(m,x): #Bessel Integral
def f(theta):
return (1/mathy.pi)*mathy.cos(m*theta - x*mathy.sin(theta)) #I replaced np. with mathy. for this line
N = 1000
A = 0
a=0
b=mathy.pi
h = ((b-a)/N)
for k in range(1,N,2):
A += 4*f(a + h*k)
for k in range(2,N,2):
A += 2*f(a + h*k)
Idx = (h/3)*(A + f(a)+f(b))
return Idx
def I(lmda,r): #Intensity
k = (mathy.pi/lmda)
return ((J(1,k*r))/(k*r))**2
wavelength = .5 # microm meters
I0 = 1
points = 500
sepration = 0.2
Intensity = np.empty([points,points],np.float)
for i in range(points):
y = sepration*i
for j in range(points):
x = sepration*j
r = np.sqrt((x)**2+(y)**2)
if r < 0.000000000001:
Intensity[i,j]= 0.5 #this is the lim as r -> 0, I -> 0.5
else:
Intensity[i,j] = I0*I(wavelength,r)
plt.imshow(Intensity,vmax=0.01,cmap='hot')
plt.gray()
plt.show()
Your code seems to run fine - if I reduce
N
down to 100 (from 1000) and the image size (points
) down to 50 (from 500). I get the following image after around 4s of execution time:Here's what we get when profiling your code with cProfile:
So it looks like most of your execution time is spent in
f
. You could possibly optimise this function, or alternatively try running your code using PyPy. PyPy is excellent at optimising this sort of thing. You need to install their version of numpy (see http://pypy.readthedocs.org/en/latest/getting-started.html#). But PyPy completes your original code in 40s on my machine (with the plotting stuff removed).EDIT
I haven't got plotlib installed in PyPy on my system, so I replaced your plt calls at the end with
And created a separate program that I execute with normal Python containing:
This produces the image below with the following modifications to your code: