Don't try to render PDBs, part II: remove "molecule" from "struct image".

14 files changed, 104 insertions, 60 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index e696fcbf..ef9338c5 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -23,7 +23,7 @@ process_hkl_LDADD = @LIBS@
 
 indexamajig_SOURCES = indexamajig.c hdf5-file.c utils.c cell.c image.c \
                       intensities.c peaks.c index.c filters.c \
-                      diffraction.c detector.c sfac.c dirax.c
+                      diffraction.c detector.c sfac.c dirax.c reflections.c
 indexamajig_LDADD = @LIBS@
 if HAVE_OPENCL
 indexamajig_SOURCES += diffraction-gpu.c cl-utils.c
diff --git a/src/Makefile.in b/src/Makefile.in
index 87fb7e4c..4b59a13d 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -69,15 +69,16 @@ hdfsee_OBJECTS = $(am_hdfsee_OBJECTS)
 hdfsee_DEPENDENCIES =
 am__indexamajig_SOURCES_DIST = indexamajig.c hdf5-file.c utils.c \
 	cell.c image.c intensities.c peaks.c index.c filters.c \
-	diffraction.c detector.c sfac.c dirax.c diffraction-gpu.c \
-	cl-utils.c
+	diffraction.c detector.c sfac.c dirax.c reflections.c \
+	diffraction-gpu.c cl-utils.c
 @HAVE_OPENCL_TRUE@am__objects_1 = diffraction-gpu.$(OBJEXT) \
 @HAVE_OPENCL_TRUE@	cl-utils.$(OBJEXT)
 am_indexamajig_OBJECTS = indexamajig.$(OBJEXT) hdf5-file.$(OBJEXT) \
 	utils.$(OBJEXT) cell.$(OBJEXT) image.$(OBJEXT) \
 	intensities.$(OBJEXT) peaks.$(OBJEXT) index.$(OBJEXT) \
 	filters.$(OBJEXT) diffraction.$(OBJEXT) detector.$(OBJEXT) \
-	sfac.$(OBJEXT) dirax.$(OBJEXT) $(am__objects_1)
+	sfac.$(OBJEXT) dirax.$(OBJEXT) reflections.$(OBJEXT) \
+	$(am__objects_1)
 indexamajig_OBJECTS = $(am_indexamajig_OBJECTS)
 indexamajig_DEPENDENCIES =
 am__pattern_sim_SOURCES_DIST = pattern_sim.c diffraction.c utils.c \
@@ -231,7 +232,7 @@ process_hkl_SOURCES = process_hkl.c sfac.c statistics.c cell.c utils.c \
 process_hkl_LDADD = @LIBS@
 indexamajig_SOURCES = indexamajig.c hdf5-file.c utils.c cell.c image.c \
 	intensities.c peaks.c index.c filters.c diffraction.c \
-	detector.c sfac.c dirax.c $(am__append_3)
+	detector.c sfac.c dirax.c reflections.c $(am__append_3)
 indexamajig_LDADD = @LIBS@
 @HAVE_GTK_TRUE@hdfsee_SOURCES = hdfsee.c displaywindow.c render.c hdf5-file.c utils.c image.c \
 @HAVE_GTK_TRUE@                 filters.c
diff --git a/src/cell.c b/src/cell.c
index 721756c9..42f34dfa 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -615,3 +615,45 @@ double resolution(UnitCell *cell, signed int h, signed int k, signed int l)
 
 	return oneoverd / 2;
 }
+
+
+UnitCell *load_cell_from_pdb(const char *filename)
+{
+	FILE *fh;
+	char *rval;
+	UnitCell *cell = NULL;
+
+	fh = fopen(filename, "r");
+
+	do {
+
+		char line[1024];
+
+		rval = fgets(line, 1023, fh);
+
+		if ( strncmp(line, "CRYST1", 6) == 0 ) {
+
+			float a, b, c, al, be, ga;
+			int r;
+
+			r = sscanf(line+7, "%f %f %f %f %f %f",
+			           &a, &b, &c, &al, &be, &ga);
+			if ( r != 6 ) {
+				ERROR("Couldn't understand CRYST1 line\n");
+				return NULL;
+			}
+
+			cell = cell_new_from_parameters(a*1e-10,
+			                                     b*1e-10, c*1e-10,
+	                                                     deg2rad(al),
+	                                                     deg2rad(be),
+	                                                     deg2rad(ga));
+
+		}
+
+	} while ( rval != NULL );
+
+	fclose(fh);
+
+	return cell;
+}
diff --git a/src/cell.h b/src/cell.h
index 5a534268..924cb0ef 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -77,4 +77,6 @@ extern void cell_print(UnitCell *cell);
 
 extern UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose);
 
