# DNA diffraction pattern

In an important paper in 1953 [1], Rosalind Franklin published the X-ray diffraction pattern of DNA from calf thymus which displays a characteristic 'X' shape of diffraction spots indicative of a helical structure.

The diffraction pattern of a uniform, continuous helix consists of a series of "layer lines" of spacing $1/p$ in reciprocal space where $p$ is the helix pitch (the height of one complete turn of the helix, measured parallel to its axis). The intensity distribution along the $n$th layer line is proportional to square of the $n$th Bessel function, $J_n(2\pi r R)$, where $r$ is the radius of the helix and $R$ the radial co-ordinate in reciprocal space.

Consider the diffraction pattern of a helix with $p = 34$ Å and $r=10$ Å. The code listing below produces an SVG image of the diffraction pattern of a helix.

import numpy as np
from scipy.special import jn
import pylab

# Vertical range of the diffraction pattern: plot nlayer line layers above and
# below the centre horizontal
nlayers = 5
ymin, ymax = -nlayers, nlayers

# Horizontal range of the diffraction pattern, x = 2pi.r.R
xmin, xmax = -10, 10
npts = 4000
x = np.linspace(xmin, xmax, npts)

# Diffraction pattern along each line layer: |Jn(x)|^2
# for n = 0, 1, ..., nlayers-1
layers = np.array([jn(i, x)**2 for i in range(nlayers)])

# Obtain the indexes of the maxima in each layer
maxi = [(np.diff(np.sign(np.diff(layers[i,:]))) < 0).nonzero()[0] + 1
for i in range(nlayers)]

# Create the SVG image, using circles of different radii for diffraction spots
svg_name='eg8-dna-diffraction.svg'
canvas_width = canvas_height = 500
fo = open(svg_name, 'w')
print("""<?xml version="1.0" encoding="utf-8"?>
<svg xmlns="http://www.w3.org/2000/svg"
width="{}" height="{}" style="background: {}">""".format(
canvas_width, canvas_height, '#ffffff'), file=fo)

def svg_circle(r, cx, cy):
""" Return the SVG mark up for a circle of radius r centred at (cx,cy). """
return r'<circle r="{}" cx="{}" cy="{}"/>'.format(r, cx, cy)

# For each spot in each layer, draw a circle on the canvas. The circle radius
# is the scaled value of the diffraction intensity maximum, with a ceiling
# value of spot_max_radius because the centre spots are very intense
for i in range(nlayers):
for j in maxi[i]:
sx = (x[j] - xmin)/(xmax-xmin) * canvas_width
sy = (i - ymin)/(ymax - ymin) * canvas_height
if i:
# The pattern is symmetric about the centre horizontal:
# duplicate the layers with i > 0
sy = canvas_height - sy

• The two-dimensional array, layers holds the diffraction intensity in each line layer, calculated as the square of a Bessel function.
• For plotting the pattern, we need to find the indexes of the maxima in the layers array: this line of code finds these maxima by determining where the differences between neighbouring items go from positive to negative.
• The (x,y) coordinates in the reciprocal space of the diffraction pattern are mapped onto the canvas coordinates, (sx,sy).