From 17aec9c77e9b75510e43b0fb073df75f66849761 Mon Sep 17 00:00:00 2001
From: Alexandra Tolstikova <alexandra.tolstikova@desy.de>
Date: Thu, 16 Nov 2017 13:46:23 +0100
Subject: pattern_sim: read spectrum from file

---
 libcrystfel/src/image.h |  6 ++++--
 libcrystfel/src/utils.h |  2 ++
 src/diffraction.c       | 34 ++++++++++++++++++++++++++++++++++
 src/diffraction.h       |  3 +++
 src/pattern_sim.c       | 17 +++++++++++++++++
 5 files changed, 60 insertions(+), 2 deletions(-)

diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h
index dffcdb3a..29ce8ae1 100644
--- a/libcrystfel/src/image.h
+++ b/libcrystfel/src/image.h
@@ -60,6 +60,7 @@ struct imagefile_field_list;
  * @SPECTRUM_TOPHAT: A top hat distribution of wavelengths
  * @SPECTRUM_SASE: A simulated SASE spectrum
  * @SPECTRUM_TWOCOLOUR: A spectrum containing two peaks
+ * @SPECTRUM_FROMFILE: An arbitrary spectrum read from input file
  *
  * A %SpectrumType represents a type of X-ray energy spectrum to use for
  * generating simulated data.
@@ -67,7 +68,8 @@ struct imagefile_field_list;
 typedef enum {
 	SPECTRUM_TOPHAT,
 	SPECTRUM_SASE,
-	SPECTRUM_TWOCOLOUR
+	SPECTRUM_TWOCOLOUR,
+	SPECTRUM_FROMFILE
 } SpectrumType;
 
 
@@ -124,7 +126,7 @@ typedef struct _imagefeaturelist ImageFeatureList;
 /* Structure describing a wavelength sample from a spectrum */
 struct sample
 {
-	double k;
+	double k; /* 1/m */
 	double weight;
 };
 
diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h
index a759ff15..2c756ad7 100644
--- a/libcrystfel/src/utils.h
+++ b/libcrystfel/src/utils.h
@@ -185,6 +185,8 @@ static inline int within_tolerance(double a, double b, double percent)
 /* Photon energy (eV) to wavelength (m) */
 #define ph_eV_to_lambda(a) ph_en_to_lambda(eV_to_J(a))
 
+/* Photon energy (eV) to k (1/m) */
+#define ph_eV_to_k(a) ((a)*ELECTRON_CHARGE/PLANCK/C_VACUO)
 
 
 /* ------------------------------ Message macros ---------------------------- */
diff --git a/src/diffraction.c b/src/diffraction.c
index 4cceba0e..97750f58 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -608,6 +608,40 @@ struct sample *generate_tophat(struct image *image)
 }
 
 
+struct sample *generate_spectrum_fromfile(struct image *image,
+					  char *spectrum_fn)
+{
+	struct sample *spectrum;
+	int i;
+	double k, w;
+	double w_sum = 0;
+
+	spectrum = malloc(image->nsamples * sizeof(struct sample));
+	if ( spectrum == NULL ) return NULL;
+
+	FILE *f;
+	f = fopen(spectrum_fn, "r");
+
+	int nsamples = 0;
+	for ( i=0; i<image->nsamples; i++ ) {
+		if (fscanf(f, "%lf %lf", &k, &w) != EOF) {
+			spectrum[i].k = ph_eV_to_k(k);
+			spectrum[i].weight = w;
+			w_sum += w;
+			nsamples += 1;
+		} else break;
+	}
+
+	for ( i=0; i<nsamples; i++ ) {
+		spectrum[i].weight /= w_sum;
+	}
+
+	image->spectrum_size = nsamples;
+
+	return spectrum;
+}
+
+
 struct sample *generate_SASE(struct image *image, gsl_rng *rng)
 {
 	struct sample *spectrum;
diff --git a/src/diffraction.h b/src/diffraction.h
index ca4e3bde..0f95f83c 100644
--- a/src/diffraction.h
+++ b/src/diffraction.h
@@ -60,4 +60,7 @@ extern struct sample *generate_SASE(struct image *image, gsl_rng *rng);
 
 extern struct sample *generate_twocolour(struct image *image);
 
+extern struct sample *generate_spectrum_fromfile(struct image *image,
+					         char *spectrum_fn);
+
 #endif	/* DIFFRACTION_H */
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 6c9eb880..d47c4b60 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -368,6 +368,7 @@ int main(int argc, char *argv[])
 	char *outfile = NULL;
 	char *geometry = NULL;
 	char *spectrum_str = NULL;
+	char *spectrum_fn = NULL;
 	GradientMethod grad;
 	SpectrumType spectrum_type;
 	int ndone = 0;    /* Number of simulations done (images or not) */
@@ -429,6 +430,7 @@ int main(int argc, char *argv[])
 		{"photon-energy",      1, NULL,                9},
 		{"nphotons",           1, NULL,               10},
 		{"beam-radius",        1, NULL,               11},
+		{"spectrum-file",      1, NULL,               12},
 
 		{0, 0, NULL, 0}
 	};
@@ -590,6 +592,9 @@ int main(int argc, char *argv[])
 			}
 			break;
 
+			case 12 :
+			spectrum_fn = strdup(optarg);
+			break;
 
 			case 0 :
 			break;
@@ -695,6 +700,8 @@ int main(int argc, char *argv[])
                     strcasecmp(spectrum_str, "twocolours") == 0 ||
                     strcasecmp(spectrum_str, "twocolors") == 0) {
 		spectrum_type = SPECTRUM_TWOCOLOUR;
+	} else if ( strcasecmp(spectrum_str, "fromfile") == 0) {
+		spectrum_type = SPECTRUM_FROMFILE;
 	} else {
 		ERROR("Unrecognised spectrum type '%s'\n", spectrum_str);
 		return 1;
@@ -855,6 +862,11 @@ int main(int argc, char *argv[])
 		STATUS("                 X-ray spectrum: two colour, "
 		       "separation %.5f %%\n", image.bw*100.0);
 		break;
+
+		case SPECTRUM_FROMFILE:
+		STATUS("                 X-ray spectrum: from %s\n",
+			spectrum_fn);
+		break;
 	}
 	if ( random_size ) {
 		STATUS("                   Crystal size: random, between "
@@ -977,6 +989,11 @@ int main(int argc, char *argv[])
                         image.spectrum = generate_twocolour(&image);
                         break;
 
+                        case SPECTRUM_FROMFILE :
+                        image.spectrum = generate_spectrum_fromfile(&image,
+								   spectrum_fn);
+                        break;
+
 		}
 
 		cell_get_parameters(cell, &a, &b, &c, &d, &d, &d);
-- 
cgit v1.2.3