Density plot with python making a Diffraction pattern with Bessel Integrals but it wont stop running

Question

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()

Answer 1

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:

$ python -m cProfile -s time bessel.py | head -n 10
         361598 function calls (359660 primitive calls) in 3.470 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
   252399    2.250    0.000    2.250    0.000 bessel.py:24(f)
     2499    0.821    0.000    3.079    0.001 bessel.py:23(J)
        1    0.027    0.027    3.472    3.472 bessel.py:15(<module>)
     2499    0.015    0.000    3.093    0.001 bessel.py:45(I)
        1    0.013    0.013    0.013    0.013 backend_macosx.py:1(<module>)

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

#plt.imshow(Intensity,vmax=0.01,cmap='hot')
#plt.gray()
#plt.show()
np.savetxt('bessel.txt', Intensity)

And created a separate program that I execute with normal Python containing:

import numpy as np
import pylab as plt

Intensity = np.loadtxt('bessel.txt')

plt.imshow(Intensity,vmax=0.01,cmap='hot')
plt.gray()
plt.show()

This produces the image below with the following modifications to your code:

sepration = 0.008 # decrease separation

Intensity = np.empty([points,points],np.float)

for i in range(points):
    y = sepration*(i - points/2) # centre on image
    for j in range(points):
        x = sepration*(j - points/2)