+extern UnitCell *load_cell_from_pdb(const char *filename);
+
 #endif	/* CELL_H */
diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c
index cc7b9e5b..741f7b69 100644
--- a/src/diffraction-gpu.c
+++ b/src/diffraction-gpu.c
@@ -107,7 +107,7 @@ static void check_sinc_lut(struct gpu_context *gctx, int n)
 
 
 void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
-                         int na, int nb, int nc, int no_sfac)
+                         int na, int nb, int nc, UnitCell *ucell)
 {
 	cl_int err;
 	double ax, ay, az;
@@ -130,9 +130,7 @@ void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
 		return;
 	}
 
-	cell_get_cartesian(image->molecule->cell, &ax, &ay, &az,
-		                                  &bx, &by, &bz,
-		                                  &cx, &cy, &cz);
+	cell_get_cartesian(ucell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
 	cell[0] = ax;  cell[1] = ay;  cell[2] = az;
 	cell[3] = bx;  cell[4] = by;  cell[5] = bz;
 	cell[6] = cx;  cell[7] = cy;  cell[8] = cz;
diff --git a/src/diffraction-gpu.h b/src/diffraction-gpu.h
index b23dab63..a9e46b66 100644
--- a/src/diffraction-gpu.h
+++ b/src/diffraction-gpu.h
@@ -24,7 +24,7 @@ struct gpu_context;
 #if HAVE_OPENCL
 
 extern void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
-                                int na, int nb, int nc);
+                                int na, int nb, int nc, UnitCell *ucell);
 extern struct gpu_context *setup_gpu(int no_sfac, struct image *image,
                                      double *intensities);
 extern void cleanup_gpu(struct gpu_context *gctx);
@@ -32,7 +32,7 @@ extern void cleanup_gpu(struct gpu_context *gctx);
 #else
 
 static void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
