In [1]:
import matplotlib.pyplot as plt
import numpy as np
import datetime
import scipy.interpolate
import scipy.optimize
%matplotlib notebook

Scattering from an X-ray beam coming perpendicular to the screen (from your eyes)¶

This is like dropping identical stones in a pond at identical times. The ripples created would look like the images below. Each blue point is a scatterer (stone dropped in the pond).

When scattering occurs from more than one point, constructive and destructive interference of the waves occurs.

The only adjustment in real life is that the ripples emanate in 3D (i.e. are spheres rather than circles).

In [2]:
fig, ax = plt.subplots(1, 3)

a = np.linspace(-3, 3, 601)
b = np.linspace(-3, 3, 601)
center=[0, 0]
A, B = np.meshgrid(a, b)
# Distance from center
dist = np.sqrt((A - center[0]) ** 2 + (B - center[1]) ** 2)*10
phase = np.sin(dist)

ax[0].contourf(A, B, phase, cmap='RdGy', levels=100)  # You can choose a different colormap
ax[0].scatter([0], [0])



x = np.linspace(-2.5, 3.5, 601)
y = np.linspace(-3, 3, 601)
X, Y = np.meshgrid(x, y)
# Distance from center
dist = np.sqrt((X - center[0]) ** 2 + (Y - center[1]) ** 2)*10
phase = np.sin(dist)



center1=[1, 0]
# Distance from center
dist1 = np.sqrt((X - center1[0]) ** 2 + (Y - center1[1]) ** 2)*10
phase1 = np.sin(dist1)

sum1 = phase + phase1
diff = phase-phase1


center2=[0.5, 0.5]
# Distance from center
dist2 = np.sqrt((X - center2[0]) ** 2 + (Y - center2[1]) ** 2)*10
phase2 = np.sin(dist2)

center3=[0.5, -0.5]
# Distance from center
dist3 = np.sqrt((X - center3[0]) ** 2 + (Y - center3[1]) ** 2)*10
phase3 = np.sin(dist3)

sum3 = phase + phase1 + phase2 + phase3
ax[1].contourf(X, Y, sum1, cmap='RdGy', levels=100)  # You can choose a different colormap

ax[2].contourf(X, Y, sum3, cmap='RdGy', levels=100)  # You can choose a different colormap

ax[1].set_xticks([])
ax[1].set_yticks([])
ax[2].set_xticks([])
ax[2].set_yticks([])
ax[1].axis('square')
ax[2].axis('square')
ax[0].set_xticks([])
ax[0].set_yticks([])
ax[0].axis('square')

ax[1].scatter([0, 1], [0, 0], marker='o')
ax[2].scatter([0, 1, 0.5, 0.5], [0, 0, -0.5, 0.5], marker='o')

center=[0, 0]

plt.show()

#plt.savefig('Diffrac images\\Pond.png', dpi=300)
No description has been provided for this image

Now try it with a larger array of scatterers:¶

(This will take 30-60 seconds, use the printouts if you want to track the progress)

In [3]:
start = -10
start2 = -10
mult=20
atomsx = []
atomsy = []

for a in range(start, 11):
    #print(a)         # printouts if you want to track progress
    #print(datetime.datetime.now())      #printouts if you want to track progress
    for b in range(start2, 11):
        atomsx.append(a)
        atomsy.append(b)
        center = [a, b]   # center of the scatterer
        x = np.linspace(-100, 100, 1001)
        y = np.linspace(-100, 100, 1001)
        X, Y = np.meshgrid(x, y)
        # Distance from center
        dist = np.sqrt(((X - center[0])*mult) ** 2 + ((Y - center[1])*mult) ** 2)
        phase = np.cos(dist)
        if a == start and b == start2:
            phasetot = phase
        else:
            phasetot += phase
In [4]:
fig, ax = plt.subplots()




# Plot the colormap
ax.contourf(X, Y, phasetot, cmap='RdGy', levels=100)
ax.scatter(atomsx, atomsy, s=0.1)



ax.set_xticks([])
ax.set_yticks([])


ax.set_aspect('equal', adjustable='box')
ax.set_xlim(-100, 100)
ax.set_ylim(-100, 100)

#plt.savefig('Diffrac images\\Pond 2.png', dpi=300)
No description has been provided for this image
Out[4]:
(-100.0, 100.0)

Rather neatly, constructive interference occurs only in specific directions (this depends on the wavelength of light used).¶

Now consider an X-ray beam that is travelling within the plane of the screen¶

