aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorThomas White <taw@bitwiz.org.uk>2015-05-10 23:32:59 +0200
committerThomas White <taw@physics.org>2015-05-19 13:57:51 +0200
commit0f07feaec6e2c96769ded00c7c31ed61bf7e019d (patch)
tree07aaaf7bb0bc7c4ae7cda922c5c931d20cdf97b3 /src
parentd5c4f2c4f68ac0893b7f37d8835cd01c849360e8 (diff)
Refine B factor as well
Diffstat (limited to 'src')
-rw-r--r--src/post-refinement.c41
1 files changed, 31 insertions, 10 deletions
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 3008b7af..57d7c5d9 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -174,8 +174,8 @@ static double volume_fraction(double rlow, double rhigh, double pr,
}
-/* Return the gradient of pG wrt parameter 'k' given the current status
- * of the crystal. */
+/* Return the gradient of pGexp(-Bs^2) wrt parameter 'k' given the current
+ * status of the crystal. */
double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
{
double glow, ghigh;
@@ -184,11 +184,20 @@ double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
double R = crystal_get_profile_radius(cr);
double osf = crystal_get_osf(cr);
double gr;
+ signed int hi, ki, li;
+ double s, B;
+ get_indices(refl, &hi, &ki, &li);
+ s = resolution(crystal_get_cell(cr), hi, ki, li);
+ B = crystal_get_Bfac(cr);
get_partial(refl, &rlow, &rhigh, &p);
if ( k == GPARAM_OSF ) {
- return p;
+ return p*exp(-B*s*s);
+ }
+
+ if ( k == GPARAM_BFAC ) {
+ return -s*s*osf*p*exp(-B*s*s);
}
if ( k == GPARAM_R ) {
@@ -283,6 +292,12 @@ static void apply_shift(Crystal *cr, int k, double shift)
crystal_set_profile_radius(cr, t);
break;
+ case GPARAM_BFAC :
+ t = crystal_get_Bfac(cr);
+ t += shift;
+ crystal_set_Bfac(cr, t);
+ break;
+
case GPARAM_OSF :
if ( isnan(shift) ) {
ERROR("Refusing nan shift of OSF\n");
@@ -335,7 +350,7 @@ static double pr_iterate(Crystal *cr, const RefList *full,
const int verbose = 0;
int num_params = 0;
enum gparam rv[32];
- double G;
+ double G, B;
*n_filtered = 0;
@@ -356,7 +371,7 @@ static double pr_iterate(Crystal *cr, const RefList *full,
/* If we are scaling, refine scale factors (duh) */
if ( !no_scale ) {
rv[num_params++] = GPARAM_OSF;
- //rv[num_params++] = GPARAM_BFAC; FIXME: Next
+ rv[num_params++] = GPARAM_BFAC;
}
reflections = crystal_get_reflections(cr);
@@ -365,6 +380,7 @@ static double pr_iterate(Crystal *cr, const RefList *full,
v = gsl_vector_calloc(num_params);
G = crystal_get_osf(cr);
+ B = crystal_get_Bfac(cr);
/* Construct the equations, one per reflection in this image */
for ( refl = first_refl(reflections, &iter);
@@ -376,7 +392,7 @@ static double pr_iterate(Crystal *cr, const RefList *full,
double w;
double I_partial;
int k;
- double p, L;
+ double p, L, s;
Reflection *match;
double gradients[num_params];
@@ -393,6 +409,7 @@ static double pr_iterate(Crystal *cr, const RefList *full,
p = get_partiality(refl);
L = get_lorentz(refl);
+ s = resolution(crystal_get_cell(cr), ha, ka, la);
/* Actual measurement of this reflection from this pattern? */
I_partial = get_intensity(refl);
@@ -429,7 +446,7 @@ static double pr_iterate(Crystal *cr, const RefList *full,
}
- delta_I = I_partial - (G * p * I_full / L);
+ delta_I = I_partial - G*exp(-B*s*s)*p*I_full / L;
v_c = w * delta_I * gradients[k];
v_curr = gsl_vector_get(v, k);
gsl_vector_set(v, k, v_curr + v_c);
@@ -477,6 +494,10 @@ static double pr_iterate(Crystal *cr, const RefList *full,
static double guide_dev(Crystal *cr, const RefList *full)
{
double dev = 0.0;
+ double G, B;
+
+ G = crystal_get_osf(cr);
+ B = crystal_get_Bfac(cr);
/* For each reflection */
Reflection *refl;
@@ -486,7 +507,7 @@ static double guide_dev(Crystal *cr, const RefList *full)
refl != NULL;
refl = next_refl(refl, iter) ) {
- double G, p, L;
+ double p, L, s;
signed int h, k, l;
Reflection *full_version;
double I_full, I_partial;
@@ -503,16 +524,16 @@ static double guide_dev(Crystal *cr, const RefList *full)
* scale_intensities() might not yet have been called, so the
* full version may not have been calculated yet. */
- G = crystal_get_osf(cr);
p = get_partiality(refl);
L = get_lorentz(refl);
I_partial = get_intensity(refl);
I_full = get_intensity(full_version);
+ s = resolution(crystal_get_cell(cr), h, k, l);
//STATUS("%3i %3i %3i %5.2f %5.2f %5.2f %5.2f %5.2f\n",
// h, k, l, G, p, I_partial, I_full,
// I_partial - p*G*I_full);
- dev += pow(I_partial - G*p*I_full/L, 2.0);
+ dev += pow(I_partial - G*exp(-B*s*s)*p*I_full/L, 2.0);
}