-                                int na, int nb, int nc)
+                                int na, int nb, int nc, UnitCell *ucell)
 {
 	/* Do nothing */
 	ERROR("This copy of CrystFEL was not compiled with OpenCL support.\n");
diff --git a/src/diffraction.c b/src/diffraction.c
index 01624be8..c47920f5 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -171,7 +171,7 @@ 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, int do_water)
+                     double *intensities, UnitCell *cell, int do_water)
 {
 	unsigned int xs, ys;
 	double ax, ay, az;
@@ -179,11 +179,7 @@ void get_diffraction(struct image *image, int na, int nb, int nc,
 	double cx, cy, cz;
 	float k, klow, bwstep;
 
-	if ( image->molecule == NULL ) return;
-
-	cell_get_cartesian(image->molecule->cell, &ax, &ay, &az,
-		                                  &bx, &by, &bz,
-		                                  &cx, &cy, &cz);
+	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
 
 	/* Allocate (and zero) the "diffraction array" */
 	image->data = calloc(image->width * image->height, sizeof(float));
diff --git a/src/diffraction.h b/src/diffraction.h
index dd0795fb..58ae10c3 100644
--- a/src/diffraction.h
+++ b/src/diffraction.h
@@ -20,7 +20,7 @@
 #include "cell.h"
 
 extern void get_diffraction(struct image *image, int na, int nb, int nc,
-                            double *intensities, int do_water);
+                            double *intensities, 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/get_hkl.c b/src/get_hkl.c
index 29ff0786..5e5331b4 100644
--- a/src/get_hkl.c
+++ b/src/get_hkl.c
@@ -110,8 +110,6 @@ int main(int argc, char *argv[])
 	char *output = NULL;
 	unsigned int *counts;
 	signed int h, k, l;
-	FILE *fh;
-	char *rval;
 
 	/* Long options */
 	const struct option longopts[] = {
diff --git a/src/image.h b/src/image.h
index 5f2cf3d0..2885b734 100644
--- a/src/image.h
+++ b/src/image.h
@@ -72,7 +72,6 @@ struct image {
 
 	float			*data;
 	double			*twotheta;
-	struct molecule		*molecule;
 	UnitCell		*indexed_cell;
 	UnitCell		*candidate_cells[MAX_CELL_CANDIDATES];
 	int			ncells;
diff --git a/src/index.c b/src/index.c
index daa16b27..f2b8ffb9 100644
--- a/src/index.c
+++ b/src/index.c
@@ -57,8 +57,8 @@ static void write_drx(struct image *image)
 }
 
 
-void index_pattern(struct image *image, IndexingMethod indm, int no_match,
-                   int verbose)
+void index_pattern(struct image *image, UnitCell *cell, IndexingMethod indm,
+                   int no_match, int verbose)
 {
 	int i;
 	int nc = 0;
@@ -115,8 +115,7 @@ void index_pattern(struct image *image, IndexingMethod indm, int no_match,
 			STATUS("--------------------\n");
 		}
 
-		new_cell = match_cell(image->candidate_cells[i],
-			              image->molecule->cell, verbose);
+		new_cell = match_cell(image->candidate_cells[i], cell, verbose);
 		image->indexed_cell = new_cell;
 		if ( new_cell != NULL ) {
 			STATUS("Matched on attempt %i.\n", i);
diff --git a/src/index.h b/src/index.h
index 26ce1c44..b30db837 100644
--- a/src/index.h
+++ b/src/index.h
@@ -18,6 +18,9 @@
 #endif
 
 
+#include "cell.h"
+
+
 typedef enum {
 	INDEXING_NONE,
 	INDEXING_DIRAX,
@@ -25,7 +28,7 @@ typedef enum {
 } IndexingMethod;
 
 
-extern void index_pattern(struct image *image, IndexingMethod indm,
-                          int no_match, int verbose);
+extern void index_pattern(struct image *image, UnitCell *cell,
+                          IndexingMethod indm, int no_match, int verbose);
 
 #endif	/* INDEX_H */
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 5a47254f..e0871fe7 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -32,6 +32,7 @@
 #include "detector.h"
 #include "sfac.h"
 #include "filters.h"
+#include "reflections.h"
 
 
 static void show_help(const char *s)
@@ -77,6 +78,8 @@ static void show_help(const char *s)
 "      --no-match          Don't attempt to match the indexed cell to the\n"
 "                           model, just proceed with the one generated by the\n"
 "                           auto-indexing procedure.\n"
+"     --intensities=<file>  Specify file containing reflection intensities\n"
+"                            to use.\n"
 );
 }
 
@@ -138,30 +141,22 @@ static struct image *get_simage(struct image *template, int alternate)
 
 	image->lambda = ph_en_to_lambda(eV_to_J(1.8e3));
 
-	image->molecule = copy_molecule(template->molecule);
-	free(image->molecule->cell);
-	image->molecule->cell = cell_new_from_cell(template->indexed_cell);
-
 	return image;
 }
 
 
-static void simulate_and_write(struct image *simage,
-                               struct gpu_context **gctx)
+static void simulate_and_write(struct image *simage, struct gpu_context **gctx,
+                               double *intensities, UnitCell *cell)
 {
 	/* Set up GPU if necessary */
 	if ( (gctx != NULL) && (*gctx == NULL) ) {
-		*gctx = setup_gpu(0, simage, simage->molecule);
+		*gctx = setup_gpu(0, simage, intensities);
 	}
 
 	if ( (gctx != NULL) && (*gctx != NULL) ) {
-		get_diffraction_gpu(*gctx, simage, 24, 24, 40);
+		get_diffraction_gpu(*gctx, simage, 24, 24, 40, cell);
 	} else {
-		get_diffraction(simage, 24, 24, 40, 0, 0);
-	}
-	if ( simage->molecule == NULL ) {
-		ERROR("Couldn't open molecule.pdb\n");
-		return;
+		get_diffraction(simage, 24, 24, 40, intensities, cell, 0);
 	}
 	record_image(simage, 0);
 
@@ -193,6 +188,9 @@ int main(int argc, char *argv[])
 	IndexingMethod indm;
 	char *indm_str = NULL;
 	struct image image;
+	UnitCell *cell;
+	double *intensities = NULL;
+	char *intfile = NULL;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -210,6 +208,7 @@ int main(int argc, char *argv[])
 		{"no-match",           0, &config_nomatch,     1},
 		{"verbose",            0, &config_verbose,     1},
 		{"alternate",          0, &config_alternate,   1},
+		{"intensities",        0, NULL,               'q'},
 		{0, 0, NULL, 0}
 	};
 
@@ -232,6 +231,11 @@ int main(int argc, char *argv[])
 			break;
 		}
 
+		case 'q' : {
+			intfile = strdup(optarg);
+			break;
+		}
+
 		case 0 : {
 			break;
 		}
@@ -257,6 +261,12 @@ int main(int argc, char *argv[])
 		return 1;
 	}
 
+	if ( intfile != NULL ) {
+		intensities = read_reflections(intfile);
+	} else {
+		intensities = NULL;
+	}
+
 	if ( indm_str == NULL ) {
 		STATUS("You didn't specify an indexing method, so I won't"
 		       " try to index anything.\n"
@@ -273,10 +283,9 @@ int main(int argc, char *argv[])
 	}
 	free(indm_str);
 
-	image.molecule = load_molecule();
-	if ( image.molecule == NULL ) {
-		ERROR("You must provide molecule.pdb in the working "
-		      "directory.\n");
+	cell = load_cell_from_pdb("molecule.pdb");
+	if ( cell == NULL ) {
+		ERROR("Couldn't read unit cell (from molecule.pdb)\n");
 		return 1;
 	}
 
@@ -354,7 +363,7 @@ int main(int argc, char *argv[])
 
 		/* Calculate orientation matrix (by magic) */
 		if ( config_writedrx || (indm != INDEXING_NONE) ) {
-			index_pattern(&image, indm, config_nomatch,
+			index_pattern(&image, cell, indm, config_nomatch,
 			              config_verbose);
 		}
 
@@ -376,16 +385,17 @@ int main(int argc, char *argv[])
 		/* Simulate if requested */
 		if ( config_simulate ) {
 			if ( config_gpu ) {
-				simulate_and_write(simage, &gctx);
+				simulate_and_write(simage, &gctx, intensities,
+				                   cell);
 			} else {
-				simulate_and_write(simage, NULL);
+				simulate_and_write(simage, NULL, intensities,
+				                   cell);
 			}
 		}
 
 		/* Finished with alternate image */
 		if ( simage->twotheta != NULL ) free(simage->twotheta);
 		if ( simage->data != NULL ) free(simage->data);
-		free_molecule(simage->molecule);
 		free(simage);
 
 		/* Only free cell if found */
@@ -402,7 +412,6 @@ done:
 	} while ( rval != NULL );
 
 	fclose(fh);
-	free_molecule(image.molecule);
 
 	STATUS("There were %i images.\n", n_images);
 	STATUS("%i hits were found.\n", n_hits);
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 4fd7da66..2ff12226 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -162,6 +162,7 @@ int main(int argc, char *argv[])
 	int number = 1;   /* Number used for filename of image */
 	int n_images = 1; /* Generate one image by default */
 	int done = 0;
+	UnitCell *cell;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -240,7 +241,7 @@ int main(int argc, char *argv[])
 	image.width = 1024;
 	image.height = 1024;
 	image.lambda = ph_en_to_lambda(eV_to_J(1790.0));  /* Wavelength */
-	image.molecule = load_molecule();
+	cell = load_cell_from_pdb("molecule.pdb");
 
 	#include "geometry-lcls.tmp"
 
@@ -281,7 +282,7 @@ int main(int argc, char *argv[])
 		image.data = NULL;
 		image.twotheta = NULL;
 
-		cell_get_parameters(image.molecule->cell, &a, &b, &c, &d, &d, &d);
+		cell_get_parameters(cell, &a, &b, &c, &d, &d, &d);
 		STATUS("Particle size = %i x %i x %i (=%5.2f x %5.2f x %5.2f nm)\n",
 	               na, nb, nc, na*a/1.0e-9, nb*b/1.0e-9, nc*c/1.0e-9);
 
@@ -290,15 +291,11 @@ int main(int argc, char *argv[])
 				gctx = setup_gpu(config_nosfac, &image,
 				                 intensities);
 			}
-			get_diffraction_gpu(gctx, &image, na, nb, nc);
+			get_diffraction_gpu(gctx, &image, na, nb, nc, cell);
 		} else {
-			get_diffraction(&image, na, nb, nc, intensities,
+			get_diffraction(&image, na, nb, nc, intensities, cell,
 			                !config_nowater);
 		}
-		if ( image.molecule == NULL ) {
-			ERROR("Couldn't open molecule.pdb\n");
-			return 1;
-		}
 		if ( image.data == NULL ) {
 			ERROR("Diffraction calculation failed.\n");
 			goto skip;
@@ -307,7 +304,7 @@ int main(int argc, char *argv[])
 		record_image(&image, !config_nonoise);
 
 		if ( config_nearbragg ) {
-			output_intensities(&image, image.molecule->cell);
+			output_intensities(&image, cell);
 		}
 
 		if ( config_powder ) {
@@ -359,8 +356,8 @@ skip:
 
 	free(image.det.panels);
 	free(powder);
-	free(image.molecule->reflections);
-	free(image.molecule);
+	free(cell);
+	free(intensities);
 
 	return 0;
 }