aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/facetron.c1
-rw-r--r--src/geometry.c52
-rw-r--r--src/image.h3
-rw-r--r--src/post-refinement.c53
-rw-r--r--src/post-refinement.h5
5 files changed, 78 insertions, 36 deletions
diff --git a/src/facetron.c b/src/facetron.c
index f561b554..2bac7511 100644
--- a/src/facetron.c
+++ b/src/facetron.c
@@ -467,6 +467,7 @@ int main(int argc, char *argv[])
images[i].det = det;
images[i].beam = beam;
images[i].osf = 1.0;
+ images[i].profile_radius = 0.005e9;
/* Muppet proofing */
images[i].data = NULL;
diff --git a/src/geometry.c b/src/geometry.c
index 0ef1f1c3..664817e0 100644
--- a/src/geometry.c
+++ b/src/geometry.c
@@ -103,26 +103,17 @@ static double excitation_error(double xl, double yl, double zl,
static double partiality(double r1, double r2, double r)
{
double q1, q2;
- double p, p1, p2;
+ double p1, p2;
/* Calculate degrees of penetration */
q1 = (r1 + r)/(2.0*r);
q2 = (r2 + r)/(2.0*r);
- /* Clamp */
- if ( q1 > 1.0 ) q1 = 1.0;
- if ( q1 < 0.0 ) q1 = 0.0;
- if ( q2 > 1.0 ) q2 = 1.0;
- if ( q2 < 0.0 ) q2 = 0.0;
-
/* Convert to partiality */
p1 = 3.0*pow(q1,2.0) - 2.0*pow(q1,3.0);
p2 = 3.0*pow(q2,2.0) - 2.0*pow(q2,3.0);
- /* Input values may have been backwards */
- p = fabs(p2 - p1);
-
- return p;
+ return p2 - p1;
}
@@ -143,8 +134,6 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
double klow, kcen, khigh; /* Wavenumber */
/* Bounding sphere for the shape transform approximation */
const double profile_cutoff = 0.02e9; /* 0.02 nm^-1 */
- /* Actual radius of the profile */
- const double profile_radius = 0.005e9;
cpeaks = malloc(sizeof(struct cpeak)*MAX_CPEAKS);
if ( cpeaks == NULL ) {
@@ -178,6 +167,8 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
double xda, yda; /* Position on detector */
int close, inside;
double part; /* Partiality */
+ int clamp_low = 0;
+ int clamp_high = 0;
/* Ignore central beam */
if ( (h==0) && (k==0) && (l==0) ) continue;
@@ -216,21 +207,36 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
/* If the "lower" Ewald sphere is a long way away, use the
* position at which the Ewald sphere would just touch the
* reflection. */
- if ( rlow > profile_cutoff ) rlow = profile_cutoff;
- if ( rlow < -profile_cutoff ) rlow = -profile_cutoff;
+ if ( rlow < -profile_cutoff ) {
+ rlow = -profile_cutoff;
+ clamp_low = -1;
+ }
+ if ( rlow > +profile_cutoff ) {
+ rlow = +profile_cutoff;
+ clamp_low = +1;
+ }
+ /* Likewise the "higher" Ewald sphere */
+ if ( rhigh < -profile_cutoff ) {
+ rhigh = -profile_cutoff;
+ clamp_high = -1;
+ }
+ if ( rhigh > +profile_cutoff ) {
+ rhigh = +profile_cutoff;
+ clamp_high = +1;
+ }
+ /* The six possible combinations of clamp_{low,high} (including
+ * zero) correspond to the six situations in Table 3 of Rossmann
+ * et al. (1979). */
- /* As above, other side. */
- if ( rhigh > profile_cutoff ) rhigh = profile_cutoff;
- if ( rhigh < -profile_cutoff ) rhigh = -profile_cutoff;
+ /* Calculate partiality and reject if too small */
+ part = partiality(rlow, rhigh, image->profile_radius);
+ if ( part < 0.1 ) continue;
/* Locate peak on detector. */
p = locate_peak(xl, yl, zl, kcen, image->det, &xda, &yda);
if ( p == -1 ) continue;
- part = partiality(rlow, rhigh, profile_radius);
-
- if ( part < 0.1 ) continue;
-
+ /* Add peak to list */
cpeaks[np].h = h;
cpeaks[np].k = k;
cpeaks[np].l = l;
@@ -239,6 +245,8 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
cpeaks[np].r1 = rlow;
cpeaks[np].r2 = rhigh;
cpeaks[np].p = part;
+ cpeaks[np].clamp1 = clamp_low;
+ cpeaks[np].clamp2 = clamp_high;
np++;
if ( output ) {
diff --git a/src/image.h b/src/image.h
index ae2f3179..0d2f6c2e 100644
--- a/src/image.h
+++ b/src/image.h
@@ -66,6 +66,8 @@ struct cpeak
double r1; /* First excitation error */
double r2; /* Second excitation error */
double p; /* Partiality */
+ int clamp1; /* Clamp status for r1 */
+ int clamp2; /* Clamp status for r2 */
/* Location in image */
int x;
@@ -102,6 +104,7 @@ struct image {
int f0_available; /* 0 if f0 wasn't available
* from the input. */
double osf; /* Overall scaling factor */
+ double profile_radius; /* Radius of reflection */
int width;
int height;
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 5ca0dc98..e2fa6d50 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -28,12 +28,33 @@
#include "geometry.h"
+/* Returns dp/dr at "r" */
+static double partiality_gradient(double r, double profile_radius)
+{
+ double q, dpdq, dqdr;
+
+ /* Calculate degree of penetration */
+ q = (r + profile_radius)/(2.0*profile_radius);
+
+ /* dp/dq */
+ dpdq = 6.0*(q-pow(q, 2.0));
+
+ /* dq/dr */
+ dqdr = 1.0 / (2.0*profile_radius);
+
+ return dpdq * dqdr;
+}
+
+
/* Return the gradient of parameter 'k' given the current status of 'image'. */
double gradient(struct image *image, int k,
- struct cpeak spot, double I_partial)
+ struct cpeak spot, double I_partial, double r)
{
double ds;
double nom, den;
+ double r1g = 0.0;
+ double r2g = 0.0;
+ double g = 0.0;
ds = 2.0 * resolution(image->indexed_cell, spot.h, spot.k, spot.l);
@@ -43,9 +64,20 @@ double gradient(struct image *image, int k,
return I_partial;
case REF_DIV :
- nom = sqrt(2.0) * ds * sin(image->div);
- den = sqrt(1.0 - cos(image->div));
- return nom/den;
+ /* Calculate the gradient of r1 and r2 wrt divergence */
+ if ( spot.clamp1 == 0 ) {
+ nom = sqrt(2.0) * ds * sin(image->div/2.0);
+ den = sqrt(1.0 - cos(image->div/2.0));
+ r1g = nom/den;
+ g += r1g * partiality_gradient(spot.r1, r);
+ }
+ if ( spot.clamp2 == 0 ) {
+ nom = sqrt(2.0) * ds * sin(image->div/2.0);
+ den = sqrt(1.0 - cos(image->div/2.0));
+ r2g = nom/den;
+ g += r2g * partiality_gradient(spot.r2, r);
+ }
+ return g;
}
@@ -171,7 +203,7 @@ double pr_iterate(struct image *image, double *i_full, const char *sym,
for ( k=0; k<NUM_PARAMS; k++ ) {
int g;
- double v_c;
+ double v_c, gr;
for ( g=0; g<NUM_PARAMS; g++ ) {
@@ -179,16 +211,19 @@ double pr_iterate(struct image *image, double *i_full, const char *sym,
M_curr = gsl_matrix_get(M, g, k);
- M_c = gradient(image, g, spots[h], I_partial)
- * gradient(image, k, spots[h], I_partial);
+ M_c = gradient(image, g, spots[h], I_partial,
+ image->profile_radius)
+ * gradient(image, k, spots[h], I_partial,
+ image->profile_radius);
M_c *= pow(I_full, 2.0);
gsl_matrix_set(M, g, k, M_curr + M_c);
}
- v_c = delta_I * I_full * gradient(image, k, spots[h],
- I_partial);
+ gr = gradient(image, k, spots[h], I_partial,
+ image->profile_radius);
+ v_c = delta_I * I_full * gr;
gsl_vector_set(v, k, v_c);
}
diff --git a/src/post-refinement.h b/src/post-refinement.h
index 1ef7a1fd..8dea914d 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -30,11 +30,6 @@ enum {
NUM_PARAMS
};
-
-/* Return the gradient of parameter 'k' given the current status of 'image'. */
-double gradient(struct image *image, int k, struct cpeak spot,
- double I_partial);
-
/* Apply the given shift to the 'k'th parameter of 'image'. */
void apply_shift(struct image *image, int k, double shift);
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-09 02:39:45 +0000