Convert pattern_sim to DataTemplate

3 files changed, 91 insertions, 174 deletions

diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c
index 419b5957..74f03ebf 100644
--- a/src/diffraction-gpu.c
+++ b/src/diffraction-gpu.c
@@ -38,6 +38,7 @@
 #include <string.h>
 #include <complex.h>
 
+#define CL_TARGET_OPENCL_VERSION 220
 #ifdef HAVE_CL_CL_H
 #include <CL/cl.h>
 #else
@@ -187,18 +188,18 @@ static int do_panels(struct gpu_context *gctx, struct image *image,
 	if ( set_arg_float(gctx, 2, weight) ) return 1;
 
 	/* Iterate over panels */
-	for ( i=0; i<image->det->n_panels; i++ ) {
+	for ( i=0; i<image->detgeom->n_panels; i++ ) {
 
 		size_t dims[2];
 		size_t ldims[2];
-		struct panel *p;
+		struct detgeom_panel *p;
 		cl_mem diff;
 		size_t diff_size;
 		float *diff_ptr;
 		int fs, ss;
 		cl_int err;
 
-		p = &image->det->panels[i];
+		p = &image->detgeom->panels[i];
 
 		/* Buffer for the results of this panel */
 		diff_size = p->w * p->h * sizeof(cl_float);
@@ -220,8 +221,8 @@ static int do_panels(struct gpu_context *gctx, struct image *image,
 		if ( set_arg_float(gctx, 9, p->ssx) ) return 1;
 		if ( set_arg_float(gctx, 10, p->ssy) ) return 1;
 		if ( set_arg_float(gctx, 11, p->ssz) ) return 1;
-		if ( set_arg_float(gctx, 12, p->res) ) return 1;
-		if ( set_arg_float(gctx, 13, p->clen) ) return 1;
+		if ( set_arg_float(gctx, 12, 1.0/p->pixel_pitch) ) return 1;
+		if ( set_arg_float(gctx, 13, p->cnz*p->pixel_pitch) ) return 1;
 
 		dims[0] = p->w * sampling;
 		dims[1] = p->h * sampling;
@@ -325,13 +326,13 @@ int get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
 	if ( set_arg_mem(gctx, 19, gctx->sinc_luts[nc-1]) ) return 1;
 
 	/* Allocate memory for the result */
-	image->dp = malloc(image->det->n_panels * sizeof(float *));
+	image->dp = malloc(image->detgeom->n_panels * sizeof(float *));
 	if ( image->dp == NULL ) {
 		ERROR("Couldn't allocate memory for result.\n");
 		return 1;
 	}
-	for ( i=0; i<image->det->n_panels; i++ ) {
-		struct panel *p = &image->det->panels[i];
+	for ( i=0; i<image->detgeom->n_panels; i++ ) {
+		struct detgeom_panel *p = &image->detgeom->panels[i];
 		image->dp[i] = calloc(p->w * p->h, sizeof(float));
 		if ( image->dp[i] == NULL ) {
 			ERROR("Couldn't allocate memory for panel %i\n", i);
diff --git a/src/diffraction.c b/src/diffraction.c
index 8762b0af..4f17a27d 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -376,7 +376,7 @@ static void diffraction_panel(struct image *image, const double *intensities,
 	int fs, ss;
 	const int nxs = 4;
 	const int nys = 4;
-	struct panel *p = &image->det->panels[pn];
+	struct detgeom_panel *p = &image->detgeom->panels[pn];
 
 	weight /= nxs*nys;
 
@@ -386,17 +386,22 @@ static void diffraction_panel(struct image *image, const double *intensities,
 		int idx;
 		double f_lattice, I_lattice;
 		double I_molecule;
-		struct rvec q;
 		int xs, ys;
 		float xo, yo;
 
 		for ( xs=0; xs<nxs; xs++ ) {
 		for ( ys=0; ys<nys; ys++ ) {
 
+			double qv[3];
+			struct rvec q;
+
 			xo = (1.0/nxs) * xs;
 			yo = (1.0/nys) * ys;
 
-			q = get_q_for_panel(p, fs+xo, ss+yo, NULL, k);
+			detgeom_transform_coords(p, fs+xo, ss+yo,
+			                         image->lambda, qv);
+
+			q.u = qv[0]; q.v = qv[1]; q.w = qv[2];
 
 			f_lattice = lattice_factor(q, ax, ay, az,
 					           bx, by, bz,
@@ -435,7 +440,7 @@ static void diffraction_at_k(struct image *image, const double *intensities,
 {
 	int i;
 
-	for ( i=0; i<image->det->n_panels; i++ ) {
+	for ( i=0; i<image->detgeom->n_panels; i++ ) {
 		diffraction_panel(image, intensities, phases, flags, cell, m,
 		                  sym, k, ax, ay, az, bx, by, bz, cx, cy, cz,
 		                  lut_a, lut_b, lut_c, i, weight);
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 916257e3..4d74589b 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -46,12 +46,12 @@
 #include <cell.h>
 #include <cell-utils.h>
 #include <utils.h>
-#include <detector.h>
 #include <peaks.h>
 #include <symmetry.h>
 #include <reflist.h>
 #include <reflist-utils.h>
 #include <stream.h>
+#include <datatemplate.h>
 
 #include "diffraction.h"
 #include "diffraction-gpu.h"
@@ -117,7 +117,8 @@ static void show_help(const char *s)
 }
 
 
-static void record_panel(struct panel *p, float *dp, int do_poisson,
+static void record_panel(struct detgeom_panel *p, float *dp,
+                         int do_poisson,
                          gsl_rng *rng, double ph_per_e, double background,
 			 double lambda,
                          int *n_neg1, int *n_inf1, int *n_nan1,
@@ -144,16 +145,16 @@ static void record_panel(struct panel *p, float *dp, int do_poisson,
 		if ( isnan(intensity) ) (*n_nan1)++;
 
 		/* Area of one pixel */
-		pix_area = pow(1.0/p->res, 2.0);
-		Lsq = pow(p->clen, 2.0);
+		pix_area = pow(p->pixel_pitch, 2.0);
+		Lsq = pow(p->cnz*p->pixel_pitch, 2.0);
 
 		/* Calculate distance from crystal to pixel */
 		xs = fs*p->fsx + ss*p->ssx;
 		ys = ss*p->fsy + ss*p->ssy;
-		rx = (xs + p->cnx) / p->res;
-		ry = (ys + p->cny) / p->res;
+		rx = (xs + p->cnx) * p->pixel_pitch;
+		ry = (ys + p->cny) * p->pixel_pitch;
 		dsq = pow(rx, 2.0) + pow(ry, 2.0);
-		twotheta = atan2(sqrt(dsq), p->clen);
+		twotheta = atan2(sqrt(dsq), p->cnz*p->pixel_pitch);
 
 		/* Area of pixel as seen from crystal (approximate) */
 		proj_area = pix_area * cos(twotheta);
@@ -216,11 +217,9 @@ void record_image(struct image *image, int do_poisson, double background,
 	       "Total energy = %5.3f microJ\n",
 	       nphotons, energy_density/1e7, total_energy*1e6);
 
-	fill_in_adu(image);
+	for ( pn=0; pn<image->detgeom->n_panels; pn++ ) {
 
-	for ( pn=0; pn<image->det->n_panels; pn++ ) {
-
-		record_panel(&image->det->panels[pn], image->dp[pn],
+		record_panel(&image->detgeom->panels[pn], image->dp[pn],
 		             do_poisson, rng, ph_per_e, background,
 			     image->lambda,
 		             &n_neg1, &n_inf1, &n_nan1,
@@ -364,28 +363,8 @@ static int random_ncells(double len, double min_m, double max_m)
 }
 
 
-static void fixup_geom(struct detector *det)
-{
-	int i;
-
-	for ( i=0; i<det->n_panels; i++ ) {
-		det->panels[i].clen += det->panels[i].coffset;
-	}
-}
-
-
-static int geom_contains_references(struct detector *det)
-{
-	int i;
-
-	for ( i=0; i<det->n_panels; i++ ) {
-		if ( det->panels[i].clen_from != NULL ) return 1;
-	}
-	return 0;
-}
-
-
-static void add_metadata(const char *filename, struct quaternion q,
+static void add_metadata(const char *filename,
+                         struct quaternion q,
                          UnitCell *cell)
 {
 	hid_t fh, gh, sh, dh;
@@ -488,9 +467,10 @@ static void add_metadata(const char *filename, struct quaternion q,
 int main(int argc, char *argv[])
 {
 	int c;
-	struct image image;
+	struct image *image;
+	DataTemplate *dtempl;
 	struct gpu_context *gctx = NULL;
-	struct image powder;
+	struct image *powder;
 	char *intfile = NULL;
 	double *intensities;
 	char *rval;
@@ -536,8 +516,6 @@ int main(int argc, char *argv[])
 	double photon_energy = 9000.0;
 	double sase_spike_width = 2e6;  /* m^-1 */
 	double twocol_sep = 8e6;  /* m^-1 */
-	struct beam_params beam;
-	int i;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -802,7 +780,7 @@ int main(int argc, char *argv[])
 			ERROR("You cannot use -r and --template together.\n");
 			return 1;
 		}
-		st = open_stream_for_read(template_file);
+		st = stream_open_for_read(template_file);
 		if ( st == NULL ) {
 			ERROR("Failed to open stream.\n");
 			return 1;
@@ -836,26 +814,13 @@ int main(int argc, char *argv[])
 		ERROR("You need to specify a geometry file with --geometry\n");
 		return 1;
 	}
-	image.beam = &beam;
-	image.det = get_detector_geometry(geometry, image.beam);
-	if ( image.det == NULL ) {
-		ERROR("Failed to read detector geometry from '%s'\n", geometry);
+	dtempl = data_template_new_from_file(geometry);
+	if ( dtempl == NULL ) {
+		ERROR("Failed to read detector geometry from '%s'\n",
+		      geometry);
 		return 1;
 	}
 	free(geometry);
-	if ( (beam.photon_energy > 0.0) && (beam.photon_energy_from == NULL) ) {
-		ERROR("WARNING: An explicit photon energy was found in the "
-		      "geometry file.  It will be ignored!\n");
-		ERROR("The value given on the command line "
-		      "(with --photon-energy) will be used instead.\n");
-	}
-	if ( geom_contains_references(image.det) ) {
-		ERROR("Geometry file contains a reference to an HDF5 location"
-		      " for the camera length.  Change it to a numerical value "
-		      " and try again.\n");
-		return 1;
-	}
-	fixup_geom(image.det);
 
 	if ( spectrum_str == NULL ) {
 		STATUS("You didn't specify a spectrum type, so"
@@ -952,28 +917,6 @@ int main(int argc, char *argv[])
 
 	}
 
-	image.lambda = ph_en_to_lambda(eV_to_J(photon_energy));
-	image.bw = bandwidth;
-
-	/* Initialise stuff */
-	image.filename = NULL;
-	image.features = NULL;
-	image.bad = NULL;
-	image.dp = malloc(image.det->n_panels*sizeof(float *));
-	if ( image.dp == NULL ) {
-		ERROR("Failed to allocate data\n");
-		return 1;
-	}
-	for ( i=0; i<image.det->n_panels; i++ ) {
-		image.dp[i] = calloc(image.det->panels[i].w *
-		                     image.det->panels[i].h,
-		                     sizeof(float));
-		if ( image.dp[i] == NULL ) {
-			ERROR("Failed to allocate data\n");
-			return 1;
-		}
-	}
-
 	rng = gsl_rng_alloc(gsl_rng_mt19937);
 	if ( config_random ) {
 		FILE *fh;
@@ -987,30 +930,16 @@ int main(int argc, char *argv[])
 		fclose(fh);
 	}
 
-	powder.det = image.det;
-	powder.beam = NULL;
-	powder.lambda = 0.0;
-	powder.spectrum = NULL;
-	powder.dp = malloc(image.det->n_panels*sizeof(float *));
-	if ( powder.dp == NULL ) {
-		ERROR("Failed to allocate powder data\n");
-		return 1;
-	}
-	for ( i=0; i<image.det->n_panels; i++ ) {
-		powder.dp[i] = calloc(image.det->panels[i].w *
-		                      image.det->panels[i].h,
-		                      sizeof(float));
-		if ( powder.dp[i] == NULL ) {
-			ERROR("Failed to allocate powder data\n");
-			return 1;
-		}
-	}
+	image = image_create_for_simulation(dtempl);
+	powder = image_create_for_simulation(dtempl);
 
 	/* Splurge a few useful numbers */
 	STATUS("Simulation parameters:\n");
-	STATUS("                  Photon energy: %.2f eV (wavelength %.5f A)\n",
-	       photon_energy, image.lambda*1e10);
-	STATUS("              Number of photons: %.0f (%.2f mJ)\n", nphotons,
+	STATUS("                     Wavelength: %.5f A (photon energy %.2f eV)\n",
+	       image->lambda*1e10,
+	       ph_lambda_to_eV(image->lambda));
+	STATUS("              Number of photons: %.0f (%.2f mJ)\n",
+	       nphotons,
 	       eV_to_J(photon_energy)*nphotons*1e3);
 	STATUS("                Beam divergence: not simulated\n");
 	STATUS("                    Beam radius: %.2f microns\n",
@@ -1021,22 +950,22 @@ int main(int argc, char *argv[])
 
 		case SPECTRUM_TOPHAT:
 		STATUS("                 X-ray spectrum: top hat, "
-		       "width %.5f %%\n", image.bw*100.0);
+		       "width %.5f %%\n", image->bw*100.0);
 		break;
 
 		case SPECTRUM_GAUSSIAN:
 		STATUS("                 X-ray spectrum: Gaussian, "
-		       "bandwidth %.5f %%\n", image.bw*100.0);
+		       "bandwidth %.5f %%\n", image->bw*100.0);
 		break;
 
 		case SPECTRUM_SASE:
 		STATUS("                 X-ray spectrum: SASE, "
-		       "bandwidth %.5f %%\n", image.bw*100.0);
+		       "bandwidth %.5f %%\n", image->bw*100.0);
 		break;
 
 		case SPECTRUM_TWOCOLOUR:
 		STATUS("                 X-ray spectrum: two colour, "
-		       "separation %.5f %%\n", image.bw*100.0);
+		       "separation %.5f %%\n", image->bw*100.0);
 		break;
 
 		case SPECTRUM_FROMFILE:
@@ -1066,12 +995,13 @@ int main(int argc, char *argv[])
 		int err = 0;
 		int pi;
 
-		for ( pi=0; pi<image.det->n_panels; pi++ ) {
+		/* Reset image data to zero for new pattern */
+		for ( pi=0; pi<image->detgeom->n_panels; pi++ ) {
 			long j;
-			long np = image.det->panels[pi].w
-			            * image.det->panels[pi].h;
+			long np = image->detgeom->panels[pi].w
+			            * image->detgeom->panels[pi].h;
 			for ( j=0; j<np; j++ ) {
-				image.dp[pi][j] = 0.0;
+				image->dp[pi][j] = 0.0;
 			}
 		}
 
@@ -1115,70 +1045,54 @@ int main(int argc, char *argv[])
 
 		} else {
 
-			struct image image;
-			int i;
+			struct image *templ_image;
 			Crystal *cr;
 
-			image.det = NULL;
-
 			/* Get data from next chunk */
-			if ( read_chunk(st, &image) ) break;
+			templ_image = stream_read_chunk(st, dtempl,
+			                                STREAM_CRYSTALS);
+			if ( templ_image == NULL ) break;
+			if ( templ_image->n_crystals == 0 ) continue;
 
-			free(image.filename);
-			image_feature_list_free(image.features);
+			cr = templ_image->crystals[0];
+			cell = cell_new_from_cell(crystal_get_cell(cr));
 
-			if ( image.n_crystals == 0 ) continue;
-
-			cr = image.crystals[0];
-			cell = crystal_get_cell(cr);
-
-			if ( image.n_crystals > 1 ) {
+			if ( templ_image->n_crystals > 1 ) {
 				ERROR("Using the first crystal only.\n");
 			}
 
-			for ( i=1; i<image.n_crystals; i++ ) {
-
-				Crystal *cr = image.crystals[i];
-				cell = crystal_get_cell(cr);
-
-				reflist_free(crystal_get_reflections(cr));
-				cell_free(crystal_get_cell(cr));
-				crystal_free(cr);
-
-			}
-
-			free(image.crystals);
+			image_free(templ_image);
 
 		}
 
 		switch ( spectrum_type ) {
 
 			case SPECTRUM_TOPHAT :
-			image.spectrum = spectrum_generate_tophat(image.lambda,
-			                                          image.bw);
+			image->spectrum = spectrum_generate_tophat(image->lambda,
+			                                           image->bw);
 			break;
 
 			case SPECTRUM_GAUSSIAN :
-			image.spectrum = spectrum_generate_gaussian(image.lambda,
-			                                            image.bw);
+			image->spectrum = spectrum_generate_gaussian(image->lambda,
+			                                             image->bw);
 			break;
 
 
 			case SPECTRUM_SASE :
-			image.spectrum = spectrum_generate_sase(image.lambda,
-			                                        image.bw,
-			                                        sase_spike_width,
-			                                        rng);
+			image->spectrum = spectrum_generate_sase(image->lambda,
+			                                         image->bw,
+			                                         sase_spike_width,
+			                                         rng);
 			break;
 
-                        case SPECTRUM_TWOCOLOUR :
-                        image.spectrum = spectrum_generate_twocolour(image.lambda,
-                                                                     image.bw,
-                                                                     twocol_sep);
+			case SPECTRUM_TWOCOLOUR :
+			image->spectrum = spectrum_generate_twocolour(image->lambda,
+			                                              image->bw,
+			                                              twocol_sep);
                         break;
 
-                        case SPECTRUM_FROMFILE :
-                        image.spectrum = spectrum_load(spectrum_fn);
+			case SPECTRUM_FROMFILE :
+			image->spectrum = spectrum_load(spectrum_fn);
                         break;
 
 		}
@@ -1195,12 +1109,12 @@ int main(int argc, char *argv[])
 				                 intensities, flags, sym_str,
 				                 gpu_dev);
 			}
-			err = get_diffraction_gpu(gctx, &image, na, nb, nc,
+			err = get_diffraction_gpu(gctx, image, na, nb, nc,
 			                          cell, no_fringes, flat,
 			                          nsamples);
 
 		} else {
-			get_diffraction(&image, na, nb, nc, intensities, phases,
+			get_diffraction(image, na, nb, nc, intensities, phases,
 			                flags, cell, grad, sym, no_fringes, flat,
 			                nsamples);
 		}
@@ -1210,28 +1124,28 @@ int main(int argc, char *argv[])
 			goto skip;
 		}
 
-		record_image(&image, !config_nonoise, background, rng,
+		record_image(image, !config_nonoise, background, rng,
 		             beam_radius, nphotons);
 
 		if ( powder_fn != NULL ) {
 
 			int pn;
 
-			for ( pn=0; pn<image.det->n_panels; pn++ ) {
+			for ( pn=0; pn<image->detgeom->n_panels; pn++ ) {
 
 				size_t w, i;
 
-				w = image.det->panels[pn].w
-				  * image.det->panels[pn].h;
+				w = image->detgeom->panels[pn].w
+				  * image->detgeom->panels[pn].h;
 
 				for ( i=0; i<w; i++ ) {
-					powder.dp[pn][i] += (double)image.dp[pn][i];
+					powder->dp[pn][i] += (double)image->dp[pn][i];
 				}
 
 			}
 
 			if ( !(ndone % 10) ) {
-				hdf5_write_image(powder_fn, &powder, NULL);
+				image_write(powder, dtempl, powder_fn);
 			}
 		}
 
@@ -1249,7 +1163,9 @@ int main(int argc, char *argv[])
 			number++;
 
 			/* Write the output file */
-			hdf5_write_image(filename, &image, NULL);
+			image_write(image, dtempl, filename);
+
+			/* Add some pattern_sim-specific metadata */
 			add_metadata(filename, orientation, cell);
 
 		}
@@ -1264,21 +1180,16 @@ skip:
 	} while ( !done );
 
 	if ( powder_fn != NULL ) {
-		hdf5_write_image(powder_fn, &powder, NULL);
+		image_write(powder, dtempl, powder_fn);
 	}
 
 	if ( gctx != NULL ) {
 		cleanup_gpu(gctx);
 	}
 
-	for ( i=0; i<image.det->n_panels; i++ ) {
-		free(image.dp[i]);
-		free(powder.dp[i]);
-	}
-	free(image.dp);
-	free(powder.dp);
+	image_free(image);
+	image_free(powder);
 	free(intfile);
-	free(image.det);
 	cell_free(input_cell);
 	free(intensities);
 	free(outfile);