This code plots the diffraction pattern produced by a point-like source of light in a telescope

telesc_diffr.py

import numpy as np
import scipy.special as sp
from pylab import imshow,show,gray,hot
import matplotlib.pyplot as plt
from timeit import default_timer as timer

def func(x,m,theta):
    return np.cos(m*theta-x*np.sin(theta))

wl = 500E-9     # wavelength, in meters
size = 2.0      # image size, in micrometers
N = 500         # plot resolution
G = 5E12        # Gigantic constant

start = timer()
dists = np.array([[np.sqrt(float((N/2-i)**2+(N/2-j)**2)) for i in range(N)] for j in range(N)])
dists = 1E-6*np.sqrt(size)*dists/np.max(dists)     # distances matrix

image = np.zeros([N,N],float)
for i in range(N):
    for j in range(N):
        if abs(i-N/2) < 2 and abs(j-N/2)<2:     # This deals with the limit r -> 0
            image[i,j] = 0.5**2
        else:
            image[i,j] = (sp.jv(1,G*2.*np.pi*wl*dists[i,j])/(G*2.*np.pi*wl*dists[i,j]))**2
time = timer()-start
print "Calculation time = %f seconds" % time

imshow(image,vmax=0.01)
gray()
show()