More refinement work...

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

4 files changed, 23 insertions, 24 deletions

diff --git a/src/iprtest.c b/src/iprtest.c
index cb2fb10..7f9541d 100644
--- a/src/iprtest.c
+++ b/src/iprtest.c
@@ -117,8 +117,8 @@ int main(int argc, char *argv[]) {
 	
 	ctx->image = malloc(sizeof(ImageRecord));
 	ctx->image->image = NULL;
-	ctx->image->tilt = 20.0;
-	ctx->image->omega = 0.0;
+	ctx->image->tilt = 10.0;
+	ctx->image->omega = 90.0;
 	ctx->image->slop = 0.0;
 	ctx->image->fmode = FORMULATION_PIXELSIZE;
 	ctx->image->pixel_size = 0.6e8;
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;
 	
diff --git a/src/utils.c b/src/utils.c
index 8fc3cf6..1e82abb 100644
--- a/src/utils.c
+++ b/src/utils.c
@@ -96,12 +96,14 @@ void chomp(char *s) {
 	
 }
 
-void matrix_vector_show(const gsl_matrix *m, const gsl_vector *v) {
+void matrix_vector_show(const gsl_matrix *m, const gsl_vector *v, const char *prefix) {
 
 	int i, j;
 	
+	if ( prefix == NULL ) prefix = "";
+	
 	for ( i=0; i<m->size1; i++ ) {
-		printf("[ ");
+		printf("%s[ ", prefix);
 		for ( j=0; j<m->size2; j++ ) {
 			printf("%12.8f ", gsl_matrix_get(m, i, j));
 		}
diff --git a/src/utils.h b/src/utils.h
index 47e283a..804f1af 100644
--- a/src/utils.h
+++ b/src/utils.h
@@ -30,7 +30,7 @@ extern double lambda(double voltage);
 extern double distance3d(double x1, double y1, double z1, double x2, double y2, double z2);
 extern size_t skipspace(const char *s);
 extern void chomp(char *s);
-extern void matrix_vector_show(const gsl_matrix *m, const gsl_vector *v);
+extern void matrix_vector_show(const gsl_matrix *m, const gsl_vector *v, const char *prefix);
 
 #define rad2deg(a) ((a)*180/M_PI)
 #define deg2rad(a) ((a)*M_PI/180)