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 --- src/hrs-scaling.c | 57 +++------------------------------------------------ src/partialator.c | 2 ++ src/post-refinement.c | 56 ++++++++++++++++++++++++++++++++------------------ 3 files changed, 41 insertions(+), 74 deletions(-) (limited to 'src') 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 + * 2010-2014 Thomas White * * 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 */ -- cgit v1.2.3