aboutsummaryrefslogtreecommitdiff
path: root/src/diffraction.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/diffraction.c')
-rw-r--r--src/diffraction.c106
1 files changed, 74 insertions, 32 deletions
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");
}
}