Fixed eFOM

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

1 files changed, 20 insertions, 10 deletions

diff --git a/src/ipr.c b/src/ipr.c
index 982d8c3..aa55bd5 100644
--- a/src/ipr.c
+++ b/src/ipr.c
@@ -219,23 +219,33 @@ static double ipr_efom(ReflectionContext *rctx, Basis *basis) {
 		if ( cur->type == REFLECTION_NORMAL ) {
 			
 			/* Can this basis approximately account for this reflection? */
-			double n;
-			double p, q, r, u, v, w, d, e, f;
+			double det;
+			double a11, a12, a13, a21, a22, a23, a31, a32, a33;
 			double h, k, l;
 			
 			/* Set up the coordinate transform from hkl to xyz */
-			p = basis->a.x;  q = basis->a.y;  r = basis->a.z;
-			u = basis->b.x;  v = basis->b.y;  w = basis->b.z;
-			d = basis->c.x;  e = basis->c.y;  f = basis->c.z;
+			a11 = basis->a.x;  a12 = basis->a.y;  a13 = basis->a.z;
+			a21 = basis->b.x;  a22 = basis->b.y;  a23 = basis->b.z;
+			a31 = basis->c.x;  a32 = basis->c.y;  a33 = basis->c.z;
 			
 			/* Invert the matrix to get hkl from xyz */
-			n = p*(v*f - w*e) - q*(u*f - w*d) + r*(u*e - d*v);
-			h = (p*(v*f-w*e)*cur->x + -q*(u*f-w*d)*cur->y + r*(u*e-v*d)*cur->z) / n;
-			k = (-u*(q*f-r*e)*cur->x + v*(p*f-r*d)*cur->y + -w*(p*e-q*d)*cur->z) / n;
-			l = (d*(q*w-r*v)*cur->x + -e*(p*w-r*u)*cur->y + f*(p*v-q*u)*cur->z) / n;
+			det = a11*(a22*a33 - a23*a32) - a12*(a21*a33 - a23*a31) + a13*(a21*a32 - a22*a31);
+			h = ((a22*a33-a23*a32)*cur->x + (a23*a31-a21*a33)*cur->y + (a21*a32-a22*a31)*cur->z)/ det;
+			k = ((a13*a32-a12*a33)*cur->x + (a11*a33-a13*a31)*cur->y + (a12*a31-a11*a32)*cur->z) / det;
+			l = ((a12*a23-a13*a22)*cur->x + (a13*a21-a11*a23)*cur->y + (a11*a22-a12*a21)*cur->z) / det;
+			
+			//printf("%12.0f %12.0f %12.0f ----> %3.2f %3.2f %3.2f", cur->x, cur->y, cur->z, h, k, l);
+			//printf("%f %f     %f %f     %f %f\n", cur->x, h*basis->a.x + k*basis->b.x + l*basis->c.x,
+			//					cur->y, h*basis->a.y + k*basis->b.y + l*basis->c.y,
+			//					cur->z, h*basis->a.z + k*basis->b.z + l*basis->c.z);
 			
 			/* Calculate the deviations in terms of |a|, |b| and |c| */
+			h = fabs(h);  k = fabs(k);  l = fabs(l);
 			h -= floor(h);  k -= floor(k);  l -= floor(l);
+			if ( h == 1.0 ) h = 0.0;
+			if ( k == 1.0 ) k = 0.0;
+			if ( l == 1.0 ) l = 0.0;
+			//printf(" ----> %3.2f %3.2f %3.2f\n", h, k, l);
 			if ( h*h + k*k + l*l <= 0.1*0.1*0.1 ) n_indexed++;
 			
 			n_counted++;
@@ -274,7 +284,7 @@ static Basis *ipr_choose_initial_basis(ControlContext *ctx) {
 			if ( best_basis ) free(best_basis);
 			best_efom = efom;
 			best_basis = basis;
-			printf("IP:    %3i: eFOM=%f\n", i, efom);
+			printf("IP:    %3i: eFOM = %f %%\n", i, efom*100);
 		
 		}