aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--libcrystfel/src/geometry.c42
-rw-r--r--libcrystfel/src/geometry.h8
-rw-r--r--src/hrs-scaling.c57
-rw-r--r--src/partialator.c2
-rw-r--r--src/post-refinement.c56
-rw-r--r--tests/pr_l_gradient_check.c12
-rw-r--r--tests/pr_p_gradient_check.c12
-rw-r--r--tests/pr_pl_gradient_check.c12
8 files changed, 107 insertions, 94 deletions
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index 0c0a09b8..357f64f8 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -3,11 +3,11 @@
*
* Geometry of diffraction
*
- * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2013 Thomas White <taw@physics.org>
+ * 2010-2014 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -34,6 +34,7 @@
#include <stdlib.h>
#include <assert.h>
#include <fenv.h>
+#include <gsl/gsl_sf_erf.h>
#include "utils.h"
#include "cell.h"
@@ -102,24 +103,44 @@ static signed int locate_peak(double x, double y, double z, double k,
}
-static double partiality(double rlow, double rhigh, double r)
+static double partiality(PartialityModel pmodel, double rlow, double rhigh,
+ double r)
{
double qlow, qhigh;
double plow, phigh;
+ const double ng = 2.6;
/* Calculate degrees of penetration */
qlow = (rlow + r)/(2.0*r);
qhigh = (rhigh + r)/(2.0*r);
/* Convert to partiality */
- 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);
+ switch ( pmodel ) {
+
+ default:
+ case PMODEL_UNITY:
+ plow = 1.0;
+ phigh = 0.0;
+ break;
+
+ 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;
+
+ 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;
}
static Reflection *check_reflection(struct image *image, Crystal *cryst,
+ PartialityModel pmodel,
signed int h, signed int k, signed int l,
double xl, double yl, double zl)
{
@@ -204,7 +225,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
}
/* Calculate partiality */
- part = partiality(rlow, rhigh, pr);
+ part = partiality(pmodel, rlow, rhigh, pr);
/* Add peak to list */
refl = reflection_new(h, k, l);
@@ -298,7 +319,8 @@ RefList *find_intersections(struct image *image, Crystal *cryst)
yl = h*asy + k*bsy + l*csy;
zl = h*asz + k*bsz + l*csz;
- refl = check_reflection(image, cryst, h, k, l, xl, yl, zl);
+ refl = check_reflection(image, cryst, PMODEL_SPHERE,
+ h, k, l, xl, yl, zl);
if ( refl != NULL ) {
add_refl_to_list(refl, reflections);
@@ -401,9 +423,6 @@ void update_partialities_2(Crystal *cryst, PartialityModel pmodel,
int n = 0;
if ( pmodel == PMODEL_UNITY ) {
- /* It isn't strictly necessary to set the partialities to 1,
- * because the scaling stuff will just a correction factor of
- * 1 anyway. This is just to help things make sense. */
set_unity_partialities(cryst);
return;
}
@@ -430,7 +449,8 @@ void update_partialities_2(Crystal *cryst, PartialityModel pmodel,
yl = h*asy + k*bsy + l*csy;
zl = h*asz + k*bsz + l*csz;
- vals = check_reflection(image, cryst, h, k, l, xl, yl, zl);
+ vals = check_reflection(image, cryst, pmodel,
+ h, k, l, xl, yl, zl);
if ( vals == NULL ) {
diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
index fa7ed53b..7188f6ed 100644
--- a/libcrystfel/src/geometry.h
+++ b/libcrystfel/src/geometry.h
@@ -3,12 +3,12 @@
*
* Geometry of diffraction
*
- * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
+ * Copyright © 2013-2014 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
* Copyright © 2012 Richard Kirian
*
* Authors:
- * 2010-2013 Thomas White <taw@physics.org>
+ * 2010-2014 Thomas White <taw@physics.org>
* 2012 Richard Kirian
*
* This file is part of CrystFEL.
@@ -49,6 +49,7 @@ extern "C" {
* @PMODEL_SPHERE : Intersection of sphere with excited volume of reciprocal
* space.
* @PMODEL_UNITY : Set all all partialities and Lorentz factors to 1.
+ * @PMODEL_GAUSSIAN : Gaussian profiles in 3D
*
* A %PartialityModel describes a geometrical model which can be used to
* calculate spot partialities and Lorentz correction factors.
@@ -57,6 +58,7 @@ typedef enum {
PMODEL_SPHERE,
PMODEL_UNITY,
+ PMODEL_GAUSSIAN,
} PartialityModel;
diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index 76d1ff75..ac9091ed 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -119,24 +119,7 @@ static void run_scale_job(void *vwargs, int cookie)
Ih = get_intensity(r);
}
- /* If you change this, be sure to also change
- * run_merge_job() and run_esd_job(). */
- switch ( wargs->pmodel ) {
-
- case PMODEL_UNITY :
- corr = 1.0;
- break;
-
- case PMODEL_SPHERE :
- corr = get_partiality(refl) * get_lorentz(refl);
- break;
-
- default :
- ERROR("Unrecognised partiality model!\n");
- abort();
- break;
-
- }
+ corr = get_partiality(refl) * get_lorentz(refl);
Ihl = get_intensity(refl) / corr;
@@ -285,24 +268,7 @@ static void run_merge_job(void *vwargs, int cookie)
}
- /* If you change this, be sure to also change
- * run_scale_job() and run_esd_job(). */
- switch ( wargs->pmodel ) {
-
- case PMODEL_UNITY :
- corr = 1.0;
- break;
-
- case PMODEL_SPHERE :
- corr = get_partiality(refl) * get_lorentz(refl);
- break;
-
- default :
- ERROR("Unrecognised partiality model!\n");
- abort();
- break;
-
- }
+ corr = get_partiality(refl) * get_lorentz(refl);
Ihl = get_intensity(refl) / corr;
@@ -421,24 +387,7 @@ static void run_esd_job(void *vwargs, int cookie)
num = get_temp1(f);
- /* If you change this, be sure to also change
- * run_scale_job() and run_merge_job(). */
- switch ( wargs->pmodel ) {
-
- case PMODEL_UNITY :
- corr = 1.0;
- break;
-
- case PMODEL_SPHERE :
- corr = get_partiality(refl) * get_lorentz(refl);
- break;;
-
- default :
- ERROR("Unrecognised partiality model!\n");
- abort();
- break;
-
- }
+ corr = get_partiality(refl) * get_lorentz(refl);
Ih = get_intensity(f);
Ihl = G * get_intensity(refl) / corr;
diff --git a/src/partialator.c b/src/partialator.c
index bc8704b0..c1dedad8 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -468,6 +468,8 @@ int main(int argc, char *argv[])
pmodel = PMODEL_SPHERE;
} else if ( strcmp(pmodel_str, "unity") == 0 ) {
pmodel = PMODEL_UNITY;
+ } else if ( strcmp(pmodel_str, "gaussian") == 0 ) {
+ pmodel = PMODEL_GAUSSIAN;
} else {
ERROR("Unknown partiality model '%s'.\n", pmodel_str);
return 1;
diff --git a/src/post-refinement.c b/src/post-refinement.c
index d5d8a4d1..8fd29783 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -3,11 +3,11 @@
*
* Post refinement
*
- * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2013 Thomas White <taw@physics.org>
+ * 2010-2014 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -52,39 +52,55 @@
#define MAX_CYCLES (10)
-/* Returns dp/dr at "r" */
-static double partiality_gradient(double r, double profile_radius)
+static double dpdq(double r, double profile_radius, PartialityModel pmodel)
{
- double q, dpdq, dqdr;
+ double q;
+ double ng = 3.0;
/* Calculate degree of penetration */
q = (r + profile_radius)/(2.0*profile_radius);
/* dp/dq */
- dpdq = 6.0*(q-pow(q, 2.0));
+ switch ( pmodel ) {
+
+ default:
+ case PMODEL_UNITY:
+ return 0.0;
+
+ case PMODEL_SPHERE:
+ return 6.0*(q-pow(q, 2.0));
+
+ case PMODEL_GAUSSIAN:
+ /* The flat sphere model is close enough */
+ return 6.0*(q-pow(q, 2.0));
+
+ }
+}
+
+
+/* Returns dp/dr at "r" */
+static double partiality_gradient(double r, double profile_radius,
+ PartialityModel pmodel)
+{
+ double dqdr;
/* dq/dr */
dqdr = 1.0 / (2.0*profile_radius);
- return dpdq * dqdr;
+ return dpdq(r, profile_radius, pmodel) * dqdr;
}
/* Returns dp/drad at "r" */
-static double partiality_rgradient(double r, double profile_radius)
+static double partiality_rgradient(double r, double profile_radius,
+ PartialityModel pmodel)
{
- double q, dpdq, dqdrad;
-
- /* Calculate degree of penetration */
- q = (r + profile_radius)/(2.0*profile_radius);
-
- /* dp/dq */
- dpdq = 6.0*(q-pow(q, 2.0));
+ double dqdrad;
/* dq/drad */
dqdrad = -0.5 * r * pow(profile_radius, -2.0);
- return dpdq * dqdrad;
+ return dpdq(r, profile_radius, pmodel) * dqdrad;
}
@@ -144,12 +160,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);
+ glow = partiality_gradient(rlow, r, pmodel);
} else {
glow = 0.0;
}
if ( clamp_high == 0 ) {
- ghigh = partiality_gradient(rhigh, r);
+ ghigh = partiality_gradient(rhigh, r, pmodel);
} else {
ghigh = 0.0;
}
@@ -163,8 +179,8 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
return (ds*glow + ds*ghigh) / 2.0;
case REF_R :
- gr = partiality_rgradient(rlow, r);
- gr -= partiality_rgradient(rhigh, r);
+ gr = partiality_rgradient(rlow, r, pmodel);
+ gr -= partiality_rgradient(rhigh, r, pmodel);
return gr;
/* Cell parameters and orientation */
diff --git a/tests/pr_l_gradient_check.c b/tests/pr_l_gradient_check.c
index 4b3894ba..72d353f9 100644
--- a/tests/pr_l_gradient_check.c
+++ b/tests/pr_l_gradient_check.c
@@ -279,7 +279,6 @@ int main(int argc, char *argv[])
Crystal *cr;
struct quaternion orientation;
int i;
- const PartialityModel pmodel = PMODEL_SPHERE;
int fail = 0;
int quiet = 0;
int plot = 0;
@@ -337,10 +336,19 @@ int main(int argc, char *argv[])
rng = gsl_rng_alloc(gsl_rng_mt19937);
- for ( i=0; i<1; i++ ) {
+ for ( i=0; i<2; i++ ) {
UnitCell *rot;
double val;
+ PartialityModel pmodel;
+
+ if ( i == 0 ) {
+ pmodel = PMODEL_SPHERE;
+ STATUS("Testing flat sphere model:\n");
+ } else {
+ pmodel = PMODEL_GAUSSIAN;
+ STATUS("Testing Gaussian model:\n");
+ }
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 c0582036..c3b55cfa 100644
--- a/tests/pr_p_gradient_check.c
+++ b/tests/pr_p_gradient_check.c
@@ -392,7 +392,6 @@ int main(int argc, char *argv[])
Crystal *cr;
struct quaternion orientation;
int i;
- const PartialityModel pmodel = PMODEL_SPHERE;
int fail = 0;
int quiet = 0;
int plot = 0;
@@ -450,10 +449,19 @@ int main(int argc, char *argv[])
rng = gsl_rng_alloc(gsl_rng_mt19937);
- for ( i=0; i<1; i++ ) {
+ for ( i=0; i<2; i++ ) {
UnitCell *rot;
double val;
+ PartialityModel pmodel;
+
+ if ( i == 0 ) {
+ pmodel = PMODEL_SPHERE;
+ STATUS("Testing flat sphere model:\n");
+ } else {
+ pmodel = PMODEL_GAUSSIAN;
+ STATUS("Testing Gaussian model:\n");
+ }
orientation = random_quaternion(rng);
rot = cell_rotate(cell, orientation);
diff --git a/tests/pr_pl_gradient_check.c b/tests/pr_pl_gradient_check.c
index 785d9973..ddbc1840 100644
--- a/tests/pr_pl_gradient_check.c
+++ b/tests/pr_pl_gradient_check.c
@@ -394,7 +394,6 @@ int main(int argc, char *argv[])
Crystal *cr;
struct quaternion orientation;
int i;
- const PartialityModel pmodel = PMODEL_SPHERE;
int fail = 0;
int quiet = 0;
int plot = 0;
@@ -452,10 +451,19 @@ int main(int argc, char *argv[])
rng = gsl_rng_alloc(gsl_rng_mt19937);
- for ( i=0; i<1; i++ ) {
+ for ( i=0; i<2; i++ ) {
UnitCell *rot;
double val;
+ PartialityModel pmodel;
+
+ if ( i == 0 ) {
+ pmodel = PMODEL_SPHERE;
+ STATUS("Testing flat sphere model:\n");
+ } else {
+ pmodel = PMODEL_GAUSSIAN;
+ STATUS("Testing Gaussian model:\n");
+ }
orientation = random_quaternion(rng);
rot = cell_rotate(cell, orientation);
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-28 06:42:23 +0000