Handle upper and lower CCDs correctly

5 files changed, 40 insertions, 12 deletions

diff --git a/src/detector.h b/src/detector.h
index dff41b43..e441b453 100644
--- a/src/detector.h
+++ b/src/detector.h
@@ -18,6 +18,12 @@
 
 #include "image.h"
 
+/* Position of central beam for upper and lower CCDs */
+#define UPPER_CX (492.8)
+#define UPPER_CY (437.6)
+#define LOWER_CX (494.0)
+#define LOWER_CY (772.1)
+
 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 82b70df2..b5b61f1c 100644
--- a/src/ewald.c
+++ b/src/ewald.c
@@ -18,6 +18,7 @@
 #include "utils.h"
 #include "cell.h"
 #include "ewald.h"
+#include "detector.h"
 
 
 static struct rvec quat_rot(struct rvec q, struct quaternion z)
@@ -73,8 +74,15 @@ void get_ewald(struct image *image)
 		struct rvec q;
 
 		/* Calculate q vectors for Ewald sphere */
-		rx = ((double)x - image->x_centre) / image->resolution;
-		ry = ((double)y - image->y_centre) / image->resolution;
+		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;
+		}
 		r = sqrt(pow(rx, 2.0) + pow(ry, 2.0));
 
 		twothetax = atan2(rx, image->camera_len);
diff --git a/src/hdf5-file.c b/src/hdf5-file.c
index 32bded7a..0ff98913 100644
--- a/src/hdf5-file.c
+++ b/src/hdf5-file.c
@@ -258,8 +258,8 @@ int hdf5_read(struct hdfile *f, struct image *image)
 	image->width = f->ny;
 
 	/* FIXME: The following are basically made up... */
-	image->x_centre = image->width/2 - 19;
-	image->y_centre = image->height/2 - 75;
+	image->x_centre = image->width/2;
+	image->y_centre = image->height/2;
 	image->lambda = ph_en_to_lambda(eV_to_J(1793.0));
 	image->fmode = FORMULATION_CLEN;
 	image->camera_len = 75.0e-3;  /* 75 mm camera length */
diff --git a/src/index.c b/src/index.c
index c8d5dbe6..b90a3c0b 100644
--- a/src/index.c
+++ b/src/index.c
@@ -24,8 +24,11 @@
 #include "utils.h"
 #include "peaks.h"
 #include "dirax.h"
+#include "sfac.h"
+#include "detector.h"
 
 
+/* x,y in pixels relative to central beam */
 static int map_position(struct image *image, double x, double y,
                         double *rx, double *ry, double *rz)
 {
@@ -35,13 +38,8 @@ static int map_position(struct image *image, double x, double y,
 	/* Angular description of reflection */
 	double theta, psi, k;
 
-	x -= image->x_centre;
-	y -= image->y_centre;
 	k = 1.0 / image->lambda;
 
-	/* FIXME: Don't process lower CCD for now */
-	if ( y < 0 ) return 0;
-
 	if ( image->fmode == FORMULATION_CLEN ) {
 
 		/* Convert pixels to metres */
@@ -81,14 +79,22 @@ void index_pattern(struct image *image, int no_index, int dump_peaks,
 	/* Perform 'fine' peak search */
 	search_peaks(image, dump_peaks);
 
-	/* Map positions to 3D. FIXME: Handle lower detector */
+	/* Map positions to 3D */
 	for ( i=0; i<image_feature_count(image->features); i++ ) {
 
 		struct imagefeature *f;
 
 		f = image_get_feature(image->features, i);
-		map_position(image, f->x, f->y, &f->rx, &f->ry, &f->rz);
 
+		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);
+		}
 	}
 
 	if ( use_dirax ) {
diff --git a/src/peaks.c b/src/peaks.c
index 984d4d3b..60849d0e 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -26,6 +26,14 @@
 
 #define PEAK_WINDOW_SIZE (10)
 
+static int in_streak(int x, int y)
+{
+	if ( (y>512) && (y<768) && (abs(x-493)<15) ) return 1;
+	if ( (y>768) && (abs(x-480)<15) ) return 1;
+	return 0;
+}
+
+
 struct peak {
 	int x;
 	int y;
@@ -199,7 +207,7 @@ void search_peaks(struct image *image, int dump_peaks)
 		if ( data[x+width*y] < 800 ) continue;
 
 		/* Ignore streak */
-		if ( abs(x-image->x_centre) < 15 ) continue;
+		if ( in_streak(x, y) ) continue;
 
 		/* Get gradients */
 		dx1 = data[x+width*y] - data[(x+1)+width*y];