aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2010-11-22 16:23:22 +0100
committerThomas White <taw@physics.org>2012-02-22 15:27:06 +0100
commitdb4b4370b759dff84ec0737dcefa1d61364135c0 (patch)
tree1326e68d2699772f2a49fd17b12ca193d2f84907 /src
parentdc16cc44c7ec1c775c55b813a7ea9ec17b16f656 (diff)
Add refinement of unit cell and orientation
Diffstat (limited to 'src')
-rw-r--r--src/post-refinement.c112
-rw-r--r--src/post-refinement.h9
2 files changed, 104 insertions, 17 deletions
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 5fa5a4c0..33f6780b 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -25,6 +25,7 @@
#include "peaks.h"
#include "symmetry.h"
#include "geometry.h"
+#include "cell.h"
/* Returns dp/dr at "r" */
@@ -67,34 +68,43 @@ static double partiality_rgradient(double r, double profile_radius)
double gradient(struct image *image, int k,
struct cpeak spot, double I_partial, double r)
{
- double ds;
+ double ds, tt, azi;
double nom, den;
- double r1g = 0.0;
- double r2g = 0.0;
double g = 0.0;
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+ double xl, yl, zl;
+ cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+ xl = spot.h*asx + spot.k*bsx + spot.l*csx;
+ yl = spot.h*asy + spot.k*bsy + spot.l*csy;
+ zl = spot.h*asz + spot.k*bsz + spot.l*csz;
ds = 2.0 * resolution(image->indexed_cell, spot.h, spot.k, spot.l);
+ tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+k);
+ azi = angle_between(1.0, 0.0, 0.0, xl, yl, 0.0);
+ /* Calculate the gradient of partiality wrt excitation error. */
+ if ( spot.clamp1 == 0 ) {
+ g += partiality_gradient(spot.r1, r);
+ }
+ if ( spot.clamp2 == 0 ) {
+ g += partiality_gradient(spot.r2, r);
+ }
+
+ /* For many gradients, just multiply the above number by the gradient
+ * of excitation error wrt whatever. */
switch ( k ) {
case REF_SCALE :
return I_partial;
case REF_DIV :
- /* 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;
+ nom = sqrt(2.0) * ds * sin(image->div/2.0);
+ den = sqrt(1.0 - cos(image->div/2.0));
+ return (nom/den) * g;
case REF_R :
if ( spot.clamp1 == 0 ) {
@@ -105,6 +115,26 @@ double gradient(struct image *image, int k,
}
return g;
+ /* Cell parameters and orientation */
+ case REF_ASX :
+ return spot.h * pow(sin(tt), -1.0) * cos(azi) * g;
+ case REF_BSX :
+ return spot.k * pow(sin(tt), -1.0) * cos(azi) * g;
+ case REF_CSX :
+ return spot.l * pow(sin(tt), -1.0) * cos(azi) * g;
+ case REF_ASY :
+ return spot.h * pow(sin(tt), -1.0) * sin(azi) * g;
+ case REF_BSY :
+ return spot.k * pow(sin(tt), -1.0) * sin(azi) * g;
+ case REF_CSY :
+ return spot.l * pow(sin(tt), -1.0) * sin(azi) * g;
+ case REF_ASZ :
+ return spot.h * pow(cos(tt), -1.0) * g;
+ case REF_BSZ :
+ return spot.k * pow(cos(tt), -1.0) * g;
+ case REF_CSZ :
+ return spot.l * pow(cos(tt), -1.0) * g;
+
}
ERROR("No gradient defined for parameter %i\n", k);
@@ -112,6 +142,42 @@ double gradient(struct image *image, int k,
}
+static void apply_cell_shift(UnitCell *cell, int k, double shift)
+{
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+ double a, b, c, al, be, ga;
+
+ cell_get_reciprocal(cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+ switch ( k )
+ {
+ case REF_ASX : asx += shift; break;
+ case REF_ASY : asy += shift; break;
+ case REF_ASZ : asz += shift; break;
+ case REF_BSX : bsx += shift; break;
+ case REF_BSY : bsy += shift; break;
+ case REF_BSZ : bsz += shift; break;
+ case REF_CSX : csx += shift; break;
+ case REF_CSY : csy += shift; break;
+ case REF_CSZ : csz += shift; break;
+ }
+
+ cell_set_reciprocal(cell, asx, asy, asz,
+ bsx, bsy, bsz,
+ csx, csy, csz);
+
+ if ( k == REF_CSZ ) {
+ cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga);
+ STATUS("New cell: %5.2f %5.2f %5.2f nm %5.2f %5.2f %5.2f deg\n",
+ a/1.0e-9, b/1.0e-9, c/1.0e-9,
+ rad2deg(al), rad2deg(be), rad2deg(ga));
+ }
+}
+
+
/* Apply the given shift to the 'k'th parameter of 'image'. */
void apply_shift(struct image *image, int k, double shift)
{
@@ -131,6 +197,18 @@ void apply_shift(struct image *image, int k, double shift)
image->profile_radius += shift;
break;
+ case REF_ASX :
+ case REF_ASY :
+ case REF_ASZ :
+ case REF_BSX :
+ case REF_BSY :
+ case REF_BSZ :
+ case REF_CSX :
+ case REF_CSY :
+ case REF_CSZ :
+ apply_cell_shift(image->indexed_cell, k, shift);
+ break;
+
default :
ERROR("No shift defined for parameter %i\n", k);
abort();
diff --git a/src/post-refinement.h b/src/post-refinement.h
index 660530e9..983ecc28 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -28,6 +28,15 @@ enum {
REF_SCALE,
REF_DIV,
REF_R,
+ REF_ASX,
+ REF_BSX,
+ REF_CSX,
+ REF_ASY,
+ REF_BSY,
+ REF_CSY,
+ REF_ASZ,
+ REF_BSZ,
+ REF_CSZ,
NUM_PARAMS
};