aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--src/diffraction-gpu.c19
-rw-r--r--src/diffraction.c11
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);