Use Spectrum API for simulation

13 files changed, 223 insertions, 479 deletions

diff --git a/libcrystfel/src/hdf5-file.c b/libcrystfel/src/hdf5-file.c
index ff96adb1..f88265a7 100644
--- a/libcrystfel/src/hdf5-file.c
+++ b/libcrystfel/src/hdf5-file.c
@@ -901,29 +901,32 @@ static void write_photon_energy(hid_t fh, double eV, const char *ph_en_loc)
 }
 
 
-static void write_spectrum(hid_t fh, struct sample *spectrum, int spectrum_size,
-                           int nsamples)
+static void write_spectrum(hid_t fh, Spectrum *s)
 {
 	herr_t r;
 	double *arr;
 	int i;
-	hsize_t size1d[1];
 	hid_t sh, dh, ph;
+	double kmin, kmax, step;
+	const hsize_t n = 1024;
 
 	ph = H5Pcreate(H5P_LINK_CREATE);
 	H5Pset_create_intermediate_group(ph, 1);
 
-	arr = malloc(spectrum_size*sizeof(double));
+	arr = malloc(n*sizeof(double));
 	if ( arr == NULL ) {
 		ERROR("Failed to allocate memory for spectrum.\n");
 		return;
 	}
-	for ( i=0; i<spectrum_size; i++ ) {
-		arr[i] = 1.0e10/spectrum[i].k;
+
+	/* Save the wavelength values */
+	spectrum_get_range(s, &kmin, &kmax);
+	step = (kmax-kmin)/n;
+	for ( i=0; i<n; i++ ) {
+		arr[i] = 1.0e10/(kmin+i*step);
 	}
 
-	size1d[0] = spectrum_size;
-	sh = H5Screate_simple(1, size1d, NULL);
+	sh = H5Screate_simple(1, &n, NULL);
 
 	dh = H5Dcreate2(fh, "/spectrum/wavelengths_A", H5T_NATIVE_DOUBLE,
 	                sh, ph, H5S_ALL, H5P_DEFAULT);
@@ -939,46 +942,27 @@ static void write_spectrum(hid_t fh, struct sample *spectrum, int spectrum_size,
 	}
 	H5Dclose(dh);
 
-	for ( i=0; i<spectrum_size; i++ ) {
-		arr[i] = spectrum[i].weight;
+	/* Save the probability density values */
+	for ( i=0; i<n; i++ ) {
+		arr[i] = spectrum_get_density_at_k(s, kmin+i*step);
 	}
 
-	dh = H5Dcreate2(fh, "/spectrum/weights", H5T_NATIVE_DOUBLE, sh,
+	dh = H5Dcreate2(fh, "/spectrum/pdf", H5T_NATIVE_DOUBLE, sh,
 		        H5P_DEFAULT, H5S_ALL, H5P_DEFAULT);
 	if ( dh < 0 ) {
-		ERROR("Failed to create dataset for spectrum weights.\n");
+		ERROR("Failed to create dataset for spectrum p.d.f.\n");
 		return;
 	}
 	r = H5Dwrite(dh, H5T_NATIVE_DOUBLE, H5S_ALL,
 		     H5S_ALL, H5P_DEFAULT, arr);
 	if ( r < 0 ) {
-		ERROR("Failed to write spectrum weights.\n");
-		return;
-	}
-
-	H5Dclose(dh);
-	free(arr);
-
-	size1d[0] = 1;
-	sh = H5Screate_simple(1, size1d, NULL);
-
-	dh = H5Dcreate2(fh, "/spectrum/number_of_samples", H5T_NATIVE_INT, sh,
-		        ph, H5S_ALL, H5P_DEFAULT);
-	if ( dh < 0 ) {
-		ERROR("Failed to create dataset for number of spectrum "
-		      "samples.\n");
-		return;
-	}
-
-	r = H5Dwrite(dh, H5T_NATIVE_INT, H5S_ALL,
-		     H5S_ALL, H5P_DEFAULT, &nsamples);
-	if ( r < 0 ) {
-		ERROR("Failed to write number of spectrum samples.\n");
+		ERROR("Failed to write spectrum p.d.f.\n");
 		return;
 	}
 
 	H5Dclose(dh);
 	H5Pclose(ph);
+	free(arr);
 }
 
 
@@ -1025,10 +1009,8 @@ int hdf5_write_image(const char *filename, const struct image *image,
 
 	write_photon_energy(fh, ph_lambda_to_eV(image->lambda), ph_en_loc);
 
-	if ( image->spectrum0 != NULL && image->spectrum_size > 0 ) {
-
-		write_spectrum(fh, image->spectrum0, image->spectrum_size,
-		              image->nsamples);
+	if ( image->spectrum != NULL ) {
+		write_spectrum(fh, image->spectrum);
 	}
 
 	H5Fclose(fh);
diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h
index 56f07132..706b06b1 100644
--- a/libcrystfel/src/image.h
+++ b/libcrystfel/src/image.h
@@ -3,11 +3,11 @@
  *
  * Handle images and image features
  *
- * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2019 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2009-2018 Thomas White <taw@physics.org>
+ *   2009-2019 Thomas White <taw@physics.org>
  *   2014      Valerio Mariani
  *
  *
@@ -54,6 +54,7 @@ struct imagefile_field_list;
 #include "crystal.h"
 #include "index.h"
 #include "events.h"
+#include "spectrum.h"
 
 /**
  * \file image.h
@@ -103,23 +104,6 @@ enum imagefile_type
 typedef struct _imagefeaturelist ImageFeatureList;
 
 
-struct spectrum
-{
-	int n;
-	double *ks; /* 1/m */
-	double *weights;
-};
-
-
-/** Structure describing a wavelength sample from a spectrum.
- * \deprecated Use struct spectrum instead. */
-struct sample
-{
-	double k;         /**< Wavevector in m^-1 */
-	double weight;    /**< Relative weight */
-};
-
-
 struct beam_params
 {
 	double photon_energy;       /**< eV per photon */
@@ -182,17 +166,8 @@ struct image
 	/** Monotonically increasing serial number for this image */
 	int                     serial;
 
-	/** Spectrum information (new format) */
-	struct spectrum *spectrum;
-
-	/** \name Spectrum information (old format)
-	 * @{
-	 * Array of samples, number of samples, and size of the array (may be
-	 * larger than nsamples) */
-	struct sample           *spectrum0;
-	int                     nsamples;
-	int                     spectrum_size;
-	/** @} */
+	/** Spectrum information */
+	Spectrum               *spectrum;
 
 	/** Wavelength of the incident radiation, in metres */
 	double                  lambda;
diff --git a/libcrystfel/src/spectrum.c b/libcrystfel/src/spectrum.c
index f1ada734..e083962c 100644
--- a/libcrystfel/src/spectrum.c
+++ b/libcrystfel/src/spectrum.c
@@ -451,3 +451,88 @@ Spectrum *spectrum_load(const char *filename)
 	fclose(fh);
 	return s;
 }
+
+
+/**
+ * \param wavelength Wavelength in metres
+ * \param bandwidth Bandwidth as a fraction
+ *
+ * Generates a top-hat spectrum centered on 'wavelength', where the width of the
+ * flat top is bandwidth/wavelength
+ *
+ * \returns A newly-allocated \ref Spectrum, or NULL on error.
+ */
+Spectrum *spectrum_generate_tophat(double wavelength, double bandwidth)
+{
+	Spectrum *s;
+	double kvals[2];
+	double samp[2];
+	double kcen;
+
+	s = spectrum_new();
+	if ( s == NULL ) return NULL;
+
+	kcen = 1.0/wavelength;
+	kvals[0] = kcen - kcen*bandwidth/2.0;
+	kvals[1] = kcen + kcen*bandwidth/2.0;
+	samp[0] = 1.0;
+	samp[1] = 1.0;
+	spectrum_set_pdf(s, kvals, samp, 2);
+	return s;
+}
+
+
+/**
+ * \param wavelength Wavelength in metres
+ * \param bandwidth Bandwidth as a fraction
+ *
+ * Generates a Gaussian spectrum centered on 'wavelength', where the standard
+ * deviation of the Gaussian is bandwidth/wavelength
+ *
+ * \returns A newly-allocated \ref Spectrum, or NULL on error.
+ */
+Spectrum *spectrum_generate_gaussian(double wavelength, double bandwidth)
+{
+	Spectrum *s;
+	struct gaussian g;
+
+	s = spectrum_new();
+	if ( s == NULL ) return NULL;
+
+	g.kcen = 1.0/wavelength;
+	g.sigma = bandwidth/wavelength;
+	g.height = 1;
+	spectrum_set_gaussians(s, &g, 1);
+
+	return s;
+}
+
+
+/**
+ * \param wavelength Wavelength in metres
+ * \param bandwidth Bandwidth as a fraction
+ *
+ * Generates a top-hat spectrum centered on 'wavelength', where the width of the
+ * flat top is bandwidth/wavelength
+ *
+ * \returns A newly-allocated \ref Spectrum, or NULL on error.
+ */
+Spectrum *spectrum_generate_sase(double wavelength, double bandwidth,
+                                 double spike_width, gsl_rng *rng)
+{
+}
+
+
+/**
+ * \param wavelength Wavelength in metres
+ * \param bandwidth Bandwidth as a fraction
+ *
+ * Generates a top-hat spectrum centered on 'wavelength', where the width of the
+ * flat top is bandwidth/wavelength
+ *
+ * \returns A newly-allocated \ref Spectrum, or NULL on error.
+ */
+Spectrum *spectrum_generate_twocolour(double wavelength, double bandwidth,
+                                      double separation)
+{
+}
diff --git a/libcrystfel/src/spectrum.h b/libcrystfel/src/spectrum.h
index d31318fd..8298a4cd 100644
--- a/libcrystfel/src/spectrum.h
+++ b/libcrystfel/src/spectrum.h
@@ -33,6 +33,7 @@
 #include <config.h>
 #endif
 
+#include <gsl/gsl_rng.h>
 
 /**
  * \file spectrum.h
@@ -78,6 +79,13 @@ extern void spectrum_set_pdf(Spectrum *s, double *kcens, double *heights,
 extern void spectrum_get_range(Spectrum *s, double *kmin, double *kmax);
 extern double spectrum_get_density_at_k(Spectrum *s, double k);
 
+/* Generation of spectra */
+extern Spectrum *spectrum_generate_tophat(double wavelength, double bandwidth);
+extern Spectrum *spectrum_generate_gaussian(double wavelength, double bandwidth);
+extern Spectrum *spectrum_generate_sase(double wavelength, double bandwidth,
+                                        double spike_width, gsl_rng *rng);
+extern Spectrum *spectrum_generate_twocolour(double wavelength, double bandwidth,
+                                             double separation);
 
 #ifdef __cplusplus
 }
diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c
index 9cbfdf33..e2bce0cc 100644
--- a/src/diffraction-gpu.c
+++ b/src/diffraction-gpu.c
@@ -283,7 +283,7 @@ static int do_panels(struct gpu_context *gctx, struct image *image,
 
 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 no_fringes, int flat, int n_samples)
 {
 	double ax, ay, az;
 	double bx, by, bz;
@@ -294,6 +294,8 @@ int get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
 	int n_neg = 0;
 	int n_nan = 0;
 	int i;
+	double kmin, kmax, step;
+	double tot = 0.0;
 
 	if ( gctx == NULL ) {
 		ERROR("GPU setup failed.\n");
@@ -338,23 +340,25 @@ int get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
 		}
 	}
 
-	double tot = 0.0;
-	for ( i=0; i<image->nsamples; i++ ) {
+	spectrum_get_range(image->spectrum, &kmin, &kmax);
+	step = (kmax-kmin)/n_samples;
+	for ( i=0; i<=n_samples; i++ ) {
+
+		double k = kmin + i*step;
+		double prob;
 
-		printf("%.3f eV, weight = %.5f\n",
-		       ph_lambda_to_eV(1.0/image->spectrum0[i].k),
-		       image->spectrum0[i].weight);
+		/* Probability = p.d.f. times step width */
+		prob = step * spectrum_get_density_at_k(image->spectrum, k);
 
-		err = do_panels(gctx, image, image->spectrum0[i].k,
-		                image->spectrum0[i].weight,
-		                &n_inf, &n_neg, &n_nan);
+		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;
 
-		tot += image->spectrum0[i].weight;
+		tot += prob;
 
 	}
-	printf("total weight = %f\n", tot);
+	STATUS("Total weight = %f\n", tot);
 
 	if ( n_neg + n_inf + n_nan ) {
 		ERROR("WARNING: The GPU calculation produced %i negative"
diff --git a/src/diffraction-gpu.h b/src/diffraction-gpu.h
index bef12a8e..a4d16583 100644
--- a/src/diffraction-gpu.h
+++ b/src/diffraction-gpu.h
@@ -42,7 +42,7 @@ struct gpu_context;
 
 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 no_fringes, int flat, int n_samples);
 extern struct gpu_context *setup_gpu(int no_sfac,
                                      const double *intensities,
                                      const unsigned char *flags,
@@ -53,7 +53,7 @@ extern void cleanup_gpu(struct gpu_context *gctx);
 
 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 no_fringes, int flat, int n_samples)
 {
 	/* Do nothing */
 	ERROR("This copy of CrystFEL was not compiled with OpenCL support.\n");
diff --git a/src/diffraction.c b/src/diffraction.c
index 4dc21c59..bcbc62e8 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -443,317 +443,11 @@ static void diffraction_at_k(struct image *image, const double *intensities,
 }
 
 
-static void normalise_sampled_spectrum(struct sample *spec, int n, int nsamp)
-{
-	int i;
-	double total_w = 0.0;
-	double samp_w = 0.0;
-
-	for ( i=0; i<n; i++ ) {
-		total_w += spec[i].weight;
-	}
-	STATUS("%i energies in spectrum, %i samples requested.\n", n, nsamp);
-
-	for ( i=0; i<nsamp; i++ ) {
-		samp_w += spec[i].weight;
-	}
-
-	if ( samp_w < 0.8 * total_w ) {
-		ERROR("WARNING: Only %.1f %% of the calculated spectrum "
-		      "intensity was sampled.\n", 100.0 * samp_w / total_w);
-		ERROR("I've increased the weighting of the sampled points to "
-		      "keep the final intensity constant, but the spectrum "
-		      "might be very far from accurate.  Consider increasing "
-		      "the number of spectrum samples.\n");
-	} else {
-		STATUS("%.1f %% of calculated spectrum intensity sampled.\n",
-		       100.0 * samp_w / total_w);
-		STATUS("Normalised to keep total intensity constant.\n");
-	}
-
-	for ( i=0; i<n; i++ ) {
-		spec[i].weight /= samp_w;
-	}
-}
-
-
-static int compare_samples(const void *a, const void *b)
-{
-	struct sample *sample1 = (struct sample *)a;
-	struct sample *sample2 = (struct sample *)b;
-	if ( sample1->weight < sample2->weight ) {
-		return 1;
-	}
-	return -1;
-}
-
-
-static struct sample *get_gaussian_spectrum(double eV_cen, double eV_step,
-                                            double sigma, int spec_size,
-                                            double eV_start)
-{
-	struct sample *spectrum;
-	int i;
-	double eV;
-
-	spectrum = malloc(spec_size * sizeof(struct sample));
-	if ( spectrum == NULL ) return NULL;
-
-	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);
-		spectrum[i].weight = exp(-(pow(eV - eV_cen, 2.0)
-		                         / (2.0*sigma*sigma)));
-		eV += eV_step;
-	}
-
-	return spectrum;
-}
-
-
-static int add_sase_noise(struct sample *spectrum, int nsteps, gsl_rng *rng)
-{
-	struct sample *noise;
-	int i, j;
-	double *gaussianNoise;
-	int shiftLim = 5;
-	double noise_mean = 0.0;
-	double noise_sigma = 1.0;
-
-	if ( shiftLim > nsteps ) shiftLim = nsteps;
-
-	noise = malloc(nsteps * sizeof(struct sample));
-	if ( noise == NULL ) return 1;
-
-	gaussianNoise = malloc(3 * nsteps * sizeof(double));
-	if ( gaussianNoise == NULL ) {
-		free(noise);
-		return 1;
-	}
-
-	/* Generate Gaussian noise of length of spectrum
-	 * (replicate on both ends for circshift below) */
-	for ( i=0; i<nsteps; i++) {
-
-		noise[i].weight = 0.0;
-
-		/* Gaussian noise with mean = 0, std = 1 */
-		gaussianNoise[i] = gaussian_noise(rng, noise_mean, noise_sigma);
-		gaussianNoise[i+nsteps] = gaussianNoise[i];
-		gaussianNoise[i+2*nsteps] = gaussianNoise[i];
-	}
-
-	/* Sum Gaussian noise by circshifting by +/- shiftLim */
-	for ( i=nsteps; i<2*nsteps; i++ ) {
-		for ( j=-shiftLim; j<=shiftLim; j++ ) {
-			noise[i-nsteps].weight += gaussianNoise[i+j];
-		}
-	}
-
-	/* Normalize the number of circshift sum */
-	for ( i=0; i<nsteps; i++) {
-		noise[i].weight = noise[i].weight/(2*shiftLim+1);
-	}
-
-	/* Noise amplitude should have a 2 x Gaussian distribution */
-	for ( i=0; i<nsteps; i++ ) {
-		noise[i].weight = 2.0 * spectrum[i].weight * noise[i].weight;
-	}
-
-	/* Add noise to spectrum */
-	for ( i=0; i<nsteps; i++ ) {
-
-		spectrum[i].weight += noise[i].weight;
-
-		/* The final spectrum can not be negative */
-		if ( spectrum[i].weight < 0.0 ) spectrum[i].weight = 0.0;
-
-	}
-
-	return 0;
-}
-
-
-struct sample *generate_tophat(struct image *image)
-{
-	struct sample *spectrum;
-	int i;
-	double k, k_step;
-
-	double halfwidth = image->bw * image->lambda / 2.0;  /* m */
-	double mink = 1.0/(image->lambda + halfwidth);
-	double maxk = 1.0/(image->lambda - halfwidth);
-
-	spectrum = malloc(image->nsamples * sizeof(struct sample));
-	if ( spectrum == NULL ) return NULL;
-
-	k = mink;
-	k_step = (maxk-mink)/(image->nsamples-1);
-	for ( i=0; i<image->nsamples; i++ ) {
-		spectrum[i].k = k;
-		spectrum[i].weight = 1.0/(double)image->nsamples;
-		k += k_step;
-	}
-
-	image->spectrum_size = image->nsamples;
-	/* No need to call normalise_sampled_spectrum() in this case */
-
-	return spectrum;
-}
-
-
-struct sample *generate_spectrum_fromfile(struct image *image,
-					  char *spectrum_fn)
-{
-	struct sample *spectrum;
-	int i;
-	double k, w;
-	double w_sum = 0;
-
-	spectrum = malloc(image->nsamples * sizeof(struct sample));
-	if ( spectrum == NULL ) return NULL;
-
-	FILE *f;
-	f = fopen(spectrum_fn, "r");
-
-	int nsamples = 0;
-	for ( i=0; i<image->nsamples; i++ ) {
-		if (fscanf(f, "%lf %lf", &k, &w) != EOF) {
-			spectrum[i].k = ph_eV_to_k(k);
-			spectrum[i].weight = w;
-			w_sum += w;
-			nsamples += 1;
-		} else break;
-	}
-
-	for ( i=0; i<nsamples; i++ ) {
-		spectrum[i].weight /= w_sum;
-	}
-
-	image->spectrum_size = nsamples;
-
-	return spectrum;
-}
-
-
-struct sample *generate_SASE(struct image *image, gsl_rng *rng)
-{
-	struct sample *spectrum;
-	const int spec_size = 1024;
-	double eV_cen;  /* Central photon energy for this spectrum */
-	const double jitter_sigma_eV = 8.0;
-
-	/* Central wavelength jitters with Gaussian distribution */
-	eV_cen = gaussian_noise(rng, ph_lambda_to_eV(image->lambda),
-	                        jitter_sigma_eV);
-
-	/* Convert FWHM to standard deviation.  Note that bandwidth is taken
-	 * here to be "delta E over E" (E = photon energy), not the bandwidth in
-	 * terms of wavelength (as it is everywhere else), but the difference
-	 * should be very small */
-	double sigma = (image->bw*eV_cen) / (2.0*sqrt(2.0*log(2.0)));
-
-	/* 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_step = 6.0*sigma/(spec_size-1);
-
-	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);
-
-	/* 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);
-
-	normalise_sampled_spectrum(spectrum, spec_size, image->nsamples);
-
-	image->spectrum_size = spec_size;
-	return spectrum;
-}
-
-
-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);
-
-	/* Convert FWHM to standard deviation.  Note that bandwidth is taken
-	 * here to be "delta E over E" (E = photon energy), not the bandwidth in
-	 * terms of wavelength (as it is everywhere else), but the difference
-	 * should be very small */
-	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);
-
-	normalise_sampled_spectrum(spectrum, spec_size, image->nsamples);
-
-	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,
                      GradientMethod m, const SymOpList *sym,
-                     int no_fringes, int flat)
+                     int no_fringes, int flat, int n_samples)
 {
 	double ax, ay, az;
 	double bx, by, bz;
@@ -762,6 +456,7 @@ void get_diffraction(struct image *image, int na, int nb, int nc,
 	double *lut_b;
 	double *lut_c;
 	int i;
+	double kmin, kmax, step;
 
 	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
 
@@ -769,17 +464,23 @@ void get_diffraction(struct image *image, int na, int nb, int nc,
 	lut_b = get_sinc_lut(nb, no_fringes, flat);
 	lut_c = get_sinc_lut(nc, no_fringes, flat);
 
-	for ( i=0; i<image->nsamples; i++ ) {
+	spectrum_get_range(image->spectrum, &kmin, &kmax);
+	step = (kmax-kmin)/(n_samples+1);
+	STATUS("%i samples\n", n_samples);
+	for ( i=1; i<=n_samples; i++ ) {
 
-		printf("%.1f eV, weight = %.5f\n",
-		       ph_lambda_to_eV(1.0/image->spectrum0[i].k),
-		       image->spectrum0[i].weight);
+		double k = kmin + i*step;
+		double prob;
+
+		/* Probability = p.d.f. times step width */
+		prob = step * spectrum_get_density_at_k(image->spectrum, k);
+
+		STATUS("Wavelength: %e m, weight = %.5f\n", 1.0/k, prob);
 
 		diffraction_at_k(image, intensities, phases,
-		                flags, cell, m, sym, image->spectrum0[i].k,
+		                flags, cell, m, sym, k,
 		                ax, ay, az, bx, by, bz, cx, cy, cz,
-		                lut_a, lut_b, lut_c, image->spectrum0[i].weight);
-
+		                lut_a, lut_b, lut_c, prob);
 
 	}
 
diff --git a/src/diffraction.h b/src/diffraction.h
index 0f95f83c..16b0386f 100644
--- a/src/diffraction.h
+++ b/src/diffraction.h
@@ -52,7 +52,7 @@ 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 no_fringes, int flat, int n_samples);
 
 extern struct sample *generate_tophat(struct image *image);
 
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 681d5d31..e4568966 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -250,66 +250,6 @@ static void add_geom_beam_stuff_to_field_list(struct imagefile_field_list *copym
 }
 
 
-static struct spectrum *read_spectrum_fromfile(char *fn)
-{
-	FILE *f;
-	struct spectrum *s;
-	int i;
-	double k, w;
-	double w_sum = 0;
-
-	f = fopen(fn, "r");
-	if ( f == NULL ) {
-		ERROR("Couldn't open '%s'\n", fn);
-		return NULL;
-	}
-
-	s = malloc(sizeof(struct spectrum));
-	if ( s == NULL ) return NULL;
-
-	if ( fscanf(f, "%d", &s->n) == EOF ) {
-		return NULL;
-	}
-
-	if ( s->n <= 0 ) {
-		return NULL;
-	}
-
-	s->ks = malloc(s->n * sizeof(double));
-	if ( s->ks == NULL ) {
-		ERROR("Failed to allocate spectrum!\n");
-		return NULL;
-	}
-
-	s->weights = malloc(s->n * sizeof(double));
-	if ( s->weights == NULL ) {
-		ERROR("Failed to allocate spectrum!\n");
-		return NULL;
-	}
-
-	for ( i=0; i<s->n; i++ ) {
-		if ( fscanf(f, "%lf %lf", &k, &w) != EOF ) {
-			s->ks[i] = ph_eV_to_k(k);
-			s->weights[i] = w;
-			w_sum += w;
-		} else {
-			break;
-		}
-	}
-
-	if ( i < s->n - 1 ) {
-		ERROR("Failed to read %d lines from %s\n", s->n, fn);
-		return NULL;
-	}
-
-	for ( i=0; i<s->n; i++ ) {
-		s->weights[i] /= w_sum;
-	}
-
-	return s;
-}
-
-
 int main(int argc, char *argv[])
 {
 	int c;
@@ -1164,13 +1104,12 @@ int main(int argc, char *argv[])
 
 	/* Load spectrum from file if given */
 	if ( spectrum_fn != NULL ) {
-		iargs.spectrum = read_spectrum_fromfile(spectrum_fn);
+		iargs.spectrum = spectrum_load(spectrum_fn);
 		if ( iargs.spectrum == NULL ) {
 			ERROR("Couldn't read spectrum (from %s)\n", spectrum_fn);
 			return 1;
 		}
 		free(spectrum_fn);
-		STATUS("Read %d lines from %s\n", iargs.spectrum->n, spectrum_fn);
 	} else {
 		iargs.spectrum = NULL;
 	}
diff --git a/src/partial_sim.c b/src/partial_sim.c
index 92f9703c..0ea45ad4 100644
--- a/src/partial_sim.c
+++ b/src/partial_sim.c
@@ -841,8 +841,7 @@ int main(int argc, char *argv[])
 	image.crystals = NULL;
 	image.n_crystals = 0;
 	image.indexed_by = INDEXING_SIMULATION;
-	image.spectrum_size = 0;
-	image.spectrum0 = NULL;
+	image.spectrum = NULL;
 	image.serial = 0;
 	image.event = NULL;
 
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 5ee53514..d13baff8 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -60,6 +60,7 @@
 
 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 */
@@ -102,7 +103,12 @@ static void show_help(const char *s)
 "     --template=<file>     Take orientations from stream <file>.\n"
 "     --no-fringes          Exclude the side maxima of Bragg peaks.\n"
 "     --flat                Make Bragg peaks flat.\n"
-"     --beam-bandwidth      Beam bandwidth as a fraction. Default 1%%.\n"
+"     --beam-bandwidth      Beam bandwidth (standard deviation of wavelength as\n"
+"                            a fraction of wavelenth).  Default 0.001 (1%%)\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"
 "     --photon-energy       Photon energy in eV.  Default 9000.\n"
 "     --nphotons            Number of photons per X-ray pulse.  Default 1e12.\n"
 "     --beam-radius         Radius of beam in metres (default 1e-6).\n"
@@ -401,6 +407,8 @@ int main(int argc, char *argv[])
 	double beam_radius = 1e-6;  /* metres */
 	double bandwidth = 0.01;
 	double photon_energy = 9000.0;
+	double sase_spike_width = 2e6;  /* m^-1 */
+	double twocol_sep = 8e6;  /* m^-1 */
 	struct beam_params beam;
 	int i;
 
@@ -438,6 +446,8 @@ int main(int argc, char *argv[])
 		{"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}
 	};
@@ -603,6 +613,30 @@ int main(int argc, char *argv[])
 			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;
 
@@ -700,6 +734,8 @@ int main(int argc, char *argv[])
 		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 ||
@@ -789,7 +825,6 @@ int main(int argc, char *argv[])
 
 	image.lambda = ph_en_to_lambda(eV_to_J(photon_energy));
 	image.bw = bandwidth;
-	image.nsamples = nsamples;
 
 	/* Initialise stuff */
 	image.filename = NULL;
@@ -826,7 +861,7 @@ int main(int argc, char *argv[])
 	powder.det = image.det;
 	powder.beam = NULL;
 	powder.lambda = 0.0;
-	powder.spectrum0 = NULL;
+	powder.spectrum = NULL;
 	powder.dp = malloc(image.det->n_panels*sizeof(float *));
 	if ( powder.dp == NULL ) {
 		ERROR("Failed to allocate powder data\n");
@@ -860,6 +895,11 @@ int main(int argc, char *argv[])
 		       "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);
@@ -985,20 +1025,31 @@ int main(int argc, char *argv[])
 		switch ( spectrum_type ) {
 
 			case SPECTRUM_TOPHAT :
-			image.spectrum0 = generate_tophat(&image);
+			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.spectrum0 = generate_SASE(&image, rng);
+			image.spectrum = spectrum_generate_sase(image.lambda,
+			                                        image.bw,
+			                                        sase_spike_width,
+			                                        rng);
 			break;
 
                         case SPECTRUM_TWOCOLOUR :
-                        image.spectrum0 = generate_twocolour(&image);
+                        image.spectrum = spectrum_generate_twocolour(image.lambda,
+                                                                     image.bw,
+                                                                     twocol_sep);
                         break;
 
                         case SPECTRUM_FROMFILE :
-                        image.spectrum0 = generate_spectrum_fromfile(&image,
-								   spectrum_fn);
+                        image.spectrum = spectrum_load(spectrum_fn);
                         break;
 
 		}
@@ -1016,11 +1067,13 @@ int main(int argc, char *argv[])
 				                 gpu_dev);
 			}
 			err = get_diffraction_gpu(gctx, &image, na, nb, nc,
-			                          cell, no_fringes, flat);
+			                          cell, no_fringes, flat,
+			                          nsamples);
 
 		} else {
 			get_diffraction(&image, na, nb, nc, intensities, phases,
-			                flags, cell, grad, sym, no_fringes, flat);
+			                flags, cell, grad, sym, no_fringes, flat,
+			                nsamples);
 		}
 		if ( err ) {
 			ERROR("Diffraction calculation failed.\n");
diff --git a/src/process_image.h b/src/process_image.h
index 2a43d11d..5939ae4d 100644
--- a/src/process_image.h
+++ b/src/process_image.h
@@ -112,7 +112,7 @@ struct index_args
 	struct taketwo_options taketwo_opts;
 	struct xgandalf_options xgandalf_opts;
 	struct felix_options felix_opts;
-	struct spectrum *spectrum;
+	Spectrum *spectrum;
 	signed int wait_for_file; /* -1 means wait forever */
 };
 
diff --git a/tests/gpu_sim_check.c b/tests/gpu_sim_check.c
index f1d7866f..c9508a7c 100644
--- a/tests/gpu_sim_check.c
+++ b/tests/gpu_sim_check.c
@@ -163,21 +163,19 @@ int main(int argc, char *argv[])
 	gpu_image.det = det;
 	cpu_image.beam = NULL;
 	gpu_image.beam = NULL;
-	cpu_image.spectrum0 = NULL;
-	gpu_image.spectrum0 = NULL;
 
 	cpu_image.lambda = ph_en_to_lambda(eV_to_J(6000));
 	gpu_image.lambda = ph_en_to_lambda(eV_to_J(6000));
 	cpu_image.bw = 1.0 / 100.0;
 	gpu_image.bw = 1.0 / 100.0;
 
-	cpu_image.nsamples = 10;
-	gpu_image.nsamples = 10;
-	cpu_image.spectrum0 = generate_tophat(&cpu_image);
-	gpu_image.spectrum0 = generate_tophat(&gpu_image);
+	cpu_image.spectrum = spectrum_generate_tophat(cpu_image.lambda,
+	                                              cpu_image.bw);
+	gpu_image.spectrum = spectrum_generate_tophat(gpu_image.lambda,
+	                                              gpu_image.bw);
 
 	start = get_hires_seconds();
-	if ( get_diffraction_gpu(gctx, &gpu_image, 8, 8, 8, cell, 0, 0) ) {
+	if ( get_diffraction_gpu(gctx, &gpu_image, 8, 8, 8, cell, 0, 0, 10) ) {
 		return 1;
 	}
 	end = get_hires_seconds();
@@ -201,7 +199,7 @@ int main(int argc, char *argv[])
 
 	start = get_hires_seconds();
 	get_diffraction(&cpu_image, 8, 8, 8, NULL, NULL, NULL, cell,
-	                GRADIENT_MOSAIC, sym, 0, 0);
+	                GRADIENT_MOSAIC, sym, 0, 0, 10);
 	end = get_hires_seconds();
 	cpu_time = end - start;