aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/Makefile.am4
-rw-r--r--src/detector.c2
-rw-r--r--src/diffraction.c106
-rw-r--r--src/diffraction.h3
-rw-r--r--src/ewald.c196
-rw-r--r--src/ewald.h23
-rw-r--r--src/image.h2
-rw-r--r--src/indexamajig.c3
-rw-r--r--src/intensities.c3
-rw-r--r--src/pattern_sim.c6
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);
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-04 19:42:36 +0000