Merge the two reflection_add routines for different formulations

git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@95 bf6ca9ba-c028-0410-8290-897cf20841d1

1 files changed, 29 insertions, 68 deletions

diff --git a/src/reflections.c b/src/reflections.c
index d53af10..a1a8267 100644
--- a/src/reflections.c
+++ b/src/reflections.c
@@ -93,85 +93,46 @@ Reflection *reflection_add(ReflectionContext *reflectionctx, double x, double y,
 	
 }
 
-/* x and y in metres (in real space!) */
-void reflection_add_from_dp(ControlContext *ctx, double x, double y, double tilt_degrees, double intensity) {
+/* x and y in pixels, measured from centre of image */
+void reflection_add_from_dp(ControlContext *ctx, double x, double y, ImageRecord imagerecord, double intensity) {
 
-	double nx, ny, nz;
-
-	/* Real space */
-	double L;
+	/* "Input" space */
 	double d;
-
-	/* Shared */
-	double theta;
-	double psi;
-	
-	/* Reciprocal space */
-	double r;
-	double tilt;
-	double omega;
-	double x_temp, y_temp, z_temp;
 	
-	tilt = 2*M_PI*(tilt_degrees/360);	/* Convert to Radians */
-	omega = 2*M_PI*(ctx->omega/360);	/* Likewise */
-	
-	d = sqrt((x*x) + (y*y));
-	L = ctx->camera_length;
-	theta = atan2(d, L);
-	psi = atan2(y, x);
-	
-	r = 1/ctx->lambda;
-	x_temp = r*sin(theta)*cos(psi);
-	y_temp = -r*sin(theta)*sin(psi);	/* Minus sign to define axes as y going upwards */
-	z_temp = r- r*cos(theta);
-	
-	/* Apply the rotations...
-		First: rotate image clockwise until tilt axis is aligned horizontally. */
-	nx = x_temp*cos(omega) + y_temp*-sin(omega);
-	ny = x_temp*sin(omega) + y_temp*cos(omega);
-	nz = z_temp;
-	
-	/* Now, tilt about the x-axis ANTICLOCKWISE around +x, i.e. the "wrong" way.
-		This is because the crystal is rotated in the experiment, not the Ewald sphere. */
-	x_temp = nx; y_temp = ny; z_temp = nz;	
-	nx = x_temp;
-	ny = cos(tilt)*y_temp - sin(tilt)*z_temp;
-	nz = -sin(tilt)*y_temp - cos(tilt)*z_temp;
-	
-	reflection_add(ctx->reflectionctx, nx, ny, nz, intensity, REFLECTION_NORMAL);
-
-}
-
-/* Alternative formulation for when the input data is already in reciprocal space units 
-	x and y in metres^-1 (in reciprocal space) */
-void reflection_add_from_reciprocal(ControlContext *ctx, double x, double y, double tilt_degrees, double intensity) {
-
-	double nx, ny, nz;
-	
-	/* Input (reciprocal space) */
-	double dr;
-	
-	/* Shared */
-	double theta;
-	double psi;
-	double r;
+	/* Angular description of reflection */
+	double theta, psi, k;
 	
 	/* Reciprocal space */
 	double tilt;
 	double omega;
-	double x_temp, y_temp, z_temp;
 	
-	tilt = M_PI*(tilt_degrees/180);	/* Convert to Radians */
-	omega = M_PI*(ctx->omega/180);	/* Likewise */
+	double x_temp, y_temp, z_temp;
+	double nx, ny, nz;
 	
-	dr = sqrt((x*x) + (y*y));
-	r = 1/ctx->lambda;
-	theta = atan2(dr, r);
+	tilt = 2*M_PI*(imagerecord.tilt/360);	/* Convert to Radians */
+	omega = 2*M_PI*(imagerecord.omega/360);	/* Likewise */
+	k = 1/imagerecord.lambda;
+	
+	/* Calculate an angular description of the reflection */
+	if ( ctx->fmode == FORMULATION_CLEN ) {
+		x /= imagerecord.resolution;
+		y /= imagerecord.resolution;	/* Convert pixels to metres */
+		d = sqrt((x*x) + (y*y));
+		theta = atan2(d, imagerecord.camera_len);
+	} else if ( ctx->fmode == FORMULATION_PIXELSIZE ) {
+		x *= imagerecord.pixel_size;
+		y *= imagerecord.pixel_size;	/* Convert pixels to metres^-1 */
+		d = sqrt((x*x) + (y*y));
+		theta = atan2(d, k);
+	} else {
+		fprintf(stderr, "Unrecognised formulation mode in reflection_add_from_dp\n");
+		return;
+	}
 	psi = atan2(y, x);
 	
-	x_temp = r*sin(theta)*cos(psi);
-	y_temp = -r*sin(theta)*sin(psi);	/* Minus sign to define axes as y going upwards */
-	z_temp = r - r*cos(theta);
+	x_temp = k*sin(theta)*cos(psi);
+	y_temp = -k*sin(theta)*sin(psi);	/* Minus sign to define axes as y going upwards */
+	z_temp = k- k*cos(theta);
 	
 	/* Apply the rotations...
 		First: rotate image clockwise until tilt axis is aligned horizontally. */