aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2013-07-25 17:49:03 +0200
committerThomas White <taw@physics.org>2013-07-25 17:49:03 +0200
commitfebc21aa72b87d7f743eeb98c71ece50b301701f (patch)
tree849739ed2545751853cceafb5eb7f79873051b50
parent0431728cc5efd3321b1c60f97a830dd525cf04c8 (diff)
partialator: Account for Lorentz factor
-rw-r--r--src/post-refinement.c61
-rw-r--r--src/post-refinement.h4
-rw-r--r--src/scaling-report.c6
-rw-r--r--tests/pr_gradient_check.c2
4 files changed, 58 insertions, 15 deletions
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 9e4649a2..cbb4e5b5 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -87,8 +87,9 @@ static double partiality_rgradient(double r, double profile_radius)
}
-/* Return the gradient of parameter 'k' given the current status of 'image'. */
-double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
+/* Return the gradient of partiality wrt parameter 'k' given the current status
+ * of 'image'. */
+double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
{
double ds, azi;
double glow, ghigh;
@@ -200,6 +201,46 @@ double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
}
+/* Return the gradient of Lorentz factor wrt parameter 'k' given the current
+ * status of 'image'. */
+double l_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
+{
+ double ds;
+ signed int hs, ks, ls;
+
+ switch ( k ) {
+
+ /* Cell parameters do not affect Lorentz factor */
+ case REF_ASX :
+ case REF_BSX :
+ case REF_CSX :
+ case REF_ASY :
+ case REF_BSY :
+ case REF_CSY :
+ case REF_ASZ :
+ case REF_BSZ :
+ case REF_CSZ :
+ return 0.0;
+
+ /* Nor does change of radius */
+ case REF_R :
+ return 0.0;
+
+ default:
+ break;
+
+ }
+
+ assert(k == REF_DIV);
+
+ get_symmetric_indices(refl, &hs, &ks, &ls);
+
+ ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls);
+
+ return 0.0; /* FIXME! */
+}
+
+
static void apply_cell_shift(UnitCell *cell, int k, double shift)
{
double asx, asy, asz;
@@ -383,7 +424,7 @@ static double pr_iterate(Crystal *cr, const RefList *full,
double w;
double I_partial;
int k;
- double p;
+ double p, l;
Reflection *match;
double gradients[NUM_PARAMS];
const double osf = crystal_get_osf(cr);
@@ -398,21 +439,23 @@ static double pr_iterate(Crystal *cr, const RefList *full,
" is it marked as refinable?\n", ha, ka, la);
continue;
}
- I_full = osf * get_intensity(match);
+ I_full = get_intensity(match);
/* Actual measurement of this reflection from this pattern? */
- I_partial = get_intensity(refl);
+ I_partial = osf * get_intensity(refl);
p = get_partiality(refl);
+ l = get_lorentz(refl);
/* Calculate the weight for this reflection */
w = pow(get_esd_intensity(refl), 2.0);
- w += p * I_full * pow(get_esd_intensity(match), 2.0);
+ w += l * p * I_full * pow(get_esd_intensity(match), 2.0);
w = pow(w, -1.0);
/* Calculate all gradients for this reflection */
for ( k=0; k<NUM_PARAMS; k++ ) {
double gr;
- gr = gradient(cr, k, refl, pmodel);
+ gr = p_gradient(cr, k, refl, pmodel) * l;
+ gr += l_gradient(cr, k, refl, pmodel) * p;
gradients[k] = gr;
}
@@ -429,7 +472,7 @@ static double pr_iterate(Crystal *cr, const RefList *full,
if ( g > k ) continue;
M_c = gradients[g] * gradients[k];
- M_c *= w * pow(osf * I_full, 2.0);
+ M_c *= w * pow(I_full, 2.0);
M_curr = gsl_matrix_get(M, k, g);
gsl_matrix_set(M, k, g, M_curr + M_c);
@@ -437,7 +480,7 @@ static double pr_iterate(Crystal *cr, const RefList *full,
}
- delta_I = I_partial - (p * I_full);
+ delta_I = I_partial - (l * p * I_full);
v_c = w * delta_I * I_full * gradients[k];
v_curr = gsl_vector_get(v, k);
gsl_vector_set(v, k, v_curr + v_c);
diff --git a/src/post-refinement.h b/src/post-refinement.h
index 7b822938..f565a880 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -63,8 +63,8 @@ enum {
extern void pr_refine(Crystal *cr, const RefList *full, PartialityModel pmodel);
/* Exported so it can be poked by tests/pr_gradient_check */
-extern double gradient(Crystal *cr, int k, Reflection *refl,
- PartialityModel pmodel);
+extern double p_gradient(Crystal *cr, int k, Reflection *refl,
+ PartialityModel pmodel);
#endif /* POST_REFINEMENT_H */
diff --git a/src/scaling-report.c b/src/scaling-report.c
index a02a1796..1d7904a6 100644
--- a/src/scaling-report.c
+++ b/src/scaling-report.c
@@ -251,13 +251,13 @@ static void partiality_graph(cairo_t *cr, Crystal **crystals, int n,
if ( f == NULL ) continue;
if ( get_redundancy(f) < 2 ) continue;
- Ipart = get_intensity(refl);
- Ifull = crystal_get_osf(cryst) * get_intensity(f);
+ Ipart = crystal_get_osf(cryst) * get_intensity(refl);
+ Ifull = get_intensity(f);
//if ( Ifull < 10 ) continue; /* FIXME: Ugh */
pobs = Ipart/Ifull;
- pcalc = get_partiality(refl);
+ pcalc = get_lorentz(refl) * get_partiality(refl);
//STATUS("%4i %4i %4i : %9.6f %9.6f %e %e %e\n",
// h, k, l, pobs, pcalc,
diff --git a/tests/pr_gradient_check.c b/tests/pr_gradient_check.c
index 68a625a7..4600db9b 100644
--- a/tests/pr_gradient_check.c
+++ b/tests/pr_gradient_check.c
@@ -309,7 +309,7 @@ static double test_gradients(Crystal *cr, double incr_val, int refine,
grad = (grad1 + grad2) / 2.0;
i++;
- cgrad = gradient(cr, refine, refl, pmodel);
+ cgrad = p_gradient(cr, refine, refl, pmodel);
get_partial(refl, &r1, &r2, &p, &cl, &ch);