Remove pattern_sim and partial_sim

Use of these programs has been following this pattern for several years:

1. Neglect
2. Once yearly attempt by someone to use either tool
3. Discovery that it's totally broken
4. Bug report and fast bug fix
5. Go to 1.

For more discussion, see the issue referenced below.

Closes: https://gitlab.desy.de/thomas.white/crystfel/-/issues/81

11 files changed, 0 insertions, 4551 deletions

diff --git a/src/cl-utils.c b/src/cl-utils.c
deleted file mode 100644
index 1b6667c4..00000000
--- a/src/cl-utils.c
+++ /dev/null
@@ -1,354 +0,0 @@
-/*
- * cl-utils.c
- *
- * OpenCL utility functions
- *
- * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
- *                       a research centre of the Helmholtz Association.
- *
- * Authors:
- *   2010-2020 Thomas White <taw@physics.org>
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-
-#define CL_TARGET_OPENCL_VERSION 220
-#ifdef HAVE_CL_CL_H
-#include <CL/cl.h>
-#else
-#include <cl.h>
-#endif
-
-#include "utils.h"
-
-
-/* Return 1 if a GPU device is present, 0 if not, 2 on error. */
-int have_gpu_device()
-{
-	cl_uint nplat;
-	cl_platform_id platforms[8];
-	cl_context_properties prop[3];
-	cl_int err;
-	int i;
-
-	err = clGetPlatformIDs(8, platforms, &nplat);
-	if ( err != CL_SUCCESS ) return 2;
-	if ( nplat == 0 ) return 0;
-
-	/* Find a GPU platform in the list */
-	for ( i=0; i<nplat; i++ ) {
-
-		prop[0] = CL_CONTEXT_PLATFORM;
-		prop[1] = (cl_context_properties)platforms[i];
-		prop[2] = 0;
-
-		clCreateContextFromType(prop, CL_DEVICE_TYPE_GPU,
-		                        NULL, NULL, &err);
-
-		if ( err != CL_SUCCESS ) {
-			if ( err != CL_DEVICE_NOT_FOUND ) return 2;
-		} else {
-			return 1;
-		}
-	}
-
-	return 0;
-}
-
-
-const char *clError(cl_int err)
-{
-	switch ( err ) {
-
-		case CL_SUCCESS :
-		return "no error";
-
-		case CL_DEVICE_NOT_AVAILABLE :
-		return "device not available";
-
-		case CL_DEVICE_NOT_FOUND :
-		return "device not found";
-
-		case CL_INVALID_DEVICE_TYPE :
-		return "invalid device type";
-
-		case CL_INVALID_PLATFORM :
-		return "invalid platform";
-
-		case CL_INVALID_KERNEL :
-		return "invalid kernel";
-
-		case CL_INVALID_ARG_INDEX :
-		return "invalid argument index";
-
-		case CL_INVALID_ARG_VALUE :
-		return "invalid argument value";
-
-		case CL_INVALID_MEM_OBJECT :
-		return "invalid memory object";
-
-		case CL_INVALID_SAMPLER :
-		return "invalid sampler";
-
-		case CL_INVALID_ARG_SIZE :
-		return "invalid argument size";
-
-		case CL_INVALID_COMMAND_QUEUE :
-		return "invalid command queue";
-
-		case CL_INVALID_CONTEXT :
-		return "invalid context";
-
-		case CL_INVALID_VALUE :
-		return "invalid value";
-
-		case CL_INVALID_EVENT_WAIT_LIST :
-		return "invalid wait list";
-
-		case CL_MAP_FAILURE :
-		return "map failure";
-
-		case CL_MEM_OBJECT_ALLOCATION_FAILURE :
-		return "object allocation failure";
-
-		case CL_OUT_OF_HOST_MEMORY :
-		return "out of host memory";
-
-		case CL_OUT_OF_RESOURCES :
-		return "out of resources";
-
-		case CL_INVALID_KERNEL_NAME :
-		return "invalid kernel name";
-
-		case CL_INVALID_KERNEL_ARGS :
-		return "invalid kernel arguments";
-
-		case CL_INVALID_WORK_GROUP_SIZE :
-		return "invalid work group size";
-
-		case CL_IMAGE_FORMAT_NOT_SUPPORTED :
-		return "image format not supported";
-
-		case CL_INVALID_WORK_DIMENSION :
-		return "invalid work dimension";
-
-		default :
-		return "unknown error";
-	}
-}
-
-
-static char *get_device_string(cl_device_id dev, cl_device_info info)
-{
-	int r;
-	size_t size;
-	char *val;
-
-	r = clGetDeviceInfo(dev, info, 0, NULL, &size);
-	if ( r != CL_SUCCESS ) {
-		ERROR("Couldn't get device vendor size: %s\n",
-		      clError(r));
-		return NULL;
-	}
-	val = malloc(size);
-	r = clGetDeviceInfo(dev, info, size, val, NULL);
-	if ( r != CL_SUCCESS ) {
-		ERROR("Couldn't get dev vendor: %s\n", clError(r));
-		return NULL;
-	}
-
-	return val;
-}
-
-
-cl_device_id get_cl_dev(cl_context ctx, int n)
-{
-	cl_device_id *dev;
-	cl_int r;
-	size_t size;
-	int num_devs;
-
-	/* Get the required size of the array */
-	r = clGetContextInfo(ctx, CL_CONTEXT_DEVICES, 0, NULL, &size);
-	if ( r != CL_SUCCESS ) {
-		ERROR("Couldn't get array size for devices: %s\n", clError(r));
-		return 0;
-	}
-
-	dev = malloc(size);
-	r = clGetContextInfo(ctx, CL_CONTEXT_DEVICES, size, dev, NULL);
-	if ( r != CL_SUCCESS ) {
-		ERROR("Couldn't get device: %s\n", clError(r));
-		return 0;
-	}
-	num_devs = size/sizeof(cl_device_id);
-
-	if ( n >= num_devs ) {
-		ERROR("Device ID out of range\n");
-		return 0;
-	}
-
-	if ( n < 0 ) {
-
-		int i;
-
-		STATUS("Available devices:\n");
-		for ( i=0; i<num_devs; i++ ) {
-
-			char *vendor;
-			char *name;
-
-			vendor = get_device_string(dev[i], CL_DEVICE_VENDOR);
-			name = get_device_string(dev[i], CL_DEVICE_NAME);
-
-			STATUS("Device %i: %s %s\n", i, vendor, name);
-
-		}
-		n = 0;
-
-		STATUS("Using device 0.  Use --gpu-dev to choose another.\n");
-
-	} else {
-
-		char *vendor;
-		char *name;
-
-		vendor = get_device_string(dev[n], CL_DEVICE_VENDOR);
-		name = get_device_string(dev[n], CL_DEVICE_NAME);
-
-		STATUS("Using device %i: %s %s\n", n, vendor, name);
-
-	}
-
-	return dev[n];
-}
-
-
-static void show_build_log(cl_program prog, cl_device_id dev)
-{
-	cl_int r;
-	char log[4096];
-	size_t s;
-
-	r = clGetProgramBuildInfo(prog, dev, CL_PROGRAM_BUILD_LOG, 4096, log,
-	                          &s);
-
-	STATUS("Status: %i\n", r);
-	STATUS("%s\n", log);
-}
-
-
-cl_program load_program_from_string(const char *source_in, size_t len,
-                                    cl_context ctx, cl_device_id dev,
-                                    cl_int *err, const char *extra_cflags,
-                                    const char *insert_stuff)
-{
-	cl_program prog;
-	cl_int r;
-	char cflags[1024] = "";
-	char *insert_pos;
-	char *source;
-
-	/* Copy the code because we need to zero-terminate it */
-	source = malloc(len+1);
-	if ( source == NULL ) return 0;
-	memcpy(source, source_in, len);
-	source[len] = '\0';
-
-	if ( insert_stuff != NULL ) {
-		insert_pos = strstr(source, "INSERT_HERE");
-
-		if ( insert_pos != NULL ) {
-
-			char *source2;
-			size_t il;
-
-			source2 = malloc(strlen(source)+strlen(insert_stuff)+1);
-			if ( source2 == NULL ) {
-				free(source);
-				return 0;
-			}
-
-			il = insert_pos - source;
-			memcpy(source2, source, il);
-			memcpy(source2+il, insert_stuff, strlen(insert_stuff)+1);
-			memcpy(source2+il+strlen(insert_stuff),
-			       source+il+11, strlen(source+il+11)+1);
-			free(source);
-			source = source2;
-
-		}
-	}
-
-	prog = clCreateProgramWithSource(ctx, 1, (const char **)&source,
-	                                 NULL, err);
-	if ( *err != CL_SUCCESS ) {
-		ERROR("Couldn't load source\n");
-		return 0;
-	}
-
-	cflags[0] = '\0';
-	strncat(cflags, "-cl-no-signed-zeros ", 1023-strlen(cflags));
-	strncat(cflags, extra_cflags, 1023-strlen(cflags));
-
-	r = clBuildProgram(prog, 0, NULL, cflags, NULL, NULL);
-	if ( r != CL_SUCCESS ) {
-		ERROR("Couldn't build program\n");
-		show_build_log(prog, dev);
-		*err = r;
-		return 0;
-	}
-
-	free(source);
-	*err = CL_SUCCESS;
-	return prog;
-}
-
-
-cl_program load_program(const char *filename, cl_context ctx,
-                        cl_device_id dev, cl_int *err, const char *extra_cflags,
-                        const char *insert_stuff)
-{
-	FILE *fh;
-	char *source;
-	size_t len;
-
-	fh = fopen(filename, "r");
-	if ( fh == NULL ) {
-		ERROR("Couldn't open '%s'\n", filename);
-		*err = CL_INVALID_PROGRAM;
-		return 0;
-	}
-	source = malloc(16384);
-	if ( source == NULL ) {
-		fclose(fh);
-		return 0;
-	}
-	len = fread(source, 1, 16383, fh);
-	fclose(fh);
-
-	return load_program_from_string(source, len, ctx, dev,err, extra_cflags,
-	                                insert_stuff);
-}
diff --git a/src/cl-utils.h b/src/cl-utils.h
deleted file mode 100644
index d72ee6d5..00000000
--- a/src/cl-utils.h
+++ /dev/null
@@ -1,55 +0,0 @@
-/*
- * cl-utils.h
- *
- * OpenCL utility functions
- *
- * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
- *                       a research centre of the Helmholtz Association.
- *
- * Authors:
- *   2010-2020 Thomas White <taw@physics.org>
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-#ifndef CLUTILS_H
-#define CLUTILS_H
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#define CL_TARGET_OPENCL_VERSION 220
-#ifdef HAVE_CL_CL_H
-#include <CL/cl.h>
-#else
-#include <cl.h>
-#endif
-
-extern int have_gpu_device(void);
-extern const char *clError(cl_int err);
-extern cl_device_id get_cl_dev(cl_context ctx, int n);
-extern cl_program load_program_from_string(const char *source_in, size_t len,
-                                           cl_context ctx, cl_device_id dev,
-                                           cl_int *err, const char *extra_cflags,
-                                           const char *insert_stuff);
-extern cl_program load_program(const char *filename, cl_context ctx,
-                               cl_device_id dev, cl_int *err,
-                               const char *extra_cflags, const char *insert_stuff);
-
-
-#endif	/* CLUTILS_H */
diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c
deleted file mode 100644
index 97badc36..00000000
--- a/src/diffraction-gpu.c
+++ /dev/null
@@ -1,596 +0,0 @@
-/*
- * diffraction-gpu.c
- *
- * Calculate diffraction patterns by Fourier methods (GPU version)
- *
- * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
- *                       a research centre of the Helmholtz Association.
- *
- * Authors:
- *   2009-2020 Thomas White <taw@physics.org>
- *   2013      Alexandra Tolstikova
- *   2013-2014 Chun Hong Yoon <chun.hong.yoon@desy.de>
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <stdlib.h>
-#include <math.h>
-#include <stdio.h>
-#include <string.h>
-#include <complex.h>
-
-#define CL_TARGET_OPENCL_VERSION 220
-#ifdef HAVE_CL_CL_H
-#include <CL/cl.h>
-#else
-#include <cl.h>
-#endif
-
-#include "image.h"
-#include "utils.h"
-#include "cell.h"
-#include "diffraction.h"
-#include "cl-utils.h"
-#include "pattern_sim.h"
-
-#include "diffraction.cl.h"
-
-#define SINC_LUT_ELEMENTS (4096)
-
-
-struct gpu_context
-{
-	cl_context ctx;
-	cl_command_queue cq;
-	cl_program prog;
-	cl_kernel kern;
-	cl_mem intensities;
-	cl_mem flags;
-
-	/* Array of sinc LUTs */
-	cl_mem *sinc_luts;
-	cl_float **sinc_lut_ptrs;
-	int max_sinc_lut;  /* Number of LUTs, i.e. one greater than the maximum
-	                    * index.  This equals the highest allowable "n". */
-};
-
-
-static void check_sinc_lut(struct gpu_context *gctx, int n,
-                           int no_fringes, int flat)
-{
-	cl_int err;
-	cl_image_format fmt;
-	int i;
-
-	if ( n > gctx->max_sinc_lut ) {
-
-		gctx->sinc_luts = realloc(gctx->sinc_luts,
-		                          n*sizeof(*gctx->sinc_luts));
-		gctx->sinc_lut_ptrs = realloc(gctx->sinc_lut_ptrs,
-		                              n*sizeof(*gctx->sinc_lut_ptrs));
-
-		for ( i=gctx->max_sinc_lut; i<n; i++ ) {
-			gctx->sinc_lut_ptrs[i] = NULL;
-		}
-
-		gctx->max_sinc_lut = n;
-	}
-
-	if ( gctx->sinc_lut_ptrs[n-1] != NULL ) return;
-
-	/* Create a new sinc LUT */
-	gctx->sinc_lut_ptrs[n-1] = malloc(SINC_LUT_ELEMENTS*sizeof(cl_float));
-	gctx->sinc_lut_ptrs[n-1][0] = n;
-	if ( n == 1 ) {
-		for ( i=1; i<SINC_LUT_ELEMENTS; i++ ) {
-			gctx->sinc_lut_ptrs[n-1][i] = 1.0;
-		}
-	} else {
-		for ( i=1; i<SINC_LUT_ELEMENTS; i++ ) {
-			double x, val;
-			x = (double)i/SINC_LUT_ELEMENTS;
-			if ( (flat || no_fringes) && (x > 1.0/n) && (1.0-x > 1.0/n) ) {
-				val = 0.0;
-			} else if ( flat ) {
-				val = n;
-			} else {
-				val = fabs(sin(M_PI*n*x)/sin(M_PI*x));
-			}
-			gctx->sinc_lut_ptrs[n-1][i] = val;
-		}
-	}
-
-	fmt.image_channel_order = CL_INTENSITY;
-	fmt.image_channel_data_type = CL_FLOAT;
-
-	gctx->sinc_luts[n-1] = clCreateImage2D(gctx->ctx,
-	                                CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
-	                                &fmt, SINC_LUT_ELEMENTS, 1, 0,
-	                                gctx->sinc_lut_ptrs[n-1], &err);
-	if ( err != CL_SUCCESS ) {
-		ERROR("Couldn't create LUT for %i\n", n);
-		return;
-	}
-}
-
-
-static int set_arg_float(struct gpu_context *gctx, int idx, float val)
-{
-	cl_int err;
-	err = clSetKernelArg(gctx->kern, idx, sizeof(cl_float), &val);
-	if ( err != CL_SUCCESS ) {
-		ERROR("Couldn't set kernel argument %i: %s\n",
-		      idx, clError(err));
-		return 1;
-	}
-
-	return 0;
-}
-
-
-static int set_arg_int(struct gpu_context *gctx, int idx, int val)
-{
-	cl_int err;
-
-	err = clSetKernelArg(gctx->kern, idx, sizeof(cl_int), &val);
-	if ( err != CL_SUCCESS ) {
-		ERROR("Couldn't set kernel argument %i: %s\n",
-		      idx, clError(err));
-		return 1;
-	}
-
-	return 0;
-}
-
-
-static int set_arg_mem(struct gpu_context *gctx, int idx, cl_mem val)
-{
-	cl_int err;
-
-	err = clSetKernelArg(gctx->kern, idx, sizeof(cl_mem), &val);
-	if ( err != CL_SUCCESS ) {
-		ERROR("Couldn't set kernel argument %i: %s\n",
-		      idx, clError(err));
-		return 1;
-	}
-
-	return 0;
-}
-
-
-static int do_panels(struct gpu_context *gctx, struct image *image,
-                      double k, double weight,
-                      int *n_inf, int *n_neg, int *n_nan)
-{
-	int i;
-	const int sampling = 4;  /* This, squared, number of samples / pixel */
-
-	if ( set_arg_float(gctx, 1, k) ) return 1;
-	if ( set_arg_float(gctx, 2, weight) ) return 1;
-
-	/* Iterate over panels */
-	for ( i=0; i<image->detgeom->n_panels; i++ ) {
-
-		size_t dims[2];
-		size_t ldims[2];
-		struct detgeom_panel *p;
-		cl_mem diff;
-		size_t diff_size;
-		float *diff_ptr;
-		int fs, ss;
-		cl_int err;
-
-		p = &image->detgeom->panels[i];
-
-		/* Buffer for the results of this panel */
-		diff_size = p->w * p->h * sizeof(cl_float);
-		diff = clCreateBuffer(gctx->ctx, CL_MEM_WRITE_ONLY,
-	                              diff_size, NULL, &err);
-		if ( err != CL_SUCCESS ) {
-			ERROR("Couldn't allocate diffraction memory\n");
-			return 1;
-		}
-
-		if ( set_arg_mem(gctx, 0, diff) ) return 1;
-
-		if ( set_arg_int(gctx, 3, p->w) ) return 1;
-		if ( set_arg_float(gctx, 4, p->cnx) ) return 1;
-		if ( set_arg_float(gctx, 5, p->cny) ) return 1;
-		if ( set_arg_float(gctx, 6, p->fsx) ) return 1;
-		if ( set_arg_float(gctx, 7, p->fsy) ) return 1;
-		if ( set_arg_float(gctx, 8, p->fsz) ) return 1;
-		if ( set_arg_float(gctx, 9, p->ssx) ) return 1;
-		if ( set_arg_float(gctx, 10, p->ssy) ) return 1;
-		if ( set_arg_float(gctx, 11, p->ssz) ) return 1;
-		if ( set_arg_float(gctx, 12, 1.0/p->pixel_pitch) ) return 1;
-		if ( set_arg_float(gctx, 13, p->cnz*p->pixel_pitch) ) return 1;
-
-		dims[0] = p->w * sampling;
-		dims[1] = p->h * sampling;
-
-		ldims[0] = sampling;
-		ldims[1] = sampling;
-
-		err = clSetKernelArg(gctx->kern, 20,
-		                     sampling*sampling*sizeof(cl_float), NULL);
-		if ( err != CL_SUCCESS ) {
-			ERROR("Couldn't set local memory: %s\n", clError(err));
-			return 1;
-		}
-
-		err = clEnqueueNDRangeKernel(gctx->cq, gctx->kern, 2, NULL,
-		                             dims, ldims, 0, NULL, NULL);
-		if ( err != CL_SUCCESS ) {
-			ERROR("Couldn't enqueue diffraction kernel: %s\n",
-			      clError(err));
-			return 1;
-		}
-
-		clFinish(gctx->cq);
-
-		diff_ptr = clEnqueueMapBuffer(gctx->cq, diff, CL_TRUE,
-		                              CL_MAP_READ, 0, diff_size,
-		                              0, NULL, NULL, &err);
-		if ( err != CL_SUCCESS ) {
-			ERROR("Couldn't map diffraction buffer: %s\n",
-			      clError(err));
-			return 1;
-		}
-
-		for ( ss=0; ss<p->h; ss++ ) {
-		for ( fs=0; fs<p->w; fs++ ) {
-
-			float val;
-
-			val = diff_ptr[fs + p->w*ss];
-			if ( isinf(val) ) (*n_inf)++;
-			if ( val < 0.0 ) (*n_neg)++;
-			if ( isnan(val) ) (*n_nan)++;
-
-			image->dp[i][fs + p->w*ss] += val;
-
-		}
-		}
-
-		clEnqueueUnmapMemObject(gctx->cq, diff, diff_ptr,
-		                        0, NULL, NULL);
-
-		clReleaseMemObject(diff);
-
-	}
-
-	return 0;
-}
-
-
-int get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
-                        int na, int nb, int nc, UnitCell *ucell,
-                        int no_fringes, int flat, int n_samples)
-{
-	double ax, ay, az;
-	double bx, by, bz;
-	double cx, cy, cz;
-	cl_float16 cell;
-	cl_int err;
-	int n_inf = 0;
-	int n_neg = 0;
-	int n_nan = 0;
-	int i;
-	double kmin, kmax, step, norm;
-
-	if ( gctx == NULL ) {
-		ERROR("GPU setup failed.\n");
-		return 1;
-	}
-
-	/* Ensure all required LUTs are available */
-	check_sinc_lut(gctx, na, no_fringes, flat);
-	check_sinc_lut(gctx, nb, no_fringes, flat);
-	check_sinc_lut(gctx, nc, no_fringes, flat);
-
-	/* Unit cell */
-	cell_get_cartesian(ucell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
-	cell.s[0] = ax;  cell.s[1] = ay;  cell.s[2] = az;
-	cell.s[3] = bx;  cell.s[4] = by;  cell.s[5] = bz;
-	cell.s[6] = cx;  cell.s[7] = cy;  cell.s[8] = cz;
-
-	err = clSetKernelArg(gctx->kern, 14, sizeof(cl_float16), &cell);
-	if ( err != CL_SUCCESS ) {
-		ERROR("Couldn't set unit cell: %s\n", clError(err));
-		return 1;
-	}
-
-	if ( set_arg_mem(gctx, 15, gctx->intensities) ) return 1;
-	if ( set_arg_mem(gctx, 16, gctx->flags) ) return 1;
-	if ( set_arg_mem(gctx, 17, gctx->sinc_luts[na-1]) ) return 1;
-	if ( set_arg_mem(gctx, 18, gctx->sinc_luts[nb-1]) ) return 1;
-	if ( set_arg_mem(gctx, 19, gctx->sinc_luts[nc-1]) ) return 1;
-
-	/* Allocate memory for the result */
-	image->dp = malloc(image->detgeom->n_panels * sizeof(float *));
-	if ( image->dp == NULL ) {
-		ERROR("Couldn't allocate memory for result.\n");
-		return 1;
-	}
-	for ( i=0; i<image->detgeom->n_panels; i++ ) {
-		struct detgeom_panel *p = &image->detgeom->panels[i];
-		image->dp[i] = calloc(p->w * p->h, sizeof(float));
-		if ( image->dp[i] == NULL ) {
-			ERROR("Couldn't allocate memory for panel %i\n", i);
-			return 1;
-		}
-	}
-
-	spectrum_get_range(image->spectrum, &kmin, &kmax);
-	step = (kmax-kmin)/(n_samples+1);
-
-	/* Determine normalisation factor such that weights add up to 1 after
-	 * sampling (bins must have constant width) */
-	norm = 0.0;
-	for ( i=1; i<=n_samples; i++ ) {
-		double k = kmin + i*step;
-		norm += spectrum_get_density_at_k(image->spectrum, k);
-	}
-	for ( i=1; i<=n_samples; i++ ) {
-
-		double k = kmin + i*step;
-		double prob;
-
-		/* Probability = p.d.f. times step width */
-		prob = spectrum_get_density_at_k(image->spectrum, k)/norm;
-
-		STATUS("Wavelength: %e m, weight = %.5f\n", 1.0/k, prob);
-
-		err = do_panels(gctx, image, k, prob, &n_inf, &n_neg, &n_nan);
-		if ( err ) return 1;
-
-	}
-
-	if ( n_neg + n_inf + n_nan ) {
-		ERROR("WARNING: The GPU calculation produced %i negative"
-		      " values, %i infinities and %i NaNs.\n",
-		      n_neg, n_inf, n_nan);
-	}
-
-	return 0;
-}
-
-
-/* Setup the OpenCL stuff, create buffers, load the structure factor table */
-struct gpu_context *setup_gpu(int no_sfac,
-                              const double *intensities, unsigned char *flags,
-                              const char *sym, int dev_num)
-{
-	struct gpu_context *gctx;
-	cl_uint nplat;
-	cl_platform_id platforms[8];
-	cl_context_properties prop[3];
-	cl_int err;
-	cl_device_id dev;
-	size_t intensities_size;
-	float *intensities_ptr;
-	size_t flags_size;
-	float *flags_ptr;
-	size_t maxwgsize;
-	int iplat;
-	int have_ctx = 0;
-	char cflags[512] = "";
-	char *insert_stuff = NULL;
-
-	STATUS("Setting up GPU...\n");
-
-	err = clGetPlatformIDs(8, platforms, &nplat);
-	if ( err != CL_SUCCESS ) {
-		ERROR("Couldn't get platform IDs: %i\n", err);
-		return NULL;
-	}
-	if ( nplat == 0 ) {
-		ERROR("Couldn't find at least one platform!\n");
-		return NULL;
-	}
-
-	/* Find a GPU platform in the list */
-	for ( iplat=0; iplat<nplat; iplat++ ) {
-
-		prop[0] = CL_CONTEXT_PLATFORM;
-		prop[1] = (cl_context_properties)platforms[iplat];
-		prop[2] = 0;
-
-		gctx = malloc(sizeof(*gctx));
-		gctx->ctx = clCreateContextFromType(prop, CL_DEVICE_TYPE_GPU,
-		                                    NULL, NULL, &err);
-
-		if ( err != CL_SUCCESS ) {
-			if ( err == CL_DEVICE_NOT_FOUND ) {
-				/* No GPU device in this platform */
-				continue; /* Try next platform */
-			} else {
-				ERROR("Couldn't create OpenCL context: %s (%i)\n",
-				clError(err), err);
-				free(gctx);
-				return NULL;
-			}
-		} else {
-			STATUS("Using OpenCL platform %i (%i total)\n",
-			       iplat, nplat);
-			have_ctx = 1;
-			break;
-		}
-	}
-
-	if ( !have_ctx ) {
-		ERROR("Couldn't find a GPU device in any platform.\n");
-		return NULL;
-	}
-
-	dev = get_cl_dev(gctx->ctx, dev_num);
-
-	gctx->cq = clCreateCommandQueue(gctx->ctx, dev, 0, &err);
-	if ( err != CL_SUCCESS ) {
-		ERROR("Couldn't create OpenCL command queue\n");
-		free(gctx);
-		return NULL;
-	}
-
-	/* Create a single-precision version of the scattering factors */
-	intensities_size = IDIM*IDIM*IDIM*sizeof(cl_float);
-	intensities_ptr = malloc(intensities_size);
-	if ( intensities != NULL ) {
-		int i;
-		for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
-			intensities_ptr[i] = intensities[i];
-		}
-	} else {
-		int i;
-		for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
-			intensities_ptr[i] = 100.0;  /* Does nothing */
-		}
-		strncat(cflags, "-DFLAT_INTENSITIES ", 511-strlen(cflags));
-	}
-	gctx->intensities = clCreateBuffer(gctx->ctx,
-	                             CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
-	                             intensities_size, intensities_ptr, &err);
-	if ( err != CL_SUCCESS ) {
-		ERROR("Couldn't allocate intensities memory\n");
-		free(gctx);
-		return NULL;
-	}
-	free(intensities_ptr);
-
-	if ( sym != NULL ) {
-
-		int i, n;
-		SymOpList *pg;
-		size_t islen = 0;
-
-		insert_stuff = malloc(16384);
-		if ( insert_stuff == NULL ) return NULL;
-		insert_stuff[0] = '\0';
-
-		pg = get_pointgroup(sym);
-		n = num_equivs(pg, NULL);
-		for ( i=0; i<n; i++ ) {
-
-			IntegerMatrix *op = get_symop(pg, NULL, i);
-			char line[1024];
-
-			snprintf(line, 1023,
-			         "val += lookup_flagged_intensity(intensities, "
-			         "flags, %s, %s, %s);\n\t",
-			         get_matrix_name(op, 0),
-				 get_matrix_name(op, 1),
-				 get_matrix_name(op, 2));
-
-			islen += strlen(line);
-			if ( islen > 16383 ) {
-				ERROR("Too many symmetry operators.\n");
-				return NULL;
-			}
-			strcat(insert_stuff, line);
-
-		}
-
-		free_symoplist(pg);
-
-		printf("Inserting --->%s<---\n", insert_stuff);
-
-	} else {
-		if ( intensities != NULL ) {
-			ERROR("You gave me an intensities file but no point"
-			      " group.  I'm assuming '1'.\n");
-			strncat(cflags, "-DPG1 ", 511-strlen(cflags));
-		}
-	}
-
-	/* Create a flag array */
-	flags_size = IDIM*IDIM*IDIM*sizeof(cl_float);
-	flags_ptr = malloc(flags_size);
-	if ( flags != NULL ) {
-		int i;
-		for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
-			flags_ptr[i] = flags[i];
-		}
-	} else {
-		int i;
-		for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
-			flags_ptr[i] = 1.0;
-		}
-	}
-	gctx->flags = clCreateBuffer(gctx->ctx,
-	                             CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
-	                             flags_size, flags_ptr, &err);
-	if ( err != CL_SUCCESS ) {
-		ERROR("Couldn't allocate flag buffer\n");
-		free(gctx);
-		return NULL;
-	}
-	free(flags_ptr);
-
-	gctx->prog = load_program_from_string((char *)data_diffraction_cl,
-	                                      data_diffraction_cl_len, gctx->ctx,
-	                                      dev, &err, cflags, insert_stuff);
-	if ( err != CL_SUCCESS ) {
-		free(gctx);
-		return NULL;
-	}
-
-	gctx->kern = clCreateKernel(gctx->prog, "diffraction", &err);
-	if ( err != CL_SUCCESS ) {
-		ERROR("Couldn't create kernel\n");
-		free(gctx);
-		return NULL;
-	}
-
-	gctx->max_sinc_lut = 0;
-	gctx->sinc_lut_ptrs = NULL;
-	gctx->sinc_luts = NULL;
-
-	clGetDeviceInfo(dev, CL_DEVICE_MAX_WORK_GROUP_SIZE,
-	                sizeof(size_t), &maxwgsize, NULL);
-	STATUS("Maximum work group size = %lli\n", (long long int)maxwgsize);
-
-	return gctx;
-}
-
-
-void cleanup_gpu(struct gpu_context *gctx)
-{
-	int i;
-
-	clReleaseProgram(gctx->prog);
-	clReleaseMemObject(gctx->intensities);
-
-	/* Release LUTs */
-	for ( i=1; i<=gctx->max_sinc_lut; i++ ) {
-		if ( gctx->sinc_lut_ptrs[i-1] != NULL ) {
-			clReleaseMemObject(gctx->sinc_luts[i-1]);
-			free(gctx->sinc_lut_ptrs[i-1]);
-		}
-	}
-
-	free(gctx->sinc_luts);
-	free(gctx->sinc_lut_ptrs);
-
-	clReleaseCommandQueue(gctx->cq);
-	clReleaseContext(gctx->ctx);
-	free(gctx);
-}
diff --git a/src/diffraction-gpu.h b/src/diffraction-gpu.h
deleted file mode 100644
index 1128c8b1..00000000
--- a/src/diffraction-gpu.h
+++ /dev/null
@@ -1,77 +0,0 @@
-/*
- * diffraction-gpu.h
- *
- * Calculate diffraction patterns by Fourier methods (GPU version)
- *
- * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
- *                       a research centre of the Helmholtz Association.
- *
- * Authors:
- *   2010-2019 Thomas White <taw@physics.org>
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#ifndef DIFFRACTION_GPU_H
-#define DIFFRACTION_GPU_H
-
-#include "image.h"
-#include "cell.h"
-
-struct gpu_context;
-
-#ifdef HAVE_OPENCL
-
-extern int get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
-                               int na, int nb, int nc, UnitCell *ucell,
-                               int no_fringes, int flat, int n_samples);
-extern struct gpu_context *setup_gpu(int no_sfac,
-                                     const double *intensities,
-                                     const unsigned char *flags,
-                                     const char *sym, int dev_num);
-extern void cleanup_gpu(struct gpu_context *gctx);
-
-#else
-
-static int get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
-                               int na, int nb, int nc, UnitCell *ucell,
-                               int no_fringes, int flat, int n_samples)
-{
-	/* Do nothing */
-	ERROR("This copy of CrystFEL was not compiled with OpenCL support.\n");
-	return 1;
-}
-
-static struct gpu_context *setup_gpu(int no_sfac,
-                                     const double *intensities,
-                                     const unsigned char *flags,
-                                     const char *sym, int dev_num)
-{
-	return NULL;
-}
-
-static void cleanup_gpu(struct gpu_context *gctx)
-{
-}
-
-#endif
-
-#endif	/* DIFFRACTION_GPU_H */
diff --git a/src/diffraction.c b/src/diffraction.c
deleted file mode 100644
index c0dba8aa..00000000
--- a/src/diffraction.c
+++ /dev/null
@@ -1,502 +0,0 @@
-/*
- * diffraction.c
- *
- * Calculate diffraction patterns by Fourier methods
- *
- * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
- *                       a research centre of the Helmholtz Association.
- *
- * Authors:
- *   2009-2020 Thomas White <taw@physics.org>
- *   2013-2014 Chun Hong Yoon <chun.hong.yoon@desy.de>
- *   2013      Alexandra Tolstikova
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-
-#include <stdlib.h>
-#include <math.h>
-#include <stdio.h>
-#include <string.h>
-#include <complex.h>
-#include <assert.h>
-#include <fenv.h>
-
-#include "image.h"
-#include "utils.h"
-#include "cell.h"
-#include "diffraction.h"
-#include "symmetry.h"
-#include "pattern_sim.h"
-
-
-#define SINC_LUT_ELEMENTS (4096)
-
-
-static double *get_sinc_lut(int n, int no_fringes, int flat)
-{
-	int i;
-	double *lut;
-
-	lut = malloc(SINC_LUT_ELEMENTS*sizeof(double));
-	lut[0] = n;
-	if ( n == 1 ) {
-		for ( i=1; i<SINC_LUT_ELEMENTS; i++ ) {
-			lut[i] = 1.0;
-		}
-	} else {
-		for ( i=1; i<SINC_LUT_ELEMENTS; i++ ) {
-			double x, val;
-			x = (double)i/SINC_LUT_ELEMENTS;
-			if ( (flat || no_fringes) && (x > 1.0/n) && (1.0-x > 1.0/n) ) {
-				val = 0.0;
-			} else if ( flat ) {
-				val = n;
-			} else {
-				val = fabs(sin(M_PI*n*x)/sin(M_PI*x));
-			}
-			lut[i] = val;
-		}
-	}
-
-	return lut;
-}
-
-
-static double interpolate_lut(double *lut, double val)
-{
-	double i, pos, f;
-	unsigned int low, high;
-
-	pos = SINC_LUT_ELEMENTS * modf(fabs(val), &i);
-	low = (int)pos;  /* Discard fractional part */
-	high = low + 1;
-	f = modf(pos, &i);  /* Fraction */
-	if ( high == SINC_LUT_ELEMENTS ) high = 0;
-
-	return (1.0-f)*lut[low] + f*lut[high];
-}
-
-
-static double lattice_factor(struct rvec q, double ax, double ay, double az,
-                                            double bx, double by, double bz,
-                                            double cx, double cy, double cz,
-                                            double *lut_a, double *lut_b,
-                                            double *lut_c)
-{
-	struct rvec Udotq;
-	double f1, f2, f3;
-
-	Udotq.u = ax*q.u + ay*q.v + az*q.w;
-	Udotq.v = bx*q.u + by*q.v + bz*q.w;
-	Udotq.w = cx*q.u + cy*q.v + cz*q.w;
-
-	f1 = interpolate_lut(lut_a, Udotq.u);
-	f2 = interpolate_lut(lut_b, Udotq.v);
-	f3 = interpolate_lut(lut_c, Udotq.w);
-
-	return f1 * f2 * f3;
-}
-
-
-static double sym_lookup_intensity(const double *intensities,
-                                   const unsigned char *flags,
-                                   const SymOpList *sym,
-                                   signed int h, signed int k, signed int l)
-{
-	int i;
-	double ret = 0.0;
-
-	for ( i=0; i<num_equivs(sym, NULL); i++ ) {
-
-		signed int he;
-		signed int ke;
-		signed int le;
-		double f, val;
-
-		get_equiv(sym, NULL, i, h, k, l, &he, &ke, &le);
-
-		f = (double)lookup_arr_flag(flags, he, ke, le);
-		val = lookup_arr_intensity(intensities, he, ke, le);
-
-		ret += f*val;
-
-	}
-
-	return ret;
-}
-
-
-static double sym_lookup_phase(const double *phases,
-                               const unsigned char *flags, const SymOpList *sym,
-                               signed int h, signed int k, signed int l)
-{
-	int i;
-
-	for ( i=0; i<num_equivs(sym, NULL); i++ ) {
-
-		signed int he;
-		signed int ke;
-		signed int le;
-		int f;
-
-		get_equiv(sym, NULL, i, h, k, l, &he, &ke, &le);
-
-		f = lookup_arr_flag(flags, he, ke, le);
-
-		if ( f ) return lookup_arr_phase(phases, he, ke, le);
-
-	}
-
-	return 0.0;
-}
-
-
-static double interpolate_linear(const double *ref, const unsigned char *flags,
-                                 const SymOpList *sym, float hd,
-                                 signed int k, signed int l)
-{
-	signed int h;
-	double val1, val2;
-	float f;
-
-	h = (signed int)hd;
-	if ( hd < 0.0 ) h -= 1;
-	f = hd - (float)h;
-	assert(f >= 0.0);
-
-	val1 = sym_lookup_intensity(ref, flags, sym, h, k, l);
-	val2 = sym_lookup_intensity(ref, flags, sym, h+1, k, l);
-
-	return (1.0-f)*val1 + f*val2;
-}
-
-
-static double interpolate_bilinear(const double *ref,
-                                   const unsigned char *flags,
-                                   const SymOpList *sym,
-                                   float hd, float kd, signed int l)
-{
-	signed int k;
-	double val1, val2;
-	float f;
-
-	k = (signed int)kd;
-	if ( kd < 0.0 ) k -= 1;
-	f = kd - (float)k;
-	assert(f >= 0.0);
-
-	val1 = interpolate_linear(ref, flags, sym, hd, k, l);
-	val2 = interpolate_linear(ref, flags, sym, hd, k+1, l);
-
-	return (1.0-f)*val1 + f*val2;
-}
-
-
-static double interpolate_intensity(const double *ref,
-                                    const unsigned char *flags,
-                                    const SymOpList *sym,
-                                    float hd, float kd, float ld)
-{
-	signed int l;
-	double val1, val2;
-	float f;
-
-	l = (signed int)ld;
-	if ( ld < 0.0 ) l -= 1;
-	f = ld - (float)l;
-	assert(f >= 0.0);
-
-	val1 = interpolate_bilinear(ref, flags, sym, hd, kd, l);
-	val2 = interpolate_bilinear(ref, flags, sym, hd, kd, l+1);
-
-	return (1.0-f)*val1 + f*val2;
-}
-
-
-static double complex interpolate_phased_linear(const double *ref,
-                                                const double *phases,
-                                                const unsigned char *flags,
-                                                const SymOpList *sym,
-                                                float hd,
-                                                signed int k, signed int l)
-{
-	signed int h;
-	double val1, val2;
-	float f;
-	double ph1, ph2;
-	double re1, re2, im1, im2;
-	double re, im;
-
-	h = (signed int)hd;
-	if ( hd < 0.0 ) h -= 1;
-	f = hd - (float)h;
-	assert(f >= 0.0);
-
-	val1 = sym_lookup_intensity(ref, flags, sym, h, k, l);
-	val2 = sym_lookup_intensity(ref, flags, sym, h+1, k, l);
-	ph1 = sym_lookup_phase(phases, flags, sym, h, k, l);
-	ph2 = sym_lookup_phase(phases, flags, sym, h+1, k, l);
-
-	/* Calculate real and imaginary parts */
-	re1 = val1 * cos(ph1);
-	im1 = val1 * sin(ph1);
-	re2 = val2 * cos(ph2);
-	im2 = val2 * sin(ph2);
-
-	re = (1.0-f)*re1 + f*re2;
-	im = (1.0-f)*im1 + f*im2;
-
-	return re + im*I;
-}
-
-
-static double complex interpolate_phased_bilinear(const double *ref,
-                                                  const double *phases,
-                                                  const unsigned char *flags,
-                                                  const SymOpList *sym,
-                                                  float hd, float kd,
-                                                  signed int l)
-{
-	signed int k;
-	double complex val1, val2;
-	float f;
-
-	k = (signed int)kd;
-	if ( kd < 0.0 ) k -= 1;
-	f = kd - (float)k;
-	assert(f >= 0.0);
-
-	val1 = interpolate_phased_linear(ref, phases, flags, sym, hd, k, l);
-	val2 = interpolate_phased_linear(ref, phases, flags, sym, hd, k+1, l);
-
-	return (1.0-f)*val1 + f*val2;
-}
-
-
-static double interpolate_phased_intensity(const double *ref,
-                                           const double *phases,
-                                           const unsigned char *flags,
-                                           const SymOpList *sym,
-                                           float hd, float kd, float ld)
-{
-	signed int l;
-	double complex val1, val2;
-	float f;
-
-	l = (signed int)ld;
-	if ( ld < 0.0 ) l -= 1;
-	f = ld - (float)l;
-	assert(f >= 0.0);
-
-	val1 = interpolate_phased_bilinear(ref, phases, flags, sym,
-	                                   hd, kd, l);
-	val2 = interpolate_phased_bilinear(ref, phases, flags, sym,
-	                                   hd, kd, l+1);
-
-	return cabs((1.0-f)*val1 + f*val2);
-}
-
-
-/* Look up the structure factor for the nearest Bragg condition */
-static double molecule_factor(const double *intensities, const double *phases,
-                              const unsigned char *flags, struct rvec q,
-                              double ax, double ay, double az,
-                              double bx, double by, double bz,
-                              double cx, double cy, double cz,
-                              GradientMethod m, const SymOpList *sym)
-{
-	float hd, kd, ld;
-	signed int h, k, l;
-	double r;
-
-	hd = q.u * ax + q.v * ay + q.w * az;
-	kd = q.u * bx + q.v * by + q.w * bz;
-	ld = q.u * cx + q.v * cy + q.w * cz;
-
-	/* No flags -> flat intensity distribution */
-	if ( flags == NULL ) return 100.0;
-
-	switch ( m ) {
-
-		case GRADIENT_MOSAIC :
-		fesetround(1);  /* Round to nearest */
-		h = (signed int)rint(hd);
-		k = (signed int)rint(kd);
-		l = (signed int)rint(ld);
-		if ( abs(h) > INDMAX ) r = 0.0;
-		else if ( abs(k) > INDMAX ) r = 0.0;
-		else if ( abs(l) > INDMAX ) r = 0.0;
-		else r = sym_lookup_intensity(intensities, flags, sym, h, k, l);
-		break;
-
-		case GRADIENT_INTERPOLATE :
-		r = interpolate_intensity(intensities, flags, sym, hd, kd, ld);
-		break;
-
-		case GRADIENT_PHASED :
-		r = interpolate_phased_intensity(intensities, phases, flags,
-		                                 sym, hd, kd, ld);
-		break;
-
-		default:
-		ERROR("This gradient method not implemented yet.\n");
-		exit(1);
-	}
-
-	return r;
-}
-
-
-static void diffraction_panel(struct image *image, const double *intensities,
-                              const double *phases, const unsigned char *flags,
-                              UnitCell *cell, GradientMethod m,
-                              const SymOpList *sym, double k,
-                              double ax, double ay, double az,
-                              double bx, double by, double bz,
-                              double cx, double cy, double cz,
-                              double *lut_a, double *lut_b, double *lut_c,
-                              int pn, double weight)
-{
-	int fs, ss;
-	const int nxs = 4;
-	const int nys = 4;
-	struct detgeom_panel *p = &image->detgeom->panels[pn];
-
-	weight /= nxs*nys;
-
-	for ( ss=0; ss<p->h; ss++ ) {
-	for ( fs=0; fs<p->w; fs++ ) {
-
-		int idx;
-		double f_lattice, I_lattice;
-		double I_molecule;
-		int xs, ys;
-		float xo, yo;
-
-		for ( xs=0; xs<nxs; xs++ ) {
-		for ( ys=0; ys<nys; ys++ ) {
-
-			double qv[3];
-			struct rvec q;
-
-			xo = (1.0/nxs) * xs;
-			yo = (1.0/nys) * ys;
-
-			detgeom_transform_coords(p, fs+xo, ss+yo,
-			                         1.0/k, 0.0, 0.0, qv);
-
-			q.u = qv[0]; q.v = qv[1]; q.w = qv[2];
-
-			f_lattice = lattice_factor(q, ax, ay, az,
-					           bx, by, bz,
-					           cx, cy, cz,
-					           lut_a, lut_b, lut_c);
-
-			I_molecule = molecule_factor(intensities,
-					             phases, flags, q,
-					             ax, ay, az,
-					             bx, by, bz,
-					             cx, cy, cz,
-					             m, sym);
-
-			I_lattice = pow(f_lattice, 2.0);
-
-			idx = fs + p->w*ss;
-			image->dp[pn][idx] += I_lattice * I_molecule * weight;
-
-		}
-		}
-	}
-	progress_bar(ss, p->h-1, "Calculating diffraction");
-	}
-}
-
-
-static void diffraction_at_k(struct image *image, const double *intensities,
-                             const double *phases, const unsigned char *flags,
-                             UnitCell *cell, GradientMethod m,
-                             const SymOpList *sym, double k,
-                             double ax, double ay, double az,
-                             double bx, double by, double bz,
-                             double cx, double cy, double cz,
-                             double *lut_a, double *lut_b, double *lut_c,
-                             double weight)
-{
-	int i;
-
-	for ( i=0; i<image->detgeom->n_panels; i++ ) {
-		diffraction_panel(image, intensities, phases, flags, cell, m,
-		                  sym, k, ax, ay, az, bx, by, bz, cx, cy, cz,
-		                  lut_a, lut_b, lut_c, i, weight);
-	}
-}
-
-
-void get_diffraction(struct image *image, int na, int nb, int nc,
-                     const double *intensities, const double *phases,
-                     const unsigned char *flags, UnitCell *cell,
-                     GradientMethod m, const SymOpList *sym,
-                     int no_fringes, int flat, int n_samples)
-{
-	double ax, ay, az;
-	double bx, by, bz;
-	double cx, cy, cz;
-	double *lut_a;
-	double *lut_b;
-	double *lut_c;
-	int i;
-	double kmin, kmax, step;
-	double norm = 0.0;
-
-	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
-
-	lut_a = get_sinc_lut(na, no_fringes, flat);
-	lut_b = get_sinc_lut(nb, no_fringes, flat);
-	lut_c = get_sinc_lut(nc, no_fringes, flat);
-
-	spectrum_get_range(image->spectrum, &kmin, &kmax);
-	step = (kmax-kmin)/(n_samples+1);
-
-	/* Determine normalisation factor such that weights add up to 1 after
-	 * sampling (bins must have constant width) */
-	for ( i=1; i<=n_samples; i++ ) {
-		double k = kmin + i*step;
-		norm += spectrum_get_density_at_k(image->spectrum, k);
-	}
-	for ( i=1; i<=n_samples; i++ ) {
-
-		double k = kmin + i*step;
-		double prob;
-
-		/* Probability = p.d.f. times step width */
-		prob = spectrum_get_density_at_k(image->spectrum, k)/norm;
-
-		STATUS("Wavelength: %e m, weight = %.5f\n", 1.0/k, prob);
-
-		diffraction_at_k(image, intensities, phases,
-		                 flags, cell, m, sym, k,
-		                 ax, ay, az, bx, by, bz, cx, cy, cz,
-		                 lut_a, lut_b, lut_c, prob);
-
-	}
-
-	free(lut_a);
-	free(lut_b);
-	free(lut_c);
-}
diff --git a/src/diffraction.cl.h b/src/diffraction.cl.h
deleted file mode 100644
index 1f788b67..00000000
--- a/src/diffraction.cl.h
+++ /dev/null
@@ -1,513 +0,0 @@
-/*
- * This file was generated from data/diffraction.cl
- * using the following command:
- * xxd -i data/diffraction.cl src/diffraction.cl.h
- *
- * If you have 'xxd' installed, you can run the script
- * data/gen-resources to re-create this file.
- */
-
-unsigned char data_diffraction_cl[] = {
-  0x2f, 0x2a, 0x0a, 0x20, 0x2a, 0x20, 0x64, 0x69, 0x66, 0x66, 0x72, 0x61,
-  0x63, 0x74, 0x69, 0x6f, 0x6e, 0x2e, 0x63, 0x6c, 0x0a, 0x20, 0x2a, 0x0a,
-  0x20, 0x2a, 0x20, 0x47, 0x50, 0x55, 0x20, 0x63, 0x61, 0x6c, 0x63, 0x75,
-  0x6c, 0x61, 0x74, 0x69, 0x6f, 0x6e, 0x20, 0x6b, 0x65, 0x72, 0x6e, 0x65,
-  0x6c, 0x20, 0x66, 0x6f, 0x72, 0x20, 0x74, 0x72, 0x75, 0x6e, 0x63, 0x61,
-  0x74, 0x65, 0x64, 0x20, 0x6c, 0x61, 0x74, 0x74, 0x69, 0x63, 0x65, 0x20,
-  0x64, 0x69, 0x66, 0x66, 0x72, 0x61, 0x63, 0x74, 0x69, 0x6f, 0x6e, 0x0a,
-  0x20, 0x2a, 0x0a, 0x20, 0x2a, 0x20, 0x43, 0x6f, 0x70, 0x79, 0x72, 0x69,
-  0x67, 0x68, 0x74, 0x20, 0xc2, 0xa9, 0x20, 0x32, 0x30, 0x31, 0x32, 0x2d,
-  0x32, 0x30, 0x31, 0x34, 0x20, 0x44, 0x65, 0x75, 0x74, 0x73, 0x63, 0x68,
-  0x65, 0x73, 0x20, 0x45, 0x6c, 0x65, 0x6b, 0x74, 0x72, 0x6f, 0x6e, 0x65,
-  0x6e, 0x2d, 0x53, 0x79, 0x6e, 0x63, 0x68, 0x72, 0x6f, 0x74, 0x72, 0x6f,
-  0x6e, 0x20, 0x44, 0x45, 0x53, 0x59, 0x2c, 0x0a, 0x20, 0x2a, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x61, 0x20, 0x72,
-  0x65, 0x73, 0x65, 0x61, 0x72, 0x63, 0x68, 0x20, 0x63, 0x65, 0x6e, 0x74,
-  0x72, 0x65, 0x20, 0x6f, 0x66, 0x20, 0x74, 0x68, 0x65, 0x20, 0x48, 0x65,
-  0x6c, 0x6d, 0x68, 0x6f, 0x6c, 0x74, 0x7a, 0x20, 0x41, 0x73, 0x73, 0x6f,
-  0x63, 0x69, 0x61, 0x74, 0x69, 0x6f, 0x6e, 0x2e, 0x0a, 0x20, 0x2a, 0x0a,
-  0x20, 0x2a, 0x20, 0x41, 0x75, 0x74, 0x68, 0x6f, 0x72, 0x73, 0x3a, 0x0a,
-  0x20, 0x2a, 0x20, 0x20, 0x20, 0x32, 0x30, 0x30, 0x39, 0x2d, 0x32, 0x30,
-  0x31, 0x34, 0x20, 0x54, 0x68, 0x6f, 0x6d, 0x61, 0x73, 0x20, 0x57, 0x68,
-  0x69, 0x74, 0x65, 0x20, 0x3c, 0x74, 0x61, 0x77, 0x40, 0x70, 0x68, 0x79,
-  0x73, 0x69, 0x63, 0x73, 0x2e, 0x6f, 0x72, 0x67, 0x3e, 0x0a, 0x20, 0x2a,
-  0x20, 0x20, 0x20, 0x32, 0x30, 0x31, 0x33, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x41, 0x6c, 0x65, 0x78, 0x61, 0x6e, 0x64, 0x72, 0x61, 0x20, 0x54,
-  0x6f, 0x6c, 0x73, 0x74, 0x69, 0x6b, 0x6f, 0x76, 0x61, 0x0a, 0x20, 0x2a,
-  0x0a, 0x20, 0x2a, 0x20, 0x54, 0x68, 0x69, 0x73, 0x20, 0x66, 0x69, 0x6c,
-  0x65, 0x20, 0x69, 0x73, 0x20, 0x70, 0x61, 0x72, 0x74, 0x20, 0x6f, 0x66,
-  0x20, 0x43, 0x72, 0x79, 0x73, 0x74, 0x46, 0x45, 0x4c, 0x2e, 0x0a, 0x20,
-  0x2a, 0x0a, 0x20, 0x2a, 0x20, 0x43, 0x72, 0x79, 0x73, 0x74, 0x46, 0x45,
-  0x4c, 0x20, 0x69, 0x73, 0x20, 0x66, 0x72, 0x65, 0x65, 0x20, 0x73, 0x6f,
-  0x66, 0x74, 0x77, 0x61, 0x72, 0x65, 0x3a, 0x20, 0x79, 0x6f, 0x75, 0x20,
-  0x63, 0x61, 0x6e, 0x20, 0x72, 0x65, 0x64, 0x69, 0x73, 0x74, 0x72, 0x69,
-  0x62, 0x75, 0x74, 0x65, 0x20, 0x69, 0x74, 0x20, 0x61, 0x6e, 0x64, 0x2f,
-  0x6f, 0x72, 0x20, 0x6d, 0x6f, 0x64, 0x69, 0x66, 0x79, 0x0a, 0x20, 0x2a,
-  0x20, 0x69, 0x74, 0x20, 0x75, 0x6e, 0x64, 0x65, 0x72, 0x20, 0x74, 0x68,
-  0x65, 0x20, 0x74, 0x65, 0x72, 0x6d, 0x73, 0x20, 0x6f, 0x66, 0x20, 0x74,
-  0x68, 0x65, 0x20, 0x47, 0x4e, 0x55, 0x20, 0x47, 0x65, 0x6e, 0x65, 0x72,
-  0x61, 0x6c, 0x20, 0x50, 0x75, 0x62, 0x6c, 0x69, 0x63, 0x20, 0x4c, 0x69,
-  0x63, 0x65, 0x6e, 0x73, 0x65, 0x20, 0x61, 0x73, 0x20, 0x70, 0x75, 0x62,
-  0x6c, 0x69, 0x73, 0x68, 0x65, 0x64, 0x20, 0x62, 0x79, 0x0a, 0x20, 0x2a,
-  0x20, 0x74, 0x68, 0x65, 0x20, 0x46, 0x72, 0x65, 0x65, 0x20, 0x53, 0x6f,
-  0x66, 0x74, 0x77, 0x61, 0x72, 0x65, 0x20, 0x46, 0x6f, 0x75, 0x6e, 0x64,
-  0x61, 0x74, 0x69, 0x6f, 0x6e, 0x2c, 0x20, 0x65, 0x69, 0x74, 0x68, 0x65,
-  0x72, 0x20, 0x76, 0x65, 0x72, 0x73, 0x69, 0x6f, 0x6e, 0x20, 0x33, 0x20,
-  0x6f, 0x66, 0x20, 0x74, 0x68, 0x65, 0x20, 0x4c, 0x69, 0x63, 0x65, 0x6e,
-  0x73, 0x65, 0x2c, 0x20, 0x6f, 0x72, 0x0a, 0x20, 0x2a, 0x20, 0x28, 0x61,
-  0x74, 0x20, 0x79, 0x6f, 0x75, 0x72, 0x20, 0x6f, 0x70, 0x74, 0x69, 0x6f,
-  0x6e, 0x29, 0x20, 0x61, 0x6e, 0x79, 0x20, 0x6c, 0x61, 0x74, 0x65, 0x72,
-  0x20, 0x76, 0x65, 0x72, 0x73, 0x69, 0x6f, 0x6e, 0x2e, 0x0a, 0x20, 0x2a,
-  0x0a, 0x20, 0x2a, 0x20, 0x43, 0x72, 0x79, 0x73, 0x74, 0x46, 0x45, 0x4c,
-  0x20, 0x69, 0x73, 0x20, 0x64, 0x69, 0x73, 0x74, 0x72, 0x69, 0x62, 0x75,
-  0x74, 0x65, 0x64, 0x20, 0x69, 0x6e, 0x20, 0x74, 0x68, 0x65, 0x20, 0x68,
-  0x6f, 0x70, 0x65, 0x20, 0x74, 0x68, 0x61, 0x74, 0x20, 0x69, 0x74, 0x20,
-  0x77, 0x69, 0x6c, 0x6c, 0x20, 0x62, 0x65, 0x20, 0x75, 0x73, 0x65, 0x66,
-  0x75, 0x6c, 0x2c, 0x0a, 0x20, 0x2a, 0x20, 0x62, 0x75, 0x74, 0x20, 0x57,
-  0x49, 0x54, 0x48, 0x4f, 0x55, 0x54, 0x20, 0x41, 0x4e, 0x59, 0x20, 0x57,
-  0x41, 0x52, 0x52, 0x41, 0x4e, 0x54, 0x59, 0x3b, 0x20, 0x77, 0x69, 0x74,
-  0x68, 0x6f, 0x75, 0x74, 0x20, 0x65, 0x76, 0x65, 0x6e, 0x20, 0x74, 0x68,
-  0x65, 0x20, 0x69, 0x6d, 0x70, 0x6c, 0x69, 0x65, 0x64, 0x20, 0x77, 0x61,
-  0x72, 0x72, 0x61, 0x6e, 0x74, 0x79, 0x20, 0x6f, 0x66, 0x0a, 0x20, 0x2a,
-  0x20, 0x4d, 0x45, 0x52, 0x43, 0x48, 0x41, 0x4e, 0x54, 0x41, 0x42, 0x49,
-  0x4c, 0x49, 0x54, 0x59, 0x20, 0x6f, 0x72, 0x20, 0x46, 0x49, 0x54, 0x4e,
-  0x45, 0x53, 0x53, 0x20, 0x46, 0x4f, 0x52, 0x20, 0x41, 0x20, 0x50, 0x41,
-  0x52, 0x54, 0x49, 0x43, 0x55, 0x4c, 0x41, 0x52, 0x20, 0x50, 0x55, 0x52,
-  0x50, 0x4f, 0x53, 0x45, 0x2e, 0x20, 0x20, 0x53, 0x65, 0x65, 0x20, 0x74,
-  0x68, 0x65, 0x0a, 0x20, 0x2a, 0x20, 0x47, 0x4e, 0x55, 0x20, 0x47, 0x65,
-  0x6e, 0x65, 0x72, 0x61, 0x6c, 0x20, 0x50, 0x75, 0x62, 0x6c, 0x69, 0x63,
-  0x20, 0x4c, 0x69, 0x63, 0x65, 0x6e, 0x73, 0x65, 0x20, 0x66, 0x6f, 0x72,
-  0x20, 0x6d, 0x6f, 0x72, 0x65, 0x20, 0x64, 0x65, 0x74, 0x61, 0x69, 0x6c,
-  0x73, 0x2e, 0x0a, 0x20, 0x2a, 0x0a, 0x20, 0x2a, 0x20, 0x59, 0x6f, 0x75,
-  0x20, 0x73, 0x68, 0x6f, 0x75, 0x6c, 0x64, 0x20, 0x68, 0x61, 0x76, 0x65,
-  0x20, 0x72, 0x65, 0x63, 0x65, 0x69, 0x76, 0x65, 0x64, 0x20, 0x61, 0x20,
-  0x63, 0x6f, 0x70, 0x79, 0x20, 0x6f, 0x66, 0x20, 0x74, 0x68, 0x65, 0x20,
-  0x47, 0x4e, 0x55, 0x20, 0x47, 0x65, 0x6e, 0x65, 0x72, 0x61, 0x6c, 0x20,
-  0x50, 0x75, 0x62, 0x6c, 0x69, 0x63, 0x20, 0x4c, 0x69, 0x63, 0x65, 0x6e,
-  0x73, 0x65, 0x0a, 0x20, 0x2a, 0x20, 0x61, 0x6c, 0x6f, 0x6e, 0x67, 0x20,
-  0x77, 0x69, 0x74, 0x68, 0x20, 0x43, 0x72, 0x79, 0x73, 0x74, 0x46, 0x45,
-  0x4c, 0x2e, 0x20, 0x20, 0x49, 0x66, 0x20, 0x6e, 0x6f, 0x74, 0x2c, 0x20,
-  0x73, 0x65, 0x65, 0x20, 0x3c, 0x68, 0x74, 0x74, 0x70, 0x3a, 0x2f, 0x2f,
-  0x77, 0x77, 0x77, 0x2e, 0x67, 0x6e, 0x75, 0x2e, 0x6f, 0x72, 0x67, 0x2f,
-  0x6c, 0x69, 0x63, 0x65, 0x6e, 0x73, 0x65, 0x73, 0x2f, 0x3e, 0x2e, 0x0a,
-  0x20, 0x2a, 0x0a, 0x20, 0x2a, 0x2f, 0x0a, 0x0a, 0x0a, 0x2f, 0x2a, 0x20,
-  0x4d, 0x61, 0x78, 0x6d, 0x69, 0x6d, 0x75, 0x6d, 0x20, 0x69, 0x6e, 0x64,
-  0x65, 0x78, 0x20, 0x74, 0x6f, 0x20, 0x68, 0x6f, 0x6c, 0x64, 0x20, 0x76,
-  0x61, 0x6c, 0x75, 0x65, 0x73, 0x20, 0x75, 0x70, 0x20, 0x74, 0x6f, 0x20,
-  0x28, 0x63, 0x61, 0x6e, 0x20, 0x62, 0x65, 0x20, 0x69, 0x6e, 0x63, 0x72,
-  0x65, 0x61, 0x73, 0x65, 0x64, 0x20, 0x69, 0x66, 0x20, 0x6e, 0x65, 0x63,
-  0x65, 0x73, 0x73, 0x61, 0x72, 0x79, 0x29, 0x0a, 0x20, 0x2a, 0x20, 0x57,
-  0x41, 0x52, 0x4e, 0x49, 0x4e, 0x47, 0x3a, 0x20, 0x41, 0x6c, 0x74, 0x65,
-  0x72, 0x69, 0x6e, 0x67, 0x20, 0x74, 0x68, 0x69, 0x73, 0x20, 0x76, 0x61,
-  0x6c, 0x75, 0x65, 0x20, 0x63, 0x6f, 0x6e, 0x73, 0x74, 0x69, 0x74, 0x75,
-  0x74, 0x65, 0x73, 0x20, 0x61, 0x6e, 0x20, 0x41, 0x42, 0x49, 0x20, 0x63,
-  0x68, 0x61, 0x6e, 0x67, 0x65, 0x2c, 0x20, 0x61, 0x6e, 0x64, 0x20, 0x6d,
-  0x65, 0x61, 0x6e, 0x73, 0x20, 0x79, 0x6f, 0x75, 0x20, 0x6d, 0x75, 0x73,
-  0x74, 0x0a, 0x20, 0x2a, 0x20, 0x75, 0x70, 0x64, 0x61, 0x74, 0x65, 0x20,
-  0x73, 0x72, 0x63, 0x2f, 0x70, 0x61, 0x74, 0x74, 0x65, 0x72, 0x6e, 0x5f,
-  0x73, 0x69, 0x6d, 0x2e, 0x68, 0x20, 0x74, 0x68, 0x65, 0x6e, 0x20, 0x72,
-  0x65, 0x63, 0x6f, 0x6d, 0x70, 0x69, 0x6c, 0x65, 0x20, 0x61, 0x6e, 0x64,
-  0x20, 0x72, 0x65, 0x69, 0x6e, 0x73, 0x74, 0x61, 0x6c, 0x6c, 0x20, 0x65,
-  0x76, 0x65, 0x72, 0x79, 0x74, 0x68, 0x69, 0x6e, 0x67, 0x2e, 0x20, 0x2a,
-  0x2f, 0x0a, 0x23, 0x64, 0x65, 0x66, 0x69, 0x6e, 0x65, 0x20, 0x49, 0x4e,
-  0x44, 0x4d, 0x41, 0x58, 0x20, 0x31, 0x33, 0x30, 0x0a, 0x23, 0x64, 0x65,
-  0x66, 0x69, 0x6e, 0x65, 0x20, 0x49, 0x44, 0x49, 0x4d, 0x20, 0x28, 0x49,
-  0x4e, 0x44, 0x4d, 0x41, 0x58, 0x2a, 0x32, 0x20, 0x2b, 0x31, 0x29, 0x0a,
-  0x0a, 0x23, 0x69, 0x66, 0x6e, 0x64, 0x65, 0x66, 0x20, 0x4d, 0x5f, 0x50,
-  0x49, 0x0a, 0x23, 0x64, 0x65, 0x66, 0x69, 0x6e, 0x65, 0x20, 0x4d, 0x5f,
-  0x50, 0x49, 0x20, 0x28, 0x28, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x29, 0x28,
-  0x33, 0x2e, 0x31, 0x34, 0x31, 0x35, 0x39, 0x32, 0x36, 0x35, 0x29, 0x29,
-  0x0a, 0x23, 0x65, 0x6e, 0x64, 0x69, 0x66, 0x0a, 0x0a, 0x63, 0x6f, 0x6e,
-  0x73, 0x74, 0x20, 0x73, 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x72, 0x5f, 0x74,
-  0x20, 0x73, 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x72, 0x5f, 0x61, 0x20, 0x3d,
-  0x20, 0x43, 0x4c, 0x4b, 0x5f, 0x4e, 0x4f, 0x52, 0x4d, 0x41, 0x4c, 0x49,
-  0x5a, 0x45, 0x44, 0x5f, 0x43, 0x4f, 0x4f, 0x52, 0x44, 0x53, 0x5f, 0x54,
-  0x52, 0x55, 0x45, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x7c, 0x20, 0x43,
-  0x4c, 0x4b, 0x5f, 0x41, 0x44, 0x44, 0x52, 0x45, 0x53, 0x53, 0x5f, 0x52,
-  0x45, 0x50, 0x45, 0x41, 0x54, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x7c,
-  0x20, 0x43, 0x4c, 0x4b, 0x5f, 0x46, 0x49, 0x4c, 0x54, 0x45, 0x52, 0x5f,
-  0x4c, 0x49, 0x4e, 0x45, 0x41, 0x52, 0x3b, 0x0a, 0x63, 0x6f, 0x6e, 0x73,
-  0x74, 0x20, 0x73, 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x72, 0x5f, 0x74, 0x20,
-  0x73, 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x72, 0x5f, 0x62, 0x20, 0x3d, 0x20,
-  0x43, 0x4c, 0x4b, 0x5f, 0x4e, 0x4f, 0x52, 0x4d, 0x41, 0x4c, 0x49, 0x5a,
-  0x45, 0x44, 0x5f, 0x43, 0x4f, 0x4f, 0x52, 0x44, 0x53, 0x5f, 0x54, 0x52,
-  0x55, 0x45, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x7c, 0x20, 0x43, 0x4c,
-  0x4b, 0x5f, 0x41, 0x44, 0x44, 0x52, 0x45, 0x53, 0x53, 0x5f, 0x52, 0x45,
-  0x50, 0x45, 0x41, 0x54, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x7c, 0x20,
-  0x43, 0x4c, 0x4b, 0x5f, 0x46, 0x49, 0x4c, 0x54, 0x45, 0x52, 0x5f, 0x4c,
-  0x49, 0x4e, 0x45, 0x41, 0x52, 0x3b, 0x0a, 0x63, 0x6f, 0x6e, 0x73, 0x74,
-  0x20, 0x73, 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x72, 0x5f, 0x74, 0x20, 0x73,
-  0x61, 0x6d, 0x70, 0x6c, 0x65, 0x72, 0x5f, 0x63, 0x20, 0x3d, 0x20, 0x43,
-  0x4c, 0x4b, 0x5f, 0x4e, 0x4f, 0x52, 0x4d, 0x41, 0x4c, 0x49, 0x5a, 0x45,
-  0x44, 0x5f, 0x43, 0x4f, 0x4f, 0x52, 0x44, 0x53, 0x5f, 0x54, 0x52, 0x55,
-  0x45, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x7c, 0x20, 0x43, 0x4c, 0x4b,
-  0x5f, 0x41, 0x44, 0x44, 0x52, 0x45, 0x53, 0x53, 0x5f, 0x52, 0x45, 0x50,
-  0x45, 0x41, 0x54, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x7c, 0x20, 0x43,
-  0x4c, 0x4b, 0x5f, 0x46, 0x49, 0x4c, 0x54, 0x45, 0x52, 0x5f, 0x4c, 0x49,
-  0x4e, 0x45, 0x41, 0x52, 0x3b, 0x0a, 0x0a, 0x0a, 0x66, 0x6c, 0x6f, 0x61,
-  0x74, 0x34, 0x20, 0x67, 0x65, 0x74, 0x5f, 0x71, 0x28, 0x66, 0x6c, 0x6f,
-  0x61, 0x74, 0x20, 0x66, 0x73, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74,
-  0x20, 0x73, 0x73, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x72,
-  0x65, 0x73, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x63, 0x6c,
-  0x65, 0x6e, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x6b, 0x2c,
-  0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x63, 0x6f, 0x72, 0x6e,
-  0x65, 0x72, 0x5f, 0x78, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20,
-  0x63, 0x6f, 0x72, 0x6e, 0x65, 0x72, 0x5f, 0x79, 0x2c, 0x0a, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x66,
-  0x6c, 0x6f, 0x61, 0x74, 0x20, 0x66, 0x73, 0x78, 0x2c, 0x20, 0x66, 0x6c,
-  0x6f, 0x61, 0x74, 0x20, 0x66, 0x73, 0x79, 0x2c, 0x20, 0x66, 0x6c, 0x6f,
-  0x61, 0x74, 0x20, 0x66, 0x73, 0x7a, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61,
-  0x74, 0x20, 0x73, 0x73, 0x78, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74,
-  0x20, 0x73, 0x73, 0x79, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20,
-  0x73, 0x73, 0x7a, 0x29, 0x0a, 0x7b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61,
-  0x74, 0x20, 0x72, 0x78, 0x2c, 0x20, 0x72, 0x79, 0x2c, 0x20, 0x72, 0x7a,
-  0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x61, 0x7a, 0x2c,
-  0x20, 0x74, 0x74, 0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x34,
-  0x20, 0x71, 0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x78,
-  0x73, 0x2c, 0x20, 0x79, 0x73, 0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61,
-  0x74, 0x20, 0x6b, 0x78, 0x2c, 0x20, 0x6b, 0x79, 0x2c, 0x20, 0x6b, 0x7a,
-  0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x63, 0x74, 0x74,
-  0x3b, 0x0a, 0x0a, 0x09, 0x2f, 0x2a, 0x20, 0x43, 0x61, 0x6c, 0x63, 0x75,
-  0x6c, 0x61, 0x74, 0x65, 0x20, 0x33, 0x44, 0x20, 0x70, 0x6f, 0x73, 0x69,
-  0x74, 0x69, 0x6f, 0x6e, 0x20, 0x6f, 0x66, 0x20, 0x67, 0x69, 0x76, 0x65,
-  0x6e, 0x20, 0x70, 0x6f, 0x73, 0x69, 0x74, 0x69, 0x6f, 0x6e, 0x2c, 0x20,
-  0x69, 0x6e, 0x20, 0x6d, 0x20, 0x2a, 0x2f, 0x0a, 0x09, 0x72, 0x78, 0x20,
-  0x3d, 0x20, 0x28, 0x63, 0x6f, 0x72, 0x6e, 0x65, 0x72, 0x5f, 0x78, 0x20,
-  0x2b, 0x20, 0x66, 0x73, 0x2a, 0x66, 0x73, 0x78, 0x20, 0x2b, 0x20, 0x73,
-  0x73, 0x2a, 0x73, 0x73, 0x78, 0x29, 0x20, 0x2f, 0x20, 0x72, 0x65, 0x73,
-  0x3b, 0x0a, 0x09, 0x72, 0x79, 0x20, 0x3d, 0x20, 0x28, 0x63, 0x6f, 0x72,
-  0x6e, 0x65, 0x72, 0x5f, 0x79, 0x20, 0x2b, 0x20, 0x66, 0x73, 0x2a, 0x66,
-  0x73, 0x79, 0x20, 0x2b, 0x20, 0x73, 0x73, 0x2a, 0x73, 0x73, 0x79, 0x29,
-  0x20, 0x2f, 0x20, 0x72, 0x65, 0x73, 0x3b, 0x0a, 0x09, 0x72, 0x7a, 0x20,
-  0x3d, 0x20, 0x63, 0x6c, 0x65, 0x6e, 0x20, 0x2b, 0x20, 0x28, 0x66, 0x73,
-  0x2a, 0x66, 0x73, 0x7a, 0x20, 0x2b, 0x20, 0x73, 0x73, 0x2a, 0x73, 0x73,
-  0x7a, 0x29, 0x2f, 0x72, 0x65, 0x73, 0x3b, 0x0a, 0x0a, 0x09, 0x63, 0x74,
-  0x74, 0x20, 0x3d, 0x20, 0x72, 0x7a, 0x20, 0x2f, 0x20, 0x73, 0x71, 0x72,
-  0x74, 0x28, 0x72, 0x78, 0x2a, 0x72, 0x78, 0x20, 0x2b, 0x20, 0x72, 0x79,
-  0x2a, 0x72, 0x79, 0x20, 0x2b, 0x20, 0x72, 0x7a, 0x2a, 0x72, 0x7a, 0x29,
-  0x3b, 0x20, 0x20, 0x2f, 0x2a, 0x20, 0x63, 0x6f, 0x73, 0x28, 0x32, 0x74,
-  0x68, 0x65, 0x74, 0x61, 0x29, 0x20, 0x2a, 0x2f, 0x0a, 0x09, 0x74, 0x74,
-  0x20, 0x3d, 0x20, 0x61, 0x63, 0x6f, 0x73, 0x28, 0x63, 0x74, 0x74, 0x29,
-  0x3b, 0x0a, 0x09, 0x61, 0x7a, 0x20, 0x3d, 0x20, 0x61, 0x74, 0x61, 0x6e,
-  0x32, 0x28, 0x72, 0x79, 0x2c, 0x20, 0x72, 0x78, 0x29, 0x3b, 0x0a, 0x0a,
-  0x09, 0x6b, 0x78, 0x20, 0x3d, 0x20, 0x6b, 0x2a, 0x6e, 0x61, 0x74, 0x69,
-  0x76, 0x65, 0x5f, 0x73, 0x69, 0x6e, 0x28, 0x74, 0x74, 0x29, 0x2a, 0x6e,
-  0x61, 0x74, 0x69, 0x76, 0x65, 0x5f, 0x63, 0x6f, 0x73, 0x28, 0x61, 0x7a,
-  0x29, 0x3b, 0x0a, 0x09, 0x6b, 0x79, 0x20, 0x3d, 0x20, 0x6b, 0x2a, 0x6e,
-  0x61, 0x74, 0x69, 0x76, 0x65, 0x5f, 0x73, 0x69, 0x6e, 0x28, 0x74, 0x74,
-  0x29, 0x2a, 0x6e, 0x61, 0x74, 0x69, 0x76, 0x65, 0x5f, 0x73, 0x69, 0x6e,
-  0x28, 0x61, 0x7a, 0x29, 0x3b, 0x0a, 0x09, 0x6b, 0x7a, 0x20, 0x3d, 0x20,
-  0x6b, 0x2a, 0x28, 0x63, 0x74, 0x74, 0x20, 0x2d, 0x20, 0x31, 0x2e, 0x30,
-  0x29, 0x3b, 0x0a, 0x0a, 0x09, 0x71, 0x20, 0x3d, 0x20, 0x28, 0x66, 0x6c,
-  0x6f, 0x61, 0x74, 0x34, 0x29, 0x28, 0x6b, 0x78, 0x2c, 0x20, 0x6b, 0x79,
-  0x2c, 0x20, 0x6b, 0x7a, 0x2c, 0x20, 0x30, 0x2e, 0x30, 0x29, 0x3b, 0x0a,
-  0x0a, 0x09, 0x72, 0x65, 0x74, 0x75, 0x72, 0x6e, 0x20, 0x71, 0x3b, 0x0a,
-  0x7d, 0x0a, 0x0a, 0x0a, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x6c, 0x61,
-  0x74, 0x74, 0x69, 0x63, 0x65, 0x5f, 0x66, 0x61, 0x63, 0x74, 0x6f, 0x72,
-  0x28, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x31, 0x36, 0x20, 0x63, 0x65, 0x6c,
-  0x6c, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x34, 0x20, 0x71, 0x2c,
-  0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x72, 0x65,
-  0x61, 0x64, 0x5f, 0x6f, 0x6e, 0x6c, 0x79, 0x20, 0x69, 0x6d, 0x61, 0x67,
-  0x65, 0x32, 0x64, 0x5f, 0x74, 0x20, 0x66, 0x75, 0x6e, 0x63, 0x5f, 0x61,
-  0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x72,
-  0x65, 0x61, 0x64, 0x5f, 0x6f, 0x6e, 0x6c, 0x79, 0x20, 0x69, 0x6d, 0x61,
-  0x67, 0x65, 0x32, 0x64, 0x5f, 0x74, 0x20, 0x66, 0x75, 0x6e, 0x63, 0x5f,
-  0x62, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x72, 0x65, 0x61, 0x64, 0x5f, 0x6f, 0x6e, 0x6c, 0x79, 0x20, 0x69, 0x6d,
-  0x61, 0x67, 0x65, 0x32, 0x64, 0x5f, 0x74, 0x20, 0x66, 0x75, 0x6e, 0x63,
-  0x5f, 0x63, 0x29, 0x0a, 0x7b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74,
-  0x20, 0x66, 0x31, 0x2c, 0x20, 0x66, 0x32, 0x2c, 0x20, 0x66, 0x33, 0x2c,
-  0x20, 0x76, 0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x34, 0x20,
-  0x55, 0x64, 0x6f, 0x74, 0x71, 0x3b, 0x0a, 0x0a, 0x09, 0x55, 0x64, 0x6f,
-  0x74, 0x71, 0x2e, 0x78, 0x20, 0x3d, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e,
-  0x73, 0x30, 0x2a, 0x71, 0x2e, 0x78, 0x20, 0x2b, 0x20, 0x63, 0x65, 0x6c,
-  0x6c, 0x2e, 0x73, 0x31, 0x2a, 0x71, 0x2e, 0x79, 0x20, 0x2b, 0x20, 0x63,
-  0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x32, 0x2a, 0x71, 0x2e, 0x7a, 0x3b, 0x0a,
-  0x09, 0x55, 0x64, 0x6f, 0x74, 0x71, 0x2e, 0x79, 0x20, 0x3d, 0x20, 0x63,
-  0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x33, 0x2a, 0x71, 0x2e, 0x78, 0x20, 0x2b,
-  0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x34, 0x2a, 0x71, 0x2e, 0x79,
-  0x20, 0x2b, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x35, 0x2a, 0x71,
-  0x2e, 0x7a, 0x3b, 0x0a, 0x09, 0x55, 0x64, 0x6f, 0x74, 0x71, 0x2e, 0x7a,
-  0x20, 0x3d, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x36, 0x2a, 0x71,
-  0x2e, 0x78, 0x20, 0x2b, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x37,
-  0x2a, 0x71, 0x2e, 0x79, 0x20, 0x2b, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e,
-  0x73, 0x38, 0x2a, 0x71, 0x2e, 0x7a, 0x3b, 0x0a, 0x0a, 0x09, 0x2f, 0x2a,
-  0x20, 0x4c, 0x6f, 0x6f, 0x6b, 0x20, 0x75, 0x70, 0x20, 0x76, 0x61, 0x6c,
-  0x75, 0x65, 0x73, 0x20, 0x66, 0x72, 0x6f, 0x6d, 0x20, 0x70, 0x72, 0x65,
-  0x63, 0x61, 0x6c, 0x63, 0x75, 0x6c, 0x61, 0x74, 0x65, 0x64, 0x20, 0x73,
-  0x69, 0x6e, 0x63, 0x28, 0x29, 0x20, 0x74, 0x61, 0x62, 0x6c, 0x65, 0x20,
-  0x2a, 0x2f, 0x0a, 0x09, 0x66, 0x31, 0x20, 0x3d, 0x20, 0x72, 0x65, 0x61,
-  0x64, 0x5f, 0x69, 0x6d, 0x61, 0x67, 0x65, 0x66, 0x28, 0x66, 0x75, 0x6e,
-  0x63, 0x5f, 0x61, 0x2c, 0x20, 0x73, 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x72,
-  0x5f, 0x61, 0x2c, 0x20, 0x28, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x32, 0x29,
-  0x28, 0x55, 0x64, 0x6f, 0x74, 0x71, 0x2e, 0x78, 0x2c, 0x20, 0x30, 0x2e,
-  0x30, 0x29, 0x29, 0x2e, 0x73, 0x30, 0x3b, 0x0a, 0x09, 0x66, 0x32, 0x20,
-  0x3d, 0x20, 0x72, 0x65, 0x61, 0x64, 0x5f, 0x69, 0x6d, 0x61, 0x67, 0x65,
-  0x66, 0x28, 0x66, 0x75, 0x6e, 0x63, 0x5f, 0x62, 0x2c, 0x20, 0x73, 0x61,
-  0x6d, 0x70, 0x6c, 0x65, 0x72, 0x5f, 0x62, 0x2c, 0x20, 0x28, 0x66, 0x6c,
-  0x6f, 0x61, 0x74, 0x32, 0x29, 0x28, 0x55, 0x64, 0x6f, 0x74, 0x71, 0x2e,
-  0x79, 0x2c, 0x20, 0x30, 0x2e, 0x30, 0x29, 0x29, 0x2e, 0x73, 0x30, 0x3b,
-  0x0a, 0x09, 0x66, 0x33, 0x20, 0x3d, 0x20, 0x72, 0x65, 0x61, 0x64, 0x5f,
-  0x69, 0x6d, 0x61, 0x67, 0x65, 0x66, 0x28, 0x66, 0x75, 0x6e, 0x63, 0x5f,
-  0x63, 0x2c, 0x20, 0x73, 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x72, 0x5f, 0x63,
-  0x2c, 0x20, 0x28, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x32, 0x29, 0x28, 0x55,
-  0x64, 0x6f, 0x74, 0x71, 0x2e, 0x7a, 0x2c, 0x20, 0x30, 0x2e, 0x30, 0x29,
-  0x29, 0x2e, 0x73, 0x30, 0x3b, 0x0a, 0x0a, 0x09, 0x72, 0x65, 0x74, 0x75,
-  0x72, 0x6e, 0x20, 0x66, 0x31, 0x20, 0x2a, 0x20, 0x66, 0x32, 0x20, 0x2a,
-  0x20, 0x66, 0x33, 0x3b, 0x0a, 0x7d, 0x0a, 0x0a, 0x0a, 0x66, 0x6c, 0x6f,
-  0x61, 0x74, 0x20, 0x6c, 0x6f, 0x6f, 0x6b, 0x75, 0x70, 0x5f, 0x69, 0x6e,
-  0x74, 0x65, 0x6e, 0x73, 0x69, 0x74, 0x79, 0x28, 0x67, 0x6c, 0x6f, 0x62,
-  0x61, 0x6c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x2a, 0x69, 0x6e,
-  0x74, 0x65, 0x6e, 0x73, 0x69, 0x74, 0x69, 0x65, 0x73, 0x2c, 0x0a, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x73, 0x69,
-  0x67, 0x6e, 0x65, 0x64, 0x20, 0x69, 0x6e, 0x74, 0x20, 0x68, 0x2c, 0x20,
-  0x73, 0x69, 0x67, 0x6e, 0x65, 0x64, 0x20, 0x69, 0x6e, 0x74, 0x20, 0x6b,
-  0x2c, 0x20, 0x73, 0x69, 0x67, 0x6e, 0x65, 0x64, 0x20, 0x69, 0x6e, 0x74,
-  0x20, 0x6c, 0x29, 0x0a, 0x7b, 0x0a, 0x09, 0x69, 0x6e, 0x74, 0x20, 0x69,
-  0x64, 0x78, 0x3b, 0x0a, 0x0a, 0x09, 0x2f, 0x2a, 0x20, 0x4f, 0x75, 0x74,
-  0x20, 0x6f, 0x66, 0x20, 0x72, 0x61, 0x6e, 0x67, 0x65, 0x3f, 0x20, 0x2a,
-  0x2f, 0x0a, 0x09, 0x69, 0x66, 0x20, 0x28, 0x20, 0x28, 0x61, 0x62, 0x73,
-  0x28, 0x68, 0x29, 0x20, 0x3e, 0x20, 0x49, 0x4e, 0x44, 0x4d, 0x41, 0x58,
-  0x29, 0x20, 0x7c, 0x7c, 0x20, 0x28, 0x61, 0x62, 0x73, 0x28, 0x6b, 0x29,
-  0x20, 0x3e, 0x20, 0x49, 0x4e, 0x44, 0x4d, 0x41, 0x58, 0x29, 0x20, 0x7c,
-  0x7c, 0x20, 0x28, 0x61, 0x62, 0x73, 0x28, 0x6c, 0x29, 0x20, 0x3e, 0x20,
-  0x49, 0x4e, 0x44, 0x4d, 0x41, 0x58, 0x29, 0x20, 0x29, 0x20, 0x7b, 0x0a,
-  0x09, 0x09, 0x72, 0x65, 0x74, 0x75, 0x72, 0x6e, 0x20, 0x30, 0x2e, 0x30,
-  0x3b, 0x0a, 0x09, 0x7d, 0x0a, 0x0a, 0x09, 0x68, 0x20, 0x3d, 0x20, 0x28,
-  0x68, 0x3e, 0x3d, 0x30, 0x29, 0x20, 0x3f, 0x20, 0x68, 0x20, 0x3a, 0x20,
-  0x68, 0x2b, 0x49, 0x44, 0x49, 0x4d, 0x3b, 0x0a, 0x09, 0x6b, 0x20, 0x3d,
-  0x20, 0x28, 0x6b, 0x3e, 0x3d, 0x30, 0x29, 0x20, 0x3f, 0x20, 0x6b, 0x20,
-  0x3a, 0x20, 0x6b, 0x2b, 0x49, 0x44, 0x49, 0x4d, 0x3b, 0x0a, 0x09, 0x6c,
-  0x20, 0x3d, 0x20, 0x28, 0x6c, 0x3e, 0x3d, 0x30, 0x29, 0x20, 0x3f, 0x20,
-  0x6c, 0x20, 0x3a, 0x20, 0x6c, 0x2b, 0x49, 0x44, 0x49, 0x4d, 0x3b, 0x0a,
-  0x0a, 0x09, 0x69, 0x64, 0x78, 0x20, 0x3d, 0x20, 0x68, 0x20, 0x2b, 0x20,
-  0x28, 0x49, 0x44, 0x49, 0x4d, 0x2a, 0x6b, 0x29, 0x20, 0x2b, 0x20, 0x28,
-  0x49, 0x44, 0x49, 0x4d, 0x2a, 0x49, 0x44, 0x49, 0x4d, 0x2a, 0x6c, 0x29,
-  0x3b, 0x0a, 0x0a, 0x09, 0x72, 0x65, 0x74, 0x75, 0x72, 0x6e, 0x20, 0x69,
-  0x6e, 0x74, 0x65, 0x6e, 0x73, 0x69, 0x74, 0x69, 0x65, 0x73, 0x5b, 0x69,
-  0x64, 0x78, 0x5d, 0x3b, 0x0a, 0x7d, 0x0a, 0x0a, 0x0a, 0x66, 0x6c, 0x6f,
-  0x61, 0x74, 0x20, 0x6c, 0x6f, 0x6f, 0x6b, 0x75, 0x70, 0x5f, 0x66, 0x6c,
-  0x61, 0x67, 0x67, 0x65, 0x64, 0x5f, 0x69, 0x6e, 0x74, 0x65, 0x6e, 0x73,
-  0x69, 0x74, 0x79, 0x28, 0x67, 0x6c, 0x6f, 0x62, 0x61, 0x6c, 0x20, 0x66,
-  0x6c, 0x6f, 0x61, 0x74, 0x20, 0x2a, 0x69, 0x6e, 0x74, 0x65, 0x6e, 0x73,
-  0x69, 0x74, 0x69, 0x65, 0x73, 0x2c, 0x20, 0x67, 0x6c, 0x6f, 0x62, 0x61,
-  0x6c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x2a, 0x66, 0x6c, 0x61,
-  0x67, 0x73, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x73,
-  0x69, 0x67, 0x6e, 0x65, 0x64, 0x20, 0x69, 0x6e, 0x74, 0x20, 0x68, 0x2c,
-  0x20, 0x73, 0x69, 0x67, 0x6e, 0x65, 0x64, 0x20, 0x69, 0x6e, 0x74, 0x20,
-  0x6b, 0x2c, 0x20, 0x73, 0x69, 0x67, 0x6e, 0x65, 0x64, 0x20, 0x69, 0x6e,
-  0x74, 0x20, 0x6c, 0x29, 0x0a, 0x7b, 0x0a, 0x09, 0x72, 0x65, 0x74, 0x75,
-  0x72, 0x6e, 0x20, 0x6c, 0x6f, 0x6f, 0x6b, 0x75, 0x70, 0x5f, 0x69, 0x6e,
-  0x74, 0x65, 0x6e, 0x73, 0x69, 0x74, 0x79, 0x28, 0x69, 0x6e, 0x74, 0x65,
-  0x6e, 0x73, 0x69, 0x74, 0x69, 0x65, 0x73, 0x2c, 0x20, 0x68, 0x2c, 0x20,
-  0x6b, 0x2c, 0x20, 0x6c, 0x29, 0x0a, 0x09, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x2a, 0x20, 0x6c, 0x6f, 0x6f, 0x6b, 0x75, 0x70, 0x5f,
-  0x69, 0x6e, 0x74, 0x65, 0x6e, 0x73, 0x69, 0x74, 0x79, 0x28, 0x66, 0x6c,
-  0x61, 0x67, 0x73, 0x2c, 0x20, 0x68, 0x2c, 0x20, 0x6b, 0x2c, 0x20, 0x6c,
-  0x29, 0x3b, 0x0a, 0x7d, 0x0a, 0x0a, 0x0a, 0x66, 0x6c, 0x6f, 0x61, 0x74,
-  0x20, 0x6d, 0x6f, 0x6c, 0x65, 0x63, 0x75, 0x6c, 0x65, 0x5f, 0x66, 0x61,
-  0x63, 0x74, 0x6f, 0x72, 0x28, 0x67, 0x6c, 0x6f, 0x62, 0x61, 0x6c, 0x20,
-  0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x2a, 0x69, 0x6e, 0x74, 0x65, 0x6e,
-  0x73, 0x69, 0x74, 0x69, 0x65, 0x73, 0x2c, 0x20, 0x67, 0x6c, 0x6f, 0x62,
-  0x61, 0x6c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x2a, 0x66, 0x6c,
-  0x61, 0x67, 0x73, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x31, 0x36, 0x20, 0x63,
-  0x65, 0x6c, 0x6c, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x34, 0x20,
-  0x71, 0x29, 0x0a, 0x7b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20,
-  0x68, 0x66, 0x2c, 0x20, 0x6b, 0x66, 0x2c, 0x20, 0x6c, 0x66, 0x3b, 0x0a,
-  0x09, 0x69, 0x6e, 0x74, 0x20, 0x68, 0x2c, 0x20, 0x6b, 0x2c, 0x20, 0x6c,
-  0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x76, 0x61, 0x6c,
-  0x20, 0x3d, 0x20, 0x30, 0x2e, 0x30, 0x3b, 0x0a, 0x0a, 0x09, 0x23, 0x69,
-  0x66, 0x64, 0x65, 0x66, 0x20, 0x46, 0x4c, 0x41, 0x54, 0x5f, 0x49, 0x4e,
-  0x54, 0x45, 0x4e, 0x53, 0x49, 0x54, 0x49, 0x45, 0x53, 0x0a, 0x09, 0x72,
-  0x65, 0x74, 0x75, 0x72, 0x6e, 0x20, 0x31, 0x30, 0x30, 0x2e, 0x30, 0x3b,
-  0x0a, 0x09, 0x23, 0x65, 0x6c, 0x73, 0x65, 0x0a, 0x0a, 0x09, 0x68, 0x66,
-  0x20, 0x3d, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x30, 0x2a, 0x71,
-  0x2e, 0x78, 0x20, 0x2b, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x31,
-  0x2a, 0x71, 0x2e, 0x79, 0x20, 0x2b, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e,
-  0x73, 0x32, 0x2a, 0x71, 0x2e, 0x7a, 0x3b, 0x20, 0x20, 0x2f, 0x2a, 0x20,
-  0x68, 0x20, 0x2a, 0x2f, 0x0a, 0x09, 0x6b, 0x66, 0x20, 0x3d, 0x20, 0x63,
-  0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x33, 0x2a, 0x71, 0x2e, 0x78, 0x20, 0x2b,
-  0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x34, 0x2a, 0x71, 0x2e, 0x79,
-  0x20, 0x2b, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x35, 0x2a, 0x71,
-  0x2e, 0x7a, 0x3b, 0x20, 0x20, 0x2f, 0x2a, 0x20, 0x6b, 0x20, 0x2a, 0x2f,
-  0x0a, 0x09, 0x6c, 0x66, 0x20, 0x3d, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e,
-  0x73, 0x36, 0x2a, 0x71, 0x2e, 0x78, 0x20, 0x2b, 0x20, 0x63, 0x65, 0x6c,
-  0x6c, 0x2e, 0x73, 0x37, 0x2a, 0x71, 0x2e, 0x79, 0x20, 0x2b, 0x20, 0x63,
-  0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x38, 0x2a, 0x71, 0x2e, 0x7a, 0x3b, 0x20,
-  0x20, 0x2f, 0x2a, 0x20, 0x6c, 0x20, 0x2a, 0x2f, 0x0a, 0x0a, 0x09, 0x68,
-  0x20, 0x3d, 0x20, 0x72, 0x6f, 0x75, 0x6e, 0x64, 0x28, 0x68, 0x66, 0x29,
-  0x3b, 0x0a, 0x09, 0x6b, 0x20, 0x3d, 0x20, 0x72, 0x6f, 0x75, 0x6e, 0x64,
-  0x28, 0x6b, 0x66, 0x29, 0x3b, 0x0a, 0x09, 0x6c, 0x20, 0x3d, 0x20, 0x72,
-  0x6f, 0x75, 0x6e, 0x64, 0x28, 0x6c, 0x66, 0x29, 0x3b, 0x0a, 0x0a, 0x09,
-  0x2f, 0x2a, 0x20, 0x53, 0x79, 0x6d, 0x6d, 0x65, 0x74, 0x72, 0x79, 0x20,
-  0x73, 0x74, 0x75, 0x66, 0x66, 0x20, 0x67, 0x6f, 0x65, 0x73, 0x20, 0x68,
-  0x65, 0x72, 0x65, 0x20, 0x2a, 0x2f, 0x0a, 0x09, 0x49, 0x4e, 0x53, 0x45,
-  0x52, 0x54, 0x5f, 0x48, 0x45, 0x52, 0x45, 0x0a, 0x0a, 0x09, 0x72, 0x65,
-  0x74, 0x75, 0x72, 0x6e, 0x20, 0x76, 0x61, 0x6c, 0x3b, 0x0a, 0x09, 0x23,
-  0x65, 0x6e, 0x64, 0x69, 0x66, 0x20, 0x2f, 0x2a, 0x20, 0x46, 0x4c, 0x41,
-  0x54, 0x5f, 0x49, 0x4e, 0x54, 0x45, 0x4e, 0x53, 0x49, 0x54, 0x49, 0x49,
-  0x45, 0x53, 0x20, 0x2a, 0x2f, 0x0a, 0x7d, 0x0a, 0x0a, 0x0a, 0x6b, 0x65,
-  0x72, 0x6e, 0x65, 0x6c, 0x20, 0x76, 0x6f, 0x69, 0x64, 0x20, 0x64, 0x69,
-  0x66, 0x66, 0x72, 0x61, 0x63, 0x74, 0x69, 0x6f, 0x6e, 0x28, 0x67, 0x6c,
-  0x6f, 0x62, 0x61, 0x6c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x2a,
-  0x64, 0x69, 0x66, 0x66, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20,
-  0x6b, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x77, 0x65, 0x69,
-  0x67, 0x68, 0x74, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x69, 0x6e, 0x74, 0x20, 0x77, 0x2c, 0x20,
-  0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x63, 0x6f, 0x72, 0x6e, 0x65, 0x72,
-  0x5f, 0x78, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x63, 0x6f,
-  0x72, 0x6e, 0x65, 0x72, 0x5f, 0x79, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x66, 0x6c, 0x6f, 0x61,
-  0x74, 0x20, 0x66, 0x73, 0x78, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74,
-  0x20, 0x66, 0x73, 0x79, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20,
-  0x66, 0x73, 0x7a, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x73,
-  0x73, 0x78, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x73, 0x73,
-  0x79, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x73, 0x73, 0x7a,
-  0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x72, 0x65, 0x73, 0x2c,
-  0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x63, 0x6c, 0x65, 0x6e, 0x2c,
-  0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x31, 0x36, 0x20, 0x63, 0x65, 0x6c,
-  0x6c, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x67, 0x6c, 0x6f, 0x62, 0x61, 0x6c, 0x20, 0x66, 0x6c,
-  0x6f, 0x61, 0x74, 0x20, 0x2a, 0x69, 0x6e, 0x74, 0x65, 0x6e, 0x73, 0x69,
-  0x74, 0x69, 0x65, 0x73, 0x2c, 0x20, 0x67, 0x6c, 0x6f, 0x62, 0x61, 0x6c,
-  0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x2a, 0x66, 0x6c, 0x61, 0x67,
-  0x73, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x72, 0x65, 0x61, 0x64, 0x5f, 0x6f, 0x6e, 0x6c, 0x79,
-  0x20, 0x69, 0x6d, 0x61, 0x67, 0x65, 0x32, 0x64, 0x5f, 0x74, 0x20, 0x66,
-  0x75, 0x6e, 0x63, 0x5f, 0x61, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x72, 0x65, 0x61, 0x64, 0x5f,
-  0x6f, 0x6e, 0x6c, 0x79, 0x20, 0x69, 0x6d, 0x61, 0x67, 0x65, 0x32, 0x64,
-  0x5f, 0x74, 0x20, 0x66, 0x75, 0x6e, 0x63, 0x5f, 0x62, 0x2c, 0x0a, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x72,
-  0x65, 0x61, 0x64, 0x5f, 0x6f, 0x6e, 0x6c, 0x79, 0x20, 0x69, 0x6d, 0x61,
-  0x67, 0x65, 0x32, 0x64, 0x5f, 0x74, 0x20, 0x66, 0x75, 0x6e, 0x63, 0x5f,
-  0x63, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x6c, 0x6f, 0x63, 0x61, 0x6c, 0x20, 0x66, 0x6c, 0x6f,
-  0x61, 0x74, 0x20, 0x2a, 0x74, 0x6d, 0x70, 0x29, 0x0a, 0x7b, 0x0a, 0x09,
-  0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x66, 0x73, 0x2c, 0x20, 0x73, 0x73,
-  0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x66, 0x5f, 0x6c,
-  0x61, 0x74, 0x74, 0x69, 0x63, 0x65, 0x2c, 0x20, 0x49, 0x5f, 0x6c, 0x61,
-  0x74, 0x74, 0x69, 0x63, 0x65, 0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61,
-  0x74, 0x20, 0x49, 0x5f, 0x6d, 0x6f, 0x6c, 0x65, 0x63, 0x75, 0x6c, 0x65,
-  0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x34, 0x20, 0x71, 0x3b,
-  0x0a, 0x09, 0x63, 0x6f, 0x6e, 0x73, 0x74, 0x20, 0x69, 0x6e, 0x74, 0x20,
-  0x6c, 0x73, 0x30, 0x20, 0x3d, 0x20, 0x67, 0x65, 0x74, 0x5f, 0x6c, 0x6f,
-  0x63, 0x61, 0x6c, 0x5f, 0x73, 0x69, 0x7a, 0x65, 0x28, 0x30, 0x29, 0x3b,
-  0x0a, 0x09, 0x63, 0x6f, 0x6e, 0x73, 0x74, 0x20, 0x69, 0x6e, 0x74, 0x20,
-  0x6c, 0x73, 0x31, 0x20, 0x3d, 0x20, 0x67, 0x65, 0x74, 0x5f, 0x6c, 0x6f,
-  0x63, 0x61, 0x6c, 0x5f, 0x73, 0x69, 0x7a, 0x65, 0x28, 0x31, 0x29, 0x3b,
-  0x0a, 0x09, 0x63, 0x6f, 0x6e, 0x73, 0x74, 0x20, 0x69, 0x6e, 0x74, 0x20,
-  0x6c, 0x69, 0x30, 0x20, 0x3d, 0x20, 0x67, 0x65, 0x74, 0x5f, 0x6c, 0x6f,
-  0x63, 0x61, 0x6c, 0x5f, 0x69, 0x64, 0x28, 0x30, 0x29, 0x3b, 0x0a, 0x09,
-  0x63, 0x6f, 0x6e, 0x73, 0x74, 0x20, 0x69, 0x6e, 0x74, 0x20, 0x6c, 0x69,
-  0x31, 0x20, 0x3d, 0x20, 0x67, 0x65, 0x74, 0x5f, 0x6c, 0x6f, 0x63, 0x61,
-  0x6c, 0x5f, 0x69, 0x64, 0x28, 0x31, 0x29, 0x3b, 0x0a, 0x09, 0x63, 0x6f,
-  0x6e, 0x73, 0x74, 0x20, 0x69, 0x6e, 0x74, 0x20, 0x6c, 0x73, 0x20, 0x3d,
-  0x20, 0x6c, 0x73, 0x30, 0x20, 0x2a, 0x20, 0x6c, 0x73, 0x31, 0x3b, 0x0a,
-  0x0a, 0x09, 0x2f, 0x2a, 0x20, 0x43, 0x61, 0x6c, 0x63, 0x75, 0x6c, 0x61,
-  0x74, 0x65, 0x20, 0x66, 0x72, 0x61, 0x63, 0x74, 0x69, 0x6f, 0x6e, 0x61,
-  0x6c, 0x20, 0x63, 0x6f, 0x6f, 0x72, 0x64, 0x69, 0x6e, 0x61, 0x74, 0x65,
-  0x73, 0x20, 0x69, 0x6e, 0x20, 0x66, 0x73, 0x2f, 0x73, 0x73, 0x20, 0x2a,
-  0x2f, 0x0a, 0x09, 0x66, 0x73, 0x20, 0x3d, 0x20, 0x63, 0x6f, 0x6e, 0x76,
-  0x65, 0x72, 0x74, 0x5f, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x28, 0x67, 0x65,
-  0x74, 0x5f, 0x67, 0x6c, 0x6f, 0x62, 0x61, 0x6c, 0x5f, 0x69, 0x64, 0x28,
-  0x30, 0x29, 0x29, 0x20, 0x2f, 0x20, 0x63, 0x6f, 0x6e, 0x76, 0x65, 0x72,
-  0x74, 0x5f, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x28, 0x6c, 0x73, 0x30, 0x29,
-  0x3b, 0x0a, 0x09, 0x73, 0x73, 0x20, 0x3d, 0x20, 0x63, 0x6f, 0x6e, 0x76,
-  0x65, 0x72, 0x74, 0x5f, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x28, 0x67, 0x65,
-  0x74, 0x5f, 0x67, 0x6c, 0x6f, 0x62, 0x61, 0x6c, 0x5f, 0x69, 0x64, 0x28,
-  0x31, 0x29, 0x29, 0x20, 0x2f, 0x20, 0x63, 0x6f, 0x6e, 0x76, 0x65, 0x72,
-  0x74, 0x5f, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x28, 0x6c, 0x73, 0x31, 0x29,
-  0x3b, 0x0a, 0x0a, 0x09, 0x2f, 0x2a, 0x20, 0x47, 0x65, 0x74, 0x20, 0x74,
-  0x68, 0x65, 0x20, 0x73, 0x63, 0x61, 0x74, 0x74, 0x65, 0x72, 0x69, 0x6e,
-  0x67, 0x20, 0x76, 0x65, 0x63, 0x74, 0x6f, 0x72, 0x20, 0x2a, 0x2f, 0x0a,
-  0x09, 0x71, 0x20, 0x3d, 0x20, 0x67, 0x65, 0x74, 0x5f, 0x71, 0x28, 0x66,
-  0x73, 0x2c, 0x20, 0x73, 0x73, 0x2c, 0x20, 0x72, 0x65, 0x73, 0x2c, 0x20,
-  0x63, 0x6c, 0x65, 0x6e, 0x2c, 0x20, 0x6b, 0x2c, 0x0a, 0x09, 0x20, 0x20,
-  0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x63, 0x6f, 0x72, 0x6e,
-  0x65, 0x72, 0x5f, 0x78, 0x2c, 0x20, 0x63, 0x6f, 0x72, 0x6e, 0x65, 0x72,
-  0x5f, 0x79, 0x2c, 0x20, 0x66, 0x73, 0x78, 0x2c, 0x20, 0x66, 0x73, 0x79,
-  0x2c, 0x20, 0x66, 0x73, 0x7a, 0x2c, 0x20, 0x73, 0x73, 0x78, 0x2c, 0x20,
-  0x73, 0x73, 0x79, 0x2c, 0x20, 0x73, 0x73, 0x7a, 0x29, 0x3b, 0x0a, 0x0a,
-  0x09, 0x2f, 0x2a, 0x20, 0x43, 0x61, 0x6c, 0x63, 0x75, 0x6c, 0x61, 0x74,
-  0x65, 0x20, 0x74, 0x68, 0x65, 0x20, 0x64, 0x69, 0x66, 0x66, 0x72, 0x61,
-  0x63, 0x74, 0x69, 0x6f, 0x6e, 0x20, 0x2a, 0x2f, 0x0a, 0x09, 0x66, 0x5f,
-  0x6c, 0x61, 0x74, 0x74, 0x69, 0x63, 0x65, 0x20, 0x3d, 0x20, 0x6c, 0x61,
-  0x74, 0x74, 0x69, 0x63, 0x65, 0x5f, 0x66, 0x61, 0x63, 0x74, 0x6f, 0x72,
-  0x28, 0x63, 0x65, 0x6c, 0x6c, 0x2c, 0x20, 0x71, 0x2c, 0x20, 0x66, 0x75,
-  0x6e, 0x63, 0x5f, 0x61, 0x2c, 0x20, 0x66, 0x75, 0x6e, 0x63, 0x5f, 0x62,
-  0x2c, 0x20, 0x66, 0x75, 0x6e, 0x63, 0x5f, 0x63, 0x29, 0x3b, 0x0a, 0x09,
-  0x49, 0x5f, 0x6d, 0x6f, 0x6c, 0x65, 0x63, 0x75, 0x6c, 0x65, 0x20, 0x3d,
-  0x20, 0x6d, 0x6f, 0x6c, 0x65, 0x63, 0x75, 0x6c, 0x65, 0x5f, 0x66, 0x61,
-  0x63, 0x74, 0x6f, 0x72, 0x28, 0x69, 0x6e, 0x74, 0x65, 0x6e, 0x73, 0x69,
-  0x74, 0x69, 0x65, 0x73, 0x2c, 0x20, 0x66, 0x6c, 0x61, 0x67, 0x73, 0x2c,
-  0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2c, 0x20, 0x71, 0x29, 0x3b, 0x0a, 0x09,
-  0x49, 0x5f, 0x6c, 0x61, 0x74, 0x74, 0x69, 0x63, 0x65, 0x20, 0x3d, 0x20,
-  0x70, 0x6f, 0x77, 0x28, 0x66, 0x5f, 0x6c, 0x61, 0x74, 0x74, 0x69, 0x63,
-  0x65, 0x2c, 0x20, 0x32, 0x2e, 0x30, 0x66, 0x29, 0x3b, 0x0a, 0x0a, 0x09,
-  0x74, 0x6d, 0x70, 0x5b, 0x6c, 0x69, 0x30, 0x20, 0x2b, 0x20, 0x6c, 0x73,
-  0x30, 0x2a, 0x6c, 0x69, 0x31, 0x5d, 0x20, 0x3d, 0x20, 0x49, 0x5f, 0x6d,
-  0x6f, 0x6c, 0x65, 0x63, 0x75, 0x6c, 0x65, 0x20, 0x2a, 0x20, 0x49, 0x5f,
-  0x6c, 0x61, 0x74, 0x74, 0x69, 0x63, 0x65, 0x3b, 0x0a, 0x0a, 0x09, 0x62,
-  0x61, 0x72, 0x72, 0x69, 0x65, 0x72, 0x28, 0x43, 0x4c, 0x4b, 0x5f, 0x4c,
-  0x4f, 0x43, 0x41, 0x4c, 0x5f, 0x4d, 0x45, 0x4d, 0x5f, 0x46, 0x45, 0x4e,
-  0x43, 0x45, 0x29, 0x3b, 0x0a, 0x0a, 0x09, 0x2f, 0x2a, 0x20, 0x46, 0x69,
-  0x72, 0x73, 0x74, 0x20, 0x74, 0x68, 0x72, 0x65, 0x61, 0x64, 0x20, 0x69,
-  0x6e, 0x20, 0x67, 0x72, 0x6f, 0x75, 0x70, 0x20, 0x73, 0x75, 0x6d, 0x73,
-  0x20, 0x74, 0x68, 0x65, 0x20, 0x73, 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x73,
-  0x20, 0x2a, 0x2f, 0x0a, 0x09, 0x69, 0x66, 0x20, 0x28, 0x20, 0x6c, 0x69,
-  0x30, 0x20, 0x2b, 0x20, 0x6c, 0x69, 0x31, 0x20, 0x3d, 0x3d, 0x20, 0x30,
-  0x20, 0x29, 0x20, 0x7b, 0x0a, 0x0a, 0x09, 0x09, 0x69, 0x6e, 0x74, 0x20,
-  0x69, 0x3b, 0x0a, 0x09, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x73,
-  0x75, 0x6d, 0x20, 0x3d, 0x20, 0x30, 0x2e, 0x30, 0x3b, 0x0a, 0x09, 0x09,
-  0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x76, 0x61, 0x6c, 0x3b, 0x0a, 0x09,
-  0x09, 0x69, 0x6e, 0x74, 0x20, 0x69, 0x64, 0x78, 0x3b, 0x0a, 0x0a, 0x09,
-  0x09, 0x69, 0x64, 0x78, 0x20, 0x3d, 0x20, 0x63, 0x6f, 0x6e, 0x76, 0x65,
-  0x72, 0x74, 0x5f, 0x69, 0x6e, 0x74, 0x5f, 0x72, 0x74, 0x7a, 0x28, 0x66,
-  0x73, 0x29, 0x20, 0x2b, 0x20, 0x77, 0x2a, 0x63, 0x6f, 0x6e, 0x76, 0x65,
-  0x72, 0x74, 0x5f, 0x69, 0x6e, 0x74, 0x5f, 0x72, 0x74, 0x7a, 0x28, 0x73,
-  0x73, 0x29, 0x3b, 0x0a, 0x0a, 0x09, 0x09, 0x66, 0x6f, 0x72, 0x20, 0x28,
-  0x20, 0x69, 0x3d, 0x30, 0x3b, 0x20, 0x69, 0x3c, 0x6c, 0x73, 0x3b, 0x20,
-  0x69, 0x2b, 0x2b, 0x20, 0x29, 0x20, 0x73, 0x75, 0x6d, 0x20, 0x2b, 0x3d,
-  0x20, 0x74, 0x6d, 0x70, 0x5b, 0x69, 0x5d, 0x3b, 0x0a, 0x0a, 0x09, 0x09,
-  0x76, 0x61, 0x6c, 0x20, 0x3d, 0x20, 0x77, 0x65, 0x69, 0x67, 0x68, 0x74,
-  0x20, 0x2a, 0x20, 0x73, 0x75, 0x6d, 0x20, 0x2f, 0x20, 0x63, 0x6f, 0x6e,
-  0x76, 0x65, 0x72, 0x74, 0x5f, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x28, 0x6c,
-  0x73, 0x29, 0x3b, 0x0a, 0x09, 0x09, 0x64, 0x69, 0x66, 0x66, 0x5b, 0x69,
-  0x64, 0x78, 0x5d, 0x20, 0x3d, 0x20, 0x76, 0x61, 0x6c, 0x3b, 0x0a, 0x0a,
-  0x09, 0x7d, 0x0a, 0x0a, 0x7d, 0x0a
-};
-unsigned int data_diffraction_cl_len = 6006;
diff --git a/src/diffraction.h b/src/diffraction.h
deleted file mode 100644
index 63b022c9..00000000
--- a/src/diffraction.h
+++ /dev/null
@@ -1,66 +0,0 @@
-/*
- * diffraction.h
- *
- * Calculate diffraction patterns by Fourier methods
- *
- * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
- *                       a research centre of the Helmholtz Association.
- *
- * Authors:
- *   2009-2019 Thomas White <taw@physics.org>
- *   2013-2014 Chun Hong Yoon <chun.hong.yoon@desy.de>
- *   2013      Alexandra Tolstikova
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#ifndef DIFFRACTION_H
-#define DIFFRACTION_H
-
-#include <gsl/gsl_rng.h>
-
-#include "image.h"
-#include "cell.h"
-#include "symmetry.h"
-
-
-typedef enum {
-	GRADIENT_MOSAIC,
-	GRADIENT_INTERPOLATE,
-	GRADIENT_PHASED
-} GradientMethod;
-
-extern void get_diffraction(struct image *image, int na, int nb, int nc,
-                            const double *intensities, const double *phases,
-                            const unsigned char *flags, UnitCell *cell,
-                            GradientMethod m, const SymOpList *sym,
-                            int no_fringes, int flat, int n_samples);
-
-extern struct sample *generate_tophat(struct image *image);
-
-extern struct sample *generate_SASE(struct image *image, gsl_rng *rng);
-
-extern struct sample *generate_twocolour(struct image *image);
-
-extern struct sample *generate_spectrum_fromfile(struct image *image,
-					         char *spectrum_fn);
-
-#endif	/* DIFFRACTION_H */
diff --git a/src/list_tmp.h b/src/list_tmp.h
deleted file mode 100644
index 0c1a84e4..00000000
--- a/src/list_tmp.h
+++ /dev/null
@@ -1,104 +0,0 @@
-/*
- * Template for creating indexed 3D lists of a given type, usually indexed
- * as signed h,k,l values where -INDMAX<={h,k,l}<=+INDMAX.
- *
- * These are used, for example, for:
- *  - a list of 'double complex' values for storing structure factors,
- *  - a list of 'double' values for storing reflection intensities,
- *  - a list of 'unsigned int' values for counts of some sort.
- *
- * When LABEL and TYPE are #defined appropriately, including this header
- * defines functions such as:
- *  - new_list_<LABEL>(), for creating a new list of the given type,
- *  - set_<LABEL>(), for setting a value in a list,
- *  - lookup_<LABEL>(), for retrieving values from a list,
- * ..and so on.
- *
- * See src/utils.h for more information.
- *
- */
-
-/*
- * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
- *                       a research centre of the Helmholtz Association.
- *
- * Authors:
- *   2009-2012 Thomas White <taw@physics.org>
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-#include <stdio.h>
-
-#define ERROR_T(...) fprintf(stderr, __VA_ARGS__)
-
-
-static inline void LABEL(set_arr)(TYPE *ref, signed int h,
-                              signed int k, signed int l,
-                              TYPE i)
-{
-	int idx;
-
-	if ( (abs(h) > INDMAX) || (abs(k) > INDMAX) || (abs(l) > INDMAX) ) {
-		ERROR_T("\nReflection %i %i %i is out of range!\n", h, k, l);
-		ERROR_T("You need to re-configure INDMAX and re-run.\n");
-		exit(1);
-	}
-
-	if ( h < 0 ) h += IDIM;
-	if ( k < 0 ) k += IDIM;
-	if ( l < 0 ) l += IDIM;
-
-	idx = h + (IDIM*k) + (IDIM*IDIM*l);
-	ref[idx] = i;
-}
-
-
-static inline TYPE LABEL(lookup_arr)(const TYPE *ref, signed int h,
-                                 signed int k, signed int l)
-{
-	int idx;
-
-	if ( (abs(h) > INDMAX) || (abs(k) > INDMAX) || (abs(l) > INDMAX) ) {
-		ERROR_T("\nReflection %i %i %i is out of range!\n", h, k, l);
-		ERROR_T("You need to re-configure INDMAX and re-run.\n");
-		exit(1);
-	}
-
-	if ( h < 0 ) h += IDIM;
-	if ( k < 0 ) k += IDIM;
-	if ( l < 0 ) l += IDIM;
-
-	idx = h + (IDIM*k) + (IDIM*IDIM*l);
-	return ref[idx];
-}
-
-
-static inline TYPE *LABEL(new_arr)(void)
-{
-	TYPE *r;
-	size_t r_size;
-	r_size = IDIM*IDIM*IDIM*sizeof(TYPE);
-	r = malloc(r_size);
-	memset(r, 0, r_size);
-	return r;
-}
-
-
-#undef LABEL
-#undef TYPE
-#undef ERROR_T
diff --git a/src/partial_sim.c b/src/partial_sim.c
deleted file mode 100644
index f5588ceb..00000000
--- a/src/partial_sim.c
+++ /dev/null
@@ -1,1034 +0,0 @@
-/*
- * partial_sim.c
- *
- * Generate partials for testing scaling
- *
- * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
- *                       a research centre of the Helmholtz Association.
- *
- * Authors:
- *   2011-2020 Thomas White <taw@physics.org>
- *   2014      Valerio Mariani
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <stdarg.h>
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-#include <unistd.h>
-#include <getopt.h>
-#include <assert.h>
-#include <pthread.h>
-#include <gsl/gsl_rng.h>
-
-#include <image.h>
-#include <utils.h>
-#include <reflist-utils.h>
-#include <symmetry.h>
-#include <geometry.h>
-#include <stream.h>
-#include <thread-pool.h>
-#include <cell-utils.h>
-
-#include "version.h"
-
-
-/* Number of bins for partiality graph */
-#define NBINS 50
-
-
-static void mess_up_cell(Crystal *cr, double cnoise, gsl_rng *rng)
-{
-	double ax, ay, az;
-	double bx, by, bz;
-	double cx, cy, cz;
-	UnitCell *cell = crystal_get_cell(cr);
-
-	//STATUS("Real:\n");
-	//cell_print(cell);
-
-	cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
-	ax = flat_noise(rng, ax, cnoise*fabs(ax)/100.0);
-	ay = flat_noise(rng, ay, cnoise*fabs(ay)/100.0);
-	az = flat_noise(rng, az, cnoise*fabs(az)/100.0);
-	bx = flat_noise(rng, bx, cnoise*fabs(bx)/100.0);
-	by = flat_noise(rng, by, cnoise*fabs(by)/100.0);
-	bz = flat_noise(rng, bz, cnoise*fabs(bz)/100.0);
-	cx = flat_noise(rng, cx, cnoise*fabs(cx)/100.0);
-	cy = flat_noise(rng, cy, cnoise*fabs(cy)/100.0);
-	cz = flat_noise(rng, cz, cnoise*fabs(cz)/100.0);
-	cell_set_reciprocal(cell, ax, ay, az, bx, by, bz, cx, cy, cz);
-
-	//STATUS("Changed:\n");
-	//cell_print(cell);
-}
-
-
-/* For each reflection in "partial", fill in what the intensity would be
- * according to "full" */
-static void calculate_partials(Crystal *cr,
-                               RefList *full, const SymOpList *sym,
-                               int random_intensities,
-                               pthread_rwlock_t *full_lock,
-                               unsigned long int *n_ref, double *p_hist,
-                               double *p_max, double max_q, double full_stddev,
-                               double noise_stddev, gsl_rng *rng,
-                               UnitCell *template_cell, RefList *template_reflist)
-{
-	Reflection *refl;
-	RefListIterator *iter;
-
-	for ( refl = first_refl(crystal_get_reflections(cr), &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) )
-	{
-		signed int h, k, l;
-		Reflection *rfull;
-		double L, p, Ip, If;
-		int bin;
-		double res;
-
-		get_indices(refl, &h, &k, &l);
-		get_asymm(sym, h, k, l, &h, &k, &l);
-		p = get_partiality(refl);
-		L = get_lorentz(refl);
-
-		pthread_rwlock_rdlock(full_lock);
-		rfull = find_refl(full, h, k, l);
-		pthread_rwlock_unlock(full_lock);
-
-		if ( rfull == NULL ) {
-			if ( random_intensities ) {
-
-				pthread_rwlock_wrlock(full_lock);
-
-				/* In the gap between the unlock and the wrlock,
-				 * the reflection might have been created by
-				 * another thread.  So, we must check again */
-				rfull = find_refl(full, h, k, l);
-				if ( rfull == NULL ) {
-					rfull = add_refl(full, h, k, l);
-					If = fabs(gaussian_noise(rng, 0.0,
-					                         full_stddev));
-					set_intensity(rfull, If);
-					set_redundancy(rfull, 1);
-				} else {
-					If = get_intensity(rfull);
-				}
-				pthread_rwlock_unlock(full_lock);
-
-			} else {
-				set_redundancy(refl, 0);
-				If = 0.0;
-			}
-		} else {
-			If = get_intensity(rfull);
-			if ( random_intensities ) {
-				lock_reflection(rfull);
-				int red = get_redundancy(rfull);
-				set_redundancy(rfull, red+1);
-				unlock_reflection(rfull);
-			}
-		}
-
-		Ip = crystal_get_osf(cr) * L * p * If;
-
-		res = resolution(crystal_get_cell(cr), h, k, l);
-		bin = NBINS*2.0*res/max_q;
-		if ( (bin < NBINS) && (bin>=0) ) {
-			p_hist[bin] += p;
-			n_ref[bin]++;
-			if ( p > p_max[bin] ) p_max[bin] = p;
-		} else {
-			STATUS("Reflection out of histogram range: %e %i %f\n",
-			       res, bin,  p);
-		}
-
-		Ip = gaussian_noise(rng, Ip, noise_stddev);
-
-		set_intensity(refl, Ip);
-		set_esd_intensity(refl, noise_stddev);
-	}
-}
-
-
-static void draw_and_write_image(struct image *image,
-                                 const DataTemplate *dtempl,
-                                 RefList *reflections,
-                                 gsl_rng *rng, double background)
-{
-	Reflection *refl;
-	RefListIterator *iter;
-	int i;
-
-	image->dp = malloc(image->detgeom->n_panels*sizeof(float *));
-	if ( image->dp == NULL ) {
-		ERROR("Failed to allocate data\n");
-		return;
-	}
-	for ( i=0; i<image->detgeom->n_panels; i++ ) {
-
-		int j;
-		struct detgeom_panel *p = &image->detgeom->panels[i];
-
-		image->dp[i] = calloc(p->w * p->h, sizeof(float));
-		if ( image->dp[i] == NULL ) {
-			ERROR("Failed to allocate data\n");
-			return;
-		}
-		for ( j=0; j<p->w*p->h; j++ ) {
-			image->dp[i][j] = poisson_noise(rng, background);
-		}
-
-	}
-
-	for ( refl = first_refl(reflections, &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) )
-	{
-		double Ip;
-		double dfs, dss;
-		int fs, ss;
-		struct detgeom_panel *p;
-		signed int pn;
-
-		Ip = get_intensity(refl);
-
-		get_detector_pos(refl, &dfs, &dss);
-		pn = get_panel_number(refl);
-		assert(pn < image->detgeom->n_panels);
-		p = &image->detgeom->panels[pn];
-
-		/* Explicit rounding, downwards */
-		fs = dfs;  ss = dss;
-		assert(fs >= 0);
-		assert(ss >= 0);
-		assert(fs < p->w);
-		assert(ss < p->h);
-
-		image->dp[pn][fs + p->w*ss] += Ip;
-	}
-
-	image_write(image, dtempl, image->filename);
-
-	for ( i=0; i<image->detgeom->n_panels; i++ ) {
-		free(image->dp[i]);
-	}
-	free(image->dp);
-}
-
-
-static void show_help(const char *s)
-{
-	printf("Syntax: %s [options]\n\n", s);
-	printf(
-"Generate a stream containing partials from a reflection list.\n"
-"\n"
-" -h, --help               Display this help message.\n"
-"     --version            Print CrystFEL version number and exit.\n"
-"\n"
-"You need to provide the following basic options:\n"
-" -i, --input=<file>       Read reflections from <file>.\n"
-"                           Default: generate random ones instead (see -r).\n"
-" -o, --output=<file>      Write partials in stream format to <file>.\n"
-"     --images=<prefix>    Write images to <prefix>NNN.h5.\n"
-" -g. --geometry=<file>    Get detector geometry from file.\n"
-" -p, --pdb=<file>         PDB file from which to get the unit cell.\n"
-"\n"
-" -y, --symmetry=<sym>     Symmetry of the input reflection list.\n"
-" -n <n>                   Simulate <n> patterns.  Default: 2.\n"
-" -r, --save-random=<file> Save randomly generated intensities to file.\n"
-"     --pgraph=<file>      Save a histogram of partiality values to file.\n"
-" -c, --cnoise=<val>       Amount of reciprocal space cell noise, in percent.\n"
-"     --osf-stddev=<val>   Standard deviation of the scaling factors.\n"
-"     --full-stddev=<val>  Standard deviation of the randomly\n"
-"                           generated full intensities, if not using -i.\n"
-"     --noise-stddev=<val> Set the standard deviation of the noise.\n"
-"     --background=<val>   Background level in photons.  Default 3000.\n"
-"     --profile-radius     Reciprocal space reflection profile radius in m^-1.\n"
-"                           Default 0.001e9 m^-1\n"
-"     --really-random      Be non-deterministic.\n"
-"\n"
-);
-}
-
-
-struct partial_sim_queue_args
-{
-	RefList *full;
-	pthread_rwlock_t full_lock;
-	const DataTemplate *dtempl;
-
-	int n_done;
-	int n_started;
-	int n_to_do;
-
-	SymOpList *sym;
-	int random_intensities;
-	UnitCell *cell;
-	double cnoise;
-	double osf_stddev;
-	double full_stddev;
-	double noise_stddev;
-	double background;
-	double profile_radius;
-
-	double max_q;
-
-	char *image_prefix;
-
-	/* The overall histogram */
-	double p_hist[NBINS];
-	unsigned long int n_ref[NBINS];
-	double p_max[NBINS];
-
-	Stream *stream;
-	gsl_rng **rngs;
-	Stream *template_stream;
-};
-
-
-struct partial_sim_worker_args
-{
-	struct partial_sim_queue_args *qargs;
-	struct image *image;
-
-	UnitCell *template_cell;
-	RefList *template_reflist;
-
-	/* Histogram for this image */
-	double p_hist[NBINS];
-	unsigned long int n_ref[NBINS];
-	double p_max[NBINS];
-
-	int n;
-};
-
-
-static void *create_job(void *vqargs)
-{
-	struct partial_sim_worker_args *wargs;
-	struct partial_sim_queue_args *qargs = vqargs;
-
-	/* All done already? */
-	if ( qargs->n_started == qargs->n_to_do ) return NULL;
-
-	wargs = malloc(sizeof(struct partial_sim_worker_args));
-
-	wargs->qargs = qargs;
-
-	if ( qargs->template_stream != NULL ) {
-
-		struct image *image;
-
-		image = stream_read_chunk(qargs->template_stream,
-		                          STREAM_REFLECTIONS);
-		if ( image == NULL ) {
-			ERROR("Failed to read template chunk!\n");
-			free(wargs);
-			return NULL;
-		}
-		if ( image->n_crystals != 1 ) {
-			ERROR("Template stream must have exactly one crystal "
-			      "per frame.\n");
-			free(wargs);
-			return NULL;
-		}
-
-		wargs->template_cell = crystal_get_cell(image->crystals[0]);
-		wargs->template_reflist = crystal_get_reflections(image->crystals[0]);
-
-		image_free(image);
-
-	} else {
-		wargs->template_cell = NULL;
-		wargs->template_reflist = NULL;
-	}
-
-	qargs->n_started++;
-	wargs->n = qargs->n_started;
-
-	return wargs;
-}
-
-
-static void run_job(void *vwargs, int cookie)
-{
-	struct partial_sim_worker_args *wargs = vwargs;
-	struct partial_sim_queue_args *qargs = wargs->qargs;
-	int i;
-	Crystal *cr;
-	double osf;
-	struct image *image;
-
-	image = image_create_for_simulation(qargs->dtempl);
-	if ( image == NULL ) {
-		ERROR("Failed to create image.\n");
-		return;
-	}
-
-	cr = crystal_new();
-	if ( cr == NULL ) {
-		ERROR("Failed to create crystal.\n");
-		return;
-	}
-	crystal_set_image(cr, image);
-	image_add_crystal(image, cr);
-
-	do {
-		osf = gaussian_noise(qargs->rngs[cookie], 1.0,
-		                     qargs->osf_stddev);
-	} while ( osf <= 0.0 );
-	crystal_set_osf(cr, osf);
-	crystal_set_mosaicity(cr, 0.0);
-	crystal_set_profile_radius(cr, qargs->profile_radius);
-
-	if ( wargs->template_cell == NULL ) {
-		/* Set up a random orientation */
-		struct quaternion orientation;
-		char tmp[128];
-		orientation = random_quaternion(qargs->rngs[cookie]);
-		crystal_set_cell(cr, cell_rotate(qargs->cell, orientation));
-		snprintf(tmp, 127, "quaternion = %f %f %f %f",
-		         orientation.w, orientation.x, orientation.y, orientation.z);
-		crystal_add_notes(cr, tmp);
-	} else {
-		crystal_set_cell(cr, wargs->template_cell);
-	}
-
-	image->filename = malloc(256);
-	if ( image->filename == NULL ) {
-		ERROR("Failed to allocate filename.\n");
-		return;
-	}
-	if ( qargs->image_prefix != NULL ) {
-		snprintf(image->filename, 255, "%s%i.h5",
-		         qargs->image_prefix, wargs->n);
-	} else {
-		snprintf(image->filename, 255, "dummy-%i.h5", wargs->n);
-	}
-
-	if ( wargs->template_reflist == NULL ) {
-		RefList *reflections;
-		reflections = predict_to_res(cr, qargs->max_q);
-		crystal_set_reflections(cr, reflections);
-		calculate_partialities(cr, PMODEL_XSPHERE);
-	} else {
-		crystal_set_reflections(cr, wargs->template_reflist);
-		update_predictions(cr);
-		calculate_partialities(cr, PMODEL_XSPHERE);
-	}
-
-	for ( i=0; i<NBINS; i++ ) {
-		wargs->n_ref[i] = 0;
-		wargs->p_hist[i] = 0.0;
-		wargs->p_max[i] = 0.0;
-	}
-
-	calculate_partials(cr, qargs->full,
-	                   qargs->sym, qargs->random_intensities,
-	                   &qargs->full_lock,
-	                   wargs->n_ref, wargs->p_hist, wargs->p_max,
-	                   qargs->max_q, qargs->full_stddev,
-	                   qargs->noise_stddev, qargs->rngs[cookie],
-	                   wargs->template_cell, wargs->template_reflist);
-
-	if ( qargs->image_prefix != NULL ) {
-		draw_and_write_image(image, qargs->dtempl,
-		                     crystal_get_reflections(cr),
-		                     qargs->rngs[cookie],
-		                     qargs->background);
-	}
-
-	/* Give a slightly incorrect cell in the stream */
-	mess_up_cell(cr, qargs->cnoise, qargs->rngs[cookie]);
-
-	wargs->image = image;
-}
-
-
-static void finalise_job(void *vqargs, void *vwargs)
-{
-	struct partial_sim_worker_args *wargs = vwargs;
-	struct partial_sim_queue_args *qargs = vqargs;
-	int i;
-	int ret;
-
-	ret = stream_write_chunk(qargs->stream, wargs->image,
-	                         STREAM_REFLECTIONS);
-	if ( ret != 0 ) {
-		ERROR("WARNING: error writing stream file.\n");
-	}
-
-	for ( i=0; i<NBINS; i++ ) {
-		qargs->n_ref[i] += wargs->n_ref[i];
-		qargs->p_hist[i] += wargs->p_hist[i];
-		if ( wargs->p_max[i] > qargs->p_max[i] ) {
-			qargs->p_max[i] = wargs->p_max[i];
-		}
-	}
-
-	qargs->n_done++;
-	progress_bar(qargs->n_done, qargs->n_to_do, "Simulating");
-
-	image_free(wargs->image);
-}
-
-
-int main(int argc, char *argv[])
-{
-	int c;
-	char *input_file = NULL;
-	char *output_file = NULL;
-	char *geomfile = NULL;
-	char *cellfile = NULL;
-	DataTemplate *dtempl;
-	RefList *full = NULL;
-	char *sym_str = NULL;
-	SymOpList *sym;
-	UnitCell *cell = NULL;
-	Stream *stream;
-	int n = 2;
-	int random_intensities = 0;
-	char *save_file = NULL;
-	struct partial_sim_queue_args qargs;
-	int n_threads = 1;
-	char *rval;
-	int i;
-	char *phist_file = NULL;
-	int config_random = 0;
-	char *image_prefix = NULL;
-	Stream *template_stream = NULL;
-	char *template = NULL;
-	struct image *test_image;
-
-	/* Default simulation parameters */
-	double profile_radius = 0.001e9;
-	double osf_stddev = 2.0;
-	double full_stddev = 1000.0;
-	double noise_stddev = 20.0;
-	double background = 3000.0;
-	double cnoise = 0.0;
-
-	/* Long options */
-	const struct option longopts[] = {
-		{"help",               0, NULL,               'h'},
-		{"version",            0, NULL,               'v'},
-		{"beam",               1, NULL,               'b'},
-		{"output",             1, NULL,               'o'},
-		{"input",              1, NULL,               'i'},
-		{"pdb",                1, NULL,               'p'},
-		{"geometry",           1, NULL,               'g'},
-		{"symmetry",           1, NULL,               'y'},
-		{"save-random",        1, NULL,               'r'},
-		{"cnoise",             1, NULL,               'c'},
-
-		{"pgraph",             1, NULL,                2},
-		{"osf-stddev",         1, NULL,                3},
-		{"full-stddev",        1, NULL,                4},
-		{"noise-stddev",       1, NULL,                5},
-		{"images",             1, NULL,                6},
-		{"background",         1, NULL,                7},
-		{"beam-divergence",    1, NULL,                8},
-		{"beam-bandwidth",     1, NULL,                9},
-		{"profile-radius",     1, NULL,               10},
-		{"photon-energy",      1, NULL,               11},
-		{"template-stream",    1, NULL,               12},
-
-		{"really-random",      0, &config_random,      1},
-
-		{0, 0, NULL, 0}
-	};
-
-	/* Short options */
-	while ((c = getopt_long(argc, argv, "hi:o:p:g:y:n:r:j:c:vb:",
-	                        longopts, NULL)) != -1)
-	{
-		switch (c) {
-
-			case 'h' :
-			show_help(argv[0]);
-			return 0;
-
-			case 'v' :
-			printf("CrystFEL: %s\n",
-			       crystfel_version_string());
-			printf("%s\n",
-			       crystfel_licence_string());
-			return 0;
-
-			case 'b' :
-			ERROR("WARNING: This version of CrystFEL no longer "
-			      "uses beam files.  Please remove the beam file "
-			      "from your partial_sim command line.\n");
-			return 1;
-
-			case 'i' :
-			input_file = strdup(optarg);
-			break;
-
-			case 'o' :
-			output_file = strdup(optarg);
-			break;
-
-			case 'p' :
-			cellfile = strdup(optarg);
-			break;
-
-			case 'g' :
-			geomfile = strdup(optarg);
-			break;
-
-			case 'y' :
-			sym_str = strdup(optarg);
-			break;
-
-			case 'n' :
-			n = atoi(optarg);
-			break;
-
-			case 'r' :
-			save_file = strdup(optarg);
-			break;
-
-			case 'j' :
-			n_threads = atoi(optarg);
-			break;
-
-			case 'c' :
-			cnoise = strtod(optarg, &rval);
-			if ( *rval != '\0' ) {
-				ERROR("Invalid cell noise value.\n");
-				return 1;
-			}
-			break;
-
-			case 2 :
-			phist_file = strdup(optarg);
-			break;
-
-			case 3 :
-			osf_stddev = strtod(optarg, &rval);
-			if ( *rval != '\0' ) {
-				ERROR("Invalid OSF standard deviation.\n");
-				return 1;
-			}
-			if ( osf_stddev < 0.0 ) {
-				ERROR("Invalid OSF standard deviation.");
-				ERROR(" (must be positive).\n");
-				return 1;
-			}
-			break;
-
-			case 4 :
-			full_stddev = strtod(optarg, &rval);
-			if ( *rval != '\0' ) {
-				ERROR("Invalid full standard deviation.\n");
-				return 1;
-			}
-			if ( full_stddev < 0.0 ) {
-				ERROR("Invalid full standard deviation.");
-				ERROR(" (must be positive).\n");
-				return 1;
-			}
-			break;
-
-			case 5 :
-			noise_stddev = strtod(optarg, &rval);
-			if ( *rval != '\0' ) {
-				ERROR("Invalid noise standard deviation.\n");
-				return 1;
-			}
-			if ( noise_stddev < 0.0 ) {
-				ERROR("Invalid noise standard deviation.");
-				ERROR(" (must be positive).\n");
-				return 1;
-			}
-			break;
-
-			case 6 :
-			image_prefix = strdup(optarg);
-			break;
-
-			case 7 :
-			background = strtod(optarg, &rval);
-			if ( *rval != '\0' ) {
-				ERROR("Invalid background level.\n");
-				return 1;
-			}
-			if ( background < 0.0 ) {
-				ERROR("Background level must be positive.\n");
-				return 1;
-			}
-			break;
-
-			case 8 :
-			ERROR("--beam-divergence is no longer used.\n");
-			ERROR("The 'xsphere' partiality model does not take divergence into account.\n");
-			return 1;
-
-			case 9 :
-			ERROR("--beam-bandwidth is no longer used.\n");
-			ERROR("Set the bandwidth in the geometry file instead.\n");
-			return 1;
-
-			case 10 :
-			profile_radius = strtod(optarg, &rval);
-			if ( *rval != '\0' ) {
-				ERROR("Invalid profile radius.\n");
-				return 1;
-			}
-			if ( profile_radius < 0.0 ) {
-				ERROR("Profile radius must be positive.\n");
-				return 1;
-			}
-			break;
-
-			case 11 :
-			ERROR("--photon-energy is no longer used.\n");
-			ERROR("Set the photon energy in the geometry file instead.\n");
-			return 1;
-
-			case 12 :
-			template = strdup(optarg);
-			break;
-
-			case 0 :
-			break;
-
-			case '?' :
-			return 1;
-
-			default :
-			ERROR("Unhandled option '%c'\n", c);
-			break;
-
-		}
-	}
-
-	if ( n_threads < 1 ) {
-		ERROR("Invalid number of threads.\n");
-		return 1;
-	}
-
-	if ( (n_threads > 1) && (image_prefix != NULL) ) {
-		ERROR("Option \"--images\" is incompatible with \"-j\".\n");
-		return 1;
-	}
-
-	/* Load cell */
-	if ( cellfile == NULL ) {
-		ERROR("You need to give a PDB file with the unit cell.\n");
-		return 1;
-	}
-	cell = load_cell_from_file(cellfile);
-	if ( cell == NULL ) {
-		ERROR("Failed to get cell from '%s'\n", cellfile);
-		return 1;
-	}
-	free(cellfile);
-
-	if ( !cell_is_sensible(cell) ) {
-		ERROR("Invalid unit cell parameters:\n");
-		cell_print(cell);
-		return 1;
-	}
-
-	/* Load geometry */
-	if ( geomfile == NULL ) {
-		ERROR("You need to give a geometry file.\n");
-		return 1;
-	}
-	dtempl = data_template_new_from_file(geomfile);
-	if ( dtempl == NULL ) {
-		ERROR("Failed to read geometry from '%s'\n", geomfile);
-		return 1;
-	}
-
-	if ( save_file == NULL ) save_file = strdup("partial_sim.hkl");
-
-	/* Load (full) reflections */
-	if ( input_file != NULL ) {
-
-		RefList *as;
-		char *sym_str_fromfile = NULL;
-
-		full = read_reflections_2(input_file, &sym_str_fromfile);
-		if ( full == NULL ) {
-			ERROR("Failed to read reflections from '%s'\n",
-			      input_file);
-			return 1;
-		}
-
-		/* If we don't have a point group yet, and if the file provides
-		 * one, use the one from the file */
-		if ( (sym_str == NULL) && (sym_str_fromfile != NULL) ) {
-			sym_str = sym_str_fromfile;
-			STATUS("Using symmetry from reflection file: %s\n",
-			       sym_str);
-		}
-
-		/* If we still don't have a point group, use "1" */
-		if ( sym_str == NULL ) sym_str = strdup("1");
-
-		pointgroup_warning(sym_str);
-		sym = get_pointgroup(sym_str);
-
-		if ( check_list_symmetry(full, sym) ) {
-			ERROR("The input reflection list does not appear to"
-			      " have symmetry %s\n", symmetry_name(sym));
-			if ( cell_get_lattice_type(cell) == L_MONOCLINIC ) {
-				ERROR("You may need to specify the unique axis "
-				      "in your point group.  The default is "
-				      "unique axis c.\n");
-				ERROR("See 'man crystfel' for more details.\n");
-			}
-			return 1;
-		}
-
-		as = asymmetric_indices(full, sym);
-		reflist_free(full);
-		full = as;
-
-	} else {
-		random_intensities = 1;
-		if ( sym_str == NULL ) sym_str = strdup("1");
-		sym = get_pointgroup(sym_str);
-	}
-
-	if ( n < 1 ) {
-		ERROR("Number of patterns must be at least 1.\n");
-		return 1;
-	}
-
-	if ( output_file == NULL ) {
-		ERROR("You must give a filename for the output.\n");
-		return 1;
-	}
-	stream = stream_open_for_write(output_file, dtempl);
-	if ( stream == NULL ) {
-		ERROR("Couldn't open output file '%s'\n", output_file);
-		return 1;
-	}
-	free(output_file);
-
-	stream_write_geometry_file(stream, geomfile);
-	stream_write_commandline_args(stream, argc, argv);
-
-	if ( template != NULL ) {
-		template_stream = stream_open_for_read(template);
-		if ( template_stream == NULL ) {
-			ERROR("Couldn't open template stream '%s'\n", template);
-			return 1;
-		}
-	}
-
-	test_image = image_create_for_simulation(dtempl);
-	if ( test_image == NULL ) {
-		ERROR("Could not create image structure.\n");
-		ERROR("Does the geometry file contain references?\n");
-		return 1;
-	}
-
-	STATUS("Simulation parameters:\n");
-	STATUS("                     Wavelength: %.5f A (photon energy %.2f eV)\n",
-	       test_image->lambda*1e10,
-	       ph_lambda_to_eV(test_image->lambda));
-	STATUS("                Beam divergence: 0 (not modelled)\n");
-	STATUS("                 Beam bandwidth: %.5f %%\n",
-	       test_image->bw*100.0);
-	STATUS("Reciprocal space profile radius: %e m^-1\n",
-	       profile_radius);
-	if ( image_prefix != NULL ) {
-		STATUS("                     Background: %.2f detector units\n",
-		       background);
-	} else {
-		STATUS("                     Background: none (no image "
-		       "output)\n");
-	}
-	STATUS("               Partiality model: xsphere (hardcoded)\n");
-	STATUS("       Noise standard deviation: %.2f detector units\n",
-	       noise_stddev);
-	if ( random_intensities ) {
-		STATUS("               Full intensities: randomly generated: "
-		       "abs(Gaussian(sigma=%.2f)), symmetry %s\n",
-		       full_stddev, sym_str);
-	} else {
-		STATUS("               Full intensities: from %s (symmetry %s)\n",
-		       input_file, sym_str);
-	}
-	STATUS("   Max error in cell components: %.2f %%\n", cnoise);
-	STATUS("Scale factor standard deviation: %.2f\n", osf_stddev);
-	if ( template_stream != NULL ) {
-		STATUS("Crystal orientations and reflections to use from %s\n",
-		       template);
-	}
-
-	if ( random_intensities ) {
-		full = reflist_new();
-	}
-
-	qargs.full = full;
-	pthread_rwlock_init(&qargs.full_lock, NULL);
-	qargs.n_to_do = n;
-	qargs.n_done = 0;
-	qargs.n_started = 0;
-	qargs.sym = sym;
-	qargs.dtempl = dtempl;
-	qargs.random_intensities = random_intensities;
-	qargs.cell = cell;
-	qargs.stream = stream;
-	qargs.cnoise = cnoise;
-	qargs.osf_stddev = osf_stddev;
-	qargs.full_stddev = full_stddev;
-	qargs.noise_stddev = noise_stddev;
-	qargs.background = background;
-	qargs.max_q = detgeom_max_resolution(test_image->detgeom,
-	                                     test_image->lambda);
-	qargs.image_prefix = image_prefix;
-	qargs.profile_radius = profile_radius;
-	qargs.template_stream = template_stream;
-
-	qargs.rngs = malloc(n_threads * sizeof(gsl_rng *));
-	if ( qargs.rngs == NULL ) {
-		ERROR("Failed to allocate RNGs\n");
-		return 1;
-	}
-
-	if ( config_random ) {
-
-		FILE *fh;
-
-		fh = fopen("/dev/urandom", "r");
-		if ( fh == NULL ) {
-			ERROR("Failed to open /dev/urandom.  Try again without"
-			      " --really-random.\n");
-			free(qargs.rngs);
-			return 1;
-		}
-
-		for ( i=0; i<n_threads; i++ ) {
-
-			unsigned long int seed;
-			qargs.rngs[i] = gsl_rng_alloc(gsl_rng_mt19937);
-
-			if ( fread(&seed, sizeof(seed), 1, fh) == 1 ) {
-				gsl_rng_set(qargs.rngs[i], seed);
-			} else {
-				ERROR("Failed to seed RNG %i\n", i);
-			}
-
-		}
-
-		fclose(fh);
-
-	} else {
-		gsl_rng *rng_for_seeds;
-		rng_for_seeds = gsl_rng_alloc(gsl_rng_mt19937);
-		for ( i=0; i<n_threads; i++ ) {
-			qargs.rngs[i] = gsl_rng_alloc(gsl_rng_mt19937);
-			gsl_rng_set(qargs.rngs[i], gsl_rng_get(rng_for_seeds));
-		}
-		gsl_rng_free(rng_for_seeds);
-	}
-
-	for ( i=0; i<NBINS; i++ ) {
-		qargs.n_ref[i] = 0;
-		qargs.p_hist[i] = 0.0;
-		qargs.p_max[i] = 0.0;
-	}
-
-	run_threads(n_threads, run_job, create_job, finalise_job,
-	            &qargs, n, 0, 0, 0);
-
-	if ( random_intensities ) {
-		STATUS("Writing full intensities to %s\n", save_file);
-		write_reflist_2(save_file, full, sym);
-	}
-
-	if ( phist_file != NULL ) {
-
-		FILE *fh;
-
-		fh = fopen(phist_file, "w");
-
-		if ( fh != NULL ) {
-
-			double overall_max = 0.0;
-			double overall_mean = 0.0;
-			long long int overall_total = 0;
-
-			for ( i=0; i<NBINS; i++ ) {
-
-				double rcen;
-
-				if ( qargs.p_max[i] > overall_max ) {
-					overall_max = qargs.p_max[i];
-				}
-
-				overall_mean += qargs.p_hist[i];
-				overall_total += qargs.n_ref[i];
-
-				rcen = i/(double)NBINS*qargs.max_q
-					  + qargs.max_q/(2.0*NBINS);
-				fprintf(fh, "%.2f %7lu %.3f %.3f\n", rcen/1.0e9,
-					qargs.n_ref[i],
-					qargs.p_hist[i]/qargs.n_ref[i],
-					qargs.p_max[i]);
-
-			}
-
-			fclose(fh);
-
-			overall_mean /= overall_total;
-
-			STATUS("Overall max partiality = %.2f\n",
-			       overall_max);
-			STATUS("Overall mean partiality = %.2f\n", overall_mean);
-			STATUS("Total number of reflections = %lli\n",
-			       overall_total);
-
-		} else {
-			ERROR("Failed to open file '%s' for writing.\n",
-			      phist_file);
-		}
-
-	}
-
-	for ( i=0; i<n_threads; i++ ) {
-		gsl_rng_free(qargs.rngs[i]);
-	}
-	free(qargs.rngs);
-	pthread_rwlock_destroy(&qargs.full_lock);
-	stream_close(stream);
-	cell_free(cell);
-	data_template_free(dtempl);
-	free_symoplist(sym);
-	reflist_free(full);
-	free(save_file);
-	free(geomfile);
-	free(input_file);
-
-	return 0;
-}
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
deleted file mode 100644
index 9769a0d4..00000000
--- a/src/pattern_sim.c
+++ /dev/null
@@ -1,1183 +0,0 @@
-/*
- * pattern_sim.c
- *
- * Simulate diffraction patterns from small crystals
- *
- * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
- *                       a research centre of the Helmholtz Association.
- *
- * Authors:
- *   2009-2020 Thomas White <taw@physics.org>
- *   2013-2014 Chun Hong Yoon <chun.hong.yoon@desy.de>
- *   2014      Valerio Mariani
- *   2013      Alexandra Tolstikova
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <stdarg.h>
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-#include <unistd.h>
-#include <getopt.h>
-#include <hdf5.h>
-
-#include <image.h>
-#include <cell.h>
-#include <cell-utils.h>
-#include <utils.h>
-#include <peaks.h>
-#include <symmetry.h>
-#include <reflist.h>
-#include <reflist-utils.h>
-#include <stream.h>
-#include <datatemplate.h>
-
-#include "diffraction.h"
-#include "diffraction-gpu.h"
-#include "pattern_sim.h"
-#include "version.h"
-
-
-enum spectrum_type {
-	SPECTRUM_TOPHAT,    /**< A top hat distribution */
-	SPECTRUM_GAUSSIAN,  /**< A Gaussian distribution */
-	SPECTRUM_SASE,      /**< A simulated SASE spectrum */
-	SPECTRUM_TWOCOLOUR, /**< A spectrum containing two peaks */
-	SPECTRUM_FROMFILE   /**< An arbitrary spectrum read from a file */
-};
-
-static void show_help(const char *s)
-{
-	printf("Syntax: %s [options]\n\n", s);
-	printf(
-"Simulate diffraction patterns from small crystals probed with femtosecond\n"
-"pulses of X-rays from a free electron laser.\n"
-"\n"
-" -h, --help                Display this help message.\n"
-"     --version             Print CrystFEL version number and exit.\n"
-"\n"
-" -p, --pdb=<file>          File from which to get the unit cell.\n"
-"                            (The actual Bragg intensities come from the\n"
-"                            intensities file)\n"
-"     --gpu                 Use the GPU to speed up the calculation.\n"
-"     --gpu-dev=<n>         Use GPU device <n>.  Omit this option to see the\n"
-"                            available devices.\n"
-" -g, --geometry=<file>     Get detector geometry from file.\n"
-" -n, --number=<N>          Generate N images.  Default 1.\n"
-"     --no-images           Do not output any HDF5 files.\n"
-" -o, --output=<filename>   Output HDF5 filename.  Default: sim.h5.\n"
-" -r, --random-orientation  Use randomly generated orientations.\n"
-"     --powder=<file>       Write a summed pattern of all images simulated by\n"
-"                            this invocation as the given filename.\n"
-" -i, --intensities=<file>  Specify file containing reflection intensities\n"
-"                            (and phases) to use.\n"
-" -y, --symmetry=<sym>      The symmetry of the intensities file.\n"
-" -t, --gradients=<method>  Select method for calculation of shape transforms\n"
-"     --really-random       Seed the random number generator with /dev/urandom.\n"
-"     --min-size=<s>        Minimum crystal size in nm.\n"
-"     --max-size=<s>        Maximum crystal size in nm.\n"
-"     --no-noise            Do not calculate Poisson noise.\n"
-" -s, --sample-spectrum=<N> Use N samples from spectrum. Default 3.\n"
-" -x, --spectrum=<type>     Type of spectrum to simulate.\n"
-"     --background=<N>      Add N photons of Poisson background (default 0).\n"
-"     --template=<file>     Take orientations from stream <file>.\n"
-"     --no-fringes          Exclude the side maxima of Bragg peaks.\n"
-"     --flat                Make Bragg peaks flat.\n"
-"     --sase-spike-width    SASE spike width (standard deviation in m^-1)\n"
-"                            Default 2e6 m^-1\n"
-"     --twocol-separation   Separation between peaks in two-colour mode in m^-1\n"
-"                            Default 8e6 m^-1\n"
-"     --nphotons            Number of photons per X-ray pulse.  Default 1e12.\n"
-"     --beam-radius         Radius of beam in metres (default 1e-6).\n"
-);
-}
-
-
-static void record_panel(struct detgeom_panel *p, float *dp,
-                         int do_poisson,
-                         gsl_rng *rng, double ph_per_e, double background,
-			 double lambda,
-                         int *n_neg1, int *n_inf1, int *n_nan1,
-                         int *n_neg2, int *n_inf2, int *n_nan2,
-			 double *max_tt)
-{
-	int fs, ss;
-
-	for ( ss=0; ss<p->h; ss++ ) {
-	for ( fs=0; fs<p->w; fs++ ) {
-
-		double counts;
-		double intensity, sa;
-		double pix_area, Lsq;
-		double xs, ys, rx, ry;
-		double dsq, proj_area;
-		float dval;
-		double twotheta;
-
-		intensity = (double)dp[fs + p->w*ss];
-		if ( isinf(intensity) ) (*n_inf1)++;
-		if ( intensity < 0.0 ) (*n_neg1)++;
-		if ( isnan(intensity) ) (*n_nan1)++;
-
-		/* Area of one pixel */
-		pix_area = pow(p->pixel_pitch, 2.0);
-		Lsq = pow(p->cnz*p->pixel_pitch, 2.0);
-
-		/* Calculate distance from crystal to pixel */
-		xs = fs*p->fsx + ss*p->ssx;
-		ys = ss*p->fsy + ss*p->ssy;
-		rx = (xs + p->cnx) * p->pixel_pitch;
-		ry = (ys + p->cny) * p->pixel_pitch;
-		dsq = pow(rx, 2.0) + pow(ry, 2.0);
-		twotheta = atan2(sqrt(dsq), p->cnz*p->pixel_pitch);
-
-		/* Area of pixel as seen from crystal (approximate) */
-		proj_area = pix_area * cos(twotheta);
-
-		/* Projected area of pixel divided by distance squared */
-		sa = proj_area / (dsq + Lsq);
-
-		if ( do_poisson ) {
-			counts = poisson_noise(rng, intensity * ph_per_e * sa);
-		} else {
-			counts = intensity * ph_per_e * sa;
-		}
-
-		/* Number of photons in pixel */
-		dval = counts + poisson_noise(rng, background);
-
-		/* Convert to ADU */
-		dval *= p->adu_per_photon;
-
-		/* Saturation */
-		if ( dval > p->max_adu ) dval = p->max_adu;
-
-		dp[fs + p->w*ss] = dval;
-
-		/* Sanity checks */
-		if ( isinf(dp[fs + p->w*ss]) ) n_inf2++;
-		if ( isnan(dp[fs + p->w*ss]) ) n_nan2++;
-		if ( dp[fs + p->w*ss] < 0.0 ) n_neg2++;
-
-		if ( twotheta > *max_tt ) *max_tt = twotheta;
-
-	}
-	}
-}
-
-
-void record_image(struct image *image, int do_poisson, double background,
-                  gsl_rng *rng, double beam_radius, double nphotons)
-{
-	double total_energy, energy_density;
-	double ph_per_e;
-	double area;
-	double max_tt = 0.0;
-	int n_inf1 = 0;
-	int n_neg1 = 0;
-	int n_nan1 = 0;
-	int n_inf2 = 0;
-	int n_neg2 = 0;
-	int n_nan2 = 0;
-	int pn;
-
-	/* How many photons are scattered per electron? */
-	area = M_PI*pow(beam_radius, 2.0);
-	total_energy = nphotons * ph_lambda_to_en(image->lambda);
-	energy_density = total_energy / area;
-	ph_per_e = (nphotons /area) * pow(THOMSON_LENGTH, 2.0);
-	STATUS("Fluence = %8.2e photons, "
-	       "Energy density = %5.3f kJ/cm^2, "
-	       "Total energy = %5.3f microJ\n",
-	       nphotons, energy_density/1e7, total_energy*1e6);
-
-	for ( pn=0; pn<image->detgeom->n_panels; pn++ ) {
-
-		record_panel(&image->detgeom->panels[pn], image->dp[pn],
-		             do_poisson, rng, ph_per_e, background,
-			     image->lambda,
-		             &n_neg1, &n_inf1, &n_nan1,
-		             &n_neg2, &n_inf2, &n_nan2, &max_tt);
-	}
-
-	STATUS("Max 2theta = %.2f deg, min d = %.2f nm\n",
-	        rad2deg(max_tt), (image->lambda/(2.0*sin(max_tt/2.0)))/1e-9);
-
-	STATUS("Halve the d values to get the voxel size for a synthesis.\n");
-
-	if ( n_neg1 + n_inf1 + n_nan1 + n_neg2 + n_inf2 + n_nan2 ) {
-		ERROR("WARNING: The raw calculation produced %i negative"
-		      " values, %i infinities and %i NaNs.\n",
-		      n_neg1, n_inf1, n_nan1);
-		ERROR("WARNING: After processing, there were %i negative"
-		      " values, %i infinities and %i NaNs.\n",
-		      n_neg2, n_inf2, n_nan2);
-	}
-}
-
-
-static double *intensities_from_list(RefList *list, SymOpList *sym)
-{
-	Reflection *refl;
-	RefListIterator *iter;
-	double *out = new_arr_intensity();
-	SymOpMask *m = new_symopmask(sym);
-	int neq = num_equivs(sym, NULL);
-
-	for ( refl = first_refl(list, &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) ) {
-
-		signed int h, k, l;
-		int eps;
-		double intensity = get_intensity(refl);
-
-		get_indices(refl, &h, &k, &l);
-		special_position(sym, m, h, k, l);
-		eps = neq / num_equivs(sym, m);
-
-		set_arr_intensity(out, h, k, l, intensity / eps);
-
-	}
-
-	return out;
-}
-
-
-static double *phases_from_list(RefList *list)
-{
-	Reflection *refl;
-	RefListIterator *iter;
-	double *out = new_arr_phase();
-
-	for ( refl = first_refl(list, &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) ) {
-
-		signed int h, k, l;
-		double phase = get_phase(refl, NULL);
-
-		get_indices(refl, &h, &k, &l);
-
-		set_arr_phase(out, h, k, l, phase);
-
-	}
-
-	return out;
-
-}
-
-
-static unsigned char *flags_from_list(RefList *list)
-{
-	Reflection *refl;
-	RefListIterator *iter;
-	unsigned char *out = new_arr_flag();
-
-	for ( refl = first_refl(list, &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) ) {
-
-		signed int h, k, l;
-
-		get_indices(refl, &h, &k, &l);
-
-		set_arr_flag(out, h, k, l, 1);
-
-	}
-
-	return out;
-
-}
-
-
-static struct quaternion read_quaternion()
-{
-	do {
-
-		int r;
-		float w, x, y, z;
-		char line[1024];
-		char *rval;
-
-		printf("Please input quaternion: w x y z\n");
-		rval = fgets(line, 1023, stdin);
-		if ( rval == NULL ) return invalid_quaternion();
-		chomp(line);
-
-		r = sscanf(line, "%f %f %f %f", &w, &x, &y, &z);
-		if ( r == 4 ) {
-
-			struct quaternion quat;
-
-			quat.w = w;
-			quat.x = x;
-			quat.y = y;
-			quat.z = z;
-
-			return quat;
-
-		} else {
-			ERROR("Invalid rotation '%s'\n", line);
-		}
-
-	} while ( 1 );
-}
-
-
-static int random_ncells(double len, double min_m, double max_m)
-{
-	int min_cells, max_cells, wid;
-
-	min_cells = min_m / len;
-	max_cells = max_m / len;
-	wid = max_cells - min_cells;
-
-	return wid*random()/RAND_MAX + min_cells;
-}
-
-
-static void add_metadata(const char *filename,
-                         struct quaternion q,
-                         UnitCell *cell)
-{
-	hid_t fh, gh, sh, dh;
-	herr_t r;
-	hsize_t size[1];
-	float data[4];
-	double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz;
-
-	cell_get_reciprocal(cell, &asx, &asy, &asz,
-	                          &bsx, &bsy, &bsz,
-	                          &csx, &csy, &csz);
-
-	fh = H5Fopen(filename, H5F_ACC_RDWR, H5P_DEFAULT);
-	if ( fh < 0 ) {
-		ERROR("Failed to open file to add metadata.\n");
-		return;
-	}
-
-	gh = H5Gcreate2(fh, "pattern_sim", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-	if ( gh < 0 ) {
-		ERROR("Failed to create metadata group.\n");
-		H5Fclose(fh);
-		return;
-	}
-
-	size[0] = 4;
-	sh = H5Screate_simple(1, size, NULL);
-
-	dh = H5Dcreate2(gh, "quaternion", H5T_NATIVE_FLOAT, sh,
-	                H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-	if ( dh < 0 ) {
-		ERROR("Couldn't create dataset\n");
-		H5Fclose(fh);
-		return;
-	}
-
-	data[0] = q.w;
-	data[1] = q.x;
-	data[2] = q.y;
-	data[3] = q.z;
-	r = H5Dwrite(dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, data);
-	if ( r < 0 ) {
-		ERROR("Failed to write quaternion to file\n");
-	}
-	H5Dclose(dh);
-
-	size[0] = 3;
-	sh = H5Screate_simple(1, size, NULL);
-
-	dh = H5Dcreate2(gh, "astar", H5T_NATIVE_FLOAT, sh,
-	                H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-	if ( dh < 0 ) {
-		ERROR("Couldn't create dataset\n");
-		H5Fclose(fh);
-		return;
-	}
-	data[0] = asx;
-	data[1] = asy;
-	data[2] = asz;
-	r = H5Dwrite(dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, data);
-	if ( r < 0 ) {
-		ERROR("Failed to write astar to file.\n");
-	}
-	H5Dclose(dh);
-
-	dh = H5Dcreate2(gh, "bstar", H5T_NATIVE_FLOAT, sh,
-	                H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-	if ( dh < 0 ) {
-		ERROR("Couldn't create dataset\n");
-		H5Fclose(fh);
-		return;
-	}
-	data[0] = bsx;
-	data[1] = bsy;
-	data[2] = bsz;
-	r = H5Dwrite(dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, data);
-	if ( r < 0 ) {
-		ERROR("Failed to write bstar to file.\n");
-	}
-	H5Dclose(dh);
-
-	dh = H5Dcreate2(gh, "cstar", H5T_NATIVE_FLOAT, sh,
-	                H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-	if ( dh < 0 ) {
-		ERROR("Couldn't create dataset\n");
-		H5Fclose(fh);
-		return;
-	}
-	data[0] = csx;
-	data[1] = csy;
-	data[2] = csz;
-	r = H5Dwrite(dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, data);
-	if ( r < 0 ) {
-		ERROR("Failed to write cstar to file.\n");
-	}
-	H5Dclose(dh);
-
-	H5Gclose(gh);
-	H5Fclose(fh);
-}
-
-
-int main(int argc, char *argv[])
-{
-	int argn;
-	struct image *image;
-	DataTemplate *dtempl;
-	struct gpu_context *gctx = NULL;
-	struct image *powder;
-	char *intfile = NULL;
-	double *intensities;
-	char *rval;
-	double *phases;
-	unsigned char *flags;
-	int config_randomquat = 0;
-	int config_noimages = 0;
-	int config_nonoise = 0;
-	int config_nosfac = 0;
-	int config_gpu = 0;
-	int config_random = 0;
-	char *powder_fn = NULL;
-	char *cell_filename = NULL;
-	char *grad_str = NULL;
-	char *outfile = NULL;
-	char *geometry = NULL;
-	char *spectrum_str = NULL;
-	char *spectrum_fn = NULL;
-	GradientMethod grad;
-	enum spectrum_type spectrum_type;
-	int ndone = 0;    /* Number of simulations done (images or not) */
-	int number = 1;   /* Number used for filename of image */
-	int n_images = 1; /* Generate one image by default */
-	int done = 0;
-	UnitCell *input_cell;
-	struct quaternion orientation;
-	int gpu_dev = -1;
-	int random_size = 0;
-	double min_size = 0.0;
-	double max_size = 0.0;
-	char *sym_str = NULL;
-	SymOpList *sym;
-	int nsamples = 3;
-	gsl_rng *rng;
-	double background = 0.0;
-	char *template_file = NULL;
-	Stream *st = NULL;
-	int no_fringes = 0;
-	int flat = 0;
-	double nphotons = 1e12;
-	double beam_radius = 1e-6;  /* metres */
-	double sase_spike_width = 2e6;  /* m^-1 */
-	double twocol_sep = 8e6;  /* m^-1 */
-
-	/* Long options */
-	const struct option longopts[] = {
-		{"help",               0, NULL,               'h'},
-		{"version",            0, NULL,               'v'},
-		{"gpu",                0, &config_gpu,         1},
-		{"beam",               1, NULL,               'b'},
-		{"random-orientation", 0, NULL,               'r'},
-		{"number",             1, NULL,               'n'},
-		{"no-images",          0, &config_noimages,    1},
-		{"no-noise",           0, &config_nonoise,     1},
-		{"intensities",        1, NULL,               'i'},
-		{"symmetry",           1, NULL,               'y'},
-		{"powder",             1, NULL,               'w'},
-		{"gradients",          1, NULL,               't'},
-		{"pdb",                1, NULL,               'p'},
-		{"output",             1, NULL,               'o'},
-		{"geometry",           1, NULL,               'g'},
-		{"sample-spectrum",    1, NULL,               's'},
-		{"type-spectrum",      1, NULL,               'x'},
-		{"spectrum",           1, NULL,               'x'},
-		{"really-random",      0, &config_random,      1},
-		{"no-fringes",         0, &no_fringes,         1},
-		{"flat",               0, &flat,               1},
-
-		{"gpu-dev",            1, NULL,                2},
-		{"min-size",           1, NULL,                3},
-		{"max-size",           1, NULL,                4},
-		{"background",         1, NULL,                5},
-		{"template",           1, NULL,                6},
-		{"beam-bandwidth",     1, NULL,                7},
-		{"photon-energy",      1, NULL,                9},
-		{"nphotons",           1, NULL,               10},
-		{"beam-radius",        1, NULL,               11},
-		{"spectrum-file",      1, NULL,               12},
-		{"sase-spike-width",   1, NULL,               13},
-		{"twocol-separation",  1, NULL,               14},
-
-		{0, 0, NULL, 0}
-	};
-
-	/* Short options */
-	while ((argn = getopt_long(argc, argv, "hrn:i:t:p:o:g:y:s:x:vb:",
-	                        longopts, NULL)) != -1) {
-
-		switch (argn) {
-
-			case 'h' :
-			show_help(argv[0]);
-			return 0;
-
-			case 'v' :
-			printf("CrystFEL: %s\n",
-			       crystfel_version_string());
-			printf("%s\n",
-			       crystfel_licence_string());
-			return 0;
-
-			case 'b' :
-			ERROR("WARNING: This version of CrystFEL no longer "
-			      "uses beam files.  Please remove the beam file "
-			      "from your pattern_sim command line.\n");
-			return 1;
-
-			case 'r' :
-			config_randomquat = 1;
-			break;
-
-			case 'n' :
-			n_images = strtol(optarg, &rval, 10);
-			if ( *rval != '\0' ) {
-				ERROR("Invalid number of images.\n");
-				return 1;
-			}
-			break;
-
-			case 'i' :
-			intfile = strdup(optarg);
-			break;
-
-			case 't' :
-			grad_str = strdup(optarg);
-			break;
-
-			case 'p' :
-			cell_filename = strdup(optarg);
-			break;
-
-			case 'o' :
-			outfile = strdup(optarg);
-			break;
-
-			case 'w' :
-			powder_fn = strdup(optarg);
-			break;
-
-			case 'g' :
-			geometry = strdup(optarg);
-			break;
-
-			case 'y' :
-			sym_str = strdup(optarg);
-			break;
-
-			case 's' :
-			nsamples = strtol(optarg, &rval, 10);
-			if ( *rval != '\0' ) {
-				ERROR("Invalid number of spectrum samples.\n");
-				return 1;
-			}
-			break;
-
-			case 'x' :
-			spectrum_str = strdup(optarg);
-			break;
-
-			case 2 :
-			gpu_dev = atoi(optarg);
-			break;
-
-			case 3 :
-			min_size = strtod(optarg, &rval);
-			if ( *rval != '\0' ) {
-				ERROR("Invalid minimum size.\n");
-				return 1;
-			}
-			min_size /= 1e9;
-			random_size++;
-			break;
-
-			case 4 :
-			max_size = strtod(optarg, &rval);
-			if ( *rval != '\0' ) {
-				ERROR("Invalid maximum size.\n");
-				return 1;
-			}
-			max_size /= 1e9;
-			random_size++;
-			break;
-
-			case 5 :
-			background = strtod(optarg, &rval);
-			if ( *rval != '\0' ) {
-				ERROR("Invalid background level.\n");
-				return 1;
-			}
-			break;
-
-			case 6 :
-			template_file = strdup(optarg);
-			break;
-
-			case 7 :
-			ERROR("--beam-bandwidth is no longer used.\n");
-			ERROR("Set the bandwidth in the geometry file instead.\n");
-			return 1;
-
-			case 9 :
-			ERROR("--photon-energy is no longer used.\n");
-			ERROR("Set the photon energy in the geometry file instead.\n");
-			return 1;
-
-			case 10 :
-			nphotons = strtod(optarg, &rval);
-			if ( *rval != '\0' ) {
-				ERROR("Invalid number of photons.\n");
-				return 1;
-			}
-			if ( nphotons < 0.0 ) {
-				ERROR("Number of photons must be positive.\n");
-				return 1;
-			}
-			break;
-
-			case 11 :
-			beam_radius = strtod(optarg, &rval);
-			if ( *rval != '\0' ) {
-				ERROR("Invalid beam radius.\n");
-				return 1;
-			}
-			if ( beam_radius < 0.0 ) {
-				ERROR("Beam radius must be positive.\n");
-				return 1;
-			}
-			break;
-
-			case 12 :
-			spectrum_fn = strdup(optarg);
-			break;
-
-			case 13 :
-			sase_spike_width = strtod(optarg, &rval);
-			if ( *rval != '\0' ) {
-				ERROR("Invalid SASE spike width.\n");
-				return 1;
-			}
-			if ( beam_radius < 0.0 ) {
-				ERROR("SASE spike width must be positive.\n");
-				return 1;
-			}
-			break;
-
-			case 14 :
-			twocol_sep = strtod(optarg, &rval);
-			if ( *rval != '\0' ) {
-				ERROR("Invalid two colour separation.\n");
-				return 1;
-			}
-			if ( beam_radius < 0.0 ) {
-				ERROR("Two-colour separation must be positive.\n");
-				return 1;
-			}
-			break;
-
-			case 0 :
-			break;
-
-			case '?' :
-			break;
-
-			default :
-			ERROR("Unhandled option '%c'\n", argn);
-			break;
-
-		}
-
-	}
-
-	if ( random_size == 1 ) {
-		ERROR("You must specify both --min-size and --max-size.\n");
-		return 1;
-	}
-
-	if ( cell_filename == NULL ) {
-		cell_filename = strdup("molecule.pdb");
-	}
-
-	if ( outfile == NULL ) {
-		if ( n_images == 1 ) {
-			outfile = strdup("sim.h5");
-		} else {
-			outfile = strdup("sim");
-		}
-	}
-
-	if ( template_file != NULL ) {
-		if ( config_randomquat ) {
-			ERROR("You cannot use -r and --template together.\n");
-			return 1;
-		}
-		st = stream_open_for_read(template_file);
-		if ( st == NULL ) {
-			ERROR("Failed to open stream.\n");
-			return 1;
-		}
-		free(template_file);
-	}
-
-	if ( grad_str == NULL ) {
-		STATUS("You didn't specify a gradient calculation method, so"
-		       " I'm using the 'mosaic' method, which is fastest.\n");
-		grad = GRADIENT_MOSAIC;
-	} else if ( strcmp(grad_str, "mosaic") == 0 ) {
-		grad = GRADIENT_MOSAIC;
-	} else if ( strcmp(grad_str, "interpolate") == 0) {
-		grad = GRADIENT_INTERPOLATE;
-	} else if ( strcmp(grad_str, "phased") == 0) {
-		grad = GRADIENT_PHASED;
-	} else {
-		ERROR("Unrecognised gradient method '%s'\n", grad_str);
-		return 1;
-	}
-	free(grad_str);
-
-	if ( config_gpu && (grad != GRADIENT_MOSAIC) ) {
-		ERROR("Only the mosaic method can be used for gradients when"
-		      "calculating on the GPU.\n");
-		return 1;
-	}
-
-	if ( geometry == NULL ) {
-		ERROR("You need to specify a geometry file with --geometry\n");
-		return 1;
-	}
-	dtempl = data_template_new_from_file(geometry);
-	if ( dtempl == NULL ) {
-		ERROR("Failed to read detector geometry from '%s'\n",
-		      geometry);
-		return 1;
-	}
-	free(geometry);
-
-	if ( spectrum_str == NULL ) {
-		STATUS("You didn't specify a spectrum type, so"
-		       " I'm using a 'tophat' spectrum.\n");
-		spectrum_type = SPECTRUM_TOPHAT;
-	} else if ( strcasecmp(spectrum_str, "tophat") == 0) {
-		spectrum_type = SPECTRUM_TOPHAT;
-	} else if ( strcasecmp(spectrum_str, "gaussian") == 0) {
-		spectrum_type = SPECTRUM_GAUSSIAN;
-	} else if ( strcasecmp(spectrum_str, "sase") == 0) {
-		spectrum_type = SPECTRUM_SASE;
-	} else if ( strcasecmp(spectrum_str, "twocolour") == 0 ||
-                    strcasecmp(spectrum_str, "twocolor") == 0 ||
-                    strcasecmp(spectrum_str, "twocolours") == 0 ||
-                    strcasecmp(spectrum_str, "twocolors") == 0) {
-		spectrum_type = SPECTRUM_TWOCOLOUR;
-	} else if ( strcasecmp(spectrum_str, "fromfile") == 0) {
-		spectrum_type = SPECTRUM_FROMFILE;
-	} else {
-		ERROR("Unrecognised spectrum type '%s'\n", spectrum_str);
-		return 1;
-	}
-	free(spectrum_str);
-
-	/* Load unit cell */
-	input_cell = load_cell_from_file(cell_filename);
-	if ( input_cell == NULL ) {
-		exit(1);
-	}
-
-	if ( intfile == NULL ) {
-
-		/* Gentle reminder */
-		STATUS("You didn't specify the file containing the ");
-		STATUS("reflection intensities (with --intensities).\n");
-		STATUS("I'll simulate a flat intensity distribution.\n");
-		intensities = NULL;
-		phases = NULL;
-		flags = NULL;
-
-		if ( sym_str == NULL ) sym_str = strdup("1");
-		pointgroup_warning(sym_str);
-		sym = get_pointgroup(sym_str);
-
-	} else {
-
-		RefList *reflections;
-		char *sym_str_fromfile = NULL;
-
-		reflections = read_reflections_2(intfile, &sym_str_fromfile);
-		if ( reflections == NULL ) {
-			ERROR("Problem reading input file %s\n", intfile);
-			return 1;
-		}
-
-		/* If we don't have a point group yet, and if the file provides
-		 * one, use the one from the file */
-		if ( (sym_str == NULL) && (sym_str_fromfile != NULL) ) {
-			sym_str = sym_str_fromfile;
-			STATUS("Using symmetry from reflection file: %s\n",
-			       sym_str);
-		}
-
-		/* If we still don't have a point group, use "1" */
-		if ( sym_str == NULL ) sym_str = strdup("1");
-
-		pointgroup_warning(sym_str);
-		sym = get_pointgroup(sym_str);
-
-		if ( grad == GRADIENT_PHASED ) {
-			phases = phases_from_list(reflections);
-		} else {
-			phases = NULL;
-		}
-		intensities = intensities_from_list(reflections, sym);
-		flags = flags_from_list(reflections);
-
-		/* Check that the intensities have the correct symmetry */
-		if ( check_list_symmetry(reflections, sym) ) {
-			ERROR("The input reflection list does not appear to"
-			      " have symmetry %s\n", symmetry_name(sym));
-			if ( cell_get_lattice_type(input_cell) == L_MONOCLINIC )
-			{
-				ERROR("You may need to specify the unique axis "
-				      "in your point group.  The default is "
-				      "unique axis c.\n");
-				ERROR("See 'man crystfel' for more details.\n");
-			}
-			return 1;
-		}
-
-		reflist_free(reflections);
-
-	}
-
-	rng = gsl_rng_alloc(gsl_rng_mt19937);
-	if ( config_random ) {
-		FILE *fh;
-		unsigned long int seed;
-		fh = fopen("/dev/urandom", "r");
-		if ( fread(&seed, sizeof(seed), 1, fh) == 1 ) {
-			gsl_rng_set(rng, seed);
-		} else {
-			ERROR("Failed to seed random number generator\n");
-		}
-		fclose(fh);
-	}
-
-	image = image_create_for_simulation(dtempl);
-	powder = image_create_for_simulation(dtempl);
-
-	/* Splurge a few useful numbers */
-	STATUS("Simulation parameters:\n");
-	STATUS("                     Wavelength: %.5f A (photon energy %.2f eV)\n",
-	       image->lambda*1e10,
-	       ph_lambda_to_eV(image->lambda));
-	STATUS("              Number of photons: %.0f (%.2f mJ)\n",
-	       nphotons, eV_to_J(ph_lambda_to_eV(image->lambda))*nphotons*1e3);
-	STATUS("                Beam divergence: not simulated\n");
-	STATUS("                    Beam radius: %.2f microns\n",
-	      beam_radius*1e6);
-	STATUS("                     Background: %.2f photons\n", background);
-
-	switch ( spectrum_type ) {
-
-		case SPECTRUM_TOPHAT:
-		STATUS("                 X-ray spectrum: top hat, "
-		       "width %.5f %%\n", image->bw*100.0);
-		break;
-
-		case SPECTRUM_GAUSSIAN:
-		STATUS("                 X-ray spectrum: Gaussian, "
-		       "bandwidth %.5f %%\n", image->bw*100.0);
-		break;
-
-		case SPECTRUM_SASE:
-		STATUS("                 X-ray spectrum: SASE, "
-		       "bandwidth %.5f %%\n", image->bw*100.0);
-		break;
-
-		case SPECTRUM_TWOCOLOUR:
-		STATUS("                 X-ray spectrum: two colour, "
-		       "bandwidth %.5f %%, separation %.5e m^-1\n",
-		       image->bw*100.0, twocol_sep);
-		break;
-
-		case SPECTRUM_FROMFILE:
-		STATUS("                 X-ray spectrum: from %s\n",
-			spectrum_fn);
-		break;
-	}
-	if ( random_size ) {
-		STATUS("                   Crystal size: random, between "
-		       "%.2f and %.2f nm along each of a, b and c\n",
-		       min_size*1e9, max_size*1e9);
-	} else {
-		STATUS("                   Crystal size: 8 unit cells along "
-		       "each of a, b and c\n");
-	}
-	if ( intfile == NULL ) {
-		STATUS("               Full intensities: all equal\n");
-	} else {
-		STATUS("               Full intensities: from %s (symmetry %s)\n",
-		       intfile, sym_str);
-	}
-	do {
-
-		int na, nb, nc;
-		double a, b, c, dis;
-		UnitCell *cell;
-		int err = 0;
-		int pi;
-
-		/* Reset image data to zero for new pattern */
-		for ( pi=0; pi<image->detgeom->n_panels; pi++ ) {
-			long j;
-			long np = image->detgeom->panels[pi].w
-			            * image->detgeom->panels[pi].h;
-			for ( j=0; j<np; j++ ) {
-				image->dp[pi][j] = 0.0;
-			}
-		}
-
-		if ( random_size ) {
-
-			double alen, blen, clen, discard;
-
-			cell_get_parameters(input_cell, &alen, &blen, &clen,
-			                    &discard, &discard, &discard);
-
-			na = random_ncells(alen, min_size, max_size);
-			nb = random_ncells(blen, min_size, max_size);
-			nc = random_ncells(clen, min_size, max_size);
-
-		} else {
-
-			na = 8;
-			nb = 8;
-			nc = 8;
-
-		}
-
-		if ( st == NULL ) {
-
-			if ( config_randomquat ) {
-				orientation = random_quaternion(rng);
-			} else {
-				orientation = read_quaternion();
-			}
-
-			STATUS("Orientation is %5.3f %5.3f %5.3f %5.3f\n",
-			       orientation.w, orientation.x,
-			       orientation.y, orientation.z);
-
-			if ( !quaternion_valid(orientation) ) {
-				ERROR("Orientation modulus is not zero!\n");
-				return 1;
-			}
-
-			cell = cell_rotate(input_cell, orientation);
-
-		} else {
-
-			struct image *templ_image;
-			Crystal *cr;
-
-			/* Get data from next chunk */
-			templ_image = stream_read_chunk(st, 0);
-			if ( templ_image == NULL ) break;
-			if ( templ_image->n_crystals == 0 ) continue;
-
-			cr = templ_image->crystals[0];
-			cell = cell_new_from_cell(crystal_get_cell(cr));
-
-			if ( templ_image->n_crystals > 1 ) {
-				ERROR("Using the first crystal only.\n");
-			}
-
-			image_free(templ_image);
-
-		}
-
-		switch ( spectrum_type ) {
-
-			case SPECTRUM_TOPHAT :
-			image->spectrum = spectrum_generate_tophat(image->lambda,
-			                                           image->bw);
-			break;
-
-			case SPECTRUM_GAUSSIAN :
-			image->spectrum = spectrum_generate_gaussian(image->lambda,
-			                                             image->bw);
-			break;
-
-
-			case SPECTRUM_SASE :
-			image->spectrum = spectrum_generate_sase(image->lambda,
-			                                         image->bw,
-			                                         sase_spike_width,
-			                                         rng);
-			break;
-
-			case SPECTRUM_TWOCOLOUR :
-			image->spectrum = spectrum_generate_twocolour(image->lambda,
-			                                              image->bw,
-			                                              twocol_sep);
-                        break;
-
-			case SPECTRUM_FROMFILE :
-			image->spectrum = spectrum_load(spectrum_fn);
-                        break;
-
-		}
-
-		cell_get_parameters(cell, &a, &b, &c, &dis, &dis, &dis);
-		STATUS("Particle size = %i x %i x %i"
-		       " ( = %5.2f x %5.2f x %5.2f nm)\n",
-	               na, nb, nc, na*a/1.0e-9, nb*b/1.0e-9, nc*c/1.0e-9);
-
-		if ( config_gpu ) {
-
-			if ( gctx == NULL ) {
-				gctx = setup_gpu(config_nosfac,
-				                 intensities, flags, sym_str,
-				                 gpu_dev);
-			}
-			err = get_diffraction_gpu(gctx, image, na, nb, nc,
-			                          cell, no_fringes, flat,
-			                          nsamples);
-
-		} else {
-			get_diffraction(image, na, nb, nc, intensities, phases,
-			                flags, cell, grad, sym, no_fringes, flat,
-			                nsamples);
-		}
-		if ( err ) {
-			ERROR("Diffraction calculation failed.\n");
-			done = 1;
-			goto skip;
-		}
-
-		record_image(image, !config_nonoise, background, rng,
-		             beam_radius, nphotons);
-
-		if ( powder_fn != NULL ) {
-
-			int pn;
-
-			for ( pn=0; pn<image->detgeom->n_panels; pn++ ) {
-
-				size_t w, i;
-
-				w = image->detgeom->panels[pn].w
-				  * image->detgeom->panels[pn].h;
-
-				for ( i=0; i<w; i++ ) {
-					powder->dp[pn][i] += (double)image->dp[pn][i];
-				}
-
-			}
-
-			if ( !(ndone % 10) ) {
-				image_write(powder, dtempl, powder_fn);
-			}
-		}
-
-		if ( !config_noimages ) {
-
-			char h5filename[1024];
-
-			if ( n_images != 1 ) {
-				snprintf(h5filename, 1023, "%s-%i.h5",
-				         outfile, number);
-			} else {
-				strncpy(h5filename, outfile, 1023);
-			}
-
-			number++;
-
-			/* Write the output file */
-			image_write(image, dtempl, h5filename);
-
-			/* Add some pattern_sim-specific metadata */
-			add_metadata(h5filename, orientation, cell);
-
-		}
-
-		cell_free(cell);
-
-skip:
-		ndone++;
-
-		if ( n_images && (ndone >= n_images) ) done = 1;
-
-	} while ( !done );
-
-	if ( powder_fn != NULL ) {
-		image_write(powder, dtempl, powder_fn);
-	}
-
-	if ( gctx != NULL ) {
-		cleanup_gpu(gctx);
-	}
-
-	image_free(image);
-	image_free(powder);
-	free(intfile);
-	cell_free(input_cell);
-	free(intensities);
-	free(outfile);
-	free(cell_filename);
-	free(sym_str);
-	free_symoplist(sym);
-
-	gsl_rng_free(rng);
-
-	return 0;
-}
diff --git a/src/pattern_sim.h b/src/pattern_sim.h
deleted file mode 100644
index f013c0fb..00000000
--- a/src/pattern_sim.h
+++ /dev/null
@@ -1,67 +0,0 @@
-/*
- * pattern_sim.h
- *
- * Simulate diffraction patterns from small crystals
- *
- * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
- *                       a research centre of the Helmholtz Association.
- *
- * Authors:
- *   2011-2021 Thomas White <taw@physics.org>
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#ifndef PATTERN_SIM_H
-#define PATTERN_SIM_H
-
-
-/* Maxmimum index to hold values up to (can be increased if necessary)
- * WARNING: Altering this value constitutes an ABI change, and means you must
- * update data/diffraction.cl then recompile and reinstall everything. */
-#define INDMAX 130
-
-/* Array size */
-#define IDIM (INDMAX*2 +1)
-#define LIST_SIZE (IDIM*IDIM*IDIM)
-
-/* Create functions for storing reflection intensities indexed as h,k,l */
-#define LABEL(x) x##_intensity
-#define TYPE double
-#include "list_tmp.h"
-
-/* As above, but for phase values */
-#define LABEL(x) x##_phase
-#define TYPE double
-#include "list_tmp.h"
-
-/* As above, but for (unsigned) integer counts */
-#define LABEL(x) x##_count
-#define TYPE unsigned int
-#include "list_tmp.h"
-
-/* As above, but for simple flags */
-#define LABEL(x) x##_flag
-#define TYPE unsigned char
-#include "list_tmp.h"
-
-
-#endif	/* PATTERN_SIM_H */