In [36]:
import numpy as np
from scipy import special
from matplotlib.pylab import plt
%matplotlib inline
In [5]:
D = 6.5 #diameter of the telescope
In [3]:
#Referring to the image in Fig 10.5: 1D PSF
nu=np.arange(1,70)*np.pi/10.
#print x
i_P = (2*special.j1(nu)/nu)**2
plt.plot(nu/np.pi,i_P)
plt.yscale('log')
ax = plt.gca()
ax.set_ylim([0.0001,1])
plt.xlabel(r'$\nu / \pi$')
plt.ylabel('log(PSF)')
plt.show()
In [4]:
# 2D PSF
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
from scipy import special
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
u = np.linspace(0.1, 7 * np.pi, 100)#radius
v = np.linspace(0, 2 * np.pi, 100)#angle
x = np.outer(u, np.sin(v))#generating polar grid
y = np.outer(u, np.cos(v))
#x, y = np.meshgrid(u*np.sin(v),u*np.cos(v),sparse=True)
z = np.outer((2*special.j1(u)/u)**2,np.ones(np.size(u)))
Z = np.log10(z)
ax.plot_surface(x/np.pi,y/np.pi,z,cmap=plt.cm.jet)
#zticks = [1e-4, 1e-3, 1e-2, 1e-1, 1e-0]
#ax.set_zticks(np.log10(zticks))
#ax.set_zticklabels(zticks)
plt.xlabel(r'$\nu / \pi$')
plt.ylabel(r'$\nu / \pi$')
#ax.set_zlim(0, 0.1)
ax.set_zlabel('log(PSF)')
plt.show()
In [9]:
# Compute the FFT of the image
OTF = np.fft.fft2(z).real
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x/np.pi,y/np.pi,OTF,cmap=plt.cm.jet)
#zticks = [1e-4, 1e-3, 1e-2, 1e-1, 1e-0]
#ax.set_zticks(np.log10(zticks))
#ax.set_zticklabels(zticks)
plt.xlabel(r'$\nu / \pi$')
plt.ylabel(r'$\nu / \pi$')
#ax.set_zlim(0, 0.1)
plt.show()
In [13]:
# Contour plot of the intensity projection on the image plane
plt.figure(figsize=(7, 6), dpi=80, facecolor='w', edgecolor='k')
ax = fig.add_subplot(111)
CS = plt.contourf(x, y, z,cmap=plt.get_cmap('gray'))
#plt.clabel(CS, inline=1, fontsize=10)
plt.title('log(PSF) of the Airy Disk')
ax.set_aspect('equal')
plt.colorbar()
plt.show()
In [16]:
# 2D PSF
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
from scipy import special
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
u = np.linspace(-7*np.pi, 7*np.pi, 100) #radius
v = np.linspace(-7*np.pi, 7*np.pi, 100) #angle
x, y = np.meshgrid(u, v, sparse=True)
z = np.outer((np.sin(u)/u)**2,(np.sin(v)/v)**2)
Z = np.log10(z)
ax.plot_surface(x/np.pi,y/np.pi,z,cmap=plt.cm.jet)
# zticks = [1e-4, 1e-3, 1e-2, 1e-1, 1e-0]
# ax.set_zticks(np.log10(zticks))
# ax.set_zticklabels(zticks)
plt.xlabel(r'$\nu / \pi$')
plt.ylabel(r'$\nu / \pi$')
ax.set_zlabel('log(PSF)')
ax.set_zlim(0, 0.1)
plt.show()
Test the Inverse Fourier Transform and its form
In [37]:
from PIL import Image
from numpy import eye
arr = (eye(200)*255).astype('uint8') # sample array
im = Image.fromarray(arr) # monochromatic image
imrgb = Image.merge('RGB', (im,im,im)) # color image
imrgb.show()
In [46]:
# generate the uniform distribution code
n = 10
a = np.ones((n, n))
plt.imshow(a)
plt.show()
In [52]:
# Linear aperture pillbox
n = 501 # number of the total pixels in this image
# rectangular
h = np.ones((n, n))
# circular mask
a, b = (n-1)/2, (n-1)/2 # center of the mask
r = 10 # radius of the circular obscuration
y,x = np.ogrid[-a:n-a, -b:n-b]
mask = x*x + y*y <= r*r
array = np.ones((n, n))
array[mask] = 255
plt.imshow(array)
plt.show()
In [66]:
# OTF is the 2d Fourier transform of the PSF, MTF is its modulus
OTF = np.fft.fft2(array)
MTF = abs(OTF) # MTF is the modulus of the OTF
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
u = np.linspace(0,500,num=501)
v = np.linspace(0,500,num=501)
x, y = np.meshgrid(u, v, sparse=True)
ax.plot_surface(x,y,MTF,cmap=plt.cm.jet)
plt.show()
In [16]:
#1D OTF
nu=np.arange(1,100)*np.pi/100.
#print x
i_P = (2*special.j1(nu)/nu)**2
OTF = np.pi**2 /2
from scipy import integrate
integration_omega = lambda omega: (2*special.j1(omega)/omega)**2 * special.j0(2*np.pi*omega*nu_i)*omega
integrate.quad(integration_omega, 0, np.inf)
plt.plot(nu/np.pi,OTF)
plt.yscale('log')
ax = plt.gca()
plt.xlabel(r'$\nu / \pi$')
plt.ylabel('MTF')
plt.show()
In [17]:
import matplotlib.pyplot as plt
t = np.arange(400)
n = np.zeros((400,), dtype=complex)
n[40:60] = np.exp(1j*np.random.uniform(0, 2*np.pi, (20,)))
s = np.fft.ifft(n)
plt.plot(t, s.real, 'b-', t, s.imag, 'r--')
plt.legend(('real', 'imaginary'))
plt.show()
In [ ]:
In [34]:
nu = np.linspace(0.1,7*np.pi,num=100)
print nu
E = np.ones(100)
PSF = abs(np.fft.ifft(E))
I = np.convolve(E**2,PSF,'same')
print len(I)
print I
plt.plot(nu, I)
plt.show()