WIP on per-panel detector geometry

9 files changed, 42 insertions, 34 deletions

diff --git a/src/detector.c b/src/detector.c
index 1bcddef1..94e11a9e 100644
--- a/src/detector.c
+++ b/src/detector.c
@@ -148,7 +148,6 @@ void record_image(struct image *image, int do_water, int do_poisson,
 	int x, y;
 	double total_energy, energy_density;
 	double ph_per_e;
-	double pix_area, Lsq;
 	double area;
 
 	/* How many photons are scattered per electron? */
@@ -161,10 +160,6 @@ void record_image(struct image *image, int do_water, int do_poisson,
 	       "Total energy = %5.3f microJ\n",
 	       FLUENCE, energy_density/1e7, total_energy*1e6);
 
-	/* Area of one pixel */
-	pix_area = pow(1.0/image->resolution, 2.0);
-	Lsq = pow(image->camera_len, 2.0);
-
 	image->hdr = malloc(image->width * image->height * sizeof(double));
 
 	for ( x=0; x<image->width; x++ ) {
@@ -174,11 +169,27 @@ void record_image(struct image *image, int do_water, int do_poisson,
 		double cf;
 		double intensity, sa, water;
 		double complex val;
+		double pix_area, Lsq;
 		double dsq, proj_area;
+		int p;
+		int found = 0;
 
 		val = image->sfacs[x + image->width*y];
 		intensity = pow(cabs(val), 2.0);
 
+		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) ) {
+				found = 1;
+			}
+		}
+		if ( !found ) {
+			ERROR("No mapping found for %i,%i\n", x, y);
+			return;
+		}
+
 		if ( do_water ) {
 
 			struct rvec q;
@@ -193,12 +204,18 @@ void record_image(struct image *image, int do_water, int do_poisson,
 
 		}
 
+		/* Area of one pixel */
+		pix_area = pow(1.0/image->det.panels[p].res, 2.0);
+		Lsq = pow(image->det.panels[p].clen, 2.0);
+
 		/* Area of pixel as seen from crystal (approximate) */
 		proj_area = pix_area * cos(image->twotheta[x + image->width*y]);
 
 		/* Calculate distance from crystal to pixel */
-		dsq = pow(((double)x - image->x_centre)/image->resolution, 2.0);
-		dsq += pow(((double)y - image->y_centre)/image->resolution, 2.0);
+		dsq = pow(((double)x - image->det.panels[p].cx)
+		                               / image->det.panels[p].res, 2.0);
+		dsq += pow(((double)y - image->det.panels[p].cy)
+		                               / image->det.panels[p].res, 2.0);
 
 		/* Projected area of pixel divided by distance squared */
 		sa = proj_area / (dsq + Lsq);
diff --git a/src/detector.h b/src/detector.h
index 30c2f44c..4a983313 100644
--- a/src/detector.h
+++ b/src/detector.h
@@ -26,8 +26,10 @@ struct 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    cx;      /* Location of centre */
 	float    cy;
+	float    clen;    /* Camera length */
+	float    res;     /* Resolution */
 };
 
 struct detector
diff --git a/src/diffraction.c b/src/diffraction.c
index fb993bf6..b71c28ce 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -155,18 +155,18 @@ struct rvec get_q(struct image *image, unsigned int xs, unsigned int ys,
 		  && (y >= image->det.panels[p].min_y)
 		  && (y <= image->det.panels[p].max_y) ) {
 			rx = ((float)xs - (sampling*image->det.panels[p].cx))
-			               / (sampling * image->resolution);
+			               / (sampling * image->det.panels[p].res);
 			ry = ((float)ys - (sampling*image->det.panels[p].cy))
-			               / (sampling * image->resolution);
+			               / (sampling * image->det.panels[p].res);
 			break;
 		}
 	}
 
 	/* Calculate q-vector for this sub-pixel */
 	r = sqrt(pow(rx, 2.0) + pow(ry, 2.0));
-	twothetax = atan2(rx, image->camera_len);
-	twothetay = atan2(ry, image->camera_len);
-	twotheta = atan2(r, image->camera_len);
+	twothetax = atan2(rx, image->det.panels[p].clen);
+	twothetay = atan2(ry, image->det.panels[p].clen);
+	twotheta = atan2(r, image->det.panels[p].clen);
 
 	if ( ttp != NULL ) *ttp = twotheta;
 
