From 687eed4650e74a86f45cd533dd94a0a2ace1d1e8 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 8 Feb 2010 10:57:09 +0100 Subject: Rework multisampling (remove Ewald module) --- src/diffraction.c | 106 +++++++++++++++++++++++++++++++++++++----------------- 1 file changed, 74 insertions(+), 32 deletions(-) (limited to 'src/diffraction.c') 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; pdet.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; xwidth; x++ ) { - for ( y=0; yheight; y++ ) { + /* Needed later for Lorentz calculation */ + image->twotheta = malloc(image->width * image->height * sizeof(double)); + + for ( xs=0; xswidth*SAMPLING; xs++ ) { + for ( ys=0; ysheight*SAMPLING; ys++ ) { double f_lattice; double complex f_molecule; struct rvec q; - double complex val; - int i; - - for ( i=0; inspheres; 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"); } } -- cgit v1.2.3