aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/detector.c2
-rw-r--r--src/diffraction.c27
-rw-r--r--src/ewald.c116
-rw-r--r--src/image.h3
-rw-r--r--src/intensities.c2
-rw-r--r--src/pattern_sim.c4
6 files changed, 125 insertions, 29 deletions
diff --git a/src/detector.c b/src/detector.c
index 2c313cdc..ccf37024 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[x + image->width*y],
+ water = water_intensity(image->qvecs[0][x + image->width*y],
ph_lambda_to_en(image->lambda),
BEAM_RADIUS, WATER_RADIUS);
intensity += water;
diff --git a/src/diffraction.c b/src/diffraction.c
index 27037839..b9430ba7 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -174,20 +174,25 @@ void get_diffraction(struct image *image, int na, int nb, int nc, int no_sfac)
double complex f_molecule;
struct rvec q;
double complex val;
+ int i;
- q = image->qvecs[x + image->width*y];
+ for ( i=0; i<image->nspheres; i++ ) {
- f_lattice = lattice_factor(q, ax,ay,az,bx,by,bz,cx,cy,cz,
- na, nb, nc);
- if ( no_sfac ) {
- f_molecule = 1.0;
- } else {
- f_molecule = molecule_factor(image->molecule, q,
- ax,ay,az,bx,by,bz,cx,cy,cz);
- }
+ 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;
+ val = f_molecule * f_lattice;
+ image->sfacs[x + image->width*y] = val;
+
+ }
}
progress_bar(x, image->width-1, "Calculating lattice factors");
diff --git a/src/ewald.c b/src/ewald.c
index 94fb4ed6..f24e9273 100644
--- a/src/ewald.c
+++ b/src/ewald.c
@@ -21,6 +21,11 @@
#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;
@@ -52,18 +57,9 @@ static struct rvec quat_rot(struct rvec q, struct quaternion z)
}
-void get_ewald(struct image *image)
+static void add_sphere(struct image *image, double k)
{
int x, y;
- double k; /* Wavenumber */
-
- k = 1/image->lambda;
-
- image->qvecs = malloc(image->width * image->height
- * sizeof(struct rvec));
-
- image->twotheta = malloc(image->width * image->height
- * sizeof(double));
for ( x=0; x<image->width; x++ ) {
for ( y=0; y<image->height; y++ ) {
@@ -73,8 +69,8 @@ void get_ewald(struct image *image)
double r;
double twothetax, twothetay, twotheta;
double qx, qy, qz;
- struct rvec q;
- int p;
+ 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++ ) {
@@ -86,23 +82,80 @@ void get_ewald(struct image *image)
/ 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;
- q.u = qx; q.v = qy; q.w = qz;
- image->qvecs[x + image->width*y] = quat_rot(q,
+ /* 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 = 0;
+ 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);
- image->twotheta[x + image->width*y] = twotheta;
+ }
+ }
if ( (x==0) && (y==(int)image->y_centre) ) {
double s;
@@ -123,4 +176,37 @@ void get_ewald(struct image *image)
}
}
+
+}
+
+
+void get_ewald(struct image *image)
+{
+ double kc; /* Wavenumber */
+ int i, kstep;
+
+ kc = 1/image->lambda; /* Centre */
+
+ image->qvecs = malloc(image->width * image->height
+ * sizeof(struct rvec *));
+
+ image->twotheta = malloc(image->width * image->height
+ * sizeof(double));
+
+ /* Create the spheres */
+ image->nspheres = SAMPLING*SAMPLING*BWSAMPLING;
+ for ( i=0; i<image->nspheres; i++ ) {
+ image->qvecs[i] = malloc(image->width * image->height
+ * sizeof(struct rvec));
+ }
+
+ for ( kstep=0; kstep<BWSAMPLING; kstep++ ) {
+
+ double k;
+
+ k = kc + (kstep-(BWSAMPLING/2))*kc*(BANDWIDTH/BWSAMPLING);
+
+ add_sphere(image, k);
+
+ }
}
diff --git a/src/image.h b/src/image.h
index 166b20df..ff87b1bc 100644
--- a/src/image.h
+++ b/src/image.h
@@ -76,7 +76,8 @@ struct image {
int *hdr; /* Actual counts */
int16_t *data; /* Integer counts after bloom */
double complex *sfacs;
- struct rvec *qvecs;
+ struct rvec **qvecs;
+ int nspheres;
double *twotheta;
struct molecule *molecule;
UnitCell *indexed_cell;
diff --git a/src/intensities.c b/src/intensities.c
index d5b3604c..edc15910 100644
--- a/src/intensities.c
+++ b/src/intensities.c
@@ -72,7 +72,7 @@ void output_intensities(struct image *image, UnitCell *cell)
int found = 0;
int j;
- q = image->qvecs[x + image->width*y];
+ q = image->qvecs[0][x + image->width*y];
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 bfc6be37..70384250 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -247,6 +247,7 @@ 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;
@@ -302,6 +303,9 @@ 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);
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-16 16:48:56 +0000