aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/partialator.c2
-rw-r--r--src/post-refinement.c135
2 files changed, 81 insertions, 56 deletions
diff --git a/src/partialator.c b/src/partialator.c
index 2ef01dba..e5083611 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -348,6 +348,8 @@ int main(int argc, char *argv[])
pmodel = PMODEL_UNITY;
} else if ( strcmp(pmodel_str, "gaussian") == 0 ) {
pmodel = PMODEL_GAUSSIAN;
+ } else if ( strcmp(pmodel_str, "thin") == 0 ) {
+ pmodel = PMODEL_THIN;
} else {
ERROR("Unknown partiality model '%s'.\n", pmodel_str);
return 1;
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 5de8a246..99ddf6b1 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -55,15 +55,23 @@
#define MAX_CYCLES (10)
-static double dpdq(double r, double profile_radius, PartialityModel pmodel)
+static double dpdq(double r, double profile_radius)
{
double q;
- double ng = 3.0;
/* Calculate degree of penetration */
q = (r + profile_radius)/(2.0*profile_radius);
- /* dp/dq */
+ return 6.0*(q-q*q);
+}
+
+
+/* Returns dp/dr at "r" */
+static double partiality_gradient(double r, double profile_radius,
+ PartialityModel pmodel)
+{
+ double dqdr; /* dq/dr */
+
switch ( pmodel ) {
default:
@@ -71,39 +79,46 @@ static double dpdq(double r, double profile_radius, PartialityModel pmodel)
return 0.0;
case PMODEL_SPHERE:
- return 6.0*(q-pow(q, 2.0));
+ dqdr = 1.0 / (2.0*profile_radius);
+ return dpdq(r, profile_radius) * dqdr;
case PMODEL_GAUSSIAN:
- /* The flat sphere model is close enough */
- return 6.0*(q-pow(q, 2.0));
+ /* FIXME: Get a proper gradient */
+ dqdr = 1.0 / (2.0*profile_radius);
+ return dpdq(r, profile_radius) * dqdr;
+
+ case PMODEL_THIN:
+ return -(2.0*r)/(profile_radius*profile_radius);
}
}
-/* Returns dp/dr at "r" */
-static double partiality_gradient(double r, double profile_radius,
- PartialityModel pmodel)
+/* Returns dp/drad at "r" */
+static double partiality_rgradient(double r, double profile_radius,
+ PartialityModel pmodel)
{
- double dqdr;
+ double dqdrad; /* dq/drad */
- /* dq/dr */
- dqdr = 1.0 / (2.0*profile_radius);
+ switch ( pmodel ) {
- return dpdq(r, profile_radius, pmodel) * dqdr;
-}
+ default:
+ case PMODEL_UNITY:
+ return 0.0;
+ case PMODEL_SPHERE:
+ dqdrad = -0.5 * r / (profile_radius * profile_radius);
+ return dpdq(r, profile_radius) * dqdrad;
-/* Returns dp/drad at "r" */
-static double partiality_rgradient(double r, double profile_radius,
- PartialityModel pmodel)
-{
- double dqdrad;
+ case PMODEL_GAUSSIAN:
+ /* FIXME: Get a proper gradient */
+ dqdrad = -0.5 * r / (profile_radius * profile_radius);
+ return dpdq(r, profile_radius) * dqdrad;
- /* dq/drad */
- dqdrad = -0.5 * r * pow(profile_radius, -2.0);
+ case PMODEL_THIN:
+ return 2.0*r*r*pow(profile_radius, -3.0);
- return dpdq(r, profile_radius, pmodel) * dqdrad;
+ }
}
@@ -173,19 +188,48 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
ghigh = 0.0;
}
+ if ( k == REF_R ) {
+ switch ( pmodel ) {
+
+ default:
+ case PMODEL_UNITY:
+ return 0.0;
+
+ case PMODEL_GAUSSIAN:
+ gr = partiality_rgradient(rlow, r, pmodel);
+ gr -= partiality_rgradient(rhigh, r, pmodel);
+ return gr;
+
+ case PMODEL_SPHERE:
+ gr = partiality_rgradient(rlow, r, pmodel);
+ gr -= partiality_rgradient(rhigh, r, pmodel);
+ return gr;
+
+ case PMODEL_THIN:
+ return 2.0*rlow*rlow/(r*r*r);
+ }
+ }
+
+ if ( k == REF_DIV ) {
+ switch ( pmodel ) {
+
+ default:
+ case PMODEL_UNITY:
+ return 0.0;
+
+ case PMODEL_GAUSSIAN:
+ case PMODEL_SPHERE:
+ return (ds*glow + ds*ghigh) / 2.0;
+
+ case PMODEL_THIN:
+ return 0.0;
+ }
+ }
+
/* For many gradients, just multiply the above number by the gradient
* of excitation error wrt whatever. */
switch ( k ) {
- case REF_DIV :
- /* Small angle approximation */
- return (ds*glow + ds*ghigh) / 2.0;
-
- case REF_R :
- gr = partiality_rgradient(rlow, r, pmodel);
- gr -= partiality_rgradient(rhigh, r, pmodel);
- return gr;
-
/* Cell parameters and orientation */
case REF_ASX :
return hs * sin(phi) * cos(azi) * (ghigh-glow);
@@ -227,37 +271,16 @@ double l_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
{
double ds;
signed int hs, ks, ls;
+ double L;
- switch ( k ) {
-
- /* Cell parameters do not affect Lorentz factor */
- case REF_ASX :
- case REF_BSX :
- case REF_CSX :
- case REF_ASY :
- case REF_BSY :
- case REF_CSY :
- case REF_ASZ :
- case REF_BSZ :
- case REF_CSZ :
- return 0.0;
-
- /* Nor does change of radius */
- case REF_R :
- return 0.0;
-
- default:
- break;
-
- }
-
- assert(k == REF_DIV);
+ if ( k != REF_DIV ) return 0.0;
get_symmetric_indices(refl, &hs, &ks, &ls);
ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls);
- return -ds*pow(get_lorentz(refl), 2.0) / LORENTZ_SCALE;
+ L = get_lorentz(refl);
+ return -ds*L*L / LORENTZ_SCALE;
}