This changes our calculations slightly. Instead of waves being scattered by the scatterers at the same time, as in our pond/ripple example, the wave has travelled different distances to reach different scatterers. As such, the starting phase of the wave emanating from each scatterer is different.

Polar coordinates have been used here in order to simplify the next step in understanding diffraction.

(The cell below takes approximately 120 seconds, use printouts if you want to track progress)

In [5]:
start = -3
start2 = -3
mult=15
atomsr = []
atomstheta = []
beamstart = [10000, 10000]
beamstartrtheta = [np.sqrt(beamstart[0]**2 + beamstart[1]**2), np.arctan(beamstart[1]/beamstart[0])]

for a in range(start, 4):
    #print(a)
    #print(datetime.datetime.now())    #printouts for tracking progress
    for b in range(start2, 4):
        atomr = np.sqrt(a**2 + b**2)
        atomsr.append(atomr)
        if a != 0:
            atomtheta = np.arctan(b/a)
            if a < 0 :
                atomtheta -= np.pi
        elif b > 0:
            atomtheta = np.pi/2
        elif b <= 0:
            atomtheta = -np.pi/2
        atomstheta.append(atomtheta)
        center = [a, b]
        r = np.linspace(0, 100, 10001)
        theta = np.linspace(0, 2*np.pi, 3601)
        R, Theta = np.meshgrid(r, theta)
        X, Y = R*np.cos(Theta), R*np.sin(Theta)
        
        # Distance travelled to scatterer
        dist1 = np.sqrt((a-beamstart[0])**2 + (b-beamstart[1])**2) * mult
        # Distance from scatterer
        dist2 = np.sqrt((X - a)**2 + (Y-b)**2) * mult
        phase = np.cos(dist1 + dist2)
        if a == start and b == start2:
            phasetot = phase
        else:
            phasetot += phase
In [6]:
fig = plt.figure()

# Change the coordinate system from scalar to polar
ax = fig.add_subplot(projection='polar')
ax.yaxis.grid(False)
ax.set_thetagrids([])
ax.set_rgrids([])

ax.contourf(Theta, R, phasetot, cmap='RdGy', levels=50)
ax.scatter(atomstheta, atomsr, s=1)
ax.plot([np.pi/4, beamstartrtheta[1]], [2, beamstartrtheta[0]], linewidth=1)
ax.set_ylim(0, 100)

#plt.savefig('Diffrac images\\Diffrac 1.png', dpi=300)
No description has been provided for this image
Out[6]:
(0.0, 100.0)

The blue square is a series of individual scatterers. The blue line shows the origin of the incident X-ray beam.¶

How does this become a diffraction pattern?¶

It is clear to see the regions of constructive and destructive interference. In an experiment, we would simply detect the photons (when they behave as a particle rather than a wave), which would give us a value at each angle of diffraction. For this simulation however, I have only calculated the wave equations.

I cannot take the displacement of the wave at one point, as even constructively interfered waves produce points at which the displacement is equal to 0 (consider a simple cosine wave, it crosses 0 regularly). We actually need the amplitude of the wavefront at each angle. To do this, I must pick a value of $\theta$ from the centre of the sample and analyse the value of the wavefront against the distance ($r$) for this value of $\theta$.

This looks like the below for the wavefront emanating directly to the right from the scatterers.

In [7]:
points = np.vstack((R.flatten(), Theta.flatten())).T  # Shape (N, 2)
values = phasetot.flatten()  # Shape (N,)

i = len(r)
x_new = points[0:i][:, 0]
y_new = values[0:i]

fig, ax = plt.subplots(2, 1)

ax[0].scatter(x_new, y_new, s=0.1)
ax[1].scatter(x_new, y_new, s=0.3)
ax[1].set_xlim(90, 100)

plt.subplots_adjust(hspace=0.3)

ax[0].set_xlabel(r'$r$')
ax[1].set_xlabel(r'$r$')
ax[1].set_ylim(-4, 4)

#plt.savefig('Diffrac images\\Diffrac 2.png', dpi=300)
No description has been provided for this image
Out[7]:
(-4.0, 4.0)

Close to the scatterers, this appears as a very chaotic interference pattern. However, the distance between sample and detector is always much greater than the distance across the sample in a diffraction experiment.¶

Close to the sample, the amplitude of the wavefront is unclear using this approach. In reality, the amplitude at these points is actually time-dependent (the phase difference is not constant). Phasors would be required to determine the amplitude, but this will not be covered.

In the limit of infinite distance between the sample and detector, the wavefront breaks down into simply one wave with a certain amplitude, frequency and phase shift. This is because at infinite distance, the direction of each scattered wave is almost identical. This means we are simply summing waves with different phases that appeared to originate from the same point. These waves all have the same frequency (the frequency of the incident beam) and wavevector (direction), and as such only the phase shift is relevant in determining interference patterns.

