diff options
| -rw-r--r-- | libcrystfel/src/spectrum.c | 58 | ||||
| -rw-r--r-- | tests/spectrum_check.c | 24 |
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; } |