More refinement work...

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

1 files changed, 16 insertions, 19 deletions

diff --git a/src/refine.c b/src/refine.c
index 14681c9..90a02b1 100644
--- a/src/refine.c
+++ b/src/refine.c
@@ -66,7 +66,7 @@ static void refine_fix_unconstrained(gsl_matrix *M) {
 			
 				int zero2 = 1;
 				
-				if ( row_altered ) break;
+				if ( row_altered ) break; /* This row is now fixed, so move on to the next */
 				
 				for ( i=0; i<M->size1; i++ ) {
 					if ( gsl_matrix_get(M, i, col) != 0.0 ) {
@@ -79,7 +79,6 @@ static void refine_fix_unconstrained(gsl_matrix *M) {
 					gsl_matrix_set(M, row, col, 1.0);
 					modified = 1;
 					row_altered = 1;
-					continue;
 				}
 			
 			}
@@ -116,6 +115,9 @@ ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image) {
 		int s;
 		double det;
 		
+		if ( axis == AXIS_X ) printf("RF: ------------------------------------------------------------------ Refining x coordinates -----\n");
+		if ( axis == AXIS_Y ) printf("RF: ------------------------------------------------------------------ Refining y coordinates -----\n");
+		
 		gsl_matrix_set_zero(M);
 		gsl_vector_set_zero(p);
 		
@@ -197,11 +199,11 @@ ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image) {
 		
 		/* Do the fitting */
 		refine_fix_unconstrained(M);
-		matrix_vector_show(M, p);
+		matrix_vector_show(M, p, "RF: ");
 		perm = gsl_permutation_alloc(M->size1);
 		gsl_linalg_LU_decomp(M, perm, &s);
 		det = gsl_linalg_LU_det(M, s);
-		printf("Determinant of M = %f\n", det);
+		printf("RF: Determinant of M = %f\n", det);
 		
 		q = gsl_vector_alloc(3);	/* This is the "answer" */
 		gsl_vector_set_zero(q);
@@ -230,31 +232,26 @@ ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image) {
 	gsl_vector_free(p);
 	gsl_vector_free(q);
 	
-	printf("a should change by %+7.5f %+7.5f nm^-1 in the image plane\n", dax/1e9, day/1e9);
-	printf("b should change by %+7.5f %+7.5f nm^-1 in the image plane\n", dbx/1e9, dby/1e9);
-	printf("c should change by %+7.5f %+7.5f nm^-1 in the image plane\n", dcx/1e9, dcy/1e9);
+	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);
+	printf("RF: c should change by %+7.5f %+7.5f nm^-1 in the image plane\n", dcx/1e9, dcy/1e9);
 	
 	/* Update the cell */
 	mapping_rotate(dax, day, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt);
-	printf("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,
+	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);
-	if ( dlx/cell->a.x < 0.1 ) cell->a.x += dlx;
-	if ( dly/cell->a.y < 0.1 ) cell->a.y += dly;
-	if ( dlz/cell->a.z < 0.1 ) cell->a.z += dlz;
+	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("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,
+	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);
-	if ( dlx/cell->b.x < 0.1 ) cell->b.x += dlx;
-	if ( dly/cell->b.y < 0.1 ) cell->b.y += dly;
-	if ( dlz/cell->b.z < 0.1 ) cell->b.z += dlz;
+	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("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,
+	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);
-	if ( dlx/cell->c.x < 0.1 ) cell->c.x += dlx;
-	if ( dly/cell->c.y < 0.1 ) cell->c.y += dly;
-	if ( dlz/cell->c.z < 0.1 ) cell->c.z += dlz;
+	cell->c.x += dlx;  cell->c.y += dly;  cell->c.z += dlz;
 	
 	return NULL;