pattern_sim: read spectrum from file

3 files changed, 54 insertions, 0 deletions

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);