diff options
Diffstat (limited to 'src')
-rw-r--r-- | src/diffraction-gpu.c | 2 | ||||
-rw-r--r-- | src/diffraction.c | 76 | ||||
-rw-r--r-- | src/diffraction.h | 2 | ||||
-rw-r--r-- | src/partial_sim.c | 71 | ||||
-rw-r--r-- | src/pattern_sim.c | 96 |
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 */ |