From 844947d19261a2c71280b9f638ad7f5e88e16c73 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 26 Nov 2009 12:20:41 +0100 Subject: Add more options, including random orientations --- src/pattern_sim.c | 165 +++++++++++++++++++++++++++++++++--------------------- src/utils.c | 50 +++++++++++++++++ src/utils.h | 15 +++++ 3 files changed, 167 insertions(+), 63 deletions(-) diff --git a/src/pattern_sim.c b/src/pattern_sim.c index ea7ed053..5e2df13b 100644 --- a/src/pattern_sim.c +++ b/src/pattern_sim.c @@ -35,35 +35,84 @@ static void show_help(const char *s) { printf("Syntax: %s\n\n", s); printf("Simulate diffraction patterns from small crystals\n"); - printf(" probed with femosecond pulses from a free electron laser\n\n"); - printf(" -h, --help Display this help message\n"); - printf(" --simulation-details Show details of the simulation\n"); + printf(" probed with femosecond pulses from a free electron laser.\n\n"); + printf(" -h, --help Display this help message\n"); + printf(" --simulation-details Show details of the simulation\n"); + printf(" --near-bragg Output h,k,l,I near Bragg conditions\n"); + printf(" -r, --random-orientation Use a randomly generated orientation\n"); + printf(" (a new orientation will be used for each image)\n"); + printf(" -n, --number= Generate N images. Default 1\n"); } -static void show_details(const char *s) +static void show_details() { - printf("%s: Simulation Details\n\n", s); - printf("Simulates diffraction patterns from small crystals\n"); - printf("probed with femtosecond pulses from a free electron laser\n\n"); + printf("This program simulates diffraction patterns from small crystals illuminated\n"); + printf("with femtosecond X-ray pulses from a free electron laser.\n\n"); + printf("Scattering Factors\n"); + printf("------------------\n"); + printf("Scattering factors\n"); +} + + +static struct quaternion read_quaternion() +{ + do { + + int r; + float w, x, y, z; + char line[1024]; + char *rval; + + printf("Please input quaternion: w x y z\n"); + rval = fgets(line, 1023, stdin); + if ( rval == NULL ) return invalid_quaternion(); + chomp(line); + + r = sscanf(line, "%f %f %f %f", &w, &x, &y, &z); + if ( r == 4 ) { + + struct quaternion quat; + + quat.w = w; + quat.x = x; + quat.y = y; + quat.z = z; + + return quat; + + } else { + fprintf(stderr, "Invalid rotation '%s'\n", line); + } + + } while ( 1 ); } int main(int argc, char *argv[]) { - int c, done; + int c; struct image image; char filename[1024]; - int number = 1; int config_simdetails = 0; + int config_nearbragg = 0; + int config_randomquat = 0; + int number = 1; /* Index for the current image */ + int n_images = 1; /* Generate one image by default */ + int done = 0; + /* Long options */ const struct option longopts[] = { - {"help", 0, NULL, 'h'}, - {"simulation-details", 0, &config_simdetails, 1}, - {0, 0, NULL, 0} + {"help", 0, NULL, 'h'}, + {"simulation-details", 0, &config_simdetails, 1}, + {"near-bragg", 0, &config_nearbragg, 1}, + {"random-orientation", 0, NULL, 'r'}, + {"number", 1, NULL, 'n'}, + {0, 0, NULL, 0} }; - while ((c = getopt_long(argc, argv, "h", longopts, NULL)) != -1) { + /* Short options */ + while ((c = getopt_long(argc, argv, "hrn:", longopts, NULL)) != -1) { switch (c) { case 'h' : { @@ -71,6 +120,16 @@ int main(int argc, char *argv[]) return 0; } + case 'r' : { + config_randomquat = 1; + break; + } + + case 'n' : { + n_images = atoi(optarg); + break; + } + case 0 : { break; } @@ -83,7 +142,7 @@ int main(int argc, char *argv[]) } if ( config_simdetails ) { - show_details(argv[0]); + show_details(); return 0; } @@ -102,64 +161,44 @@ int main(int argc, char *argv[]) /* Splurge a few useful numbers */ printf("Wavelength is %f nm\n", image.lambda/1.0e-9); -again: - - /* Read quaternion from stdin */ - done = 0; do { - int r; - float w, x, y, z; - char line[1024]; - char *rval; - - printf("Please input quaternion: w x y z\n"); - rval = fgets(line, 1023, stdin); - if ( rval == NULL ) return 0; - chomp(line); - - r = sscanf(line, "%f %f %f %f", &w, &x, &y, &z); - if ( r == 4 ) { - - printf("Rotation is: %f %f %f %f (modulus=%f)\n", - w, x, y, z, sqrtf(w*w + x*x + y*y + z*z)); - - image.orientation.w = w; - image.orientation.x = x; - image.orientation.y = y; - image.orientation.z = z; - - done = 1; - + /* Read quaternion from stdin */ + if ( config_randomquat ) { + image.orientation = random_quaternion(); } else { - fprintf(stderr, "Invalid rotation '%s'\n", line); + image.orientation = read_quaternion(); } - } while ( !done ); + if ( quaternion_valid(image.orientation) ) { + fprintf(stderr, "Orientation modulus is not zero!\n"); + return 1; + } + + /* Ensure no residual information */ + image.qvecs = NULL; + image.sfacs = NULL; + image.data = NULL; + image.twotheta = NULL; + image.hdr = NULL; - /* Ensure no residual information */ - image.qvecs = NULL; - image.sfacs = NULL; - image.data = NULL; - image.twotheta = NULL; - image.hdr = NULL; + get_diffraction(&image); + record_image(&image); - get_diffraction(&image); - record_image(&image); + snprintf(filename, 1023, "results/sim-%i.h5", number); + number++; - snprintf(filename, 1023, "results/sim-%i.h5", number); - number++; + /* Write the output file */ + hdf5_write(filename, image.data, image.width, image.height); - /* Write the output file */ - hdf5_write(filename, image.data, image.width, image.height); + /* Clean up */ + free(image.data); + free(image.qvecs); + free(image.hdr); + free(image.sfacs); + free(image.twotheta); - /* Clean up */ - free(image.data); - free(image.qvecs); - free(image.hdr); - free(image.sfacs); - free(image.twotheta); + if ( n_images && (number >= n_images) ) done = 1; - /* Do it all again */ - goto again; + } while ( !done ); } diff --git a/src/utils.c b/src/utils.c index 4f5f7da3..90720510 100644 --- a/src/utils.c +++ b/src/utils.c @@ -94,3 +94,53 @@ double poisson_noise(double expected) return k - 1; } + + +double quaternion_modulus(struct quaternion q) +{ + return sqrt(q.w*q.w + q.x*q.x + q.y*q.y + q.z*q.z); +} + + +struct quaternion normalise_quaternion(struct quaternion q) +{ + double mod; + struct quaternion r; + + mod = quaternion_modulus(q); + + r.w = q.w / mod; + r.x = q.x / mod; + r.y = q.y / mod; + r.z = q.z / mod; + + return r; +} + + +struct quaternion random_quaternion() +{ + struct quaternion q; + + q.w = (double)random()/RAND_MAX; + q.x = (double)random()/RAND_MAX; + q.y = (double)random()/RAND_MAX; + q.z = (double)random()/RAND_MAX; + normalise_quaternion(q); + + return q; +} + + +int quaternion_valid(struct quaternion q) +{ + double qmod; + + qmod = quaternion_modulus(q); + + /* Modulus = 1 to within some tolerance? + * Nasty allowance for floating-point accuracy follows... */ + if ( (qmod > 0.999) && (qmod < 1.001) ) return 1; + + return 0; +} diff --git a/src/utils.h b/src/utils.h index b4b5a58e..39988ae1 100644 --- a/src/utils.h +++ b/src/utils.h @@ -56,6 +56,11 @@ struct quaternion double z; }; +extern struct quaternion normalise_quaternion(struct quaternion q); +extern double quaternion_modulus(struct quaternion q); +extern struct quaternion random_quaternion(void); +extern int quaternion_valid(struct quaternion q); + /* --------------------------- Useful functions ----------------------------- */ @@ -67,6 +72,16 @@ extern void progress_bar(int val, int total, const char *text); extern double poisson_noise(double expected); /* Keep these ones inline, to avoid function call overhead */ +static inline struct quaternion invalid_quaternion(void) +{ + struct quaternion quat; + quat.w = 0.0; + quat.x = 0.0; + quat.y = 0.0; + quat.z = 0.0; + return quat; +} + static inline double distance(double x1, double y1, double x2, double y2) { return sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1)); -- cgit v1.2.3