aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@bitwiz.org.uk>2010-05-06 13:41:10 -0700
committerThomas White <taw@bitwiz.org.uk>2010-05-06 13:41:10 -0700
commit317d5b1c1cb6bd92383a9e1fa9790f324994279d (patch)
tree35ca39ae0544490875459f3bf8b1a65900e650e7
parent416ebc0310c4e3e5c48be409ee64aeb1dc2b3ced (diff)
Implement GRADIENT_MOSAIC and GRADIENT_INTERPOLATE
-rw-r--r--src/diffraction.c114
1 files changed, 103 insertions, 11 deletions
diff --git a/src/diffraction.c b/src/diffraction.c
index dcd6ae2c..1960f786 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -15,6 +15,7 @@
#include <stdio.h>
#include <string.h>
#include <complex.h>
+#include <assert.h>
#include "image.h"
#include "utils.h"
@@ -69,31 +70,121 @@ static double lattice_factor(struct rvec q, double ax, double ay, double az,
}
+static double interpolate_linear(const double *ref,
+ const unsigned int *counts,
+ float hd, signed int k, signed int l)
+{
+ signed int h;
+ double val1, val2;
+ float f;
+ unsigned int c1, c2;
+
+ h = (signed int)hd;
+ if ( hd < 0.0 ) h -= 1;
+ f = hd - (float)h;
+ assert(f >= 0.0);
+
+ val1 = lookup_intensity(ref, h, k, l);
+ val2 = lookup_intensity(ref, h+1, k, l);
+ c1 = lookup_count(counts, h, k, l);
+ c2 = lookup_count(counts, h+1, k, l);
+
+ if ( c1 == 0 ) {
+ ERROR("Needed intensity for %i %i %i, but don't have it.\n",
+ h, k, l);
+ return 1.0e20;
+ }
+
+ if ( c2 == 0 ) {
+ ERROR("Needed intensity for %i %i %i, but don't have it.\n",
+ h+1, k, l);
+ return 1.0e20;
+ }
+
+ val1 = val1 / (double)c1;
+ val2 = val2 / (double)c2;
+
+ return (1.0-f)*val1 + f*val2;
+}
+
+
+static double interpolate_bilinear(const double *ref,
+ const unsigned int *counts,
+ float hd, float kd, signed int l)
+{
+ signed int k;
+ double val1, val2;
+ float f;
+
+ k = (signed int)kd;
+ if ( kd < 0.0 ) k -= 1;
+ f = kd - (float)k;
+ assert(f >= 0.0);
+
+ val1 = interpolate_linear(ref, counts, hd, k, l);
+ val2 = interpolate_linear(ref, counts, hd, k+1, l);
+
+ return (1.0-f)*val1 + f*val2;
+}
+
+
+static double interpolate_intensity(const double *ref,
+ const unsigned int *counts,
+ float hd, float kd, float ld)
+{
+ signed int l;
+ double val1, val2;
+ float f;
+
+ l = (signed int)ld;
+ if ( ld < 0.0 ) l -= 1;
+ f = ld - (float)l;
+ assert(f >= 0.0);
+
+ val1 = interpolate_bilinear(ref, counts, hd, kd, l);
+ val2 = interpolate_bilinear(ref, counts, hd, kd, l+1);
+
+ return (1.0-f)*val1 + f*val2;
+}
+
+
/* Look up the structure factor for the nearest Bragg condition */
static double molecule_factor(const double *intensities,
const unsigned int *counts, struct rvec q,
double ax, double ay, double az,
double bx, double by, double bz,
- double cx, double cy, double cz)
+ double cx, double cy, double cz,
+ GradientMethod m)
{
- double hd, kd, ld;
+ float hd, kd, ld;
signed int h, k, l;
double r;
hd = q.u * ax + q.v * ay + q.w * az;
kd = q.u * bx + q.v * by + q.w * bz;
ld = q.u * cx + q.v * cy + q.w * cz;
- h = (signed int)rint(hd);
- k = (signed int)rint(kd);
- l = (signed int)rint(ld);
- if ( lookup_count(counts, h, k, l) == 0 ) {
- ERROR("Needed intensity for %i %i %i, but don't have it.\n",
- h, k, l);
- return 1.0e20;
+ switch ( m ) {
+ case GRADIENT_MOSAIC :
+ h = (signed int)rint(hd);
+ k = (signed int)rint(kd);
+ l = (signed int)rint(ld);
+ if ( lookup_count(counts, h, k, l) == 0 ) {
+ ERROR("Needed intensity for %i %i %i, but don't have it.\n",
+ h, k, l);
+ return 1.0e20;
+ }
+ r = lookup_intensity(intensities, h, k, l);
+ break;
+ case GRADIENT_INTERPOLATE :
+ r = interpolate_intensity(intensities, counts, hd, kd, ld);
+ break;
+ case GRADIENT_PHASED :
+ default:
+ ERROR("This gradient method not implemented yet.\n");
+ exit(1);
}
- r = lookup_intensity(intensities, h, k, l);
return r;
}
@@ -235,7 +326,8 @@ void get_diffraction(struct image *image, int na, int nb, int nc,
I_molecule = molecule_factor(intensities,
counts, q,
ax,ay,az,
- bx,by,bz,cx,cy,cz);
+ bx,by,bz,cx,cy,cz,
+ m);
}
I_lattice = pow(f_lattice, 2.0);