aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--libcrystfel/src/geometry.c83
-rw-r--r--libcrystfel/src/geometry.h2
-rw-r--r--src/partialator.c2
-rw-r--r--src/post-refinement.c135
-rw-r--r--tests/pr_l_gradient_check.c1
-rw-r--r--tests/pr_p_gradient_check.c9
-rw-r--r--tests/pr_pl_gradient_check.c3
7 files changed, 155 insertions, 80 deletions
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index 5cc46b2e..3586ed0b 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -103,7 +103,8 @@ static signed int locate_peak(double x, double y, double z, double k,
}
-static double partiality(PartialityModel pmodel, double rlow, double rhigh,
+static double partiality(PartialityModel pmodel,
+ double rlow, double rmid, double rhigh,
double r)
{
double qlow, qhigh;
@@ -119,23 +120,22 @@ static double partiality(PartialityModel pmodel, double rlow, double rhigh,
default:
case PMODEL_UNITY:
- plow = 1.0;
- phigh = 0.0;
- break;
+ return 1.0;
case PMODEL_SPHERE:
- plow = 3.0*pow(qlow,2.0) - 2.0*pow(qlow,3.0);
- phigh = 3.0*pow(qhigh,2.0) - 2.0*pow(qhigh,3.0);
- break;
+ plow = 3.0*qlow*qlow - 2.0*qlow*qlow*qlow;
+ phigh = 3.0*qhigh*qhigh - 2.0*qhigh*qhigh*qhigh;
+ return plow - phigh;
case PMODEL_GAUSSIAN:
plow = 0.5 * gsl_sf_erf(ng * rlow / (sqrt(2.0)*r));
phigh = 0.5 * gsl_sf_erf(ng * rhigh / (sqrt(2.0)*r));
- break;
+ return plow - phigh;
- }
+ case PMODEL_THIN:
+ return 1.0 - (rmid*rmid)/(r*r);
- return plow - phigh;
+ }
}
@@ -146,12 +146,12 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
{
const int output = 0;
double tl;
- double rlow, rhigh; /* "Excitation error" */
+ double rlow, rmid, rhigh; /* "Excitation error" */
double part; /* Partiality */
int clamp_low, clamp_high;
- double klow, khigh; /* Wavenumber */
+ double klow, kmid, khigh; /* Wavenumber */
Reflection *refl;
- double cet, cez;
+ double cet, cez; /* Centre of Ewald sphere */
double pr;
double L;
double del;
@@ -166,6 +166,7 @@ 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 */
@@ -177,14 +178,14 @@ 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 */
- if ( (signbit(rlow) == signbit(rhigh))
- && (fabs(rlow) > pr)
- && (fabs(rhigh) > pr) ) return NULL;
-
if ( unlikely(rlow < rhigh) ) {
ERROR("Reflection with rlow < rhigh!\n");
ERROR("%3i %3i %3i rlow = %e, rhigh = %e\n",
@@ -193,10 +194,41 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
return NULL;
}
+ /* Conditions for reflection to be excited at all */
+ switch ( pmodel ) {
+
+ default:
+ case PMODEL_UNITY: /* PMODEL_UNITY shouldn't end up here */
+ case PMODEL_SPHERE:
+ case PMODEL_GAUSSIAN:
+ if ( (signbit(rlow) == signbit(rhigh))
+ && (fabs(rlow) > pr)
+ && (fabs(rhigh) > pr) ) return NULL;
+ break;
+
+ case PMODEL_THIN:
+ if ( fabs(rmid) > pr ) return NULL;
+ break;
+
+ }
+
/* Lorentz factor is determined direction from the r values, before
* clamping. The multiplication by 0.01e9 to make the
* correction factor vaguely near 1. */
- L = LORENTZ_SCALE / (rlow - rhigh);
+ switch ( pmodel ) {
+
+ default:
+ case PMODEL_SPHERE:
+ case PMODEL_GAUSSIAN:
+ L = LORENTZ_SCALE / (rlow - rhigh);
+ break;
+
+ case PMODEL_UNITY: /* PMODEL_UNITY shouldn't end up here */
+ case PMODEL_THIN:
+ L = 1.0;
+ break;
+
+ }
/* If the "lower" Ewald sphere is a long way away, use the
* position at which the Ewald sphere would just touch the
@@ -225,7 +257,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
}
/* Calculate partiality */
- part = partiality(pmodel, rlow, rhigh, pr);
+ part = partiality(pmodel, rlow, rmid, rhigh, pr);
/* Add peak to list */
refl = reflection_new(h, k, l);
@@ -244,7 +276,18 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
set_detector_pos(refl, 0.0, xda, yda);
}
- set_partial(refl, rlow, rhigh, part, clamp_low, clamp_high);
+ if ( pmodel != PMODEL_THIN ) {
+ set_partial(refl, rlow, rhigh, part, clamp_low, clamp_high);
+ } else {
+ /* If we are using the TES (Thin Ewald Sphere) model, we abuse
+ * the fields as follows:
+ * rlow = the r value for the middle (only) Ewald sphere
+ * rhigh = 0.0
+ * clamp_low = 0
+ * clamp_high = +1
+ */
+ set_partial(refl, rmid, 0.0, part, 0, +1);
+ }
set_lorentz(refl, L);
set_symmetric_indices(refl, h, k, l);
set_redundancy(refl, 1);
diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
index 1f465167..8041936a 100644
--- a/libcrystfel/src/geometry.h
+++ b/libcrystfel/src/geometry.h
@@ -50,6 +50,7 @@ extern "C" {
* space.
* @PMODEL_UNITY : Set all all partialities and Lorentz factors to 1.
* @PMODEL_GAUSSIAN : Gaussian profiles in 3D
+ * @PMODEL_THIN : Thin Ewald sphere intersecting sphere
*
* A %PartialityModel describes a geometrical model which can be used to
* calculate spot partialities and Lorentz correction factors.
@@ -59,6 +60,7 @@ typedef enum {
PMODEL_SPHERE,
PMODEL_UNITY,
PMODEL_GAUSSIAN,
+ PMODEL_THIN,
} PartialityModel;
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;
}
diff --git a/tests/pr_l_gradient_check.c b/tests/pr_l_gradient_check.c
index b707629f..7342d5f8 100644
--- a/tests/pr_l_gradient_check.c
+++ b/tests/pr_l_gradient_check.c
@@ -350,6 +350,7 @@ int main(int argc, char *argv[])
pmodel = PMODEL_GAUSSIAN;
STATUS("Testing Gaussian model:\n");
}
+ /* No point testing TES model, because it has no Lorentz factor */
orientation = random_quaternion(rng);
rot = cell_rotate(cell, orientation);
diff --git a/tests/pr_p_gradient_check.c b/tests/pr_p_gradient_check.c
index 8f7a66e4..14a2b579 100644
--- a/tests/pr_p_gradient_check.c
+++ b/tests/pr_p_gradient_check.c
@@ -174,6 +174,7 @@ static Crystal *new_shifted_crystal(Crystal *cr, int refine, double incr_val)
return cr_new;
}
+
static void calc_either_side(Crystal *cr, double incr_val,
int *valid, long double *vals[3], int refine,
PartialityModel pmodel)
@@ -225,7 +226,6 @@ static void calc_either_side(Crystal *cr, double incr_val,
}
-
static double test_gradients(Crystal *cr, double incr_val, int refine,
const char *str, const char *file,
PartialityModel pmodel, int quiet, int plot)
@@ -450,7 +450,7 @@ int main(int argc, char *argv[])
rng = gsl_rng_alloc(gsl_rng_mt19937);
- for ( i=0; i<2; i++ ) {
+ for ( i=0; i<3; i++ ) {
UnitCell *rot;
double val;
@@ -459,9 +459,12 @@ int main(int argc, char *argv[])
if ( i == 0 ) {
pmodel = PMODEL_SPHERE;
STATUS("Testing flat sphere model:\n");
- } else {
+ } else if ( i == 1 ) {
pmodel = PMODEL_GAUSSIAN;
STATUS("Testing Gaussian model:\n");
+ } else {
+ pmodel = PMODEL_THIN;
+ STATUS("Testing thin Ewald sphere model:\n");
}
orientation = random_quaternion(rng);
diff --git a/tests/pr_pl_gradient_check.c b/tests/pr_pl_gradient_check.c
index f83af335..69087a4b 100644
--- a/tests/pr_pl_gradient_check.c
+++ b/tests/pr_pl_gradient_check.c
@@ -174,6 +174,7 @@ static Crystal *new_shifted_crystal(Crystal *cr, int refine, double incr_val)
return cr_new;
}
+
static void calc_either_side(Crystal *cr, double incr_val,
int *valid, long double *vals[3], int refine,
PartialityModel pmodel)
@@ -225,7 +226,6 @@ static void calc_either_side(Crystal *cr, double incr_val,
}
-
static double test_gradients(Crystal *cr, double incr_val, int refine,
const char *str, const char *file,
PartialityModel pmodel, int quiet, int plot)
@@ -465,6 +465,7 @@ int main(int argc, char *argv[])
pmodel = PMODEL_GAUSSIAN;
STATUS("Testing Gaussian model:\n");
}
+ /* No point testing TES model */
orientation = random_quaternion(rng);
rot = cell_rotate(cell, orientation);
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-28 07:08:08 +0000