Summing two cosine waves with different phases simply produces a third cosine wave with a new phase and amplitude, as demonstrated below:

Feel free to play around with the values of $a$ and $b$. If the fitted wave ever doesn't fit the summed wave perfectly, this is due to a failing of the fitting procedure.

In [8]:
x = np.linspace(-np.pi, 3*np.pi, 3601)

a = np.pi*3/4
b = -np.pi*2
y1 = np.cos(x + a)
y2 = np.cos(x + b)
ysum = y1 + y2

def cosfunc(x, shift, amplitude):
    return amplitude*np.cos(x+shift)


coefs, cov = scipy.optimize.curve_fit(cosfunc, x, ysum, p0=[3, 1])
y_fit = cosfunc(x, coefs[0], coefs[1])

fig, ax = plt.subplots()

ax.plot(x, y1, label='Wave 1')
ax.plot(x, y2, label='Wave 2')
ax.plot(x, ysum, label='Wave Sum', color='red', alpha=0.5)
ax.plot(x, y_fit, linestyle=':', label='Wave fitted', color='black')
ax.set_xlim(0, 2*np.pi)
ax.set_xticks(list(np.arange(0, 2*np.pi+0.1, np.pi/6)))
ax.set_xticklabels(['0', r'$\frac{\pi}{6}$', r'$\frac{\pi}{3}$', r'$\frac{\pi}{2}$', r'$\frac{2\pi}{3}$', r'$\frac{5\pi}{6}$', r'$\pi$', r'$\frac{7\pi}{6}$', r'$\frac{4\pi}{3}$', r'$\frac{3\pi}{2}$', r'$\frac{5\pi}{3}$', r'$\frac{11\pi}{6}$', r'$2\pi$'])
ax.set_yticks([])
ax.legend()

#plt.savefig('Diffrac images\\Diffrac 3.png', dpi=300)
No description has been provided for this image
Out[8]:
<matplotlib.legend.Legend at 0x1fc70768ad0>

Returning to the problem of 'how do I calculate a photon count from just my wave equations?'.¶

If we can fit the wavefront in each direction to a wave equation, the amplitude of the wave in the limit of infinite distance is the interference pattern!

So, to summarise, for each angle, we must produce a dataset of wave values against $r$. We must then fit this dataset to a cosine wave at large distances only, then take the amplitude of the wavefront as our expected photon count. Repeat this process for each value of $\theta$, and we have a dataset of expected photon counts against $\theta$.

Below this is performed for just the plots shown above. Plot 1 shows the full range of $r$. Plots 2 and 3 show this zoomed in on low $r$ and high $r$ values.

In [9]:
points = np.vstack((R.flatten(), Theta.flatten())).T  # Shape (N, 2)
values = phasetot.flatten()  # Shape (N,)

i = len(r)
x_new = points[0:i][:, 0]
y_new = values[0:i]

fig, ax = plt.subplots(3, 1, figsize=(5, 8))

ax[0].scatter(x_new, y_new, s=0.1)
ax[1].scatter(x_new, y_new, s=0.3)
ax[1].set_xlim(0, 30)
ax[2].scatter(x_new, y_new, s=0.3)
ax[2].set_xlim(90, 100)

plt.subplots_adjust(hspace=0.3)

ax[0].set_xlabel(r'$r$')
ax[1].set_xlabel(r'$r$')
ax[2].set_xlabel(r'$r$')
ax[2].set_ylim(-4, 4)

# Fitting just the large distance limit

def cosfunc(x, freq, shift, amplitude):
    return amplitude*np.cos(freq*x+shift)

x_new = points[i-int(len(r)/10):i][:, 0]
y_new = values[i-int(len(r)/10):i]
coefs, cov = scipy.optimize.curve_fit(cosfunc, x_new, y_new, p0=[mult, 3, 30])
x_fit = np.linspace(x_new[0], x_new[-1], 1001)
y_fit = cosfunc(x_fit, coefs[0], coefs[1], coefs[2])
ax[2].plot(x_fit, y_fit, linewidth=0.5, color='C1')
ax[0].set_xlim(0, 100)

#plt.savefig('Diffrac images\\Diffrac 4.png', dpi=300)
No description has been provided for this image
Out[9]:
(0.0, 100.0)

Now repeat this process for all values of $\theta$¶

I have set this up for datapoints every 0.1°.

We can plot these amplitudes on the outside of our scattering plot as is done below.

