Use Spectrum API for simulation

4 files changed, 118 insertions, 68 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
 }