More refinement (bored now..)

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

2 files changed, 31 insertions, 30 deletions

diff --git a/src/refine.c b/src/refine.c
index 90a02b1..4e5a3ab 100644
--- a/src/refine.c
+++ b/src/refine.c
@@ -95,9 +95,6 @@ static void refine_fix_unconstrained(gsl_matrix *M) {
 ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image) {
 	
 	ImageFeatureList *flist;
-	gsl_matrix *M;
-	gsl_vector *q;
-	gsl_vector *p;
 	int j, ns;
 	Axis axis;
 	double dax = 0.0, dbx = 0.0, dcx = 0.0;
@@ -105,19 +102,24 @@ ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image) {
 	double dlx = 0.0, dly = 0.0, dlz = 0.0;
 	
 	flist = image->features;
-	
-	M = gsl_matrix_alloc(3, 3);
-	p = gsl_vector_alloc(3);
-	
+		
 	for ( axis = AXIS_X; axis <= AXIS_Y; axis++ ) {
 	
 		gsl_permutation *perm;
 		int s;
 		double det;
+		gsl_matrix *M;
+		gsl_vector *q;
+		gsl_vector *p;
+		gsl_vector *scratch;
+		gsl_matrix *V;
+		gsl_vector *S;
 		
 		if ( axis == AXIS_X ) printf("RF: ------------------------------------------------------------------ Refining x coordinates -----\n");
 		if ( axis == AXIS_Y ) printf("RF: ------------------------------------------------------------------ Refining y coordinates -----\n");
 		
+		M = gsl_matrix_alloc(3, 3);
+		p = gsl_vector_alloc(3);
 		gsl_matrix_set_zero(M);
 		gsl_vector_set_zero(p);
 		
@@ -144,15 +146,6 @@ ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image) {
 			diy = rf->partner->y - rf->y;
 			printf("RF: Feature %3i: %3i %3i %3i dev = %+9.5f %+9.5f px ", j, h, k, l, dix, diy);
 			
-			//if ( (h<0) || (k<0) || (l<0) ) {
-			//	h = -h;
-			//	k = -k;
-			//	l = -l;
-			//	dix = -dix;
-			//	diy = -diy;
-			//	printf("\n          Fudge: %3i %3i %3i dev = %+9.5f %+9.5f px ", h, k, l, dix, diy);
-			//}
-			
 			old_x = rf->partner->x;
 			old_y = rf->partner->y;
 			rf->partner->x = dix + rf->partner->parent->x_centre;
@@ -197,17 +190,30 @@ ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image) {
 			return NULL;
 		}	
 		
-		/* Do the fitting */
+		/* Prepare for solving */
 		refine_fix_unconstrained(M);
 		matrix_vector_show(M, p, "RF: ");
+		
+		/* Calculate and display the determinant */
 		perm = gsl_permutation_alloc(M->size1);
 		gsl_linalg_LU_decomp(M, perm, &s);
 		det = gsl_linalg_LU_det(M, s);
 		printf("RF: Determinant of M = %f\n", det);
+		gsl_permutation_free(perm);		
 		
+		/* Solve the equations */
+		V = gsl_matrix_alloc(M->size2, M->size2);
+		S = gsl_vector_alloc(M->size2);
+		scratch = gsl_vector_alloc(M->size2);
+		gsl_linalg_SV_decomp(M, V, S, scratch);
 		q = gsl_vector_alloc(3);	/* This is the "answer" */
 		gsl_vector_set_zero(q);
-		gsl_linalg_LU_solve(M, perm, p, q);
+		gsl_linalg_SV_solve(M, V, S, p, q);
+		gsl_vector_free(scratch);
+		gsl_matrix_free(V);
+		gsl_vector_free(S);
+		gsl_matrix_free(M);
+		gsl_vector_free(p);
 		
 		switch ( axis ) {
 			case AXIS_X : {
@@ -224,14 +230,10 @@ ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image) {
 			}
 		}
 		
-		gsl_permutation_free(perm);
-		
+		gsl_vector_free(q);
+
 	}
 	
-	gsl_matrix_free(M);
-	gsl_vector_free(p);
-	gsl_vector_free(q);
-	
 	printf("RF: ------------------------------------------------------------------ Refinement results ---------\n");
 	printf("RF: a should change by %+7.5f %+7.5f nm^-1 in the image plane\n", dax/1e9, day/1e9);
 	printf("RF: b should change by %+7.5f %+7.5f nm^-1 in the image plane\n", dbx/1e9, dby/1e9);
@@ -239,18 +241,15 @@ ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image) {
 	
 	/* Update the cell */
 	mapping_rotate(dax, day, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt);
-	printf("RF: a changed by %+7.5f %+7.5f %+7.5f nm^-1 in reciprocal space (%+10.5f %+10.5f %+10.5f %%)\n", dlx/1e9, dly/1e9, dlz/1e9,
-										  100*dlx/cell->a.x, 100*dly/cell->a.y, 100*dlz/cell->a.z);
+	printf("RF: a changed by [ %+7.5f %+7.5f %+7.5f ] nm^-1 in reciprocal space\n", dlx/1e9, dly/1e9, dlz/1e9);
 	cell->a.x += dlx;  cell->a.y += dly; cell->a.z += dlz;
 	
 	mapping_rotate(dbx, dby, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt);
-	printf("RF: b changed by %+7.5f %+7.5f %+7.5f nm^-1 in reciprocal space (%+10.5f %+10.5f %+10.5f %%)\n", dlx/1e9, dly/1e9, dlz/1e9,
-										  100*dlx/cell->b.x, 100*dly/cell->b.y, 100*dlz/cell->b.z);
+	printf("RF: b changed by [ %+7.5f %+7.5f %+7.5f ] nm^-1 in reciprocal space\n", dlx/1e9, dly/1e9, dlz/1e9);
 	cell->b.x += dlx;  cell->b.y += dly;  cell->b.z += dlz;
 	
 	mapping_rotate(dcx, dcy, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt);
-	printf("RF: c changed by %+7.5f %+7.5f %+7.5f nm^-1 in reciprocal space (%+10.5f %+10.5f %+10.5f %%)\n", dlx/1e9, dly/1e9, dlz/1e9,
-										  100*dlx/cell->c.x, 100*dly/cell->c.y, 100*dlz/cell->c.z);
+	printf("RF: c changed by [ %+7.5f %+7.5f %+7.5f ] nm^-1 in reciprocal space\n", dlx/1e9, dly/1e9, dlz/1e9);
 	cell->c.x += dlx;  cell->c.y += dly;  cell->c.z += dlz;
 	
 	return NULL;
diff --git a/src/reflections.c b/src/reflections.c
index 8bfa7ff..58a7992 100644
--- a/src/reflections.c
+++ b/src/reflections.c
@@ -235,6 +235,8 @@ ReflectionList *reflection_list_from_cell(Basis *basis) {
 
 				double x, y, z;
 				
+				if ( h != 0 ) continue;
+				
 				x = h*basis->a.x + k*basis->b.x + l*basis->c.x;
 				y = h*basis->a.y + k*basis->b.y + l*basis->c.y;
 				z = h*basis->a.z + k*basis->b.z + l*basis->c.z;