aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--libcrystfel/src/hdf5-file.c58
-rw-r--r--libcrystfel/src/image.h35
-rw-r--r--libcrystfel/src/spectrum.c85
-rw-r--r--libcrystfel/src/spectrum.h8
-rw-r--r--src/diffraction-gpu.c26
-rw-r--r--src/diffraction-gpu.h4
-rw-r--r--src/diffraction.c329
-rw-r--r--src/diffraction.h2
-rw-r--r--src/indexamajig.c63
-rw-r--r--src/partial_sim.c3
-rw-r--r--src/pattern_sim.c73
-rw-r--r--src/process_image.h2
-rw-r--r--tests/gpu_sim_check.c14
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;
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-16 15:42:50 +0000