diff options
| -rw-r--r-- | src/diffraction-gpu.c | 19 | ||||
| -rw-r--r-- | src/diffraction.c | 11 |
2 files changed, 20 insertions, 10 deletions
diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c index e2bce0cc..4ecf0558 100644 --- a/src/diffraction-gpu.c +++ b/src/diffraction-gpu.c @@ -294,8 +294,7 @@ int get_diffraction_gpu(struct gpu_context *gctx, struct image *image, int n_neg = 0; int n_nan = 0; int i; - double kmin, kmax, step; - double tot = 0.0; + double kmin, kmax, step, norm; if ( gctx == NULL ) { ERROR("GPU setup failed.\n"); @@ -341,24 +340,28 @@ int get_diffraction_gpu(struct gpu_context *gctx, struct image *image, } spectrum_get_range(image->spectrum, &kmin, &kmax); - step = (kmax-kmin)/n_samples; - for ( i=0; i<=n_samples; i++ ) { + step = (kmax-kmin)/(n_samples+1); + + /* Determine normalisation factor such that weights add up to 1 after + * sampling (bins must have constant width) */ + for ( i=1; i<=n_samples; i++ ) { + double k = kmin + i*step; + norm += spectrum_get_density_at_k(image->spectrum, k); + } + for ( i=1; i<=n_samples; i++ ) { double k = kmin + i*step; double prob; /* Probability = p.d.f. times step width */ - prob = step * spectrum_get_density_at_k(image->spectrum, k); + prob = spectrum_get_density_at_k(image->spectrum, k)/norm; STATUS("Wavelength: %e m, weight = %.5f\n", 1.0/k, prob); err = do_panels(gctx, image, k, prob, &n_inf, &n_neg, &n_nan); if ( err ) return 1; - tot += prob; - } - STATUS("Total weight = %f\n", tot); if ( n_neg + n_inf + n_nan ) { ERROR("WARNING: The GPU calculation produced %i negative" diff --git a/src/diffraction.c b/src/diffraction.c index bcbc62e8..6f52df11 100644 --- a/src/diffraction.c +++ b/src/diffraction.c @@ -457,6 +457,7 @@ void get_diffraction(struct image *image, int na, int nb, int nc, double *lut_c; int i; double kmin, kmax, step; + double norm = 0.0; cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); @@ -466,14 +467,20 @@ void get_diffraction(struct image *image, int na, int nb, int nc, spectrum_get_range(image->spectrum, &kmin, &kmax); step = (kmax-kmin)/(n_samples+1); - STATUS("%i samples\n", n_samples); + + /* Determine normalisation factor such that weights add up to 1 after + * sampling (bins must have constant width) */ + for ( i=1; i<=n_samples; i++ ) { + double k = kmin + i*step; + norm += spectrum_get_density_at_k(image->spectrum, k); + } for ( i=1; i<=n_samples; i++ ) { double k = kmin + i*step; double prob; /* Probability = p.d.f. times step width */ - prob = step * spectrum_get_density_at_k(image->spectrum, k); + prob = spectrum_get_density_at_k(image->spectrum, k)/norm; STATUS("Wavelength: %e m, weight = %.5f\n", 1.0/k, prob); |