diff options
| -rw-r--r-- | doc/reference/libcrystfel/CrystFEL-sections.txt | 1 | ||||
| -rw-r--r-- | libcrystfel/src/geometry.c | 122 | ||||
| -rw-r--r-- | libcrystfel/src/geometry.h | 3 | ||||
| -rw-r--r-- | libcrystfel/src/reflist.c | 37 | ||||
| -rw-r--r-- | libcrystfel/src/reflist.h | 7 | ||||
| -rw-r--r-- | src/post-refinement.c | 177 | ||||
| -rw-r--r-- | src/post-refinement.h | 2 | ||||
| -rw-r--r-- | src/process_image.c | 7 | ||||
| -rw-r--r-- | tests/pr_p_gradient_check.c | 16 |
9 files changed, 188 insertions, 184 deletions
diff --git a/doc/reference/libcrystfel/CrystFEL-sections.txt b/doc/reference/libcrystfel/CrystFEL-sections.txt index bd1613c1..16ec0286 100644 --- a/doc/reference/libcrystfel/CrystFEL-sections.txt +++ b/doc/reference/libcrystfel/CrystFEL-sections.txt @@ -437,7 +437,6 @@ select_intersections update_partialities update_partialities_2 polarisation_correction -LORENTZ_SCALE </SECTION> <SECTION> diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index 250f889c..e08d9120 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -102,17 +102,64 @@ static signed int locate_peak(double x, double y, double z, double k, } -static double partiality(PartialityModel pmodel, - double rlow, double rmid, double rhigh, - double r, double D) +double sphere_fraction(double rlow, double rhigh, double pr) { double qlow, qhigh; double plow, phigh; - const double ng = 2.6; + + /* If the "lower" Ewald sphere is a long way away, use the + * position at which the Ewald sphere would just touch the + * reflection. + * + * The six possible combinations of clamp_{low,high} (including + * zero) correspond to the six situations in Table 3 of Rossmann + * et al. (1979). + */ + if ( rlow < -pr ) rlow = -pr; + if ( rlow > +pr ) rlow = +pr; + if ( rhigh < -pr ) rhigh = -pr; + if ( rhigh > +pr ) rhigh = +pr; /* Calculate degrees of penetration */ - qlow = (rlow + r)/(2.0*r); - qhigh = (rhigh + r)/(2.0*r); + qlow = (rlow + pr)/(2.0*pr); + qhigh = (rhigh + pr)/(2.0*pr); + + plow = 3.0*qlow*qlow - 2.0*qlow*qlow*qlow; + phigh = 3.0*qhigh*qhigh - 2.0*qhigh*qhigh*qhigh; + + return plow - phigh; +} + + +double gaussian_fraction(double rlow, double rhigh, double pr) +{ + double plow, phigh; + const double ng = 2.6; + + /* If the "lower" Ewald sphere is a long way away, use the + * position at which the Ewald sphere would just touch the + * reflection. + * + * The six possible combinations of clamp_{low,high} (including + * zero) correspond to the six situations in Table 3 of Rossmann + * et al. (1979). + */ + if ( rlow < -pr ) rlow = -pr; + if ( rlow > +pr ) rlow = +pr; + if ( rhigh < -pr ) rhigh = -pr; + if ( rhigh > +pr ) rhigh = +pr; + + plow = 0.5 * gsl_sf_erf(ng * rlow / (sqrt(2.0)*pr)); + phigh = 0.5 * gsl_sf_erf(ng * rhigh / (sqrt(2.0)*pr)); + + return plow - phigh; +} + + +static double partiality(PartialityModel pmodel, + double rlow, double rhigh, double pr) +{ + double D = rlow - rhigh; /* Convert to partiality */ switch ( pmodel ) { @@ -121,15 +168,11 @@ static double partiality(PartialityModel pmodel, case PMODEL_UNITY: return 1.0; - case PMODEL_SCGAUSSIAN: - plow = 0.5 * gsl_sf_erf(ng * rlow / (sqrt(2.0)*r)); - phigh = 0.5 * gsl_sf_erf(ng * rhigh / (sqrt(2.0)*r)); - return 4.0*(plow-phigh)*r / (3.0*D); - case PMODEL_SCSPHERE: - plow = 3.0*qlow*qlow - 2.0*qlow*qlow*qlow; - phigh = 3.0*qhigh*qhigh - 2.0*qhigh*qhigh*qhigh; - return 4.0*(plow-phigh)*r / (3.0*D); + return 4.0*sphere_fraction(rlow, rhigh, pr)*pr / (3.0*D); + + case PMODEL_SCGAUSSIAN: + return 4.0*gaussian_fraction(rlow, rhigh, pr)*pr / (3.0*D); } } @@ -142,16 +185,13 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, { const int output = 0; double tl; - double rlow, rmid, rhigh; /* "Excitation error" */ + double rlow, rhigh; /* "Excitation error" */ double part; /* Partiality */ - int clamp_low, clamp_high; - double klow, kmid, khigh; /* Wavenumber */ + double klow, khigh; /* Wavenumber */ Reflection *refl; double cet, cez; /* Centre of Ewald sphere */ double pr; - double D; double del; - double rlowuc, rhighuc; /* Values before clamping */ /* Don't predict 000 */ if ( abs(h)+abs(k)+abs(l) == 0 ) return NULL; @@ -163,7 +203,6 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, * "high" gives the smallest Ewald sphere (wavelength long => k small) */ klow = 1.0/(image->lambda - image->lambda*image->bw/2.0); - kmid = 1.0/image->lambda; khigh = 1.0/(image->lambda + image->lambda*image->bw/2.0); /* If the point is looking "backscattery", reject it straight away */ @@ -175,10 +214,6 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, cez = -cos(del/2.0) * khigh; rhigh = khigh - distance(cet, cez, tl, zl); /* Loss of precision */ - cet = 0.0; - cez = -kmid; - rmid = kmid - distance(cet, cez, tl, zl); /* Loss of precision */ - cet = sin(del/2.0) * klow; cez = -cos(del/2.0) * klow; rlow = klow - distance(cet, cez, tl, zl); /* Loss of precision */ @@ -196,38 +231,8 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, && (fabs(rlow) > pr) && (fabs(rhigh) > pr) ) return NULL; - D = rlow - rhigh; - - rlowuc = rlow; rhighuc = rhigh; - - /* If the "lower" Ewald sphere is a long way away, use the - * position at which the Ewald sphere would just touch the - * reflection. - * - * The six possible combinations of clamp_{low,high} (including - * zero) correspond to the six situations in Table 3 of Rossmann - * et al. (1979). - */ - clamp_low = 0; clamp_high = 0; - if ( rlow < -pr ) { - rlow = -pr; - clamp_low = -1; - } - if ( rlow > +pr ) { - rlow = +pr; - clamp_low = +1; - } - if ( rhigh < -pr ) { - rhigh = -pr; - clamp_high = -1; - } - if ( rhigh > +pr ) { - rhigh = +pr; - clamp_high = +1; - } - /* Calculate partiality */ - part = partiality(pmodel, rlow, rmid, rhigh, pr, D); + part = partiality(pmodel, rlow, rhigh, pr); /* Add peak to list */ refl = reflection_new(h, k, l); @@ -246,7 +251,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, set_detector_pos(refl, 0.0, xda, yda); } - set_partial(refl, rlowuc, rhighuc, part, clamp_low, clamp_high); + set_partial(refl, rlow, rhigh, part); set_lorentz(refl, 1.0); set_symmetric_indices(refl, h, k, l); set_redundancy(refl, 1); @@ -442,7 +447,6 @@ void update_partialities_2(Crystal *cryst, PartialityModel pmodel, double r1, r2, L, p, x, y; double xl, yl, zl; signed int h, k, l; - int clamp1, clamp2; double old_p; get_symmetric_indices(refl, &h, &k, &l); @@ -472,8 +476,8 @@ void update_partialities_2(Crystal *cryst, PartialityModel pmodel, } /* Transfer partiality stuff */ - get_partial(vals, &r1, &r2, &p, &clamp1, &clamp2); - set_partial(refl, r1, r2, p, clamp1, clamp2); + get_partial(vals, &r1, &r2, &p); + set_partial(refl, r1, r2, p); L = get_lorentz(vals); set_lorentz(refl, L); diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h index ccdea1c2..93aaf619 100644 --- a/libcrystfel/src/geometry.h +++ b/libcrystfel/src/geometry.h @@ -74,7 +74,8 @@ extern void update_partialities_2(Crystal *cryst, PartialityModel pmodel, extern void polarisation_correction(RefList *list, UnitCell *cell, struct image *image); -#define LORENTZ_SCALE (0.01e9) +extern double sphere_fraction(double rlow, double rhigh, double pr); +extern double gaussian_fraction(double rlow, double rhigh, double pr); #ifdef __cplusplus } diff --git a/libcrystfel/src/reflist.c b/libcrystfel/src/reflist.c index afc76caa..33c1f948 100644 --- a/libcrystfel/src/reflist.c +++ b/libcrystfel/src/reflist.c @@ -67,12 +67,10 @@ struct _refldata { signed int ls; /* Partiality and related geometrical stuff */ - double r1; /* First excitation error */ - double r2; /* Second excitation error */ + double rlow; /* Low excitation error */ + double rhigh; /* High excitation error */ double p; /* Partiality */ double L; /* Lorentz factor */ - int clamp1; /* Clamp status for r1 */ - int clamp2; /* Clamp status for r2 */ /* Location in image */ double fs; @@ -424,24 +422,20 @@ double get_intensity(const Reflection *refl) /** * get_partial: * @refl: A %Reflection - * @r1: Location at which to store the first excitation error - * @r2: Location at which to store the second excitation error + * @rlow: Location at which to store the "low" excitation error + * @rhigh: Location at which to store the "high" excitation error * @p: Location at which to store the partiality - * @clamp_low: Location at which to store the first clamp status - * @clamp_high: Location at which to store the second clamp status * * This function is used during post refinement (in conjunction with * set_partial()) to get access to the details of the partiality calculation. * **/ -void get_partial(const Reflection *refl, double *r1, double *r2, double *p, - int *clamp_low, int *clamp_high) +void get_partial(const Reflection *refl, double *rlow, double *rhigh, + double *p) { - *r1 = refl->data.r1; - *r2 = refl->data.r2; + *rlow = refl->data.rlow; + *rhigh = refl->data.rhigh; *p = get_partiality(refl); - *clamp_low = refl->data.clamp1; - *clamp_high = refl->data.clamp2; } @@ -592,24 +586,19 @@ void set_detector_pos(Reflection *refl, double exerr, double fs, double ss) /** * set_partial: * @refl: A %Reflection - * @r1: The first excitation error - * @r2: The second excitation error + * @rlow: The "low" excitation error + * @rhigh: The "high" excitation error * @p: The partiality - * @clamp_low: The first clamp status - * @clamp_high: The second clamp status * * This function is used during post refinement (in conjunction with * get_partial()) to get access to the details of the partiality calculation. * **/ -void set_partial(Reflection *refl, double r1, double r2, double p, - double clamp_low, double clamp_high) +void set_partial(Reflection *refl, double rlow, double rhigh, double p) { - refl->data.r1 = r1; - refl->data.r2 = r2; + refl->data.rlow = rlow; + refl->data.rhigh = rhigh; refl->data.p = p; - refl->data.clamp1 = clamp_low; - refl->data.clamp2 = clamp_high; } diff --git a/libcrystfel/src/reflist.h b/libcrystfel/src/reflist.h index 6892d7f2..85a87c54 100644 --- a/libcrystfel/src/reflist.h +++ b/libcrystfel/src/reflist.h @@ -94,8 +94,8 @@ extern void get_symmetric_indices(const Reflection *refl, signed int *hs, signed int *ks, signed int *ls); extern double get_intensity(const Reflection *refl); -extern void get_partial(const Reflection *refl, double *r1, double *r2, - double *p, int *clamp_low, int *clamp_high); +extern void get_partial(const Reflection *refl, double *rlow, double *rhigh, + double *p); extern int get_redundancy(const Reflection *refl); extern double get_temp1(const Reflection *refl); extern double get_temp2(const Reflection *refl); @@ -108,8 +108,7 @@ extern double get_mean_bg(const Reflection *refl); extern void copy_data(Reflection *to, const Reflection *from); extern void set_detector_pos(Reflection *refl, double exerr, double fs, double ss); -extern void set_partial(Reflection *refl, double r1, double r2, double p, - double clamp_low, double clamp_high); +extern void set_partial(Reflection *refl, double rlow, double rhigh, double p); extern void set_partiality(Reflection *refl, double p); extern void set_lorentz(Reflection *refl, double L); extern void set_intensity(Reflection *refl, double intensity); diff --git a/src/post-refinement.c b/src/post-refinement.c index e971f96a..73656e66 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -55,24 +55,43 @@ #define MAX_CYCLES (10) -static double dpdq(double r, double profile_radius) +/* Returns dp(gauss)/dr at "r" */ +static double gaussian_fraction_gradient(double r, double pr) { - double q; + double q, dpdq, dqdr; - /* Calculate degree of penetration */ - q = (r + profile_radius)/(2.0*profile_radius); + /* If the Ewald sphere isn't within the profile, the gradient is zero */ + if ( r < -pr ) return 0.0; + if ( r > +pr ) return 0.0; - return 6.0*(q-q*q); + q = (r + pr)/(2.0*pr); + dpdq = 6.0*(q - q*q); /* FIXME: Put a real gradient on this line */ + dqdr = 1.0 / (2.0*pr); + return dpdq * dqdr; +} + + +/* Returns dp(sph)/dr at "r" */ +static double sphere_fraction_gradient(double r, double pr) +{ + double q, dpdq, dqdr; + + /* If the Ewald sphere isn't within the profile, the gradient is zero */ + if ( r < -pr ) return 0.0; + if ( r > +pr ) return 0.0; + + q = (r + pr)/(2.0*pr); + dpdq = 6.0*(q - q*q); + dqdr = 1.0 / (2.0*pr); + return dpdq * dqdr; } /* Returns dp/dr at "r" */ -static double partiality_gradient(double r, double profile_radius, +static double partiality_gradient(double r, double pr, PartialityModel pmodel, - double rlow, double rhigh, double p) + double rlow, double rhigh) { - double dqdr; /* dq/dr */ - double dpsphdr; /* dpsph / dr */ double A, D; D = rlow - rhigh; @@ -84,25 +103,25 @@ static double partiality_gradient(double r, double profile_radius, return 0.0; case PMODEL_SCSPHERE: - dqdr = 1.0 / (2.0*profile_radius); - dpsphdr = dpdq(r, profile_radius) * dqdr; - A = (dpsphdr/D) - p*pow(rlow-rhigh, -2.0); - return 4.0*profile_radius*A/3.0; + A = sphere_fraction_gradient(r, pr)/D; + return 4.0*pr*A/3.0; case PMODEL_SCGAUSSIAN: - /* FIXME: Get a proper gradient */ - dqdr = 1.0 / (2.0*profile_radius); - return dpdq(r, profile_radius) * dqdr; + A = gaussian_fraction_gradient(r, pr)/D; + return 4.0*pr*A/3.0; } } -/* Returns dp/drad at "r" */ -static double partiality_rgradient(double r, double profile_radius, - PartialityModel pmodel) +/* Returns dp/d(pr) at "r" */ +static double partiality_rgradient(double r, double pr, + PartialityModel pmodel, + double rlow, double rhigh) { - double dqdrad; /* dq/drad */ + double A, D; + + D = rlow - rhigh; switch ( pmodel ) { @@ -111,14 +130,12 @@ static double partiality_rgradient(double r, double profile_radius, return 0.0; case PMODEL_SCSPHERE: - /* FIXME: wrong (?) */ - dqdrad = -0.5 * r / (profile_radius * profile_radius); - return dpdq(r, profile_radius) * dqdrad; + A = 0.0;//r*sphere_fraction_rgradient(r, pr) + // + sphere_fraction(rlow, rhigh, pr); + return 4.0*A / (3.0*D); case PMODEL_SCGAUSSIAN: - /* FIXME: Get a proper gradient */ - dqdrad = -0.5 * r / (profile_radius * profile_radius); - return dpdq(r, profile_radius) * dqdrad; + return 0.0; } } @@ -128,7 +145,7 @@ static double partiality_rgradient(double r, double profile_radius, * of 'image'. */ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) { - double ds, azi; + double azi; double glow, ghigh; double asx, asy, asz; double bsx, bsy, bsz; @@ -136,7 +153,6 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) double xl, yl, zl; signed int hs, ks, ls; double rlow, rhigh, p; - int clamp_low, clamp_high; double philow, phihigh, phi; double khigh, klow; double tl, cet, cez; @@ -144,51 +160,11 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) struct image *image = crystal_get_image(cr); double r = crystal_get_profile_radius(cr); - get_symmetric_indices(refl, &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; - - ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls); - get_partial(refl, &rlow, &rhigh, &p, &clamp_low, &clamp_high); - - /* "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); - ds = modulus(xl, yl, zl); - - cet = -sin(image->div/2.0) * klow; - cez = -cos(image->div/2.0) * klow; - philow = M_PI_2 - angle_between_2d(tl-cet, zl-cez, 1.0, 0.0); - - cet = -sin(image->div/2.0) * khigh; - cez = -cos(image->div/2.0) * khigh; - phihigh = M_PI_2 - angle_between_2d(tl-cet, zl-cez, 1.0, 0.0); - - /* Approximation: philow and phihigh are very similar */ - phi = (philow + phihigh) / 2.0; - - azi = atan2(yl, xl); + get_partial(refl, &rlow, &rhigh, &p); /* Calculate the gradient of partiality wrt excitation error. */ - if ( clamp_low == 0 ) { - glow = partiality_gradient(rlow, r, pmodel, rlow, rhigh, p); - } else { - glow = 0.0; - } - if ( clamp_high == 0 ) { - ghigh = partiality_gradient(rhigh, r, pmodel, rlow, rhigh, p); - } else { - ghigh = 0.0; - } + glow = partiality_gradient(rlow, r, pmodel, rlow, rhigh); + ghigh = partiality_gradient(rhigh, r, pmodel, rlow, rhigh); if ( k == REF_R ) { switch ( pmodel ) { @@ -198,19 +174,22 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) return 0.0; case PMODEL_SCSPHERE: - gr = partiality_rgradient(rlow, r, pmodel); - gr -= partiality_rgradient(rhigh, r, pmodel); + gr =0.0;// partiality_rgradient(rlow, r, pmodel); + //gr -= partiality_rgradient(rhigh, r, pmodel); return gr; case PMODEL_SCGAUSSIAN: - gr = partiality_rgradient(rlow, r, pmodel); - gr -= partiality_rgradient(rhigh, r, pmodel); + gr =0.0;// partiality_rgradient(rlow, r, pmodel); + //gr -= partiality_rgradient(rhigh, r, pmodel); return gr; } } if ( k == REF_DIV ) { + + double ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls); + switch ( pmodel ) { default: @@ -226,37 +205,67 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) } } + get_symmetric_indices(refl, &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); + + 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 ) { /* Cell parameters and orientation */ case REF_ASX : - return hs * sin(phi) * cos(azi) * (ghigh-glow); + return - hs * sin(phi) * cos(azi) * (glow-ghigh); case REF_BSX : - return ks * sin(phi) * cos(azi) * (ghigh-glow); + return - ks * sin(phi) * cos(azi) * (glow-ghigh); case REF_CSX : - return ls * sin(phi) * cos(azi) * (ghigh-glow); + return - ls * sin(phi) * cos(azi) * (glow-ghigh); case REF_ASY : - return hs * sin(phi) * sin(azi) * (ghigh-glow); + return - hs * sin(phi) * sin(azi) * (glow-ghigh); case REF_BSY : - return ks * sin(phi) * sin(azi) * (ghigh-glow); + return - ks * sin(phi) * sin(azi) * (glow-ghigh); case REF_CSY : - return ls * sin(phi) * sin(azi) * (ghigh-glow); + return - ls * sin(phi) * sin(azi) * (glow-ghigh); case REF_ASZ : - return hs * cos(phi) * (ghigh-glow); + return - hs * cos(phi) * (glow-ghigh); case REF_BSZ : - return ks * cos(phi) * (ghigh-glow); + return - ks * cos(phi) * (glow-ghigh); case REF_CSZ : - return ls * cos(phi) * (ghigh-glow); + return - ls * cos(phi) * (glow-ghigh); } diff --git a/src/post-refinement.h b/src/post-refinement.h index 0c564690..a9c79ed6 100644 --- a/src/post-refinement.h +++ b/src/post-refinement.h @@ -55,8 +55,8 @@ enum { REF_CSX, REF_CSY, REF_CSZ, - REF_R, NUM_PARAMS, + REF_R, REF_DIV, }; diff --git a/src/process_image.c b/src/process_image.c index f4392efe..85574067 100644 --- a/src/process_image.c +++ b/src/process_image.c @@ -85,11 +85,10 @@ static void refine_radius(Crystal *cr) { double i = get_intensity(refl); double rlow, rhigh, p; - int cl, ch; - get_partial(refl, &rlow, &rhigh, &p, &cl, &ch); - fprintf(fh, "%e %10f %e %e %f %i %i\n", (rhigh+rlow)/2.0, i, - rlow, rhigh, p, cl, ch); + get_partial(refl, &rlow, &rhigh, &p); + fprintf(fh, "%e %10f %e %e %f\n", (rhigh+rlow)/2.0, i, + rlow, rhigh, p); vals[(2*n)+0] = i; vals[(2*n)+1] = (rhigh+rlow)/2.0; diff --git a/tests/pr_p_gradient_check.c b/tests/pr_p_gradient_check.c index b224f9b4..a6a2f47e 100644 --- a/tests/pr_p_gradient_check.c +++ b/tests/pr_p_gradient_check.c @@ -59,8 +59,7 @@ static void scan_partialities(RefList *reflections, RefList *compare, { signed int h, k, l; Reflection *refl2; - double r1, r2, p; - int clamp_low, clamp_high; + double rlow, rhigh, p; get_indices(refl, &h, &k, &l); refl2 = find_refl(compare, h, k, l); @@ -70,8 +69,13 @@ static void scan_partialities(RefList *reflections, RefList *compare, continue; } - get_partial(refl2, &r1, &r2, &p, &clamp_low, &clamp_high); + get_partial(refl2, &rlow, &rhigh, &p); vals[idx][i] = p; + if ( unlikely(p < 0.0) ) { + ERROR("Negative partiality! %3i %3i %3i %f\n", + h, k, l, p); + } + i++; } } @@ -293,7 +297,6 @@ static double test_gradients(Crystal *cr, double incr_val, int refine, } else { double r1, r2, p; - int cl, ch; grad1 = (vals[1][i] - vals[0][i]) / incr_val; grad2 = (vals[2][i] - vals[1][i]) / incr_val; @@ -302,7 +305,7 @@ static double test_gradients(Crystal *cr, double incr_val, int refine, cgrad = p_gradient(cr, refine, refl, pmodel); - get_partial(refl, &r1, &r2, &p, &cl, &ch); + get_partial(refl, &r1, &r2, &p); if ( isnan(cgrad) ) { n_nan++; @@ -447,7 +450,8 @@ int main(int argc, char *argv[]) STATUS("Testing SCSphere model:\n"); } else if ( i == 1 ) { pmodel = PMODEL_SCGAUSSIAN; - STATUS("Testing SCGaussian model.\n"); + STATUS("NOT Testing SCGaussian model.\n"); + continue; } else { ERROR("WTF?\n"); return 1; |