aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2014-08-14 14:34:40 +0200
committerThomas White <taw@physics.org>2014-09-25 10:53:56 +0200
commita06a3f67f57de0bc85982976b9ea6d598598e014 (patch)
tree6eab60a0a40acacc27f5e07264043bfd617c859f /src
parent4ad4127092681a2190c682a736d8baecb1554328 (diff)
WIP on gradients for scsphere
Diffstat (limited to 'src')
-rw-r--r--src/post-refinement.c27
1 files changed, 21 insertions, 6 deletions
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 99ddf6b1..4c16f2d8 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -68,9 +68,14 @@ static double dpdq(double r, double profile_radius)
/* Returns dp/dr at "r" */
static double partiality_gradient(double r, double profile_radius,
- PartialityModel pmodel)
+ PartialityModel pmodel,
+ double rlow, double rhigh, double p)
{
- double dqdr; /* dq/dr */
+ double dqdr; /* dq/dr */
+ double dpsphdr; /* dpsph / dr */
+ double A, D;
+
+ D = rlow - rhigh;
switch ( pmodel ) {
@@ -90,6 +95,12 @@ static double partiality_gradient(double r, double profile_radius,
case PMODEL_THIN:
return -(2.0*r)/(profile_radius*profile_radius);
+ 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;
+
}
}
@@ -178,12 +189,12 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
/* Calculate the gradient of partiality wrt excitation error. */
if ( clamp_low == 0 ) {
- glow = partiality_gradient(rlow, r, pmodel);
+ glow = partiality_gradient(rlow, r, pmodel, rlow, rhigh, p);
} else {
glow = 0.0;
}
if ( clamp_high == 0 ) {
- ghigh = partiality_gradient(rhigh, r, pmodel);
+ ghigh = partiality_gradient(rhigh, r, pmodel, rlow, rhigh, p);
} else {
ghigh = 0.0;
}
@@ -215,14 +226,15 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
default:
case PMODEL_UNITY:
+ case PMODEL_THIN:
return 0.0;
case PMODEL_GAUSSIAN:
case PMODEL_SPHERE:
return (ds*glow + ds*ghigh) / 2.0;
- case PMODEL_THIN:
- return 0.0;
+ case PMODEL_SCSPHERE:
+ return 0.0; /* FIXME */
}
}
@@ -273,6 +285,9 @@ double l_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
signed int hs, ks, ls;
double L;
+ /* L has a non-zero gradient only for div in sphere or gaussian model */
+ if ( (pmodel != PMODEL_SPHERE)
+ && (pmodel != PMODEL_GAUSSIAN) ) return 0.0;
if ( k != REF_DIV ) return 0.0;
get_symmetric_indices(refl, &hs, &ks, &ls);