Merge branch 'master' of ssh://git.bitwiz.org.uk/crystfel

9 files changed, 293 insertions, 72 deletions

diff --git a/libcrystfel/src/hdf5-file.c b/libcrystfel/src/hdf5-file.c
index ad9495b5..ed3b1e29 100644
--- a/libcrystfel/src/hdf5-file.c
+++ b/libcrystfel/src/hdf5-file.c
@@ -384,57 +384,60 @@ int hdf5_write_image(const char *filename, struct image *image)
 
 	H5Dclose(dh);
 
-	arr = malloc(image->spectrum_size*sizeof(double));
-	if ( arr == NULL ) {
-		H5Fclose(fh);
-		return 1;
-	}
-	for ( i=0; i<image->spectrum_size; i++ ) {
-		arr[i] = 1.0e10/image->spectrum[i].k;
-	}
+	if ( image->spectrum_size > 0 ) {
 
-	size[0] = image->spectrum_size;
-	sh = H5Screate_simple(1, size, NULL);
+		arr = malloc(image->spectrum_size*sizeof(double));
+		if ( arr == NULL ) {
+			H5Fclose(fh);
+			return 1;
+		}
+		for ( i=0; i<image->spectrum_size; i++ ) {
+			arr[i] = 1.0e10/image->spectrum[i].k;
+		}
 
-	dh = H5Dcreate2(gh, "spectrum_wavelengths_A", H5T_NATIVE_DOUBLE, sh,
-	                H5P_DEFAULT, H5S_ALL, H5P_DEFAULT);
-	if ( dh < 0 ) {
-		H5Fclose(fh);
-		return 1;
-	}
-	r = H5Dwrite(dh, H5T_NATIVE_DOUBLE, H5S_ALL,
-	             H5S_ALL, H5P_DEFAULT, arr);
-	H5Dclose(dh);
+		size[0] = image->spectrum_size;
+		sh = H5Screate_simple(1, size, NULL);
 
-	for ( i=0; i<image->spectrum_size; i++ ) {
-		arr[i] = image->spectrum[i].weight;
-	}
-	dh = H5Dcreate2(gh, "spectrum_weights", H5T_NATIVE_DOUBLE, sh,
-	                H5P_DEFAULT, H5S_ALL, H5P_DEFAULT);
-	if ( dh < 0 ) {
-		H5Fclose(fh);
-		return 1;
-	}
-	r = H5Dwrite(dh, H5T_NATIVE_DOUBLE, H5S_ALL,
-	             H5S_ALL, H5P_DEFAULT, arr);
+		dh = H5Dcreate2(gh, "spectrum_wavelengths_A", H5T_NATIVE_DOUBLE,
+		                sh, H5P_DEFAULT, H5S_ALL, H5P_DEFAULT);
+		if ( dh < 0 ) {
+			H5Fclose(fh);
+			return 1;
+		}
+		r = H5Dwrite(dh, H5T_NATIVE_DOUBLE, H5S_ALL,
+			     H5S_ALL, H5P_DEFAULT, arr);
+		H5Dclose(dh);
 
-	H5Dclose(dh);
-	free(arr);
+		for ( i=0; i<image->spectrum_size; i++ ) {
+			arr[i] = image->spectrum[i].weight;
+		}
+		dh = H5Dcreate2(gh, "spectrum_weights", H5T_NATIVE_DOUBLE, sh,
+			        H5P_DEFAULT, H5S_ALL, H5P_DEFAULT);
+		if ( dh < 0 ) {
+			H5Fclose(fh);
+			return 1;
+		}
+		r = H5Dwrite(dh, H5T_NATIVE_DOUBLE, H5S_ALL,
+			     H5S_ALL, H5P_DEFAULT, arr);
 
-	size[0] = 1;
-	sh = H5Screate_simple(1, size, NULL);
+		H5Dclose(dh);
+		free(arr);
 
-	dh = H5Dcreate2(gh, "number_of_samples", H5T_NATIVE_INT, sh,
-	                H5P_DEFAULT, H5S_ALL, H5P_DEFAULT);
-	if ( dh < 0 ) {
-		H5Fclose(fh);
-		return 1;
-	}
+		size[0] = 1;
+		sh = H5Screate_simple(1, size, NULL);
+
+		dh = H5Dcreate2(gh, "number_of_samples", H5T_NATIVE_INT, sh,
+			        H5P_DEFAULT, H5S_ALL, H5P_DEFAULT);
+		if ( dh < 0 ) {
+			H5Fclose(fh);
+			return 1;
+		}
 
-	r = H5Dwrite(dh, H5T_NATIVE_INT, H5S_ALL,
-	             H5S_ALL, H5P_DEFAULT, &image->nsamples);
+		r = H5Dwrite(dh, H5T_NATIVE_INT, H5S_ALL,
+			     H5S_ALL, H5P_DEFAULT, &image->nsamples);
 
-	H5Dclose(dh);
+		H5Dclose(dh);
+	}
 
 	H5Gclose(gh);
 
diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h
index 0c189cad..6b49ad91 100644
--- a/libcrystfel/src/image.h
+++ b/libcrystfel/src/image.h
@@ -46,7 +46,8 @@
 
 typedef enum {
 	SPECTRUM_TOPHAT,
-	SPECTRUM_SASE
+	SPECTRUM_SASE,
+	SPECTRUM_TWOCOLOUR
 } SpectrumType;
 
 /* Structure describing a feature in an image */
diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c
index 1b2193fa..b86399e9 100644
--- a/libcrystfel/src/integration.c
+++ b/libcrystfel/src/integration.c
@@ -3,11 +3,11 @@
  *
  * Integration of intensities
  *
- * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
  *                  a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2013 Thomas White <taw@physics.org>
+ *   2010-2014 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -1483,6 +1483,15 @@ static void integrate_prof2d(IntegrationMethod meth, Crystal *cr,
 	setup_profile_boxes(&ic, list);
 	calculate_reference_profiles(&ic);
 
+	for ( i=0; i<ic.n_reference_profiles; i++ ) {
+		if ( ic.n_profiles_in_reference[i] == 0 ) {
+			ERROR("Reference profile %i has no contributions.\n",
+			      i);
+			free_intcontext(&ic);
+			return;
+		}
+	}
+
 	for ( i=0; i<ic.n_boxes; i++ ) {
 		struct peak_box *bx;
 		bx = &ic.boxes[i];
diff --git a/scripts/Rsplit_surface.py b/scripts/Rsplit_surface.py
index f50e0fd7..ca68a92e 100644
--- a/scripts/Rsplit_surface.py
+++ b/scripts/Rsplit_surface.py
@@ -59,7 +59,7 @@ for i in steps :
 	filename='rsplit'+str(i)+'.dat'
 
 	file = open(filename, 'r')
-	
+
 	lines = array(file.readlines())
 	N = lines.size
 
@@ -80,24 +80,25 @@ Z = griddata(x,y,z,X,Y)
 
 cm = ax.pcolormesh(X, Y, Z, cmap = plt.get_cmap('spectral'), vmin=0, vmax=100)
 
-plt.title(r"$R_{split}$ surface for %s data $N_0$ = %d" % (data_name, N0), fontsize=16)
+plt.title(r"$R_{split}$ surface for %s" % data_name, fontsize=16)
 
-ax.set_xlabel(r"1/d / nm^-1", fontsize=16)
+ax.set_xlabel(r"Resolution (1/d) / $nm^{-1}$", fontsize=16)
 ax.set_ylabel(r"Number of patterns", fontsize=16)
 ypoints = [50000,40000,30000,20000,10000,9000,8000,7000,6000,5000,4000,3000,2000,1000,900,800,700,600,500,400,300,200,100]
 temp = []
 for ypoint in ypoints : temp.append(N0/ypoint)
 ax.set_yticks(log2(temp)+1)
 labels = []
-for ypoint in ypoints : 
+for ypoint in ypoints :
 	if ypoint in [50000,10000,5000,1000,500] : labels.append(str(ypoint))
 	else : labels.append("")
-	
+
 ax.set_yticklabels(labels)
 
 ax.set_ylim(log2(temp[0])+1,log2(temp[-1])+1)
 
-fig.colorbar(cm, shrink=0.5, aspect=10)
+cb = fig.colorbar(cm, shrink=0.5, aspect=10)
+cb.set_label("Rsplit / \%")
 
 plt.savefig(data_name[:-7] + '.png')
 #plt.savefig(data_name[:-7] + '.ps')
diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c
index 1ed008e8..b5bd2e88 100644
--- a/src/diffraction-gpu.c
+++ b/src/diffraction-gpu.c
@@ -325,7 +325,7 @@ void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
 	double tot = 0.0;
 	for ( i=0; i<image->nsamples; i++ ) {
 
-		printf("%.1f eV, weight = %.5f\n",
+		printf("%.3f eV, weight = %.5f\n",
 		       ph_lambda_to_eV(1.0/image->spectrum[i].k),
 		       image->spectrum[i].weight);
 
diff --git a/src/diffraction.c b/src/diffraction.c
index 92e18521..0fa9f79b 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -434,7 +434,8 @@ static int compare_samples(const void *a, const void *b)
 
 
 static struct sample *get_gaussian_spectrum(double eV_cen, double eV_step,
-                                            double sigma, int spec_size)
+                                            double sigma, int spec_size,
+                                            double eV_start)
 {
 	struct sample *spectrum;
 	int i;
@@ -443,7 +444,12 @@ static struct sample *get_gaussian_spectrum(double eV_cen, double eV_step,
 	spectrum = malloc(spec_size * sizeof(struct sample));
 	if ( spectrum == NULL ) return NULL;
 
-	eV = eV_cen - (spec_size/2)*eV_step;
+	if (eV_start == 0) { /* eV starts at 3 sigma below the mean*/
+		eV = eV_cen - (spec_size/2)*eV_step;
+	} else {
+		eV = eV_start;
+	}
+
 	for ( i=0; i<spec_size; i++ ) {
 
 		spectrum[i].k = 1.0/ph_eV_to_lambda(eV);
@@ -568,7 +574,7 @@ struct sample *generate_SASE(struct image *image, gsl_rng *rng)
 	 * points */
 	double eV_step = 6.0*sigma/(spec_size-1);
 
-	spectrum = get_gaussian_spectrum(eV_cen, eV_step, sigma, spec_size);
+	spectrum = get_gaussian_spectrum(eV_cen, eV_step, sigma, spec_size,0);
 
 	/* Add SASE-type noise to Gaussian spectrum */
 	add_sase_noise(spectrum, spec_size, rng);
@@ -592,6 +598,70 @@ struct sample *generate_SASE(struct image *image, gsl_rng *rng)
 }
 
 
+struct sample *generate_twocolour(struct image *image)
+{
+	struct sample *spectrum;
+	struct sample *spectrum1;
+	struct sample *spectrum2;
+	int i;
+	double eV_cen1;  /* Central photon energy for first colour */
+	double eV_cen2;  /* Central photon energy for second colour */
+	double eV_cen;  /* Central photon energy for this spectrum */
+	const int spec_size = 1024;
+
+	eV_cen = ph_lambda_to_eV(image->lambda);
+
+	double halfwidth = eV_cen*image->bw/2.0; /* eV */
+
+	eV_cen1 = eV_cen - halfwidth;
+	eV_cen2 = eV_cen + halfwidth;
+
+	/* Hard-code sigma to be 1/5 of bandwidth */
+	double sigma = eV_cen*image->bw/5.0; /* eV */
+
+	/* The spectrum will be calculated to a resolution which spreads six
+	* sigmas of the original (no SASE noise) Gaussian pulse over spec_size
+	* points */
+	double eV_start = eV_cen1 - 3*sigma;
+	double eV_end = eV_cen2 + 3*sigma;
+	double eV_step = (eV_end - eV_start)/(spec_size-1);
+
+	spectrum1 = get_gaussian_spectrum(eV_cen1, eV_step, sigma, spec_size,
+	                                  eV_start);
+	spectrum2 = get_gaussian_spectrum(eV_cen2, eV_step, sigma, spec_size,
+	                                  eV_start);
+
+	spectrum = malloc(spec_size * sizeof(struct sample));
+	if ( spectrum == NULL ) return NULL;
+
+	for ( i=0; i<spec_size; i++ ) {
+		spectrum[i].weight = spectrum1[i].weight + spectrum2[i].weight;
+		spectrum[i].k = spectrum1[i].k;
+		if ( spectrum2[i].k != spectrum1[i].k ) {
+			printf("%e %e\n", spectrum1[i].k, spectrum2[i].k);
+		}
+	}
+
+	/* Normalise intensity (before taking restricted number of samples) */
+	double total_weight = 0.0;
+	for ( i=0; i<spec_size; i++ ) {
+		total_weight += spectrum[i].weight;
+	}
+
+	for ( i=0; i<spec_size; i++ ) {
+		spectrum[i].weight /= total_weight;
+	}
+
+	/* Sort samples in spectrum by weight.  Diffraction calculation will
+	* take the requested number, starting from the brightest */
+	qsort(spectrum, spec_size, sizeof(struct sample), compare_samples);
+
+	image->spectrum_size = spec_size;
+
+	return spectrum;
+}
+
+
 void get_diffraction(struct image *image, int na, int nb, int nc,
                      const double *intensities, const double *phases,
                      const unsigned char *flags, UnitCell *cell,
diff --git a/src/diffraction.h b/src/diffraction.h
index 5b67c198..e397bd3a 100644
--- a/src/diffraction.h
+++ b/src/diffraction.h
@@ -57,4 +57,6 @@ 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);
+
 #endif	/* DIFFRACTION_H */
diff --git a/src/partial_sim.c b/src/partial_sim.c
index 97e2d14d..5b05218f 100644
--- a/src/partial_sim.c
+++ b/src/partial_sim.c
@@ -169,6 +169,46 @@ static void calculate_partials(Crystal *cr,
 }
 
 
+static void draw_and_write_image(struct image *image, RefList *reflections,
+                                 gsl_rng *rng)
+{
+	Reflection *refl;
+	RefListIterator *iter;
+	const int background = 3000;
+	int i;
+
+	image->data = calloc(image->width*image->height, sizeof(float));
+
+	for ( refl = first_refl(reflections, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) )
+	{
+		double Ip;
+		double dfs, dss;
+		int fs, ss;
+
+		Ip = get_intensity(refl);
+
+		get_detector_pos(refl, &dfs, &dss);
+		fs = nearbyint(dfs);
+		ss = nearbyint(dss);
+		assert(fs >= 0);
+		assert(ss >= 0);
+		assert(fs < image->width);
+		assert(ss < image->height);
+		image->data[fs + image->width*ss] += Ip;
+
+	}
+
+	for ( i=0; i<image->width*image->height; i++ ) {
+		image->data[i] += poisson_noise(rng, background);
+	}
+
+	hdf5_write_image(image->filename, image);
+	free(image->data);
+}
+
+
 static void show_help(const char *s)
 {
 	printf("Syntax: %s [options]\n\n", s);
@@ -181,6 +221,7 @@ static void show_help(const char *s)
 " -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"
 " -b, --beam=<file>        Get beam parameters from file\n"
 " -p, --pdb=<file>         PDB file from which to get the unit cell.\n"
@@ -221,6 +262,8 @@ struct queue_args
 	struct image *template_image;
 	double max_q;
 
+	char *image_prefix;
+
 	/* The overall histogram */
 	double p_hist[NBINS];
 	unsigned long int n_ref[NBINS];
@@ -241,6 +284,8 @@ struct worker_args
 	double p_hist[NBINS];
 	unsigned long int n_ref[NBINS];
 	double p_max[NBINS];
+
+	int n;
 };
 
 
@@ -258,6 +303,7 @@ static void *create_job(void *vqargs)
 	wargs->image = *qargs->template_image;
 
 	qargs->n_started++;
+	wargs->n = qargs->n_started;
 
 	return wargs;
 }
@@ -293,7 +339,13 @@ static void run_job(void *vwargs, int cookie)
 	orientation = random_quaternion(qargs->rngs[cookie]);
 	crystal_set_cell(cr, cell_rotate(qargs->cell, orientation));
 
-	snprintf(wargs->image.filename, 255, "dummy.h5");
+	if ( qargs->image_prefix != NULL ) {
+		snprintf(wargs->image.filename, 255, "%s%i.h5",
+		         qargs->image_prefix, wargs->n);
+	} else {
+		snprintf(wargs->image.filename, 255, "dummy.h5");
+	}
+
 	reflections = find_intersections(&wargs->image, cr);
 	crystal_set_reflections(cr, reflections);
 
@@ -310,6 +362,11 @@ static void run_job(void *vwargs, int cookie)
 	                   qargs->max_q, qargs->full_stddev,
 	                   qargs->noise_stddev, qargs->rngs[cookie]);
 
+	if ( qargs->image_prefix != NULL ) {
+		draw_and_write_image(&wargs->image, reflections,
+		                     qargs->rngs[cookie]);
+	}
+
 	/* Give a slightly incorrect cell in the stream */
 	mess_up_cell(cr, qargs->cnoise, qargs->rngs[cookie]);
 
@@ -372,6 +429,7 @@ int main(int argc, char *argv[])
 	double noise_stddev = 20.0;
 	gsl_rng *rng_for_seeds;
 	int config_random = 0;
+	char *image_prefix = NULL;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -389,6 +447,7 @@ int main(int argc, char *argv[])
 		{"osf-stddev",         1, NULL,                3},
 		{"full-stddev",        1, NULL,                4},
 		{"noise-stddev",       1, NULL,                5},
+		{"images",             1, NULL,                6},
 
 		{"really-random",      0, &config_random,      1},
 
@@ -492,6 +551,10 @@ int main(int argc, char *argv[])
 			}
 			break;
 
+			case 6 :
+			image_prefix = strdup(optarg);
+			break;
+
 			case 0 :
 			break;
 
@@ -595,8 +658,8 @@ int main(int argc, char *argv[])
 	free(output_file);
 
 	image.det = det;
-	image.width = det->max_fs;
-	image.height = det->max_ss;
+	image.width = det->max_fs + 1;
+	image.height = det->max_ss + 1;
 
 	image.lambda = ph_en_to_lambda(eV_to_J(beam->photon_energy));
 	image.div = beam->divergence;
@@ -609,6 +672,7 @@ int main(int argc, char *argv[])
 	image.indexed_by = INDEXING_SIMULATION;
 	image.num_peaks = 0;
 	image.num_saturated_peaks = 0;
+	image.spectrum_size = 0;
 
 	if ( random_intensities ) {
 		full = reflist_new();
@@ -629,6 +693,7 @@ int main(int argc, char *argv[])
 	qargs.full_stddev = full_stddev;
 	qargs.noise_stddev = noise_stddev;
 	qargs.max_q = largest_q(&image);
+	qargs.image_prefix = image_prefix;
 
 	qargs.rngs = malloc(n_threads * sizeof(gsl_rng *));
 	if ( qargs.rngs == NULL ) {
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 08fc75dd..a747cd72 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -54,6 +54,7 @@
 #include "reflist.h"
 #include "reflist-utils.h"
 #include "pattern_sim.h"
+#include "stream.h"
 
 
 static void show_help(const char *s)
@@ -90,6 +91,7 @@ static void show_help(const char *s)
 " -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"
 );
 }
 
@@ -251,6 +253,8 @@ int main(int argc, char *argv[])
 	int nsamples = 3;
 	gsl_rng *rng;
 	int background = 0;
+	char *template_file = NULL;
+	Stream *st = NULL;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -276,6 +280,7 @@ int main(int argc, char *argv[])
 		{"min-size",           1, NULL,                3},
 		{"max-size",           1, NULL,                4},
 		{"background",         1, NULL,                5},
+		{"template",           1, NULL,                6},
 		{0, 0, NULL, 0}
 	};
 
@@ -375,6 +380,10 @@ int main(int argc, char *argv[])
 			}
 			break;
 
