aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2019-05-24 17:08:04 +0200
committerThomas White <taw@physics.org>2019-05-29 10:42:14 +0200
commit0e1726e01645c110e9935cc8c3e8f35291b2b82d (patch)
treeb2e2fad656dda0dbd36b2e109944bd111a57bb7b
parent581f8e8ab1b1d2833fd689540a2758bea1fc93ee (diff)
Spectrum generation
-rw-r--r--libcrystfel/src/spectrum.c58
-rw-r--r--tests/spectrum_check.c24
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;
}