Allow molecule_factor() to tell when it doesn't know the intensity

8 files changed, 40 insertions, 22 deletions

diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index 68920a55..50ca6282 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -98,12 +98,12 @@ int main(int argc, char *argv[])
 	}
 
 	cell = load_cell_from_pdb("molecule.pdb");
-	ref1 = read_reflections(afile);
+	ref1 = read_reflections(afile, NULL);
 	if ( ref1 == NULL ) {
 		ERROR("Couldn't open file '%s'\n", afile);
 		return 1;
 	}
-	ref2 = read_reflections(bfile);
+	ref2 = read_reflections(bfile, NULL);
 	if ( ref2 == NULL ) {
 		ERROR("Couldn't open file '%s'\n", bfile);
 		return 1;
diff --git a/src/diffraction.c b/src/diffraction.c
index 5ecf04f5..675aa148 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -70,10 +70,11 @@ static double lattice_factor(struct rvec q, double ax, double ay, double az,
 
 
 /* Look up the structure factor for the nearest Bragg condition */
-static double molecule_factor(double *intensities, struct rvec q,
-                                      double ax, double ay, double az,
-                                      double bx, double by, double bz,
-                                      double cx, double cy, double cz)
+static double molecule_factor(double *intensities, unsigned int *counts,
+                              struct rvec q,
+                              double ax, double ay, double az,
+                              double bx, double by, double bz,
+                              double cx, double cy, double cz)
 {
 	double hd, kd, ld;
 	signed int h, k, l;
@@ -86,8 +87,13 @@ static double molecule_factor(double *intensities, struct rvec q,
 	k = (signed int)rint(kd);
 	l = (signed int)rint(ld);
 
-	r = lookup_intensity(intensities, h, k, l);
+	if ( lookup_count(counts, h, k, l) == 0 ) {
+		ERROR("Needed intensity for %i %i %i, but don't have it.\n",
+		      h, k, l);
+		return 1.0e20;
+	}
 
+	r = lookup_intensity(intensities, h, k, l);
 	return r;
 }
 
@@ -171,7 +177,8 @@ struct rvec get_q(struct image *image, unsigned int xs, unsigned int ys,
 
 
 void get_diffraction(struct image *image, int na, int nb, int nc,
-                     double *intensities, UnitCell *cell, int do_water)
+                     double *intensities, unsigned int *counts, UnitCell *cell,
+                     int do_water)
 {
 	unsigned int xs, ys;
 	double ax, ay, az;
@@ -225,8 +232,10 @@ void get_diffraction(struct image *image, int na, int nb, int nc,
 			if ( intensities == NULL ) {
 				I_molecule = 1.0e10;
 			} else {
-				I_molecule = molecule_factor(intensities, q,
-			                            ax,ay,az,bx,by,bz,cx,cy,cz);
+				I_molecule = molecule_factor(intensities,
+				                             counts, q,
+			                                     ax,ay,az,
+			                                     bx,by,bz,cx,cy,cz);
 			}
 
 			I_lattice = pow(f_lattice, 2.0);
diff --git a/src/diffraction.h b/src/diffraction.h
index 58ae10c3..c317060f 100644
--- a/src/diffraction.h
+++ b/src/diffraction.h
@@ -20,7 +20,8 @@
 #include "cell.h"
 
 extern void get_diffraction(struct image *image, int na, int nb, int nc,
-                            double *intensities, UnitCell *cell, int do_water);
+                            double *intensities, unsigned int *counts,
+                            UnitCell *cell, int do_water);
 extern struct rvec get_q(struct image *image, unsigned int xs, unsigned int ys,
                          unsigned int sampling, float *ttp, float k);
 
diff --git a/src/indexamajig.c b/src/indexamajig.c
index c99a2b8b..c24b0e6b 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -146,7 +146,8 @@ static struct image *get_simage(struct image *template, int alternate)
 
 
 static void simulate_and_write(struct image *simage, struct gpu_context **gctx,
-                               double *intensities, UnitCell *cell)
+                               double *intensities, unsigned int *counts,
+                               UnitCell *cell)
 {
 	/* Set up GPU if necessary */
 	if ( (gctx != NULL) && (*gctx == NULL) ) {
@@ -156,7 +157,8 @@ static void simulate_and_write(struct image *simage, struct gpu_context **gctx,
 	if ( (gctx != NULL) && (*gctx != NULL) ) {
 		get_diffraction_gpu(*gctx, simage, 24, 24, 40, cell);
 	} else {
-		get_diffraction(simage, 24, 24, 40, intensities, cell, 0);
+		get_diffraction(simage, 24, 24, 40,
+		                intensities, counts, cell, 0);
 	}
 	record_image(simage, 0);
 
@@ -191,6 +193,7 @@ int main(int argc, char *argv[])
 	UnitCell *cell;
 	double *intensities = NULL;
 	char *intfile = NULL;
+	unsigned int *counts = NULL;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -262,9 +265,10 @@ int main(int argc, char *argv[])
 	}
 
 	if ( intfile != NULL ) {
-		intensities = read_reflections(intfile);
+		intensities = read_reflections(intfile, counts);
 	} else {
 		intensities = NULL;
+		counts = NULL;
 	}
 
 	if ( indm_str == NULL ) {
@@ -386,10 +390,10 @@ int main(int argc, char *argv[])
 		if ( config_simulate ) {
 			if ( config_gpu ) {
 				simulate_and_write(simage, &gctx, intensities,
-				                   cell);
+				                   counts, cell);
 			} else {
 				simulate_and_write(simage, NULL, intensities,
-				                   cell);
+				                   counts, cell);
 			}
 		}
 
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 55b92c0e..3fa1155f 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -163,6 +163,7 @@ int main(int argc, char *argv[])
 	int n_images = 1; /* Generate one image by default */
 	int done = 0;
 	UnitCell *cell;
+	unsigned int *counts;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -232,8 +233,10 @@ int main(int argc, char *argv[])
 		STATUS("reflection intensities (with --intensities).\n");
 		STATUS("I'll simulate a flat intensity distribution.\n");
 		intensities = NULL;
+		counts = NULL;
 	} else {
-		intensities = read_reflections(intfile);
+		counts = new_list_count();
+		intensities = read_reflections(intfile, counts);
 		free(intfile);
 	}
 
@@ -293,8 +296,8 @@ int main(int argc, char *argv[])
 			}
 			get_diffraction_gpu(gctx, &image, na, nb, nc, cell);
 		} else {
-			get_diffraction(&image, na, nb, nc, intensities, cell,
-			                !config_nowater);
+			get_diffraction(&image, na, nb, nc, intensities, counts,
+			                cell, !config_nowater);
 		}
 		if ( image.data == NULL ) {
 			ERROR("Diffraction calculation failed.\n");
diff --git a/src/process_hkl.c b/src/process_hkl.c
index 7b1035bc..8bc76e5e 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -268,7 +268,7 @@ int main(int argc, char *argv[])
 
 	if ( intfile != NULL ) {
 		STATUS("Comparing against '%s'\n", intfile);
-		trueref = read_reflections(intfile);
+		trueref = read_reflections(intfile, NULL);
 		free(intfile);
 	} else {
 		trueref = NULL;
diff --git a/src/reflections.c b/src/reflections.c
index deba0f8d..c943d014 100644
--- a/src/reflections.c
+++ b/src/reflections.c
@@ -84,7 +84,7 @@ void write_reflections(const char *filename, unsigned int *counts,
 }
 
 
-double *read_reflections(const char *filename)
+double *read_reflections(const char *filename, unsigned int *counts)
 {
 	double *ref;
 	FILE *fh;
@@ -110,6 +110,7 @@ double *read_reflections(const char *filename)
 		if ( r != 4 ) continue;
 
 		set_intensity(ref, h, k, l, intensity);
+		if ( counts != NULL ) set_count(counts, h, k, l, 1);
 
 	} while ( rval != NULL );
 
diff --git a/src/reflections.h b/src/reflections.h
index 14209c34..8c513a64 100644
--- a/src/reflections.h
+++ b/src/reflections.h
@@ -23,7 +23,7 @@
 extern void write_reflections(const char *filename, unsigned int *counts,
                               double *ref, int zone_axis, UnitCell *cell);
 
-extern double *read_reflections(const char *filename);
+extern double *read_reflections(const char *filename, unsigned int *counts);
 
 extern double *ideal_intensities(double complex *sfac);