diff --git a/src/hdf5-file.c b/src/hdf5-file.c
index 974760be..0a63a853 100644
--- a/src/hdf5-file.c
+++ b/src/hdf5-file.c
@@ -202,14 +202,6 @@ int hdf5_read(struct hdfile *f, struct image *image)
 		image->lambda = ph_en_to_lambda(eV_to_J(2000.0));
 	}
 
-	/* These are only used for simulation (not analysis) */
-	image->x_centre = image->width/2;
-	image->y_centre = image->height/2;
-
-	/* FIXME: The following are basically made up... */
-	image->camera_len = 67.0e-3;  /* 75 mm camera length */
-	image->resolution = 13333.3;  /* 75 micron pixel size */
-
 	return 0;
 }
 
diff --git a/src/image.h b/src/image.h
index 78194433..543d924c 100644
--- a/src/image.h
+++ b/src/image.h
@@ -88,16 +88,12 @@ struct image {
 	 * If FORMULATION_PIXELSIZE, then pixel_size only is needed.*/
 	FormulationMode		fmode;
 	double			pixel_size;
-	float			camera_len;
-	float			resolution;	/* pixels per metre */
 
 	/* Wavelength must always be given */
 	double			lambda;		/* Wavelength in m */
 
 	int			width;
 	int			height;
-	double			x_centre;
-	double			y_centre;
 
 	ImageFeatureList	*features;	/* "Experimental" features */
 	ImageFeatureList	*rflist;	/* "Predicted" features */
diff --git a/src/index.c b/src/index.c
index 8faca4db..e9f5cf44 100644
--- a/src/index.c
+++ b/src/index.c
@@ -60,10 +60,10 @@ int map_position(struct image *image, double dx, double dy,
 	if ( image->fmode == FORMULATION_CLEN ) {
 
 		/* Convert pixels to metres */
-		x /= image->resolution;
-		y /= image->resolution;	/* Convert pixels to metres */
+		x /= image->det.panels[p].res;
+		y /= image->det.panels[p].res;	/* Convert pixels to metres */
 		d = sqrt((x*x) + (y*y));
-		twotheta = atan2(d, image->camera_len);
+		twotheta = atan2(d, image->det.panels[p].clen);
 
 	} else if (image->fmode == FORMULATION_PIXELSIZE ) {
 
diff --git a/src/indexamajig.c b/src/indexamajig.c
index daf37fc9..6915802a 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -180,6 +180,9 @@ int main(int argc, char *argv[])
 		image.det.panels[0].max_y = 1023;
 		image.det.panels[0].cx = 491.9;
 		image.det.panels[0].cy = 440.7;
+		image.det.panels[0].clen = 67.0e-3;
+		image.det.panels[0].res = 13333.3;  /* 75 micron pixel size */
+
 		/* Lower panel */
 		image.det.panels[1].min_x = 0;
 		image.det.panels[1].max_x = 1023;
@@ -187,6 +190,8 @@ int main(int argc, char *argv[])
 		image.det.panels[1].max_y = 511;
 		image.det.panels[1].cx = 492.0;
 		image.det.panels[1].cy = 779.7;
+		image.det.panels[1].clen = 75.0e-3;
+		image.det.panels[1].res = 13333.3;  /* 75 micron pixel size */
 
 		STATUS("Processing '%s'\n", line);
 
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index e6e96343..e389385d 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -224,10 +224,6 @@ int main(int argc, char *argv[])
 	image.width = 1024;
 	image.height = 1024;
 	image.fmode = FORMULATION_CLEN;
-	image.x_centre = 512.5;
-	image.y_centre = 512.5;
-	image.camera_len = 0.05;  /* 5 cm (front CCD can move from 5cm-20cm) */
-	image.resolution = 13333.3; /* 75 micron pixel size */
 	image.lambda = ph_en_to_lambda(eV_to_J(2.0e3));  /* Wavelength */
 	image.molecule = load_molecule();
 
diff --git a/src/peaks.c b/src/peaks.c
index c94b5aa8..fafbe696 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -450,9 +450,9 @@ void dump_peaks(struct image *image)
 			  && (y >= image->det.panels[p].min_y)
 			  && (y <= image->det.panels[p].max_y) ) {
 				rcx = ((double)x - image->det.panels[p].cx)
-				                            / image->resolution;
+				                     / image->det.panels[p].res;
 				rcy = ((double)y - image->det.panels[p].cy)
-				                            / image->resolution;
+				                     / image->det.panels[p].res;
 				found = 1;
 			}
 		}