+			case 6 :
+			template_file = strdup(optarg);
+			break;
+
 			case 0 :
 			break;
 
@@ -406,6 +415,19 @@ int main(int argc, char *argv[])
 		}
 	}
 
+	if ( template_file != NULL ) {
+		if ( config_randomquat ) {
+			ERROR("You cannot use -r and --template together.\n");
+			return 1;
+		}
+		st = open_stream_for_read(template_file);
+		if ( st == NULL ) {
+			ERROR("Failed to open stream.\n");
+			return 1;
+		}
+		free(template_file);
+	}
+
 	if ( sym_str == NULL ) sym_str = strdup("1");
 	sym = get_pointgroup(sym_str);
 	/* sym_str is used below */
@@ -451,6 +473,11 @@ int main(int argc, char *argv[])
 		spectrum_type = SPECTRUM_TOPHAT;
 	} 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 {
 		ERROR("Unrecognised spectrum type '%s'\n", spectrum_str);
 		return 1;
@@ -578,23 +605,62 @@ int main(int argc, char *argv[])
 
 		}
 
-		/* Read quaternion from stdin */
-		if ( config_randomquat ) {
-			orientation = random_quaternion(rng);
+		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 {
-			orientation = read_quaternion();
-		}
 
-		STATUS("Orientation is %5.3f %5.3f %5.3f %5.3f\n",
-		       orientation.w, orientation.x,
-		       orientation.y, orientation.z);
+			struct image image;
+			int i;
+			Crystal *cr;
 
-		if ( !quaternion_valid(orientation) ) {
-			ERROR("Orientation modulus is not zero!\n");
-			return 1;
-		}
+			image.det = NULL;
+
+			/* Get data from next chunk */
+			if ( read_chunk(st, &image) ) break;
+
+			free(image.filename);
+			image_feature_list_free(image.features);
+
+			if ( image.n_crystals == 0 ) continue;
+
+			cr = image.crystals[0];
+			cell = crystal_get_cell(cr);
+
+			if ( image.n_crystals > 1 ) {
+				ERROR("Using the first crystal only.\n");
+			}
+
+			for ( i=1; i<image.n_crystals; i++ ) {
+
+				Crystal *cr = image.crystals[i];
+				cell = crystal_get_cell(cr);
 
-		cell = cell_rotate(input_cell, orientation);
+				reflist_free(crystal_get_reflections(cr));
+				cell_free(crystal_get_cell(cr));
+				crystal_free(cr);
+
+			}
+
+			free(image.crystals);
+
+		}
 
 		switch ( spectrum_type ) {
 
@@ -606,6 +672,10 @@ int main(int argc, char *argv[])
 			image.spectrum = generate_SASE(&image, rng);
 			break;
 
+                        case SPECTRUM_TWOCOLOUR :
+                        image.spectrum = generate_twocolour(&image);
+                        break;
+
 		}
 
 		/* Ensure no residual information */