From 52fc1532ba4a93f63b82ed4f5382d637dccac3c5 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 19 Mar 2015 17:56:37 +0100 Subject: Factorise dr/da part of gradient calculation for PR and prediction --- src/post-refinement.c | 127 ++++++++++++-------------------------------------- 1 file changed, 31 insertions(+), 96 deletions(-) (limited to 'src/post-refinement.c') diff --git a/src/post-refinement.c b/src/post-refinement.c index 745e5bd7..a532cd07 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -54,6 +54,8 @@ /* Maximum number of iterations of NLSq to do for each image per macrocycle. */ #define MAX_CYCLES (10) +/* Number of parameters to refine, in the order they appear in the enum */ +#define NUM_PARAMS (9) /* Returns dp(gauss)/dr at "r" */ static double gaussian_fraction_gradient(double r, double R) @@ -179,27 +181,17 @@ static double volume_fraction(double rlow, double rhigh, double pr, * of 'image'. */ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) { - double azi; double glow, ghigh; - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - double xl, yl, zl; - double ds; - signed int hs, ks, ls; double rlow, rhigh, p; - double philow, phihigh, phi; - double khigh, klow; - double tl, cet, cez; struct image *image = crystal_get_image(cr); double R = crystal_get_profile_radius(cr); - double D, psph; get_partial(refl, &rlow, &rhigh, &p); - if ( k == REF_R ) { + if ( k == GPARAM_R ) { double Rglow, Rghigh; + double D, psph; D = rlow - rhigh; psph = volume_fraction(rlow, rhigh, R, pmodel); @@ -215,78 +207,21 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) glow = partiality_gradient(rlow, R, pmodel, rlow, rhigh); ghigh = partiality_gradient(rhigh, R, pmodel, rlow, rhigh); - get_symmetric_indices(refl, &hs, &ks, &ls); - ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls); - - cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); - xl = hs*asx + ks*bsx + ls*csx; - yl = hs*asy + ks*bsy + ls*csy; - zl = hs*asz + ks*bsz + ls*csz; - - /* "low" gives the largest Ewald sphere (wavelength short => k large) - * "high" gives the smallest Ewald sphere (wavelength long => k small) - */ - klow = 1.0/(image->lambda - image->lambda*image->bw/2.0); - khigh = 1.0/(image->lambda + image->lambda*image->bw/2.0); - - tl = sqrt(xl*xl + yl*yl); + if ( k == GPARAM_DIV ) { - cet = -sin(image->div/2.0) * klow; - cez = -cos(image->div/2.0) * klow; - philow = angle_between_2d(tl-cet, zl-cez, 0.0, 1.0); - - cet = -sin(image->div/2.0) * khigh; - cez = -cos(image->div/2.0) * khigh; - phihigh = angle_between_2d(tl-cet, zl-cez, 0.0, 1.0); - - /* Approximation: philow and phihigh are very similar */ - phi = (philow + phihigh) / 2.0; - - azi = atan2(yl, xl); - - /* For many gradients, just multiply the above number by the gradient - * of excitation error wrt whatever. */ - switch ( k ) { + double D, psph, ds; + signed int hs, ks, ls; - /* Cell parameters and orientation */ - case REF_ASX : - return - hs * sin(phi) * cos(azi) * (glow-ghigh); - - case REF_BSX : - return - ks * sin(phi) * cos(azi) * (glow-ghigh); - - case REF_CSX : - return - ls * sin(phi) * cos(azi) * (glow-ghigh); - - case REF_ASY : - return - hs * sin(phi) * sin(azi) * (glow-ghigh); - - case REF_BSY : - return - ks * sin(phi) * sin(azi) * (glow-ghigh); - - case REF_CSY : - return - ls * sin(phi) * sin(azi) * (glow-ghigh); - - case REF_ASZ : - return - hs * cos(phi) * (glow-ghigh); - - case REF_BSZ : - return - ks * cos(phi) * (glow-ghigh); - - case REF_CSZ : - return - ls * cos(phi) * (glow-ghigh); - - case REF_DIV : D = rlow - rhigh; psph = volume_fraction(rlow, rhigh, R, pmodel); + get_symmetric_indices(refl, &hs, &ks, &ls); + ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls); + return (ds/2.0)*(glow+ghigh) - 4.0*R*psph*ds/(3.0*D*D); } - ERROR("No gradient defined for parameter %i\n", k); - abort(); + return r_gradient(crystal_get_cell(cr), k, refl, image) * (glow-ghigh); } @@ -302,15 +237,15 @@ static void apply_cell_shift(UnitCell *cell, int k, double shift) 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; + case GPARAM_ASX : asx += shift; break; + case GPARAM_ASY : asy += shift; break; + case GPARAM_ASZ : asz += shift; break; + case GPARAM_BSX : bsx += shift; break; + case GPARAM_BSY : bsy += shift; break; + case GPARAM_BSZ : bsz += shift; break; + case GPARAM_CSX : csx += shift; break; + case GPARAM_CSY : csy += shift; break; + case GPARAM_CSZ : csz += shift; break; } cell_set_reciprocal(cell, asx, asy, asz, @@ -327,7 +262,7 @@ static void apply_shift(Crystal *cr, int k, double shift) switch ( k ) { - case REF_DIV : + case GPARAM_DIV : if ( isnan(shift) ) { ERROR("NaN divergence shift\n"); } else { @@ -336,21 +271,21 @@ static void apply_shift(Crystal *cr, int k, double shift) } break; - case REF_R : + case GPARAM_R : t = crystal_get_profile_radius(cr); t += shift; crystal_set_profile_radius(cr, t); 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 : + case GPARAM_ASX : + case GPARAM_ASY : + case GPARAM_ASZ : + case GPARAM_BSX : + case GPARAM_BSY : + case GPARAM_BSZ : + case GPARAM_CSX : + case GPARAM_CSY : + case GPARAM_CSZ : apply_cell_shift(crystal_get_cell(cr), k, shift); break; -- cgit v1.2.3