Calculate reflections using reciprocal cell

1 files changed, 56 insertions, 58 deletions

diff --git a/src/relrod.c b/src/relrod.c
index df4c4df3..e344133a 100644
--- a/src/relrod.c
+++ b/src/relrod.c
@@ -62,14 +62,18 @@ void get_reflections(struct image *image, UnitCell *cell)
 	double kx, ky, kz; /* Electron wavevector ("normal" times 1/lambda) */
 	double ux, uy, uz; /* "up" vector */
 	double rx, ry, rz; /* "right" vector */
-	double ax, ay, az; /* "a" lattice parameter */
-	double bx, by, bz; /* "b" lattice parameter */
-	double cx, cy, cz; /* "c" lattice parameter */
+	double asx, asy, asz; /* "a*" lattice parameter */
+	double bsx, bsy, bsz; /* "b*" lattice parameter */
+	double csx, csy, csz; /* "c*" lattice parameter */
 	signed int h, k, l;
 
-	flist = image_feature_list_new();
-	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+	/* Get the reciprocal unit cell */
+	cell_get_reciprocal(cell, &asx, &asy, &asz,
+	                          &bsx, &bsy, &bsz,
+	                          &csx, &csy, &csz);
 
+	/* Prepare list and some parameters */
+	flist = image_feature_list_new();
 	tilt = image->tilt;
 	omega = image->omega;
 	wavenumber = 1.0/image->lambda;
@@ -96,9 +100,9 @@ void get_reflections(struct image *image, UnitCell *cell)
 		double g_sq, gn;
 
 		/* Get the coordinates of the reciprocal lattice point */
-		xl = h*ax + k*bx + l*cx;
-		yl = h*ay + k*by + l*cy;
-		zl = h*az + k*bz + l*cz;
+		xl = h*asx + k*bsx + l*csx;
+		yl = h*asy + k*bsy + l*csy;
+		zl = h*asz + k*bsz + l*csz;
 		g_sq = modulus_squared(xl, yl, zl);
 		gn = xl*nx + yl*ny + zl*nz;
 
@@ -119,6 +123,7 @@ void get_reflections(struct image *image, UnitCell *cell)
 			double gx, gy, gz;
 			double theta;
 			double x, y;
+			double dx, dy, psi;
 
 			/* Determine the intersection point */
 			xi = xl + s*nx;  yi = yl + s*ny;  zi = zl + s*nz;
@@ -130,56 +135,49 @@ void get_reflections(struct image *image, UnitCell *cell)
 					*  the sphere to the intersection */
 			theta = angle_between(-kx, -ky, -kz, gx, gy, gz);
 
-			if ( theta > 0.0 ) {
-
-				double dx, dy, psi;
-
-				/* Calculate azimuth of point in image
-				 * (anticlockwise from +x) */
-				dx = xi*rx + yi*ry + zi*rz;
-				dy = xi*ux + yi*uy + zi*uz;
-				psi = atan2(dy, dx);
-
-				/* Get image coordinates from polar
-				 * representation */
-				if ( image->fmode == FORMULATION_CLEN ) {
-					x = image->camera_len*tan(theta)*cos(psi);
-					y = image->camera_len*tan(theta)*sin(psi);
-					x *= image->resolution;
-					y *= image->resolution;
-				} else if ( image->fmode==FORMULATION_PIXELSIZE ) {
-					x = tan(theta)*cos(psi) / image->lambda;
-					y = tan(theta)*sin(psi) / image->lambda;
-					x /= image->pixel_size;
-					y /= image->pixel_size;
-				} else {
-					fprintf(stderr,
-						"Unrecognised formulation mode "
-						"in get_reflections\n");
-					return;
-				}
-
-				x += image->x_centre;
-				y += image->y_centre;
-
-				/* Sanity check */
-				if ( (x>=0) && (x<image->width)
-				  && (y>=0) && (y<image->height) ) {
-
-					/* Record the reflection
-					 *  Intensity should be multiplied by
-					 *  relrod spike function except
-					 *  reprojected reflections aren't used
-					 *  quantitatively for anything
-					 */
-					image_add_feature_reflection(flist, x,
-						y, image, 1.0);
-
-				} /* else it's outside the picture somewhere */
-
-			} /* else this is the central beam */
-
-		}
+			/* Calculate azimuth of point in image
+			 * (anticlockwise from +x) */
+			dx = xi*rx + yi*ry + zi*rz;
+			dy = xi*ux + yi*uy + zi*uz;
+			psi = atan2(dy, dx);
+
+			/* Get image coordinates from polar
+			 * representation */
+			if ( image->fmode == FORMULATION_CLEN ) {
+				x = image->camera_len*tan(theta)*cos(psi);
+				y = image->camera_len*tan(theta)*sin(psi);
+				x *= image->resolution;
+				y *= image->resolution;
+			} else if ( image->fmode==FORMULATION_PIXELSIZE ) {
+				x = tan(theta)*cos(psi) / image->lambda;
+				y = tan(theta)*sin(psi) / image->lambda;
+				x /= image->pixel_size;
+				y /= image->pixel_size;
+			} else {
+				fprintf(stderr,
+					"Unrecognised formulation mode "
+					"in get_reflections\n");
+				return;
+			}
+
+			x += image->x_centre;
+			y += image->y_centre;
+
+			/* Sanity check */
+			if ( (x>=0) && (x<image->width)
+			  && (y>=0) && (y<image->height) ) {
+
+				/* Record the reflection.
+				 * Intensity should be multiplied by relrod
+				 * spike function, except reprojected
+				 * reflections aren't used quantitatively for
+				 * anything. */
+				 printf("%i %i %i\n", h,k,l);
+				image_add_feature(flist, x, y, image, 1.0);
+
+			} /* else it's outside the picture somewhere */
+
+		} /* else reflection is not excited in this orientation */
 
 	}
 	}