Spectrum generation

2 files changed, 75 insertions, 7 deletions

diff --git a/libcrystfel/src/spectrum.c b/libcrystfel/src/spectrum.c
index 576ff948..5cb5c47d 100644
--- a/libcrystfel/src/spectrum.c
+++ b/libcrystfel/src/spectrum.c
@@ -485,7 +485,8 @@ Spectrum *spectrum_generate_tophat(double wavelength, double bandwidth)
  * \param bandwidth Bandwidth as a fraction
  *
  * Generates a Gaussian spectrum centered on 'wavelength', where the standard
- * deviation of the Gaussian is bandwidth/wavelength
+ * deviation of the Gaussian is bandwidth divided by wavelength (bandwidth
+ * fraction times k value).
  *
  * \returns A newly-allocated \ref Spectrum, or NULL on error.
  */
@@ -508,29 +509,72 @@ Spectrum *spectrum_generate_gaussian(double wavelength, double bandwidth)
 
 /**
  * \param wavelength Wavelength in metres
- * \param bandwidth Bandwidth as a fraction
+ * \param bandwidth Bandwidth as a fraction of wavelength
+ * \param spike_width The width of the SASE spikes, as a fraction of wavelength
+ * \param rng A GSL random number generator
  *
- * Generates a top-hat spectrum centered on 'wavelength', where the width of the
- * flat top is bandwidth/wavelength
+ * Generates a SASE spectrum centered on 'wavelength', with 15 spikes of width
+ * spike_width divided by wavelength (i.e. spike_width times k value).  The
+ * spikes will be distributed within a width of k values of bandwidth divided
+ * by wavelength (bandwidth times k-value).
+ *
+ * Note that CrystFEL is not an undulator simulation program.  This function
+ * just simulates a rough idea of the kind of spiky spectrum resulting from the
+ * SASE process.
  *
  * \returns A newly-allocated \ref Spectrum, or NULL on error.
  */
 Spectrum *spectrum_generate_sase(double wavelength, double bandwidth,
                                  double spike_width, gsl_rng *rng)
 {
+	int i;
+	Spectrum *s;
+	struct gaussian g[15];
+
+	s = spectrum_new();
+	if ( s == NULL ) return NULL;
+
+	for ( i=0; i<15; i++ ) {
+		g[0].kcen = 1.0/wavelength + (gsl_rng_uniform_pos(rng)-0.5) * bandwidth/wavelength;
+		g[0].sigma = spike_width/wavelength;
+		g[0].height = 1;
+	}
+
+	spectrum_set_gaussians(s, g, 15);
+
+	return s;
 }
 
 
 /**
  * \param wavelength Wavelength in metres
- * \param bandwidth Bandwidth as a fraction
+ * \param bandwidth Bandwidth as a fraction of wavelength
+ * \param separation Separation between peak centres, in m^-1
  *
- * Generates a top-hat spectrum centered on 'wavelength', where the width of the
- * flat top is bandwidth/wavelength
+ * Generates a two-colour spectrum with Gaussian peaks centered at wavenumbers
+ * 1/wavelength ± separation/2.  Each peak will have a standard deviation of
+ * bandwidth divided by wavelength (bandwidth fraction times k value).
  *
  * \returns A newly-allocated \ref Spectrum, or NULL on error.
  */
 Spectrum *spectrum_generate_twocolour(double wavelength, double bandwidth,
                                       double separation)
 {
+	Spectrum *s;
+	struct gaussian g[2];
+
+	s = spectrum_new();
+	if ( s == NULL ) return NULL;
+
+	g[0].kcen = 1.0/wavelength - separation/2.0;
+	g[0].sigma = bandwidth/wavelength;
+	g[0].height = 1;
+
+	g[1].kcen = 1.0/wavelength + separation/2.0;
+	g[1].sigma = bandwidth/wavelength;
+	g[1].height = 1;
+
+	spectrum_set_gaussians(s, g, 2);
+
+	return s;
 }
diff --git a/tests/spectrum_check.c b/tests/spectrum_check.c
index 7d37689d..e910aad3 100644
--- a/tests/spectrum_check.c
+++ b/tests/spectrum_check.c
@@ -59,12 +59,30 @@ static int check_integral(Spectrum *s, int nsamp)
 }
 
 
+static void plot_spectrum(Spectrum *s)
+{
+	double min, max, step;
+	int i;
+	const int nsamp = 100;
+
+	spectrum_get_range(s, &min, &max);
+	step = (max-min)/nsamp;
+	for ( i=0; i<=nsamp; i++ ) {
+		double x = min+i*step;
+		double y = spectrum_get_density_at_k(s, x);
+		printf("%e %e\n", x, y);
+	}
+}
+
 int main(int argc, char *argv[])
 {
 	Spectrum *s;
 	struct gaussian gauss;
+	gsl_rng *rng;
 	int r = 0;
 
+	rng = gsl_rng_alloc(gsl_rng_mt19937);
+
 	s = spectrum_new();
 	gauss.kcen = ph_eV_to_k(9000);
 	gauss.sigma = ph_eV_to_k(100);
@@ -73,5 +91,11 @@ int main(int argc, char *argv[])
 	r += check_integral(s, 100);
 	spectrum_free(s);
 
+	s = spectrum_generate_sase(ph_eV_to_k(9000), 0.01, 0.0005, rng);
+	plot_spectrum(s);
+	spectrum_free(s);
+
+	gsl_rng_free(rng);
+
 	return r;
 }