From 925610d93bbbdf590a58c5bb9fdae2af5a8fd30f Mon Sep 17 00:00:00 2001
From: Alexandra Tolstikova <alexandra.tolstikova@desy.de>
Date: Wed, 22 Nov 2017 17:46:15 +0100
Subject: indexamajig: read pink beam spectrum from file

---
 libcrystfel/src/image.h | 12 +++++++++
 src/indexamajig.c       | 70 +++++++++++++++++++++++++++++++++++++++++++++++++
 src/process_image.c     |  5 ++++
 src/process_image.h     |  1 +
 4 files changed, 88 insertions(+)

diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h
index 357a550a..1cd9f35a 100644
--- a/libcrystfel/src/image.h
+++ b/libcrystfel/src/image.h
@@ -123,6 +123,15 @@ enum imagefile_type
 /* An opaque type representing a list of image features */
 typedef struct _imagefeaturelist ImageFeatureList;
 
+
+struct spectrum
+{
+	int n;
+	double *ks; /* 1/m */
+	double *weights;
+};
+
+
 /* Structure describing a wavelength sample from a spectrum */
 struct sample
 {
@@ -213,6 +222,9 @@ struct image {
 	int                     serial;  /* Monotonically ascending serial
 	                                  * number for this image */
 
+	struct spectrum *spectrum; /* Beam spectrum for pink beam data */
+
+	// These only used in pattern_sim, to be changed to struct spectrum from above later...
 	struct sample           *spectrum0;
 	int                     nsamples; /* Number of wavelengths */
 	int                     spectrum_size;  /* Size of "spectrum" */
diff --git a/src/indexamajig.c b/src/indexamajig.c
index d4d04986..d6ebe3ea 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -64,6 +64,7 @@
 #include "integration.h"
 #include "taketwo.h"
 #include "im-sandbox.h"
+#include "image.h"
 
 
 static void show_help(const char *s)
@@ -240,6 +241,57 @@ static void add_geom_beam_stuff_to_field_list(struct imagefile_field_list *copym
 }
 
 
+static struct spectrum *read_spectrum_fromfile(char *fn)
+{
+	FILE *f;
+	f = fopen(fn, "r");
+
+	struct spectrum *s = malloc(sizeof(struct spectrum));
+
+	if ( fscanf(f, "%d", &s->n) == EOF ) {
+		return NULL;
+	}
+
+	if ( s->n <= 0 ) {
+		return NULL;
+	}
+
+	s->ks = malloc(s->n * sizeof(double));
+	if ( s->ks == NULL ) {
+		ERROR("Failed to allocate spectrum!\n");
+		return NULL;
+	}
+
+	s->weights = malloc(s->n * sizeof(double));
+	if ( s->weights == NULL ) {
+		ERROR("Failed to allocate spectrum!\n");
+		return NULL;
+	}
+
+	int i;
+	double k, w;
+	double w_sum = 0;
+	for ( i = 0; i < s->n; i++ ) {
+		if (fscanf(f, "%lf %lf", &k, &w) != EOF) {
+			s->ks[i] = ph_eV_to_k(k);
+			s->weights[i] = w;
+			w_sum += w;
+		} else break;
+	}
+
+	if ( i < s->n - 1 ) {
+		ERROR("Failed to read %d lines from %s\n", s->n, fn);
+		return NULL;
+	}
+
+	for ( i = 0; i < s->n; i++ ) {
+		s->weights[i] /= w_sum;
+	}
+
+	return s;
+}
+
+
 int main(int argc, char *argv[])
 {
 	int c;
@@ -276,6 +328,7 @@ int main(int argc, char *argv[])
 	int if_multi = 0;
 	int if_retry = 1;
 	int serial_start = 1;
+	char *spectrum_fn = NULL;
 
 	/* Defaults */
 	iargs.cell = NULL;
@@ -465,6 +518,7 @@ int main(int argc, char *argv[])
 		{"xgandalf-min-lvl",                         1, NULL, 355},
 		{"xgandalf-max-lattice-vector-length",       1, NULL, 356},
 		{"xgandalf-max-lvl",                         1, NULL, 356},
+	        {"spectrum-file",            1, NULL,        357},
 
 		{0, 0, NULL, 0}
 	};
@@ -855,6 +909,9 @@ int main(int argc, char *argv[])
 			}
 			break;
 
+			case 367:
+			spectrum_fn = strdup(optarg);
+			break;
 
 			case 0 :
 			break;
@@ -1043,6 +1100,19 @@ int main(int argc, char *argv[])
 		iargs.cell = NULL;
 	}
 
+	/* Load spectrum from file if given */
+	if ( spectrum_fn != NULL ) {
+		iargs.spectrum = read_spectrum_fromfile(spectrum_fn);
+		if ( iargs.spectrum == NULL ) {
+			ERROR("Couldn't read spectrum (from %s)\n", spectrum_fn);
+			return 1;
+		}
+		free(spectrum_fn);
+		STATUS("Read %d lines from %s\n", iargs.spectrum->n, spectrum_fn);
+	} else {
+		iargs.spectrum = NULL;
+	}
+
 	/* Parse integration diagnostic */
 	if ( int_diag != NULL ) {
 
diff --git a/src/process_image.c b/src/process_image.c
index 94c0dc26..89887436 100644
--- a/src/process_image.c
+++ b/src/process_image.c
@@ -272,6 +272,11 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 		image.bw = 0.00000001;
 	}
 
+	/* Set beam spectrum for pink beam data */
+	if ( iargs->spectrum != NULL ) {
+		image.spectrum = iargs->spectrum;
+	}
+
 	if ( image_feature_count(image.features) < iargs->min_peaks ) {
 		r = chdir(rn);
 		if ( r ) {
diff --git a/src/process_image.h b/src/process_image.h
index 2b58d7bc..14b2aae8 100644
--- a/src/process_image.h
+++ b/src/process_image.h
@@ -112,6 +112,7 @@ struct index_args
 	struct taketwo_options taketwo_opts;
     struct xgandalf_options xgandalf_opts;
 	struct felix_options felix_opts;
+	struct spectrum *spectrum;
 };
 
 
-- 
cgit v1.2.3