Two colour fixes

2 files changed, 48 insertions, 37 deletions

diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c
index 1ed008e8..b5bd2e88 100644
--- a/src/diffraction-gpu.c
+++ b/src/diffraction-gpu.c
@@ -325,7 +325,7 @@ void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
 	double tot = 0.0;
 	for ( i=0; i<image->nsamples; i++ ) {
 
-		printf("%.1f eV, weight = %.5f\n",
+		printf("%.3f eV, weight = %.5f\n",
 		       ph_lambda_to_eV(1.0/image->spectrum[i].k),
 		       image->spectrum[i].weight);
 
diff --git a/src/diffraction.c b/src/diffraction.c
index 82e5c132..0fa9f79b 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -434,7 +434,8 @@ static int compare_samples(const void *a, const void *b)
 
 
 static struct sample *get_gaussian_spectrum(double eV_cen, double eV_step,
-                                            double sigma, int spec_size, double eV_start)
+                                            double sigma, int spec_size,
+                                            double eV_start)
 {
 	struct sample *spectrum;
 	int i;
@@ -596,61 +597,71 @@ struct sample *generate_SASE(struct image *image, gsl_rng *rng)
 	return spectrum;
 }
 
+
 struct sample *generate_twocolour(struct image *image)
 {
 	struct sample *spectrum;
 	struct sample *spectrum1;
 	struct sample *spectrum2;
-        int i;
+	int i;
+	double eV_cen1;  /* Central photon energy for first colour */
+	double eV_cen2;  /* Central photon energy for second colour */
+	double eV_cen;  /* Central photon energy for this spectrum */
+	const int spec_size = 1024;
 
-        double halfwidth = image->bw/2; /* eV */
+	eV_cen = ph_lambda_to_eV(image->lambda);
 
-	const int spec_size = 1024;
-        double eV_cen1;  /* Central photon energy for first colour */
-        double eV_cen2;  /* Central photon energy for second colour */
-	eV_cen1 = ph_lambda_to_eV(image->lambda) - halfwidth;
-	eV_cen2 = ph_lambda_to_eV(image->lambda) + halfwidth;
+	double halfwidth = eV_cen*image->bw/2.0; /* eV */
+
+	eV_cen1 = eV_cen - halfwidth;
+	eV_cen2 = eV_cen + halfwidth;
 
-        /* Hard-code sigma to be 1/5 of bandwidth */
-        double sigma = image->bw/5; /* eV */
+	/* Hard-code sigma to be 1/5 of bandwidth */
+	double sigma = eV_cen*image->bw/5.0; /* eV */
 
-        /* The spectrum will be calculated to a resolution which spreads six
-         * sigmas of the original (no SASE noise) Gaussian pulse over spec_size
-         * points */
-        double eV_start = eV_cen1 - 3*sigma;
-        double eV_end = eV_cen2 + 3*sigma;
-        double eV_step = (eV_end - eV_start)/(spec_size-1);
+	/* The spectrum will be calculated to a resolution which spreads six
+	* sigmas of the original (no SASE noise) Gaussian pulse over spec_size
+	* points */
+	double eV_start = eV_cen1 - 3*sigma;
+	double eV_end = eV_cen2 + 3*sigma;
+	double eV_step = (eV_end - eV_start)/(spec_size-1);
 
-	spectrum1 = get_gaussian_spectrum(eV_cen1, eV_step, sigma, spec_size, eV_start);
-	spectrum2 = get_gaussian_spectrum(eV_cen2, eV_step, sigma, spec_size, eV_start);
+	spectrum1 = get_gaussian_spectrum(eV_cen1, eV_step, sigma, spec_size,
+	                                  eV_start);
+	spectrum2 = get_gaussian_spectrum(eV_cen2, eV_step, sigma, spec_size,
+	                                  eV_start);
 
 	spectrum = malloc(spec_size * sizeof(struct sample));
-        if ( spectrum == NULL ) return NULL;
+	if ( spectrum == NULL ) return NULL;
 
-        for ( i=0; i<spec_size; i++ ) {
-                spectrum[i].weight = spectrum1[i].weight + spectrum2[i].weight;
-                spectrum[i].k = spectrum1[i].k;
-        }
+	for ( i=0; i<spec_size; i++ ) {
+		spectrum[i].weight = spectrum1[i].weight + spectrum2[i].weight;
+		spectrum[i].k = spectrum1[i].k;
+		if ( spectrum2[i].k != spectrum1[i].k ) {
+			printf("%e %e\n", spectrum1[i].k, spectrum2[i].k);
+		}
+	}
 
-        /* Normalise intensity (before taking restricted number of samples) */
-        double total_weight = 0.0;
-        for ( i=0; i<spec_size; i++ ) {
-                total_weight += spectrum[i].weight;
-        }
+	/* Normalise intensity (before taking restricted number of samples) */
+	double total_weight = 0.0;
+	for ( i=0; i<spec_size; i++ ) {
+		total_weight += spectrum[i].weight;
+	}
 
-        for ( i=0; i<spec_size; i++ ) {
-                spectrum[i].weight /= total_weight;
-        }
+	for ( i=0; i<spec_size; i++ ) {
+		spectrum[i].weight /= total_weight;
+	}
 
-        /* Sort samples in spectrum by weight.  Diffraction calculation will
-         * take the requested number, starting from the brightest */
-        qsort(spectrum, spec_size, sizeof(struct sample), compare_samples);
+	/* Sort samples in spectrum by weight.  Diffraction calculation will
+	* take the requested number, starting from the brightest */
+	qsort(spectrum, spec_size, sizeof(struct sample), compare_samples);
 
-        image->spectrum_size = spec_size;
+	image->spectrum_size = spec_size;
 
-        return spectrum;
+	return spectrum;
 }
 
+
 void get_diffraction(struct image *image, int na, int nb, int nc,
                      const double *intensities, const double *phases,
                      const unsigned char *flags, UnitCell *cell,