diff options
| -rw-r--r-- | libcrystfel/src/geometry.c | 5 | ||||
| -rw-r--r-- | src/post-refinement.c | 8 | ||||
| -rw-r--r-- | tests/pr_p_gradient_check.c | 3 |
3 files changed, 10 insertions, 6 deletions
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index 6bd18e0e..221fe96d 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -135,6 +135,7 @@ double gaussian_fraction(double rlow, double rhigh, double R) { double plow, phigh; const double ng = 2.6; + const double sig = R/ng; /* If the "lower" Ewald sphere is a long way away, use the * position at which the Ewald sphere would just touch the @@ -149,8 +150,8 @@ double gaussian_fraction(double rlow, double rhigh, double R) if ( rhigh < -R ) rhigh = -R; if ( rhigh > +R ) rhigh = +R; - plow = 0.5 * (1.0 + gsl_sf_erf(pow(ng*rlow /R, 2.0)/sqrt(2.0))); - phigh = 0.5 * (1.0 + gsl_sf_erf(pow(ng*rhigh/R, 2.0)/sqrt(2.0))); + plow = 0.5*(1.0 + gsl_sf_erf(rlow/(sig*sqrt(2.0)))); + phigh = 0.5*(1.0 + gsl_sf_erf(rhigh/(sig*sqrt(2.0)))); return plow - phigh; } diff --git a/src/post-refinement.c b/src/post-refinement.c index 8bdc0956..4dae5096 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -59,12 +59,13 @@ static double gaussian_fraction_gradient(double r, double R) { const double ng = 2.6; + const double sig = R/ng; /* If the Ewald sphere isn't within the profile, the gradient is zero */ if ( r < -R ) return 0.0; if ( r > +R ) return 0.0; - return ng/(R*sqrt(2.0*M_PI)) * exp(-pow(r*ng/R, 2.0)/2.0); + return exp(-pow(r/sig, 2.0)/2.0) / (sig*sqrt(2.0*M_PI)); } @@ -123,11 +124,14 @@ static double sphere_fraction_rgradient(double r, double R) static double gaussian_fraction_rgradient(double r, double R) { + const double ng = 2.6; + const double sig = R/ng; + /* If the Ewald sphere isn't within the profile, the gradient is zero */ if ( r < -R ) return 0.0; if ( r > +R ) return 0.0; - return -(3.0*r/(4.0*R*R)) * (1.0 - r*r/(R*R)); + return -(ng*r/(sqrt(2.0*M_PI)*R*R))*exp(-r*r/(2.0*sig*sig)); } diff --git a/tests/pr_p_gradient_check.c b/tests/pr_p_gradient_check.c index a6a2f47e..cedfc9dc 100644 --- a/tests/pr_p_gradient_check.c +++ b/tests/pr_p_gradient_check.c @@ -450,8 +450,7 @@ int main(int argc, char *argv[]) STATUS("Testing SCSphere model:\n"); } else if ( i == 1 ) { pmodel = PMODEL_SCGAUSSIAN; - STATUS("NOT Testing SCGaussian model.\n"); - continue; + STATUS("Testing SCGaussian model.\n"); } else { ERROR("WTF?\n"); return 1; |