Improved framework for detector geometry

7 files changed, 108 insertions, 30 deletions

diff --git a/src/detector.h b/src/detector.h
index 09e4ce3b..30c2f44c 100644
--- a/src/detector.h
+++ b/src/detector.h
@@ -16,13 +16,25 @@
 #ifndef DETECTOR_H
 #define DETECTOR_H
 
+struct image;
+
 #include "image.h"
 
-/* Position of central beam for upper and lower CCDs */
-#define UPPER_CX (491.9)
-#define UPPER_CY (440.7)
-#define LOWER_CX (492.0)
-#define LOWER_CY (779.7)
+struct panel
+{
+	int      min_x;   /* Smallest x value considered to be in this panel */
+	int      max_x;   /* Largest x value considered to be in this panel */
+	int      min_y;   /* ... and so on */
+	int      max_y;
+	float    cx;     /* Location of centre */
+	float    cy;
+};
+
+struct detector
+{
+	struct panel *panels;
+	int           n_panels;
+};
 
 extern void record_image(struct image *image, int do_water, int do_poisson,
                          int do_bloom);
diff --git a/src/ewald.c b/src/ewald.c
index 19ce301e..94fb4ed6 100644
--- a/src/ewald.c
+++ b/src/ewald.c
@@ -68,20 +68,25 @@ void get_ewald(struct image *image)
 	for ( x=0; x<image->width; x++ ) {
 	for ( y=0; y<image->height; y++ ) {
 
-		double rx, ry, r;
+		double rx = 0.0;
+		double ry = 0.0;
+		double r;
 		double twothetax, twothetay, twotheta;
 		double qx, qy, qz;
 		struct rvec q;
+		int p;
 
 		/* Calculate q vectors for Ewald sphere */
-		if ( y >= 512 ) {
-			/* Top CCD */
-			rx = ((double)x - UPPER_CX) / image->resolution;
-			ry = ((double)y - UPPER_CY) / image->resolution;
-		} else {
-			/* Bottom CCD */
-			rx = ((double)x - LOWER_CX) / image->resolution;
-			ry = ((double)y - LOWER_CY) / image->resolution;
+		for ( p=0; p<image->det.n_panels; p++ ) {
+			if ( (x >= image->det.panels[p].min_x)
+			  && (x <= image->det.panels[p].max_x)
+			  && (y >= image->det.panels[p].min_y)
+			  && (y <= image->det.panels[p].max_y) ) {
+				rx = ((double)x - image->det.panels[p].cx)
+				                            / image->resolution;
+				ry = ((double)y - image->det.panels[p].cy)
+				                            / image->resolution;
+			}
 		}
 		r = sqrt(pow(rx, 2.0) + pow(ry, 2.0));
 
diff --git a/src/image.h b/src/image.h
index 59435b4e..166b20df 100644
--- a/src/image.h
+++ b/src/image.h
@@ -23,6 +23,7 @@
 
 #include "utils.h"
 #include "cell.h"
+#include "detector.h"
 
 
 /* How is the scaling of the image described? */
@@ -79,6 +80,7 @@ struct image {
 	double			*twotheta;
 	struct molecule		*molecule;
 	UnitCell		*indexed_cell;
+	struct detector		det;
 
 	struct quaternion	orientation;
 
diff --git a/src/index.c b/src/index.c
index 2e0fd027..1d737e7e 100644
--- a/src/index.c
+++ b/src/index.c
@@ -111,19 +111,30 @@ void index_pattern(struct image *image, IndexingMethod indm)
 	for ( i=0; i<image_feature_count(image->features); i++ ) {
 
 		struct imagefeature *f;
+		double rx = 0.0;
+		double ry = 0.0;
+		int p;
+		int found = 0;
 
 		f = image_get_feature(image->features, i);
 		if ( f == NULL ) continue;
 
-		if ( f->y >=512 ) {
-			/* Top half of CCD */
-			map_position(image, f->x-UPPER_CX, f->y-UPPER_CY,
-			             &f->rx, &f->ry, &f->rz);
-		} else {
-			/* Lower half of CCD */
-			map_position(image, f->x-LOWER_CX, f->y-LOWER_CY,
-			             &f->rx, &f->ry, &f->rz);
+		for ( p=0; p<image->det.n_panels; p++ ) {
+			if ( (f->x >= image->det.panels[p].min_x)
+			  && (f->x <= image->det.panels[p].max_x)
+			  && (f->y >= image->det.panels[p].min_y)
+			  && (f->y <= image->det.panels[p].max_y) ) {
+				rx = ((double)f->x - image->det.panels[p].cx);
+				ry = ((double)f->y - image->det.panels[p].cy);
+				found = 1;
+			}
 		}
+		if ( !found ) {
+			ERROR("No mapping found for %f,%f\n", f->x, f->y);
+			continue;
+		}
+
+		map_position(image, rx, ry, &f->rx, &f->ry, &f->rz);
 	}
 
 	write_drx(image);
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 67073b89..6b0ad3d9 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -166,6 +166,24 @@ int main(int argc, char *argv[])
 		image.data = NULL;
 		image.indexed_cell = NULL;
 
+		/* Set up detector configuration */
+		image.det.n_panels = 2;
+		image.det.panels = malloc(2*sizeof(struct panel));
+		/* Upper panel */
+		image.det.panels[0].min_x = 0;
+		image.det.panels[0].max_x = 1023;
+		image.det.panels[0].min_y = 512;
+		image.det.panels[0].max_y = 1023;
+		image.det.panels[0].cx = 491.9;
+		image.det.panels[0].cy = 440.7;
+		/* Lower panel */
+		image.det.panels[1].min_x = 0;
+		image.det.panels[1].max_x = 1023;
+		image.det.panels[1].min_y = 0;
+		image.det.panels[1].max_y = 511;
+		image.det.panels[1].cx = 492.0;
+		image.det.panels[1].cy = 779.7;
+
 		STATUS("Processing '%s'\n", line);
 
 		n_images++;
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 7a27ceb6..d3056643 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -220,6 +220,24 @@ int main(int argc, char *argv[])
 	image.lambda = ph_en_to_lambda(eV_to_J(2.0e3));  /* Wavelength */
 	image.molecule = load_molecule();
 
+	/* Set up detector configuration */
+	image.det.n_panels = 2;
+	image.det.panels = malloc(2*sizeof(struct panel));
+	/* Upper panel */
+	image.det.panels[0].min_x = 0;
+	image.det.panels[0].max_x = 1023;
+	image.det.panels[0].min_y = 512;
+	image.det.panels[0].max_y = 1023;
+	image.det.panels[0].cx = 491.9;
+	image.det.panels[0].cy = 440.7;
+	/* Lower panel */
+	image.det.panels[1].min_x = 0;
+	image.det.panels[1].max_x = 1023;
+	image.det.panels[1].min_y = 0;
+	image.det.panels[1].max_y = 511;
+	image.det.panels[1].cx = 492.0;
+	image.det.panels[1].cy = 779.7;
+
 	/* Splurge a few useful numbers */
 	STATUS("Wavelength is %f nm\n", image.lambda/1.0e-9);
 
diff --git a/src/peaks.c b/src/peaks.c
index 31de2a79..56fa6931 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -432,7 +432,11 @@ void dump_peaks(struct image *image)
 
 		double q, rx, ry, rz;
 		int x, y;
+		double rcx = 0.0;
+		double rcy = 0.0;
+		int found = 0;
 		struct imagefeature *f;
+		int p;
 
 		f = image_get_feature(image->features, i);
 		if ( f == NULL ) continue;
@@ -440,16 +444,24 @@ void dump_peaks(struct image *image)
 		x = f->x;
 		y = f->y;
 
-		if ( f->y >=512 ) {
-			/* Top half of CCD */
-			map_position(image, f->x-UPPER_CX, f->y-UPPER_CY,
-			             &rx, &ry, &rz);
-		} else {
-			/* Lower half of CCD */
-			map_position(image, f->x-LOWER_CX, f->y-LOWER_CY,
-			             &rx, &ry, &rz);
+		for ( p=0; p<image->det.n_panels; p++ ) {
+			if ( (x >= image->det.panels[p].min_x)
+			  && (x <= image->det.panels[p].max_x)
+			  && (y >= image->det.panels[p].min_y)
+			  && (y <= image->det.panels[p].max_y) ) {
+				rcx = ((double)x - image->det.panels[p].cx)
+				                            / image->resolution;
+				rcy = ((double)y - image->det.panels[p].cy)
+				                            / image->resolution;
+				found = 1;
+			}
+		}
+		if ( !found ) {
+			ERROR("No mapping found for %i,%i\n", x, y);
+			continue;
 		}
 
+		map_position(image, rcx, rcy, &rx, &ry, &rz);
 		q = modulus(rx, ry, rz);
 
 		printf("%i\t%i\t%f\t%i\n", x, y, q/1.0e9,