diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/Makefile.am | 5 | ||||
| -rw-r--r-- | src/diffraction.c | 2 | ||||
| -rw-r--r-- | src/image.h | 11 | ||||
| -rw-r--r-- | src/integr_sim.c | 232 | ||||
| -rw-r--r-- | src/list_tmp.h | 74 | ||||
| -rw-r--r-- | src/sfac.c | 6 | ||||
| -rw-r--r-- | src/statistics.c | 67 | ||||
| -rw-r--r-- | src/statistics.h | 22 | ||||
| -rw-r--r-- | src/utils.h | 76 |
9 files changed, 436 insertions, 59 deletions
diff --git a/src/Makefile.am b/src/Makefile.am index 3a251dfb..a10e59c0 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,4 +1,4 @@ -bin_PROGRAMS = pattern_sim #integr_sim +bin_PROGRAMS = pattern_sim integr_sim AM_CFLAGS = -Wall -g @CFLAGS@ @@ -6,5 +6,6 @@ pattern_sim_SOURCES = pattern_sim.c diffraction.c utils.c image.c cell.c \ hdf5-file.c ewald.c detector.c sfac.c pattern_sim_LDADD = @LIBS@ -integr_sim_SOURCES = integr_sim.c diffraction.c utils.c ewald.c sfac.c +integr_sim_SOURCES = integr_sim.c diffraction.c utils.c ewald.c sfac.c \ + statistics.c cell.c integr_sim_LDADD = @LIBS@ diff --git a/src/diffraction.c b/src/diffraction.c index e2d0061a..6351042d 100644 --- a/src/diffraction.c +++ b/src/diffraction.c @@ -93,7 +93,7 @@ static double complex molecule_factor(struct molecule *mol, struct threevec q, k = (signed int)rint(kd); l = (signed int)rint(ld); - r = get_integral(mol->reflections, h, k, l); + r = lookup_sfac(mol->reflections, h, k, l); return r; } diff --git a/src/image.h b/src/image.h index 5206048a..e9db7cad 100644 --- a/src/image.h +++ b/src/image.h @@ -20,6 +20,8 @@ #include <stdint.h> #include <complex.h> +#include "utils.h" + /* How is the scaling of the image described? */ typedef enum { @@ -57,15 +59,6 @@ struct threevec }; -struct quaternion -{ - double w; - double x; - double y; - double z; -}; - - /* Structure describing an image */ struct image { diff --git a/src/integr_sim.c b/src/integr_sim.c new file mode 100644 index 00000000..08788277 --- /dev/null +++ b/src/integr_sim.c @@ -0,0 +1,232 @@ +/* + * integr_sim.c + * + * Test integration of intensities + * + * (c) 2006-2009 Thomas White <thomas.white@desy.de> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <stdarg.h> +#include <stdlib.h> +#include <stdio.h> +#include <string.h> +#include <unistd.h> + +#include "cell.h" +#include "image.h" +#include "utils.h" +#include "statistics.h" +#include "sfac.h" + + +static void main_show_help(const char *s) +{ + printf("Syntax: %s\n\n", s); + printf("Test relrod integration\n\n"); + printf(" -h Display this help message\n"); +} + + +static void write_RvsQ(const char *name, double *ref, double *trueref, + unsigned int *counts, double scale, UnitCell *cell) +{ + FILE *fh; + double smax, sbracket; + signed int h, k, l; + + fh = fopen(name, "w"); + + smax = 0.0; + for ( h=-INDMAX; h<INDMAX; h++ ) { + for ( k=-INDMAX; k<INDMAX; k++ ) { + for ( l=-INDMAX; l<INDMAX; l++ ) { + double s = resolution(cell, h, k, l); + if ( (lookup_count(counts, h, k, l) > 0) && (s > smax) ) { + smax = s; + } + } + } + } + + for ( sbracket=0.0; sbracket<smax; sbracket+=smax/10.0 ) { + + double top = 0.0; + double bot = 0.0; + int n = 0; + + for ( h=-INDMAX; h<INDMAX; h++ ) { + for ( k=-INDMAX; k<INDMAX; k++ ) { + for ( l=-INDMAX; l<INDMAX; l++ ) { + + double s; + int c; + c = lookup_count(counts, h, k, l); + s = resolution(cell, h, k, l); + if ((s>=sbracket) && (s<sbracket+smax/10.0) && (c>0)) { + + double obs, calc, obsi; + + obs = lookup_intensity(ref, h, k, l); + calc = lookup_intensity(trueref, h, k, l); + + obsi = obs / (double)c; + top += fabs(obsi/scale - calc); + bot += obsi/scale; + n++; + } + + } + } + } + + fprintf(fh, "%8.5f %8.5f %i\n", sbracket+smax/20.0, top/bot, n); + + } + fclose(fh); +} + + +int main(int argc, char *argv[]) +{ + int c; + int i; + struct image image; + double *ref, *trueref; + unsigned int *counts; + size_t ref_size; + signed int h, k, l; + double R, scale; + FILE *fh; + + while ((c = getopt(argc, argv, "h")) != -1) { + + switch ( c ) { + + case 'h' : { + main_show_help(argv[0]); + return 0; + } + + } + + } + + /* Define image parameters */ + image.width = 1024; + image.height = 1024; + image.fmode = FORMULATION_CLEN; + image.x_centre = 512.5; + image.y_centre = 512.5; + image.camera_len = 0.05; /* 5 cm (front CCD can move from 5cm-20cm) */ + image.resolution = 13333.3; /* 75 micron pixel size */ + image.xray_energy = eV_to_J(2.0e3); /* 2 keV energy */ + image.lambda = ph_en_to_lambda(image.xray_energy); /* Wavelength */ + image.molecule = load_molecule(); + + /* Prepare array for integrated intensities */ + ref = new_list_intensity(); + + /* Array for sample counts */ + counts = new_list_count(); + + /* Calculate true intensities */ + get_reflections_cached(image.molecule, image.xray_energy); + /* Complex structure factors now in image.molecule->reflections */ + + for ( i=1; i<=10e3; i++ ) { + + #if 0 + image.orientation = random_quaternion(); + + /* Calculate reflections using large smax + * (rather than the actual value) */ + //get_reflections(&image, cell, 1.0e9); + + nrefl = image_feature_count(image.rflist); + for ( j=0; j<nrefl; j++ ) { + + struct imagefeature *f; + double t, s, intensity, F; + + f = image_get_feature(image.rflist, j); + + t = 100.0e-9; /* Thickness 100 nm */ + s = f->s; /* Get excitation error */ + F = structure_factor(f->h, f->k, f->l); + + /* Calculate intensity from this reflection */ + intensity = pow( F * SINC(M_PI*t*s), 2.0); + + if ( intensity < 0.1 ) continue; + + if ( (f->h == 2) && (f->k == 2) && (f->l == 2) ) { + fprintf(fh1, "%f %f\n", s, intensity); + } + + if ( (f->h == 15) && (f->k == 15) && (f->l == 15) ) { + fprintf(fh2, "%f %f\n", s, intensity); + } + + integrate_reflection(ref, f->h, f->k, f->l, intensity); + add_count(counts, f->h, f->k, f->l, 1); + } + #endif + + if ( i % 1000 == 0 ) { + + int j; + double mean_counts; + int ctot = 0; + int nmeas = 0; + double ff; + char name[64]; + + for ( j=0; j<ref_size; j++ ) { + ctot += counts[j]; + if ( counts[j] > 0 ) nmeas++; + } + mean_counts = (double)ctot/nmeas; + + ff = lookup_intensity(ref, 2, 2, 2) + / lookup_count(counts, 2, 2, 2); + + R = stat_r2(ref, trueref, counts, ref_size, &scale); + printf("%8i: R=%5.2f%% (sf=%7.4f)" + " mean meas/refl=%5.2f," + " %i reflections measured. %f\n", + i, R*100.0, scale, mean_counts, nmeas, + ff); + + /* Record graph of R against q for this N */ + snprintf(name, 63, "R_vs_q-%i.dat", i); + write_RvsQ(name, ref, trueref, counts, + scale, image.molecule->cell); + + } + + } + + fh = fopen("reflections.dat", "w"); + for ( h=-INDMAX; h<INDMAX; h++ ) { + for ( k=-INDMAX; k<INDMAX; k++ ) { + for ( l=-INDMAX; l<INDMAX; l++ ) { + int N; + N = lookup_count(counts, h, k, l); + if ( N == 0 ) continue; + double F = lookup_intensity(ref, h, k, l) / N; + fprintf(fh, "%3i %3i %3i %f\n", h, k, l, F); + } + } + } + fclose(fh); + + return 0; +} diff --git a/src/list_tmp.h b/src/list_tmp.h new file mode 100644 index 00000000..ffb5f393 --- /dev/null +++ b/src/list_tmp.h @@ -0,0 +1,74 @@ +static inline void LABEL(integrate)(TYPE *ref, signed int h, + signed int k, signed int l, + TYPE i) +{ + int idx; + + if ( (abs(h) > INDMAX) || (abs(k) > INDMAX) || (abs(l) > INDMAX) ) { + printf("\nReflection %i %i %i is out of range!\n", h, k, l); + printf("You need to re-configure INDMAX, delete the reflection" + " cache file and re-run.\n"); + exit(1); + } + + if ( h < 0 ) h += IDIM; + if ( k < 0 ) k += IDIM; + if ( l < 0 ) l += IDIM; + + idx = h + (IDIM*k) + (IDIM*IDIM*l); + ref[idx] += i; +} + + +static inline void LABEL(set)(TYPE *ref, signed int h, + signed int k, signed int l, + TYPE i) +{ + int idx; + + if ( (abs(h) > INDMAX) || (abs(k) > INDMAX) || (abs(l) > INDMAX) ) { + printf("\nReflection %i %i %i is out of range!\n", h, k, l); + printf("You need to re-configure INDMAX, delete the reflection" + " cache file and re-run.\n"); + exit(1); + } + + if ( h < 0 ) h += IDIM; + if ( k < 0 ) k += IDIM; + if ( l < 0 ) l += IDIM; + + idx = h + (IDIM*k) + (IDIM*IDIM*l); + ref[idx] = i; +} + + +static inline TYPE LABEL(lookup)(TYPE *ref, signed int h, + signed int k, signed int l) +{ + int idx; + + if ( (abs(h) > INDMAX) || (abs(k) > INDMAX) || (abs(l) > INDMAX) ) { + printf("\nReflection %i %i %i is out of range!\n", h, k, l); + printf("You need to re-configure INDMAX, delete the reflection" + " cache file and re-run.\n"); + exit(1); + } + + if ( h < 0 ) h += IDIM; + if ( k < 0 ) k += IDIM; + if ( l < 0 ) l += IDIM; + + idx = h + (IDIM*k) + (IDIM*IDIM*l); + return ref[idx]; +} + + +static inline TYPE *LABEL(new_list)(void) +{ + TYPE *r; + size_t r_size; + r_size = IDIM*IDIM*IDIM*sizeof(TYPE); + r = malloc(r_size); + memset(r, 0, r_size); + return r; +} @@ -400,7 +400,7 @@ double complex *get_reflections(struct molecule *mol, double en) &bsx, &bsy, &bsz, &csx, &csy, &csz); - reflections = reflist_new(); + reflections = new_list_sfac(); for ( h=-INDMAX; h<=INDMAX; h++ ) { for ( k=-INDMAX; k<=INDMAX; k++ ) { @@ -446,7 +446,7 @@ double complex *get_reflections(struct molecule *mol, double en) } - integrate_reflection(reflections, h, k, l, F); + set_sfac(reflections, h, k, l, F); //if ( (h!=0) || (k!=0) || (l!=0) ) { // tscat += cabs(F); //} else { @@ -477,7 +477,7 @@ void get_reflections_cached(struct molecule *mol, double en) goto calc; } - mol->reflections = reflist_new(); + mol->reflections = new_list_sfac(); r = fread(mol->reflections, sizeof(double complex), IDIM*IDIM*IDIM, fh); if ( r < IDIM*IDIM*IDIM ) { printf("Found cache file (%s), but failed to read.\n", s); diff --git a/src/statistics.c b/src/statistics.c new file mode 100644 index 00000000..adadf907 --- /dev/null +++ b/src/statistics.c @@ -0,0 +1,67 @@ +/* + * statistics.c + * + * Structure-factor statistics + * + * (c) 2007-2009 Thomas White <thomas.white@desy.de> + * + * integr_sim - Test relrod integration + * + */ + + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <math.h> +#include <stdlib.h> + + +/* By what (best-fitted) factor must the list "list" be multiplied by, + * if it were to be merged with "target"? */ +double stat_scale_intensity(double *obs, double *calc, unsigned int *c, + int size) +{ + double top = 0.0; + double bot = 0.0; + int i; + + for ( i=0; i<size; i++ ) { + + if ( c[i] > 0 ) { + double obsi; + obsi = obs[i] / (double)c[i]; + top += obsi * calc[i]; + bot += calc[i] * calc[i]; + } /* else reflection not measured so don't worry about it */ + + } + + return top/bot; +} + + +double stat_r2(double *obs, double *calc, unsigned int *c, int size, + double *scalep) +{ + double top = 0.0; + double bot = 0.0; + double scale; + int i; + scale = stat_scale_intensity(obs, calc, c, size); + *scalep = scale; + + for ( i=0; i<size; i++ ) { + + if ( c[i] > 0 ) { + double obsi; + obsi = obs[i] / (double)c[i]; + top += fabs(obsi/scale - calc[i]); + bot += obsi/scale; + } + + } /* else reflection not measured so don't worry about it */ + + return top/bot; +} diff --git a/src/statistics.h b/src/statistics.h new file mode 100644 index 00000000..bc8184ca --- /dev/null +++ b/src/statistics.h @@ -0,0 +1,22 @@ +/* + * statistics.h + * + * Structure-factor statistics + * + * (c) 2007-2009 Thomas White <thomas.white@desy.de> + * + * integr_sim - Test relrod integration + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#ifndef STATISTICS_H +#define STATISTICS_H + +double stat_r2(double *obs, double *calc, unsigned int *c, int size, + double *scalep); + +#endif /* STATISTICS_H */ diff --git a/src/utils.h b/src/utils.h index e81ea840..8e9fb6cc 100644 --- a/src/utils.h +++ b/src/utils.h @@ -22,6 +22,8 @@ #include <stdlib.h> +/* -------------------------- Fundamental constants ------------------------ */ + /* Electron charge in C */ #define ELECTRON_CHARGE (1.6021773e-19) @@ -34,9 +36,16 @@ /* Thomson scattering length (m) */ #define THOMSON_LENGTH (2.81794e-15) -/* Maxmimum index to go up to */ -#define INDMAX 40 -#define IDIM (INDMAX*2 +1) + +/* --------------------------- Useful datatypes ----------------------------- */ + +struct quaternion +{ + double w; + double x; + double y; + double z; +}; extern unsigned int biggest(signed int a, signed int b); @@ -59,6 +68,9 @@ extern void mapping_rotate(double x, double y, double z, double omega, double tilt); extern void progress_bar(int val, int total); + +/* ----------------------------- Useful macros ------------------------------ */ + #define rad2deg(a) ((a)*180/M_PI) #define deg2rad(a) ((a)*M_PI/180) @@ -77,51 +89,27 @@ extern void progress_bar(int val, int total); #define J_to_eV(a) ((a)/ELECTRON_CHARGE) -static inline void integrate_reflection(double complex *ref, signed int h, - signed int k, signed int l, - double complex i) -{ - int idx; - - if ( h < 0 ) h += IDIM; - if ( k < 0 ) k += IDIM; - if ( l < 0 ) l += IDIM; - - idx = h + (IDIM*k) + (IDIM*IDIM*l); - ref[idx] += i; -} - - -static inline double complex get_integral(double complex *ref, signed int h, - signed int k, signed int l) -{ - int idx; +/* -------------------- Indexed lists for specified types ------------------- */ - if ( (abs(h) > INDMAX) || (abs(k) > INDMAX) || (abs(l) > INDMAX) ) { - printf("\nReflection %i %i %i is out of range!\n", h, k, l); - printf("You need to re-configure INDMAX, delete the reflection" - " cache file and re-run.\n"); - exit(1); - } - - if ( h < 0 ) h += IDIM; - if ( k < 0 ) k += IDIM; - if ( l < 0 ) l += IDIM; +/* Maxmimum index to hold values up to (can be increased if necessary) */ +#define INDMAX 40 - idx = h + (IDIM*k) + (IDIM*IDIM*l); - return ref[idx]; -} +/* Array size */ +#define IDIM (INDMAX*2 +1) +/* Create functions for storing reflection intensities indexed as h,k,l */ +#define LABEL(x) x##_intensity +#define TYPE double +#include "list_tmp.h" -static inline double complex *reflist_new(void) -{ - double complex *r; - size_t r_size; - r_size = IDIM*IDIM*IDIM*sizeof(double complex); - r = malloc(r_size); - memset(r, 0, r_size); - return r; -} +/* As above, but for complex structure factors */ +#define LABEL(x) x##_sfac +#define TYPE double complex +#include "list_tmp.h" +/* As above, but for (unsigned) integer counts */ +#define LABEL(x) x##_count +#define TYPE unsigned int +#include "list_tmp.h" #endif /* UTILS_H */ |