From febc21aa72b87d7f743eeb98c71ece50b301701f Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 25 Jul 2013 17:49:03 +0200 Subject: partialator: Account for Lorentz factor --- src/post-refinement.c | 61 +++++++++++++++++++++++++++++++++++++++++++-------- src/post-refinement.h | 4 ++-- src/scaling-report.c | 6 ++--- 3 files changed, 57 insertions(+), 14 deletions(-) (limited to 'src') 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 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, -- cgit v1.2.3