aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/Makefile.am7
-rw-r--r--src/integr_sim.c232
-rw-r--r--src/process_hkl.c268
-rw-r--r--src/utils.h2
4 files changed, 273 insertions, 236 deletions
diff --git a/src/Makefile.am b/src/Makefile.am
index 47f6f002..efaf09a8 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,4 +1,4 @@
-bin_PROGRAMS = pattern_sim integr_sim
+bin_PROGRAMS = pattern_sim process_hkl
AM_CFLAGS = -Wall
@@ -6,6 +6,5 @@ pattern_sim_SOURCES = pattern_sim.c diffraction.c utils.c image.c cell.c \
hdf5-file.c ewald.c detector.c sfac.c intensities.c
pattern_sim_LDADD = @LIBS@
-integr_sim_SOURCES = integr_sim.c diffraction.c utils.c ewald.c sfac.c \
- statistics.c cell.c
-integr_sim_LDADD = @LIBS@
+process_hkl_SOURCES = process_hkl.c sfac.c statistics.c cell.c utils.c
+process_hkl_LDADD = @LIBS@
diff --git a/src/integr_sim.c b/src/integr_sim.c
deleted file mode 100644
index 08788277..00000000
--- a/src/integr_sim.c
+++ /dev/null
@@ -1,232 +0,0 @@
-/*
- * 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/process_hkl.c b/src/process_hkl.c
new file mode 100644
index 00000000..097a3eb3
--- /dev/null
+++ b/src/process_hkl.c
@@ -0,0 +1,268 @@
+/*
+ * process_hkl.c
+ *
+ * Assemble and process FEL Bragg 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 <getopt.h>
+
+#include "utils.h"
+#include "statistics.h"
+#include "sfac.h"
+
+
+static void show_help(const char *s)
+{
+ printf("Syntax: %s [options]\n\n", s);
+ printf("Assemble and process FEL Bragg intensities.\n\n");
+ printf(" -h Display this help message.\n");
+ printf(" -i <filename> Specify input filename (\"-\" for stdin).\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);
+}
+
+
+static void write_reflections(const char *filename, unsigned int *counts,
+ double *ref)
+{
+ FILE *fh;
+ signed int h, k, l;
+
+ fh = fopen(filename, "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);
+}
+
+
+static double *ideal_intensities(double complex *sfac)
+{
+ double *ref;
+ signed int h, k, l;
+
+ ref = new_list_intensity();
+
+ /* Generate ideal reflections from complex structure factors */
+ for ( h=-INDMAX; h<=INDMAX; h++ ) {
+ for ( k=-INDMAX; k<=INDMAX; k++ ) {
+ for ( l=-INDMAX; l<=INDMAX; l++ ) {
+ double complex F = lookup_sfac(sfac, h, k, l);
+ double intensity = pow(cabs(F), 2.0);
+ set_intensity(ref, h, k, l, intensity);
+ }
+ }
+ }
+
+ return ref;
+}
+
+
+static void process_reflections(double *ref, double *trueref,
+ unsigned int *counts, unsigned int n_patterns,
+ UnitCell *cell)
+{
+ int j;
+ double mean_counts;
+ int ctot = 0;
+ int nmeas = 0;
+ double ff;
+ char name[64];
+ double R, scale;
+
+ for ( j=0; j<LIST_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, LIST_SIZE, &scale);
+ STATUS("%8u: R=%5.2f%% (sf=%7.4e)"
+ " mean meas/refl=%5.2f,"
+ " %i reflections measured. %f\n",
+ n_patterns, R*100.0, scale, mean_counts, nmeas, ff);
+
+ /* Record graph of R against q for this N */
+ snprintf(name, 63, "results/R_vs_q-%u.dat", n_patterns);
+ write_RvsQ(name, ref, trueref, counts, scale, cell);
+}
+
+
+int main(int argc, char *argv[])
+{
+ int c;
+ char *filename = NULL;
+ FILE *fh;
+ unsigned int n_patterns;
+ double *ref, *trueref;
+ unsigned int *counts;
+ char *rval;
+ struct molecule *mol;
+
+ /* Long options */
+ const struct option longopts[] = {
+ {"help", 0, NULL, 'h'},
+ {"input", 1, NULL, 'i'},
+ {0, 0, NULL, 0}
+ };
+
+ /* Short options */
+ while ((c = getopt_long(argc, argv, "hi:", longopts, NULL)) != -1) {
+
+ switch (c) {
+ case 'h' : {
+ show_help(argv[0]);
+ return 0;
+ }
+
+ case 'i' : {
+ filename = strdup(optarg);
+ break;
+ }
+
+ case 0 : {
+ break;
+ }
+
+ default : {
+ return 1;
+ }
+ }
+
+ }
+
+ if ( filename == NULL ) {
+ ERROR("Please specify filename using the -i option\n");
+ return 1;
+ }
+
+ mol = load_molecule();
+ get_reflections_cached(mol, eV_to_J(2.0e3));
+
+ ref = new_list_intensity();
+ counts = new_list_count();
+ trueref = ideal_intensities(mol->reflections);
+
+ if ( strcmp(filename, "-") == 0 ) {
+ fh = stdin;
+ } else {
+ fh = fopen(filename, "r");
+ }
+ free(filename);
+ if ( fh == NULL ) {
+ ERROR("Failed to open input file\n");
+ return 1;
+ }
+
+ n_patterns = 0;
+ do {
+
+ char line[1024];
+ int h, k, l, intensity;
+ int r;
+
+ rval = fgets(line, 1023, fh);
+ if ( strncmp(line, "New pattern", 11) == 0 ) {
+ n_patterns++;
+ if ( n_patterns % 1000 == 0 ) {
+ process_reflections(ref, trueref, counts,
+ n_patterns, mol->cell);
+ }
+ }
+
+ r = sscanf(line, "%i %i %i %i", &h, &k, &l, &intensity);
+ if ( r != 4 ) continue;
+
+ integrate_intensity(ref, h, k, l, intensity);
+ integrate_count(counts, h, k, l, 1);
+
+ } while ( rval != NULL );
+
+ fclose(fh);
+
+ write_reflections("results/reflections.hkl", counts, ref);
+
+ return 0;
+}
diff --git a/src/utils.h b/src/utils.h
index 1de206bf..2f31db57 100644
--- a/src/utils.h
+++ b/src/utils.h
@@ -131,6 +131,8 @@ static inline double distance3d(double x1, double y1, double z1,
/* Array size */
#define IDIM (INDMAX*2 +1)
+#define LIST_SIZE (IDIM*IDIM*IDIM)
+
/* Create functions for storing reflection intensities indexed as h,k,l */
#define LABEL(x) x##_intensity
#define TYPE double
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-16 13:28:24 +0000