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"
              xmlns:xlink="http://www.w3.org/1999/xlink"
              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
spot_scaling, spot_max_radius = 50, 20
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
        spot_radius = min(layers[i,j]*spot_scaling, spot_max_radius)
        print(svg_circle(spot_radius, sx, sy), file=fo)
        if i:
            # The pattern is symmetric about the centre horizontal:
            # duplicate the layers with i > 0
            sy = canvas_height - sy
            print(svg_circle(spot_radius, sx, sy), file=fo)

print(r'</svg>', file=fo)

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

R. E. Franklin, R. G. Gosling, Nature 171, 740 (1953).