diff options
| -rw-r--r-- | src/post-refinement.c | 61 | ||||
| -rw-r--r-- | src/post-refinement.h | 4 | ||||
| -rw-r--r-- | src/scaling-report.c | 6 | ||||
| -rw-r--r-- | tests/pr_gradient_check.c | 2 |
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); |