In [10]:
# Flatten the X, Y, and dist arrays to make them compatible with griddata
points = np.vstack((R.flatten(), Theta.flatten())).T  # Shape (N, 2)
values = phasetot.flatten()  # Shape (N,)


def cosfunc(x, freq, shift, amplitude):
    return amplitude*np.cos(freq*x+shift)

amplitude_list = []
theta_list = []
counter = 0

for i in range(len(r), len(points)+1, len(r)):    
    counter+=1
    #if counter % 50 ==0:
        #print(counter)
        #print(datetime.datetime.now())           # printouts if you want to track progress
    theta_list.append((points[i-len(r)][1]))
    x_new = points[i-int(len(r)/10):i][:, 0]
    y_new = values[i-int(len(r)/10):i]
    coefs, cov = scipy.optimize.curve_fit(cosfunc, x_new, y_new, p0=[mult, 3, 30])
    x_fit = np.linspace(x_new[0], x_new[-1], 1001)
    y_fit = cosfunc(x_fit, coefs[0], coefs[1], coefs[2])
    amplitude_list.append(abs(coefs[2]))
In [11]:
fig = plt.figure(figsize=(8,8))

# Change the coordinate system from scalar to polar
ax = fig.add_subplot(projection='polar')
ax.yaxis.grid(False)
ax.set_thetagrids([])
ax.set_rgrids([])

ax.contourf(Theta, R, phasetot, cmap='RdGy', levels=50)
ax.scatter(atomstheta, atomsr, s=1)
ax.plot([np.pi/4, beamstartrtheta[1]], [2, 100], linewidth=1)
ax.set_ylim(0, 170)

fakelinetheta = np.linspace(0, 360, 3601)
fakeliner = np.linspace(100, 100, 3601)
ax.plot(fakelinetheta, fakeliner, linewidth=0.5, color='black')

ax.plot(theta_list, [x+100 for x in amplitude_list], linewidth=0.5, color='black')

ax.annotate(r'2$\theta = 0\degree$', xy=(-135*np.pi/180, 145.7), xytext=(-125*np.pi/180, 160), arrowprops=dict(facecolor='black', width=0.05, headwidth=5, headlength=5))
ax.annotate(r'2$\theta = 90\degree$', xy=(-45*np.pi/180, 106.5), xytext=(-65*np.pi/180, 140), arrowprops=dict(facecolor='black', width=0.05, headwidth=5, headlength=5))
ax.annotate(r'2$\theta = 90\degree$', xy=(-225*np.pi/180, 106.5), xytext=(-225*np.pi/180, 160), arrowprops=dict(facecolor='black', width=0.05, headwidth=5, headlength=5))
ax.annotate(r'2$\theta = 180\degree$', xy=(45*np.pi/180, 101), xytext=(70*np.pi/180, 120), arrowprops=dict(facecolor='black', width=0.05, headwidth=5, headlength=5))

#plt.savefig('Diffrac images\\Diffrac 5.png', dpi=300)
No description has been provided for this image
Out[11]:
Text(1.2217304763960306, 120, '2$\\theta = 180\\degree$')

In the above plot, the amplitude of the wavefront is plotted on the outside of the circle.¶

From quantum mechanics, the probability density is proportional to the square of the amplitude. This adjustment is done in the next cell.

If we unzip the plot, we get the following which is the diffraction pattern you are likely more familiar with.

There is an adjustment to $\theta$ in the below cell, as $\theta$ in polar matplotlib plots is not the same as $2\theta$ in the diffraction sense. Only run this cell once.

In [12]:
theta_list = [(x*180/np.pi)+135 for x in theta_list]
amplitude_list = [x**2 for x in amplitude_list]
In [14]:
for i in range(len(theta_list)):
    if theta_list[i] >= 360:
        theta_list[i] -= 360

xy = zip(theta_list, amplitude_list)

new_theta = [x[0] for x in sorted(xy)[0::2]]
xy = zip(theta_list, amplitude_list)
new_amplitude = [x[1] for x in sorted(xy)[0::2]]

fig, ax = plt.subplots(figsize=(7, 4))

ax.plot(new_theta, new_amplitude, linewidth=0.7)
ax.set_xlabel(r"2$\theta$ / $\degree$")
ax.set_ylabel('Expected Counts')
ax.set_xticks(list(range(0, 390, 15)))

ax.set_xlim(0, 180)
ax.set_ylim(-100, 2500)
ax.set_yticks([])

#plt.savefig('Diffrac images\\Diffrac 6.png', dpi=300)
No description has been provided for this image
In [ ]: