diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/Makefile.am | 4 | ||||
| -rw-r--r-- | src/detector.c | 2 | ||||
| -rw-r--r-- | src/diffraction.c | 106 | ||||
| -rw-r--r-- | src/diffraction.h | 3 | ||||
| -rw-r--r-- | src/ewald.c | 196 | ||||
| -rw-r--r-- | src/ewald.h | 23 | ||||
| -rw-r--r-- | src/image.h | 2 | ||||
| -rw-r--r-- | src/indexamajig.c | 3 | ||||
| -rw-r--r-- | src/intensities.c | 3 | ||||
| -rw-r--r-- | src/pattern_sim.c | 6 |
10 files changed, 81 insertions, 267 deletions
diff --git a/src/Makefile.am b/src/Makefile.am index a80ba7a2..6a40c9a0 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -8,7 +8,7 @@ AM_CFLAGS = -Wall AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" 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 \ + hdf5-file.c detector.c sfac.c intensities.c \ reflections.c pattern_sim_LDADD = @LIBS@ @@ -17,7 +17,7 @@ process_hkl_SOURCES = process_hkl.c sfac.c statistics.c cell.c utils.c \ process_hkl_LDADD = @LIBS@ indexamajig_SOURCES = indexamajig.c hdf5-file.c utils.c dirax.c cell.c image.c \ - intensities.c ewald.c peaks.c index.c filters.c \ + intensities.c peaks.c index.c filters.c \ diffraction.c detector.c sfac.c indexamajig_LDADD = @LIBS@ diff --git a/src/detector.c b/src/detector.c index ccf37024..ae7576b2 100644 --- a/src/detector.c +++ b/src/detector.c @@ -182,7 +182,7 @@ void record_image(struct image *image, int do_water, int do_poisson, if ( do_water ) { /* Add intensity contribution from water */ - water = water_intensity(image->qvecs[0][x + image->width*y], + water = water_intensity(get_q(image, x, y, 1, NULL), ph_lambda_to_en(image->lambda), BEAM_RADIUS, WATER_RADIUS); intensity += water; diff --git a/src/diffraction.c b/src/diffraction.c index b3d6bef0..bc3e685c 100644 --- a/src/diffraction.c +++ b/src/diffraction.c @@ -19,15 +19,17 @@ #include "image.h" #include "utils.h" #include "cell.h" -#include "ewald.h" #include "diffraction.h" #include "sfac.h" +#define SAMPLING (1) + + static double lattice_factor(struct rvec q, double ax, double ay, double az, - double bx, double by, double bz, - double cx, double cy, double cz, - int na, int nb, int nc) + double bx, double by, double bz, + double cx, double cy, double cz, + int na, int nb, int nc) { struct rvec Udotq; double f1, f2, f3; @@ -133,19 +135,56 @@ double water_intensity(struct rvec q, double en, } +struct rvec get_q(struct image *image, unsigned int xs, unsigned int ys, + unsigned int sampling, float *ttp) +{ + struct rvec q; + float twothetax, twothetay, twotheta, r; + float rx = 0.0; + float ry = 0.0; + int p; + + const float k = 1.0/image->lambda; + const unsigned int x = xs / sampling; + const unsigned int y = ys / sampling; /* Integer part only */ + + for ( p=0; p<image->det.n_panels; p++ ) { + if ( (x >= image->det.panels[p].min_x) + && (x <= image->det.panels[p].max_x) + && (y >= image->det.panels[p].min_y) + && (y <= image->det.panels[p].max_y) ) { + rx = ((float)xs - (sampling*image->det.panels[p].cx)) + / (sampling * image->resolution); + ry = ((float)ys - (sampling*image->det.panels[p].cy)) + / (sampling * image->resolution); + break; + } + } + + /* Calculate q-vector for this sub-pixel */ + r = sqrt(pow(rx, 2.0) + pow(ry, 2.0)); + twothetax = atan2(rx, image->camera_len); + twothetay = atan2(ry, image->camera_len); + twotheta = atan2(r, image->camera_len); + + if ( ttp != NULL ) *ttp = twotheta; + + q.u = k * sin(twothetax); + q.v = k * sin(twothetay); + q.w = k - k * cos(twotheta); + + return q; +} + + void get_diffraction(struct image *image, int na, int nb, int nc, int no_sfac) { - int x, y; + unsigned int xs, ys; double ax, ay, az; double bx, by, bz; double cx, cy, cz; double a, b, c, d; - /* Generate the array of reciprocal space vectors in image->qvecs */ - if ( image->qvecs == NULL ) { - get_ewald(image); - } - if ( image->molecule == NULL ) return; cell_get_cartesian(image->molecule->cell, &ax, &ay, &az, @@ -157,6 +196,7 @@ void get_diffraction(struct image *image, int na, int nb, int nc, int no_sfac) STATUS("Particle size = %i x %i x %i (=%5.2f x %5.2f x %5.2f nm)\n", na, nb, nc, na*a/1.0e-9, nb*b/1.0e-9, nc*c/1.0e-9); + /* Allocate (and zero) the "diffraction array" */ image->sfacs = calloc(image->width * image->height, sizeof(double complex)); @@ -167,34 +207,36 @@ void get_diffraction(struct image *image, int na, int nb, int nc, int no_sfac) } } - for ( x=0; x<image->width; x++ ) { - for ( y=0; y<image->height; y++ ) { + /* Needed later for Lorentz calculation */ + image->twotheta = malloc(image->width * image->height * sizeof(double)); + + for ( xs=0; xs<image->width*SAMPLING; xs++ ) { + for ( ys=0; ys<image->height*SAMPLING; ys++ ) { double f_lattice; double complex f_molecule; struct rvec q; - double complex val; - int i; - - for ( i=0; i<image->nspheres; i++ ) { - - q = image->qvecs[i][x + image->width*y]; - - f_lattice = lattice_factor(q, ax,ay,az,bx,by,bz,cx,cy,cz, - na, nb, nc); - if ( no_sfac ) { - f_molecule = 10000.0; - } else { - f_molecule = molecule_factor(image->molecule, q, - ax,ay,az,bx,by,bz,cx,cy,cz); - } - - val = f_molecule * f_lattice; - image->sfacs[x + image->width*y] += val; - + float twotheta; + double sw = 1.0/(SAMPLING*SAMPLING); /* Sample weight */ + + const unsigned int x = xs / SAMPLING; + const unsigned int y = ys / SAMPLING; /* Integer part only */ + + q = get_q(image, xs, ys, SAMPLING, &twotheta); + image->twotheta[x + image->width*y] = twotheta; + + f_lattice = lattice_factor(q, ax,ay,az,bx,by,bz,cx,cy,cz, + na, nb, nc); + if ( no_sfac ) { + f_molecule = 10000.0; + } else { + f_molecule = molecule_factor(image->molecule, q, + ax,ay,az,bx,by,bz,cx,cy,cz); } + image->sfacs[x + image->width*y] += sw * f_molecule * f_lattice; + } - progress_bar(x, image->width-1, "Calculating lattice factors"); + progress_bar(xs, SAMPLING*image->width-1, "Calculating lattice factors"); } } diff --git a/src/diffraction.h b/src/diffraction.h index efc59d36..36d1c21f 100644 --- a/src/diffraction.h +++ b/src/diffraction.h @@ -23,5 +23,6 @@ extern void get_diffraction(struct image *image, int na, int nb, int nc, int nosfac); extern double water_intensity(struct rvec q, double en, double beam_r, double water_r); - +extern struct rvec get_q(struct image *image, unsigned int xs, unsigned int ys, + unsigned int sampling, float *ttp); #endif /* DIFFRACTION_H */ diff --git a/src/ewald.c b/src/ewald.c deleted file mode 100644 index fca488f6..00000000 --- a/src/ewald.c +++ /dev/null @@ -1,196 +0,0 @@ -/* - * ewald.c - * - * Calculate q-vector arrays - * - * (c) 2006-2010 Thomas White <taw@physics.org> - * - * Part of CrystFEL - crystallography with a FEL - * - */ - - -#include <stdlib.h> -#include <math.h> -#include <stdio.h> - -#include "image.h" -#include "utils.h" -#include "cell.h" -#include "ewald.h" -#include "detector.h" - - -#define SAMPLING (4) -#define BWSAMPLING (10) -#define BANDWIDTH (0.05) - - -static struct rvec quat_rot(struct rvec q, struct quaternion z) -{ - struct rvec res; - double t01, t02, t03, t11, t12, t13, t22, t23, t33; - - t01 = z.w*z.x; - t02 = z.w*z.y; - t03 = z.w*z.z; - t11 = z.x*z.x; - t12 = z.x*z.y; - t13 = z.x*z.z; - t22 = z.y*z.y; - t23 = z.y*z.z; - t33 = z.z*z.z; - - res.u = (1.0 - 2.0 * (t22 + t33)) * q.u - + (2.0 * (t12 + t03)) * q.v - + (2.0 * (t13 - t02)) * q.w; - - res.v = (2.0 * (t12 - t03)) * q.u - + (1.0 - 2.0 * (t11 + t33)) * q.v - + (2.0 * (t01 + t23)) * q.w; - - res.w = (2.0 * (t02 + t13)) * q.u - + (2.0 * (t23 - t01)) * q.v - + (1.0 - 2.0 * (t11 + t22)) * q.w; - - return res; -} - - -static void add_sphere(struct image *image, double k, int soffs) -{ - int x, y; - - for ( x=0; x<image->width; x++ ) { - for ( y=0; y<image->height; y++ ) { - - double rx = 0.0; - double ry = 0.0; - double r; - double twothetax, twothetay, twotheta; - double qx, qy, qz; - struct rvec q1, q2, q3, q4; - int p, sx, sy, i; - - /* Calculate q vectors for Ewald sphere */ - for ( p=0; p<image->det.n_panels; p++ ) { - if ( (x >= image->det.panels[p].min_x) - && (x <= image->det.panels[p].max_x) - && (y >= image->det.panels[p].min_y) - && (y <= image->det.panels[p].max_y) ) { - rx = ((double)x - image->det.panels[p].cx) - / image->resolution; - ry = ((double)y - image->det.panels[p].cy) - / image->resolution; - break; - } - } - - /* Bottom left corner */ - r = sqrt(pow(rx, 2.0) + pow(ry, 2.0)); - twothetax = atan2(rx, image->camera_len); - twothetay = atan2(ry, image->camera_len); - twotheta = atan2(r, image->camera_len); - qx = k * sin(twothetax); - qy = k * sin(twothetay); - qz = k - k * cos(twotheta); - q1.u = qx; q1.v = qy; q1.w = qz; - /* 2theta value is calculated at the bottom left only */ - image->twotheta[x + image->width*y] = twotheta; - - /* Bottom right corner (using the same panel configuration!) */ - rx = ((double)(x+1) - image->det.panels[p].cx) - / image->resolution; - ry = ((double)y - image->det.panels[p].cy) - / image->resolution; - twothetax = atan2(rx, image->camera_len); - twothetay = atan2(ry, image->camera_len); - twotheta = atan2(r, image->camera_len); - qx = k * sin(twothetax); - qy = k * sin(twothetay); - qz = k - k * cos(twotheta); - q2.u = qx; q2.v = qy; q2.w = qz; - - /* Top left corner (using the same panel configuration!) */ - rx = ((double)x - image->det.panels[p].cx) - / image->resolution; - ry = ((double)(y+1) - image->det.panels[p].cy) - / image->resolution; - twothetax = atan2(rx, image->camera_len); - twothetay = atan2(ry, image->camera_len); - twotheta = atan2(r, image->camera_len); - qx = k * sin(twothetax); - qy = k * sin(twothetay); - qz = k - k * cos(twotheta); - q3.u = qx; q3.v = qy; q3.w = qz; - - /* Top right corner (using the same panel configuration!) */ - rx = ((double)(x+1) - image->det.panels[p].cx) - / image->resolution; - ry = ((double)(y+1) - image->det.panels[p].cy) - / image->resolution; - twothetax = atan2(rx, image->camera_len); - twothetay = atan2(ry, image->camera_len); - twotheta = atan2(r, image->camera_len); - qx = k * sin(twothetax); - qy = k * sin(twothetay); - qz = k - k * cos(twotheta); - q4.u = qx; q4.v = qy; q4.w = qz; - - /* Now interpolate between the values to get - * the sampling points */ - i = soffs; - for ( sx=0; sx<SAMPLING; sx++ ) { - for ( sy=0; sy<SAMPLING; sy++ ) { - - struct rvec q; - - q.u = q1.u + ((q2.u - q1.u)/SAMPLING)*sx - + ((q3.u - q1.u)/SAMPLING)*sy; - q.v = q1.v + ((q2.v - q1.v)/SAMPLING)*sx - + ((q3.v - q1.v)/SAMPLING)*sy; - q.w = q1.w + ((q2.w - q1.w)/SAMPLING)*sx - + ((q3.w - q1.w)/SAMPLING)*sy; - image->qvecs[i++][x + image->width*y] = quat_rot(q, - image->orientation); - - } - } - - } - } - -} - - -void get_ewald(struct image *image) -{ - double kc; /* Wavenumber */ - int i, kstep; - long long int mtotal = 0; - - kc = 1/image->lambda; /* Centre */ - - image->twotheta = malloc(image->width * image->height - * sizeof(double)); - - /* Create the spheres */ - image->nspheres = SAMPLING*SAMPLING*BWSAMPLING; - image->qvecs = malloc(image->nspheres * sizeof(struct rvec *)); - - for ( i=0; i<image->nspheres; i++ ) { - mtotal += image->width * image->height * sizeof(struct rvec); - image->qvecs[i] = malloc(image->width * image->height - * sizeof(struct rvec)); - } - STATUS("%i spheres, %lli Mbytes\n", image->nspheres, mtotal/(1024*1024)); - - for ( kstep=0; kstep<BWSAMPLING; kstep++ ) { - - double k; - - k = kc + (kstep-(BWSAMPLING/2))*kc*(BANDWIDTH/BWSAMPLING); - add_sphere(image, k, kstep*SAMPLING); - - } -} diff --git a/src/ewald.h b/src/ewald.h deleted file mode 100644 index 75ebb44d..00000000 --- a/src/ewald.h +++ /dev/null @@ -1,23 +0,0 @@ -/* - * ewald.h - * - * Calculate q-vector arrays - * - * (c) 2006-2010 Thomas White <taw@physics.org> - * - * Part of CrystFEL - crystallography with a FEL - * - */ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#ifndef EWALD_H -#define EWALD_H - -#include "image.h" - -extern void get_ewald(struct image *image); - -#endif /* EWALD_H */ diff --git a/src/image.h b/src/image.h index ff87b1bc..144ebe42 100644 --- a/src/image.h +++ b/src/image.h @@ -76,8 +76,6 @@ struct image { int *hdr; /* Actual counts */ int16_t *data; /* Integer counts after bloom */ double complex *sfacs; - struct rvec **qvecs; - int nspheres; double *twotheta; struct molecule *molecule; UnitCell *indexed_cell; diff --git a/src/indexamajig.c b/src/indexamajig.c index d05db543..43127e66 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -26,7 +26,6 @@ #include "hdf5-file.h" #include "index.h" #include "intensities.h" -#include "ewald.h" #include "peaks.h" #include "diffraction.h" #include "detector.h" @@ -236,7 +235,6 @@ int main(int argc, char *argv[]) /* Simulate a diffraction pattern */ image.sfacs = NULL; - image.qvecs = NULL; image.twotheta = NULL; image.hdr = NULL; @@ -245,7 +243,6 @@ int main(int argc, char *argv[]) image.orientation.x = 0.0; image.orientation.y = 0.0; image.orientation.z = 0.0; - get_ewald(&image); } diff --git a/src/intensities.c b/src/intensities.c index edc15910..8b5f51c5 100644 --- a/src/intensities.c +++ b/src/intensities.c @@ -19,6 +19,7 @@ #include "intensities.h" #include "cell.h" #include "sfac.h" +#include "diffraction.h" #define MAX_HITS (1024) @@ -72,7 +73,7 @@ void output_intensities(struct image *image, UnitCell *cell) int found = 0; int j; - q = image->qvecs[0][x + image->width*y]; + q = get_q(image, x, y, 1, NULL); hd = q.u * ax + q.v * ay + q.w * az; kd = q.u * bx + q.v * by + q.w * bz; diff --git a/src/pattern_sim.c b/src/pattern_sim.c index 70384250..5fd71fc1 100644 --- a/src/pattern_sim.c +++ b/src/pattern_sim.c @@ -247,7 +247,6 @@ int main(int argc, char *argv[]) do { int na, nb, nc; - int i; na = 8*random()/RAND_MAX + 4; nb = 8*random()/RAND_MAX + 4; @@ -270,7 +269,6 @@ int main(int argc, char *argv[]) } /* Ensure no residual information */ - image.qvecs = NULL; image.sfacs = NULL; image.data = NULL; image.twotheta = NULL; @@ -303,10 +301,6 @@ int main(int argc, char *argv[]) /* Clean up */ free(image.data); - for ( i=0; i<image.nspheres; i++ ) { - free(image.qvecs[i]); - } - free(image.qvecs); free(image.hdr); free(image.sfacs); free(image.twotheta); |