From b2a198a4a7935d4c81c0b7044d9f89e3c6932472 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 20 May 2014 17:43:18 +0200 Subject: Add Gaussian partiality model --- libcrystfel/src/geometry.c | 42 +++++++++++++++++++++++++++++++----------- libcrystfel/src/geometry.h | 8 +++++--- 2 files changed, 36 insertions(+), 14 deletions(-) (limited to 'libcrystfel') 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 + * 2010-2014 Thomas White * * This file is part of CrystFEL. * @@ -34,6 +34,7 @@ #include #include #include +#include #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 + * 2010-2014 Thomas White * 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; -- cgit v1.2.3