aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/diffraction-gpu.c142
-rw-r--r--src/diffraction.c281
-rw-r--r--src/diffraction.h12
-rw-r--r--src/pattern_sim.c63
-rw-r--r--src/pattern_sim.h2
5 files changed, 379 insertions, 121 deletions
diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c
index f41a45ec..ff4f2414 100644
--- a/src/diffraction-gpu.c
+++ b/src/diffraction-gpu.c
@@ -3,11 +3,13 @@
*
* Calculate diffraction patterns by Fourier methods (GPU version)
*
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2012 Thomas White <taw@physics.org>
+ * 2009-2014 Thomas White <taw@physics.org>
+ * 2013 Alexandra Tolstikova
+ * 2013-2014 Chun Hong Yoon <chun.hong.yoon@desy.de>
*
* This file is part of CrystFEL.
*
@@ -168,18 +170,20 @@ static int set_arg_mem(struct gpu_context *gctx, int idx, cl_mem val)
static void do_panels(struct gpu_context *gctx, struct image *image,
- int *n_inf, int *n_neg, int *n_nan, double k, double frac,
- double xo, double yo)
+ double k, double weight,
+ int *n_inf, int *n_neg, int *n_nan)
{
int i;
+ const int sampling = 4; /* This, squared, number of samples / pixel */
if ( set_arg_float(gctx, 1, k) ) return;
+ if ( set_arg_float(gctx, 2, weight) ) return;
/* Iterate over panels */
for ( i=0; i<image->det->n_panels; i++ ) {
size_t dims[2];
- size_t ldims[2] = {1, 1};
+ size_t ldims[2];
struct panel *p;
cl_mem diff;
size_t diff_size;
@@ -203,20 +207,29 @@ static void do_panels(struct gpu_context *gctx, struct image *image,
}
if ( set_arg_mem(gctx, 0, diff) ) return;
- if ( set_arg_int(gctx, 2, pan_width) ) return;
- if ( set_arg_float(gctx, 3, p->cnx) ) return;
- if ( set_arg_float(gctx, 4, p->cny) ) return;
- if ( set_arg_float(gctx, 5, p->res) ) return;
- if ( set_arg_float(gctx, 6, p->clen) ) return;
- if ( set_arg_float(gctx, 13, p->fsx) ) return;
- if ( set_arg_float(gctx, 14, p->fsy) ) return;
- if ( set_arg_float(gctx, 15, p->ssx) ) return;
- if ( set_arg_float(gctx, 16, p->ssy) ) return;
- if ( set_arg_float(gctx, 17, xo) ) return;
- if ( set_arg_float(gctx, 18, yo) ) return;
-
- dims[0] = pan_width;
- dims[1] = pan_height;
+
+ if ( set_arg_int(gctx, 3, pan_width) ) return;
+ if ( set_arg_float(gctx, 4, p->cnx) ) return;
+ if ( set_arg_float(gctx, 5, p->cny) ) return;
+ if ( set_arg_float(gctx, 6, p->fsx) ) return;
+ if ( set_arg_float(gctx, 7, p->fsy) ) return;
+ if ( set_arg_float(gctx, 8, p->ssx) ) return;
+ if ( set_arg_float(gctx, 9, p->ssy) ) return;
+ if ( set_arg_float(gctx, 10, p->res) ) return;
+ if ( set_arg_float(gctx, 11, p->clen) ) return;
+
+ dims[0] = pan_width * sampling;
+ dims[1] = pan_height * sampling;
+
+ ldims[0] = sampling;
+ ldims[1] = sampling;
+
+ err = clSetKernelArg(gctx->kern, 18,
+ sampling*sampling*sizeof(cl_float), NULL);
+ if ( err != CL_SUCCESS ) {
+ ERROR("Couldn't set local memory: %s\n", clError(err));
+ return;
+ }
err = clEnqueueNDRangeKernel(gctx->cq, gctx->kern, 2, NULL,
dims, ldims, 0, NULL, NULL);
@@ -250,7 +263,7 @@ static void do_panels(struct gpu_context *gctx, struct image *image,
tfs = p->min_fs + fs;
tss = p->min_ss + ss;
- image->data[tfs + image->width*tss] += frac*val;
+ image->data[tfs + image->width*tss] += val;
}
}
@@ -276,78 +289,54 @@ void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
int n_neg = 0;
int n_nan = 0;
int fs, ss;
- const int nkstep = 10;
- const int nxs = 4;
- const int nys = 4;
- int kstep;
- double klow, khigh, kinc, w;
- double frac;
- int xs, ys;
- int n, ntot;
+ int i;
if ( gctx == NULL ) {
ERROR("GPU setup failed.\n");
return;
}
- cell_get_cartesian(ucell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
- cell.s[0] = ax; cell.s[1] = ay; cell.s[2] = az;
- cell.s[3] = bx; cell.s[4] = by; cell.s[5] = bz;
- cell.s[6] = cx; cell.s[7] = cy; cell.s[8] = cz;
-
/* Ensure all required LUTs are available */
check_sinc_lut(gctx, na);
check_sinc_lut(gctx, nb);
check_sinc_lut(gctx, nc);
- if ( set_arg_mem(gctx, 8, gctx->intensities) ) return;
- if ( set_arg_mem(gctx, 9, gctx->sinc_luts[na-1]) ) return;
- if ( set_arg_mem(gctx, 10, gctx->sinc_luts[nb-1]) ) return;
- if ( set_arg_mem(gctx, 11, gctx->sinc_luts[nc-1]) ) return;
- if ( set_arg_mem(gctx, 12, gctx->flags) ) return;
-
/* Unit cell */
- err = clSetKernelArg(gctx->kern, 7, sizeof(cl_float16), &cell);
+ cell_get_cartesian(ucell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+ cell.s[0] = ax; cell.s[1] = ay; cell.s[2] = az;
+ cell.s[3] = bx; cell.s[4] = by; cell.s[5] = bz;
+ cell.s[6] = cx; cell.s[7] = cy; cell.s[8] = cz;
+
+ err = clSetKernelArg(gctx->kern, 12, sizeof(cl_float16), &cell);
if ( err != CL_SUCCESS ) {
ERROR("Couldn't set unit cell: %s\n", clError(err));
return;
}
+ if ( set_arg_mem(gctx, 13, gctx->intensities) ) return;
+ if ( set_arg_mem(gctx, 14, gctx->flags) ) return;
+ if ( set_arg_mem(gctx, 15, gctx->sinc_luts[na-1]) ) return;
+ if ( set_arg_mem(gctx, 16, gctx->sinc_luts[nb-1]) ) return;
+ if ( set_arg_mem(gctx, 17, gctx->sinc_luts[nc-1]) ) return;
+
/* Allocate memory for the result */
image->data = calloc(image->width * image->height, sizeof(float));
- w = image->lambda * image->bw;
- klow = 1.0 / (image->lambda + (w/2.0)); /* Smallest k */
- khigh = 1.0 / (image->lambda - (w/2.0)); /* Largest k */
- kinc = (khigh-klow) / (nkstep+1);
-
- ntot = nkstep * nxs * nys;
- frac = 1.0 / ntot;
-
- n = 0;
- for ( xs=0; xs<nxs; xs++ ) {
- for ( ys=0; ys<nys; ys++ ) {
-
- double xo, yo;
-
- xo = (1.0/nxs) * xs;
- yo = (1.0/nys) * ys;
+ double tot = 0.0;
+ for ( i=0; i<image->nsamples; i++ ) {
- for ( kstep=0; kstep<nkstep; kstep++ ) {
+ printf("%.1f eV, weight = %.5f\n",
+ ph_lambda_to_eV(1.0/image->spectrum[i].k),
+ image->spectrum[i].weight);
- double k;
+ do_panels(gctx, image, image->spectrum[i].k,
+ image->spectrum[i].weight,
+ &n_inf, &n_neg, &n_nan);
- k = klow + kstep*kinc;
+ tot += image->spectrum[i].weight;
- do_panels(gctx, image, &n_inf, &n_neg, &n_nan, k, frac,
- xo, yo);
-
- n++;
- progress_bar(n, ntot, "Simulating");
-
- }
- }
}
+ printf("total weight = %f\n", tot);
if ( n_neg + n_inf + n_nan ) {
ERROR("WARNING: The GPU calculation produced %i negative"
@@ -360,16 +349,9 @@ void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
for ( fs=0; fs<image->width; fs++ ) {
for ( ss=0; ss<image->height; ss++ ) {
- double twotheta, k;
- int idx;
-
- /* Calculate k this time round */
- k = 1.0/image->lambda;
-
- get_q(image, fs, ss, &twotheta, k);
-
- idx = fs + image->width*ss;
- image->twotheta[idx] = twotheta;
+ double twotheta;
+ get_q(image, fs, ss, &twotheta, 1.0/image->lambda);
+ image->twotheta[fs + image->width*ss] = twotheta;
}
}
@@ -437,7 +419,7 @@ struct gpu_context *setup_gpu(int no_sfac,
}
} else {
for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
- intensities_ptr[i] = 1e5;
+ intensities_ptr[i] = 100.0; /* Does nothing */
}
strncat(cflags, "-DFLAT_INTENSITIES ", 511-strlen(cflags));
}
@@ -468,6 +450,8 @@ struct gpu_context *setup_gpu(int no_sfac,
strncat(cflags, "-DPG23 ", 511-strlen(cflags));
} else if ( strcmp(sym, "m-3") == 0 ) {
strncat(cflags, "-DPGM3 ", 511-strlen(cflags));
+ } else if ( strcmp(sym, "321_H") == 0 ) {
+ strncat(cflags, "-DPG321H ", 511-strlen(cflags));
} else {
ERROR("Sorry! Point group '%s' is not currently"
" supported on the GPU."
diff --git a/src/diffraction.c b/src/diffraction.c
index 2f764987..826aac77 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -3,11 +3,13 @@
*
* Calculate diffraction patterns by Fourier methods
*
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
*
* Authors:
- * 2009-2012 Thomas White <taw@physics.org>
+ * 2009-2014 Thomas White <taw@physics.org>
+ * 2013-2014 Chun Hong Yoon <chun.hong.yoon@desy.de>
+ * 2013 Alexandra Tolstikova
*
* This file is part of CrystFEL.
*
@@ -325,7 +327,7 @@ static double molecule_factor(const double *intensities, const double *phases,
ld = q.u * cx + q.v * cy + q.w * cz;
/* No flags -> flat intensity distribution */
- if ( flags == NULL ) return 1.0e5;
+ if ( flags == NULL ) return 100.0;
switch ( m ) {
@@ -358,18 +360,250 @@ static double molecule_factor(const double *intensities, const double *phases,
}
+
+static void diffraction_at_k(struct image *image, const double *intensities,
+ const double *phases, const unsigned char *flags,
+ UnitCell *cell, GradientMethod m,
+ const SymOpList *sym, double k,
+ double ax, double ay, double az,
+ double bx, double by, double bz,
+ double cx, double cy, double cz,
+ double *lut_a, double *lut_b, double *lut_c,
+ double weight)
+{
+ unsigned int fs, ss;
+ const int nxs = 4;
+ const int nys = 4;
+
+ weight /= nxs*nys;
+
+ for ( fs=0; fs<image->width; fs++ ) {
+ for ( ss=0; ss<image->height; ss++ ) {
+
+ int idx;
+ double f_lattice, I_lattice;
+ double I_molecule;
+ struct rvec q;
+ double twotheta;
+ int xs, ys;
+ float xo, yo;
+
+ for ( xs=0; xs<nxs; xs++ ) {
+ for ( ys=0; ys<nys; ys++ ) {
+
+ xo = (1.0/nxs) * xs;
+ yo = (1.0/nys) * ys;
+
+ q = get_q(image, fs+xo, ss+yo, &twotheta, k);
+
+ f_lattice = lattice_factor(q, ax, ay, az,
+ bx, by, bz,
+ cx, cy, cz,
+ lut_a, lut_b, lut_c);
+
+ I_molecule = molecule_factor(intensities,
+ phases, flags, q,
+ ax, ay, az,
+ bx, by, bz,
+ cx, cy, cz,
+ m, sym);
+
+ I_lattice = pow(f_lattice, 2.0);
+
+ idx = fs + image->width*ss;
+ image->data[idx] += I_lattice * I_molecule * weight;
+ image->twotheta[idx] = twotheta;
+
+ }
+ }
+ }
+ progress_bar(fs, image->width-1, "Calculating diffraction");
+ }
+}
+
+
+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)
+{
+ struct sample *spectrum;
+ int i;
+ double eV;
+
+ spectrum = malloc(spec_size * sizeof(struct sample));
+ if ( spectrum == NULL ) return NULL;
+
+ eV = eV_cen - (spec_size/2)*eV_step;
+ 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)
+{
+ 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(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;
+
+ return spectrum;
+}
+
+
+struct sample *generate_SASE(struct image *image)
+{
+ struct sample *spectrum;
+ int i;
+ 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(ph_lambda_to_eV(image->lambda),
+ jitter_sigma_eV);
+
+ /* Convert FWHM to standard deviation. Note that bandwidth is taken to
+ * be "delta E over E" (E = photon energy), not the bandwidth in terms
+ * of wavelength, 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);
+
+ /* Add SASE-type noise to Gaussian spectrum */
+ add_sase_noise(spectrum, spec_size);
+
+ /* 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);
+
+ 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)
{
- unsigned int fs, ss;
double ax, ay, az;
double bx, by, bz;
double cx, cy, cz;
double *lut_a;
double *lut_b;
double *lut_c;
+ int i;
cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
@@ -383,39 +617,18 @@ void get_diffraction(struct image *image, int na, int nb, int nc,
lut_b = get_sinc_lut(nb);
lut_c = get_sinc_lut(nc);
- for ( fs=0; fs<image->width; fs++ ) {
- for ( ss=0; ss<image->height; ss++ ) {
+ for ( i=0; i<image->nsamples; i++ ) {
- int idx;
- double k;
- double f_lattice, I_lattice;
- double I_molecule;
- struct rvec q;
- double twotheta;
+ printf("%.1f eV, weight = %.5f\n",
+ ph_lambda_to_eV(1.0/image->spectrum[i].k),
+ image->spectrum[i].weight);
- /* Calculate k this time round */
- k = 1.0/image->lambda;
+ diffraction_at_k(image, intensities, phases,
+ flags, cell, m, sym, image->spectrum[i].k,
+ ax, ay, az, bx, by, bz, cx, cy, cz,
+ lut_a, lut_b, lut_c, image->spectrum[i].weight);
- q = get_q(image, fs, ss, &twotheta, k);
- f_lattice = lattice_factor(q, ax, ay, az,
- bx, by, bz,
- cx, cy, cz,
- lut_a, lut_b, lut_c);
-
- I_molecule = molecule_factor(intensities,
- phases, flags, q,
- ax,ay,az,bx,by,bz,cx,cy,cz,
- m, sym);
-
- I_lattice = pow(f_lattice, 2.0);
-
- idx = fs + image->width*ss;
- image->data[idx] = I_lattice * I_molecule;
- image->twotheta[idx] = twotheta;
-
- }
- progress_bar(fs, image->width-1, "Calculating diffraction");
}
free(lut_a);
diff --git a/src/diffraction.h b/src/diffraction.h
index b7e34a0e..a903346e 100644
--- a/src/diffraction.h
+++ b/src/diffraction.h
@@ -3,11 +3,13 @@
*
* Calculate diffraction patterns by Fourier methods
*
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
*
* Authors:
- * 2009-2012 Thomas White <taw@physics.org>
+ * 2009-2014 Thomas White <taw@physics.org>
+ * 2013-2014 Chun Hong Yoon <chun.hong.yoon@desy.de>
+ * 2013 Alexandra Tolstikova
*
* This file is part of CrystFEL.
*
@@ -48,4 +50,8 @@ extern void get_diffraction(struct image *image, int na, int nb, int nc,
const unsigned char *flags, UnitCell *cell,
GradientMethod m, const SymOpList *sym);
+extern struct sample *generate_tophat(struct image *image);
+
+extern struct sample *generate_SASE(struct image *image);
+
#endif /* DIFFRACTION_H */
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 606a173c..81a94c37 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -110,6 +110,15 @@ static void show_help(const char *s)
" --max-size=<s> Use <s> as the maximum crystal size in nm.\n"
" --min-size is also required.\n"
" --no-noise Do not calculate Poisson noise.\n"
+" -s, --sample-spectrum=<N> Use N samples from spectrum. Default 3.\n"
+" -x, --spectrum=<type> Use <type> for the calculation of spectrum.\n"
+" Choose from:\n"
+" tophat : Tophat spectrum. Bandwidth is\n"
+" taken from beam parameters.\n"
+" SASE : SASE spectrum. Random SASE pulse \n"
+" is generated from a model.\n"
+" Bandwidth is taken from beam \n"
+" parameters.\n"
);
}
@@ -253,7 +262,9 @@ int main(int argc, char *argv[])
char *outfile = NULL;
char *geometry = NULL;
char *beamfile = NULL;
+ char *spectrum_str = NULL;
GradientMethod grad;
+ SpectrumType spectrum_type;
int ndone = 0; /* Number of simulations done (images or not) */
int number = 1; /* Number used for filename of image */
int n_images = 1; /* Generate one image by default */
@@ -266,6 +277,7 @@ int main(int argc, char *argv[])
double max_size = 0.0;
char *sym_str = NULL;
SymOpList *sym;
+ int nsamples = 3;
/* Long options */
const struct option longopts[] = {
@@ -283,6 +295,9 @@ int main(int argc, char *argv[])
{"output", 1, NULL, 'o'},
{"geometry", 1, NULL, 'g'},
{"beam", 1, NULL, 'b'},
+ {"sample-spectrum", 1, NULL, 's'},
+ {"type-spectrum", 1, NULL, 'x'},
+ {"spectrum", 1, NULL, 'x'},
{"really-random", 0, &config_random, 1},
{"gpu-dev", 1, NULL, 2},
{"min-size", 1, NULL, 3},
@@ -291,7 +306,7 @@ int main(int argc, char *argv[])
};
/* Short options */
- while ((c = getopt_long(argc, argv, "hrn:i:t:p:o:g:b:y:",
+ while ((c = getopt_long(argc, argv, "hrn:i:t:p:o:g:b:y:s:x:",
longopts, NULL)) != -1) {
switch (c) {
@@ -344,6 +359,18 @@ int main(int argc, char *argv[])
sym_str = strdup(optarg);
break;
+ case 's' :
+ nsamples = strtol(optarg, &rval, 10);
+ if ( *rval != '\0' ) {
+ ERROR("Invalid number of spectrum samples.\n");
+ return 1;
+ }
+ break;
+
+ case 'x' :
+ spectrum_str = strdup(optarg);
+ break;
+
case 2 :
gpu_dev = atoi(optarg);
break;
@@ -443,6 +470,20 @@ int main(int argc, char *argv[])
return 1;
}
+ if ( spectrum_str == NULL ) {
+ STATUS("You didn't specify a spectrum type, so"
+ " I'm using a 'tophat' spectrum.\n");
+ spectrum_type = SPECTRUM_TOPHAT;
+ } else if ( strcmp(spectrum_str, "tophat") == 0) {
+ spectrum_type = SPECTRUM_TOPHAT;
+ } else if ( strcmp(spectrum_str, "SASE") == 0) {
+ spectrum_type = SPECTRUM_SASE;
+ } else {
+ ERROR("Unrecognised spectrum type '%s'\n", spectrum_str);
+ return 1;
+ }
+ free(spectrum_str);
+
if ( intfile == NULL ) {
/* Gentle reminder */
@@ -504,12 +545,13 @@ int main(int argc, char *argv[])
ERROR("Photon energy must be specified, not taken from the"
" HDF5 file. Please alter %s accordingly.\n", beamfile)
return 1;
- } else {
- double wl = ph_en_to_lambda(eV_to_J(image.beam->photon_energy));
- image.lambda = wl;
}
+
+ double wl = ph_en_to_lambda(eV_to_J(image.beam->photon_energy));
+ image.lambda = wl;
image.bw = image.beam->bandwidth;
image.div = image.beam->divergence;
+ image.nsamples = nsamples;
free(beamfile);
/* Load unit cell */
@@ -571,6 +613,18 @@ int main(int argc, char *argv[])
cell = cell_rotate(input_cell, orientation);
+ switch ( spectrum_type ) {
+
+ case SPECTRUM_TOPHAT :
+ image.spectrum = generate_tophat(&image);
+ break;
+
+ case SPECTRUM_SASE :
+ image.spectrum = generate_SASE(&image);
+ break;
+
+ }
+
/* Ensure no residual information */
image.data = NULL;
image.twotheta = NULL;
@@ -638,6 +692,7 @@ int main(int argc, char *argv[])
/* Clean up */
free(image.data);
free(image.twotheta);
+
cell_free(cell);
skip:
diff --git a/src/pattern_sim.h b/src/pattern_sim.h
index 61bd87d0..4269ea90 100644
--- a/src/pattern_sim.h
+++ b/src/pattern_sim.h
@@ -37,7 +37,7 @@
/* Maxmimum index to hold values up to (can be increased if necessary)
* WARNING: Altering this value constitutes an ABI change, and means you must
* update data/diffraction.cl then recompile and reinstall everything. */
-#define INDMAX 120
+#define INDMAX 130
/* Array size */
#define IDIM (INDMAX*2 +1)
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-30 16:31:41 +0000