diff options
| -rw-r--r-- | libcrystfel/src/geometry.c | 40 | ||||
| -rw-r--r-- | libcrystfel/src/geometry.h | 2 | ||||
| -rw-r--r-- | src/partialator.c | 2 | ||||
| -rw-r--r-- | src/post-refinement.c | 14 |
4 files changed, 51 insertions, 7 deletions
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index d8b40a16..cdd936e2 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -164,7 +164,42 @@ double gaussian_fraction(double rlow, double rhigh, double R) } +static double random_partiality(signed int h, signed int k, signed int l, + int serial) +{ + gsl_rng *rng = gsl_rng_alloc(gsl_rng_mt19937); + unsigned long int seed; + double p; + int i; + + gsl_rng_set(rng, serial); + for ( i=0; i<abs(h)+1; i++ ) { + seed = gsl_rng_get(rng); + } + if ( h < 0 ) seed = gsl_rng_get(rng); + gsl_rng_set(rng, seed); + + for ( i=0; i<abs(k)+1; i++ ) { + seed = gsl_rng_get(rng); + } + if ( k < 0 ) seed = gsl_rng_get(rng); + gsl_rng_set(rng, seed); + + for ( i=0; i<abs(l)+1; i++ ) { + seed = gsl_rng_get(rng); + } + if ( l < 0 ) seed = gsl_rng_get(rng); + gsl_rng_set(rng, seed); + + p = gsl_rng_uniform(rng); + gsl_rng_free(rng); + return p; +} + + static double partiality(PartialityModel pmodel, + signed int h, signed int k, signed int l, + int serial, double rlow, double rhigh, double pr) { double D = rlow - rhigh; @@ -182,6 +217,9 @@ static double partiality(PartialityModel pmodel, case PMODEL_SCGAUSSIAN: return 4.0*gaussian_fraction(rlow, rhigh, pr)*pr / (3.0*D); + case PMODEL_RANDOM: + return random_partiality(h, k, l, serial); + } } @@ -234,7 +272,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, && (fabs(rhigh) > pr) ) return NULL; /* Calculate partiality */ - part = partiality(pmodel, rlow, rhigh, pr); + part = partiality(pmodel, h, k, l, image->serial, rlow, rhigh, pr); if ( isnan(part) ) { ERROR("Assigning NAN partiality!\n"); diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h index 474755b2..6b035175 100644 --- a/libcrystfel/src/geometry.h +++ b/libcrystfel/src/geometry.h @@ -49,6 +49,7 @@ extern "C" { * @PMODEL_UNITY : Set all all partialities and Lorentz factors to 1. * @PMODEL_SCSPHERE : Sphere model with source coverage factor included * @PMODEL_SCGAUSSIAN : Gaussian model with source coverage factor included + * @PMODEL_RANDOM : Randomly assigned partialities * * A %PartialityModel describes a geometrical model which can be used to * calculate spot partialities and Lorentz correction factors. @@ -58,6 +59,7 @@ typedef enum { PMODEL_UNITY, PMODEL_SCSPHERE, PMODEL_SCGAUSSIAN, + PMODEL_RANDOM, } PartialityModel; diff --git a/src/partialator.c b/src/partialator.c index 70c82960..d539f784 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -530,6 +530,8 @@ int main(int argc, char *argv[]) pmodel = PMODEL_SCGAUSSIAN; } else if ( strcmp(pmodel_str, "scsphere") == 0 ) { pmodel = PMODEL_SCSPHERE; + } else if ( strcmp(pmodel_str, "random") == 0 ) { + pmodel = PMODEL_RANDOM; } else { ERROR("Unknown partiality model '%s'.\n", pmodel_str); return 1; diff --git a/src/post-refinement.c b/src/post-refinement.c index 9054ed89..12a29cb2 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -147,10 +147,11 @@ static double volume_fraction_rgradient(double r, double pr, case PMODEL_SCGAUSSIAN : return gaussian_fraction_rgradient(r, pr); - } - ERROR("No pmodel in volume_fraction_rgradient!\n"); - return 1.0; + default : + ERROR("No pmodel in volume_fraction_rgradient!\n"); + return 1.0; + } } @@ -167,10 +168,11 @@ static double volume_fraction(double rlow, double rhigh, double pr, case PMODEL_SCGAUSSIAN : return gaussian_fraction(rlow, rhigh, pr); - } - ERROR("No pmodel in volume_fraction!\n"); - return 1.0; + default : + ERROR("No pmodel in volume_fraction!\n"); + return 1.0; + } } |