diff options
author | Thomas White <taw@physics.org> | 2009-11-25 16:19:05 +0100 |
---|---|---|
committer | Thomas White <taw@physics.org> | 2009-11-25 16:19:05 +0100 |
commit | 76cbb192f4abdd5f5c280cee964357c64364c783 (patch) | |
tree | ee643b3122cc168c9c6cdcbeb03ffeb8bfed69e7 /src/integr_sim.c | |
parent | 36addbc39e3d1a9da88959b6c07af9438402b016 (diff) |
Introduce integr_sim
Diffstat (limited to 'src/integr_sim.c')
-rw-r--r-- | src/integr_sim.c | 232 |
1 files changed, 232 insertions, 0 deletions
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; +} |