WIP on SCGaussian gradients

2 files changed, 66 insertions, 19 deletions

diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index e08d9120..6bd18e0e 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -131,7 +131,7 @@ double sphere_fraction(double rlow, double rhigh, double pr)
 }
 
 
-double gaussian_fraction(double rlow, double rhigh, double pr)
+double gaussian_fraction(double rlow, double rhigh, double R)
 {
 	double plow, phigh;
 	const double ng = 2.6;
@@ -144,13 +144,13 @@ double gaussian_fraction(double rlow, double rhigh, double pr)
 	 * 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;
+	if ( rlow < -R ) rlow = -R;
+	if ( rlow > +R ) rlow = +R;
+	if ( rhigh < -R ) rhigh = -R;
+	if ( rhigh > +R ) rhigh = +R;
 
-	plow = 0.5 * gsl_sf_erf(ng * rlow / (sqrt(2.0)*pr));
-	phigh = 0.5 * gsl_sf_erf(ng * rhigh / (sqrt(2.0)*pr));
+	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)));
 
 	return plow - phigh;
 }
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 5b7759a1..8bdc0956 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -56,18 +56,15 @@
 
 
 /* Returns dp(gauss)/dr at "r" */
-static double gaussian_fraction_gradient(double r, double pr)
+static double gaussian_fraction_gradient(double r, double R)
 {
-	double q, dpdq, dqdr;
+	const double ng = 2.6;
 
 	/* 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;
+	if ( r < -R ) return 0.0;
+	if ( r > +R ) return 0.0;
 
-	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;
+	return ng/(R*sqrt(2.0*M_PI)) * exp(-pow(r*ng/R, 2.0)/2.0);
 }
 
 
@@ -124,6 +121,56 @@ static double sphere_fraction_rgradient(double r, double R)
 }
 
 
+static double gaussian_fraction_rgradient(double r, double R)
+{
+	/* 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));
+}
+
+
+static double volume_fraction_rgradient(double r, double pr,
+                                       PartialityModel pmodel)
+{
+	switch ( pmodel )
+	{
+		case PMODEL_UNITY :
+		return 1.0;
+
+		case PMODEL_SCSPHERE :
+		return sphere_fraction_rgradient(r, pr);
+
+		case PMODEL_SCGAUSSIAN :
+		return gaussian_fraction_rgradient(r, pr);
+	}
+
+	ERROR("No pmodel in volume_fraction_rgradient!\n");
+	return 1.0;
+}
+
+
+static double volume_fraction(double rlow, double rhigh, double pr,
+                              PartialityModel pmodel)
+{
+	switch ( pmodel )
+	{
+		case PMODEL_UNITY :
+		return 1.0;
+
+		case PMODEL_SCSPHERE :
+		return sphere_fraction(rlow, rhigh, pr);
+
+		case PMODEL_SCGAUSSIAN :
+		return gaussian_fraction(rlow, rhigh, pr);
+	}
+
+	ERROR("No pmodel in volume_fraction!\n");
+	return 1.0;
+}
+
+
 /* 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)
@@ -151,10 +198,10 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
 		double Rglow, Rghigh;
 
 		D = rlow - rhigh;
-		psph = sphere_fraction(rlow, rhigh, R);
+		psph = volume_fraction(rlow, rhigh, R, pmodel);
 
-		Rglow = sphere_fraction_rgradient(rlow, R);
-		Rghigh = sphere_fraction_rgradient(rhigh, R);
+		Rglow = volume_fraction_rgradient(rlow, R, pmodel);
+		Rghigh = volume_fraction_rgradient(rhigh, R, pmodel);
 
 		return 4.0*psph/(3.0*D) + (4.0*R/(3.0*D))*(Rglow - Rghigh);
 
@@ -229,7 +276,7 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
 
 		case REF_DIV :
 		D = rlow - rhigh;
-		psph = sphere_fraction(rlow, rhigh, R);
+		psph = volume_fraction(rlow, rhigh, R, pmodel);
 		return (ds/2.0)*(glow+ghigh) - 4.0*R*psph*ds/(3.0*D*D);
 
 	}