aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/diffraction-gpu.c2
-rw-r--r--src/diffraction.c76
-rw-r--r--src/diffraction.h2
-rw-r--r--src/partial_sim.c71
-rw-r--r--src/pattern_sim.c96
5 files changed, 227 insertions, 20 deletions
diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c
index 1ed008e8..b5bd2e88 100644
--- a/src/diffraction-gpu.c
+++ b/src/diffraction-gpu.c
@@ -325,7 +325,7 @@ void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
double tot = 0.0;
for ( i=0; i<image->nsamples; i++ ) {
- printf("%.1f eV, weight = %.5f\n",
+ printf("%.3f eV, weight = %.5f\n",
ph_lambda_to_eV(1.0/image->spectrum[i].k),
image->spectrum[i].weight);
diff --git a/src/diffraction.c b/src/diffraction.c
index 92e18521..0fa9f79b 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -434,7 +434,8 @@ static int compare_samples(const void *a, const void *b)
static struct sample *get_gaussian_spectrum(double eV_cen, double eV_step,
- double sigma, int spec_size)
+ double sigma, int spec_size,
+ double eV_start)
{
struct sample *spectrum;
int i;
@@ -443,7 +444,12 @@ static struct sample *get_gaussian_spectrum(double eV_cen, double eV_step,
spectrum = malloc(spec_size * sizeof(struct sample));
if ( spectrum == NULL ) return NULL;
- eV = eV_cen - (spec_size/2)*eV_step;
+ 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);
@@ -568,7 +574,7 @@ struct sample *generate_SASE(struct image *image, gsl_rng *rng)
* points */
double eV_step = 6.0*sigma/(spec_size-1);
- spectrum = get_gaussian_spectrum(eV_cen, eV_step, sigma, spec_size);
+ 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);
@@ -592,6 +598,70 @@ struct sample *generate_SASE(struct image *image, gsl_rng *rng)
}
+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);
+
+ 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);
+
+ 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,
diff --git a/src/diffraction.h b/src/diffraction.h
index 5b67c198..e397bd3a 100644
--- a/src/diffraction.h
+++ b/src/diffraction.h
@@ -57,4 +57,6 @@ extern struct sample *generate_tophat(struct image *image);
extern struct sample *generate_SASE(struct image *image, gsl_rng *rng);
+extern struct sample *generate_twocolour(struct image *image);
+
#endif /* DIFFRACTION_H */
diff --git a/src/partial_sim.c b/src/partial_sim.c
index 97e2d14d..5b05218f 100644
--- a/src/partial_sim.c
+++ b/src/partial_sim.c
@@ -169,6 +169,46 @@ static void calculate_partials(Crystal *cr,
}
+static void draw_and_write_image(struct image *image, RefList *reflections,
+ gsl_rng *rng)
+{
+ Reflection *refl;
+ RefListIterator *iter;
+ const int background = 3000;
+ int i;
+
+ image->data = calloc(image->width*image->height, sizeof(float));
+
+ for ( refl = first_refl(reflections, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) )
+ {
+ double Ip;
+ double dfs, dss;
+ int fs, ss;
+
+ Ip = get_intensity(refl);
+
+ get_detector_pos(refl, &dfs, &dss);
+ fs = nearbyint(dfs);
+ ss = nearbyint(dss);
+ assert(fs >= 0);
+ assert(ss >= 0);
+ assert(fs < image->width);
+ assert(ss < image->height);
+ image->data[fs + image->width*ss] += Ip;
+
+ }
+
+ for ( i=0; i<image->width*image->height; i++ ) {
+ image->data[i] += poisson_noise(rng, background);
+ }
+
+ hdf5_write_image(image->filename, image);
+ free(image->data);
+}
+
+
static void show_help(const char *s)
{
printf("Syntax: %s [options]\n\n", s);
@@ -181,6 +221,7 @@ static void show_help(const char *s)
" -i, --input=<file> Read reflections from <file>.\n"
" Default: generate random ones instead (see -r).\n"
" -o, --output=<file> Write partials in stream format to <file>.\n"
+" --images=<prefix> Write images to <prefix>NNN.h5.\n"
" -g. --geometry=<file> Get detector geometry from file.\n"
" -b, --beam=<file> Get beam parameters from file\n"
" -p, --pdb=<file> PDB file from which to get the unit cell.\n"
@@ -221,6 +262,8 @@ struct queue_args
struct image *template_image;
double max_q;
+ char *image_prefix;
+
/* The overall histogram */
double p_hist[NBINS];
unsigned long int n_ref[NBINS];
@@ -241,6 +284,8 @@ struct worker_args
double p_hist[NBINS];
unsigned long int n_ref[NBINS];
double p_max[NBINS];
+
+ int n;
};
@@ -258,6 +303,7 @@ static void *create_job(void *vqargs)
wargs->image = *qargs->template_image;
qargs->n_started++;
+ wargs->n = qargs->n_started;
return wargs;
}
@@ -293,7 +339,13 @@ static void run_job(void *vwargs, int cookie)
orientation = random_quaternion(qargs->rngs[cookie]);
crystal_set_cell(cr, cell_rotate(qargs->cell, orientation));
- snprintf(wargs->image.filename, 255, "dummy.h5");
+ if ( qargs->image_prefix != NULL ) {
+ snprintf(wargs->image.filename, 255, "%s%i.h5",
+ qargs->image_prefix, wargs->n);
+ } else {
+ snprintf(wargs->image.filename, 255, "dummy.h5");
+ }
+
reflections = find_intersections(&wargs->image, cr);
crystal_set_reflections(cr, reflections);
@@ -310,6 +362,11 @@ static void run_job(void *vwargs, int cookie)
qargs->max_q, qargs->full_stddev,
qargs->noise_stddev, qargs->rngs[cookie]);
+ if ( qargs->image_prefix != NULL ) {
+ draw_and_write_image(&wargs->image, reflections,
+ qargs->rngs[cookie]);
+ }
+
/* Give a slightly incorrect cell in the stream */
mess_up_cell(cr, qargs->cnoise, qargs->rngs[cookie]);
@@ -372,6 +429,7 @@ int main(int argc, char *argv[])
double noise_stddev = 20.0;
gsl_rng *rng_for_seeds;
int config_random = 0;
+ char *image_prefix = NULL;
/* Long options */
const struct option longopts[] = {
@@ -389,6 +447,7 @@ int main(int argc, char *argv[])
{"osf-stddev", 1, NULL, 3},
{"full-stddev", 1, NULL, 4},
{"noise-stddev", 1, NULL, 5},
+ {"images", 1, NULL, 6},
{"really-random", 0, &config_random, 1},
@@ -492,6 +551,10 @@ int main(int argc, char *argv[])
}
break;
+ case 6 :
+ image_prefix = strdup(optarg);
+ break;
+
case 0 :
break;
@@ -595,8 +658,8 @@ int main(int argc, char *argv[])
free(output_file);
image.det = det;
- image.width = det->max_fs;
- image.height = det->max_ss;
+ image.width = det->max_fs + 1;
+ image.height = det->max_ss + 1;
image.lambda = ph_en_to_lambda(eV_to_J(beam->photon_energy));
image.div = beam->divergence;
@@ -609,6 +672,7 @@ int main(int argc, char *argv[])
image.indexed_by = INDEXING_SIMULATION;
image.num_peaks = 0;
image.num_saturated_peaks = 0;
+ image.spectrum_size = 0;
if ( random_intensities ) {
full = reflist_new();
@@ -629,6 +693,7 @@ int main(int argc, char *argv[])
qargs.full_stddev = full_stddev;
qargs.noise_stddev = noise_stddev;
qargs.max_q = largest_q(&image);
+ qargs.image_prefix = image_prefix;
qargs.rngs = malloc(n_threads * sizeof(gsl_rng *));
if ( qargs.rngs == NULL ) {
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 08fc75dd..a747cd72 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -54,6 +54,7 @@
#include "reflist.h"
#include "reflist-utils.h"
#include "pattern_sim.h"
+#include "stream.h"
static void show_help(const char *s)
@@ -90,6 +91,7 @@ static void show_help(const char *s)
" -s, --sample-spectrum=<N> Use N samples from spectrum. Default 3.\n"
" -x, --spectrum=<type> Type of spectrum to simulate.\n"
" --background=<N> Add N photons of Poisson background (default 0).\n"
+" --template=<file> Take orientations from stream <file>.\n"
);
}
@@ -251,6 +253,8 @@ int main(int argc, char *argv[])
int nsamples = 3;
gsl_rng *rng;
int background = 0;
+ char *template_file = NULL;
+ Stream *st = NULL;
/* Long options */
const struct option longopts[] = {
@@ -276,6 +280,7 @@ int main(int argc, char *argv[])
{"min-size", 1, NULL, 3},
{"max-size", 1, NULL, 4},
{"background", 1, NULL, 5},
+ {"template", 1, NULL, 6},
{0, 0, NULL, 0}
};
@@ -375,6 +380,10 @@ int main(int argc, char *argv[])
}
break;
+ case 6 :
+ template_file = strdup(optarg);
+ break;
+
case 0 :
break;
@@ -406,6 +415,19 @@ int main(int argc, char *argv[])
}
}
+ if ( template_file != NULL ) {
+ if ( config_randomquat ) {
+ ERROR("You cannot use -r and --template together.\n");
+ return 1;
+ }
+ st = open_stream_for_read(template_file);
+ if ( st == NULL ) {
+ ERROR("Failed to open stream.\n");
+ return 1;
+ }
+ free(template_file);
+ }
+
if ( sym_str == NULL ) sym_str = strdup("1");
sym = get_pointgroup(sym_str);
/* sym_str is used below */
@@ -451,6 +473,11 @@ int main(int argc, char *argv[])
spectrum_type = SPECTRUM_TOPHAT;
} else if ( strcasecmp(spectrum_str, "sase") == 0) {
spectrum_type = SPECTRUM_SASE;
+ } else if ( strcasecmp(spectrum_str, "twocolour") == 0 ||
+ strcasecmp(spectrum_str, "twocolor") == 0 ||
+ strcasecmp(spectrum_str, "twocolours") == 0 ||
+ strcasecmp(spectrum_str, "twocolors") == 0) {
+ spectrum_type = SPECTRUM_TWOCOLOUR;
} else {
ERROR("Unrecognised spectrum type '%s'\n", spectrum_str);
return 1;
@@ -578,23 +605,62 @@ int main(int argc, char *argv[])
}
- /* Read quaternion from stdin */
- if ( config_randomquat ) {
- orientation = random_quaternion(rng);
+ if ( st == NULL ) {
+
+ if ( config_randomquat ) {
+ orientation = random_quaternion(rng);
+ } else {
+ orientation = read_quaternion();
+ }
+
+ STATUS("Orientation is %5.3f %5.3f %5.3f %5.3f\n",
+ orientation.w, orientation.x,
+ orientation.y, orientation.z);
+
+ if ( !quaternion_valid(orientation) ) {
+ ERROR("Orientation modulus is not zero!\n");
+ return 1;
+ }
+
+ cell = cell_rotate(input_cell, orientation);
+
} else {
- orientation = read_quaternion();
- }
- STATUS("Orientation is %5.3f %5.3f %5.3f %5.3f\n",
- orientation.w, orientation.x,
- orientation.y, orientation.z);
+ struct image image;
+ int i;
+ Crystal *cr;
- if ( !quaternion_valid(orientation) ) {
- ERROR("Orientation modulus is not zero!\n");
- return 1;
- }
+ image.det = NULL;
+
+ /* Get data from next chunk */
+ if ( read_chunk(st, &image) ) break;
+
+ free(image.filename);
+ image_feature_list_free(image.features);
+
+ if ( image.n_crystals == 0 ) continue;
+
+ cr = image.crystals[0];
+ cell = crystal_get_cell(cr);
+
+ if ( image.n_crystals > 1 ) {
+ ERROR("Using the first crystal only.\n");
+ }
+
+ for ( i=1; i<image.n_crystals; i++ ) {
+
+ Crystal *cr = image.crystals[i];
+ cell = crystal_get_cell(cr);
- cell = cell_rotate(input_cell, orientation);
+ reflist_free(crystal_get_reflections(cr));
+ cell_free(crystal_get_cell(cr));
+ crystal_free(cr);
+
+ }
+
+ free(image.crystals);
+
+ }
switch ( spectrum_type ) {
@@ -606,6 +672,10 @@ int main(int argc, char *argv[])
image.spectrum = generate_SASE(&image, rng);
break;
+ case SPECTRUM_TWOCOLOUR :
+ image.spectrum = generate_twocolour(&image);
+ break;
+
}
/* Ensure no residual information */
generated by cgit v1.2.3 (git 2.25.1) at 2024-09-11 18:10:28 +0000