Restore working refinement...

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

3 files changed, 22 insertions, 24 deletions

diff --git a/src/refine.c b/src/refine.c
index 987fba1..5de8df3 100644
--- a/src/refine.c
+++ b/src/refine.c
@@ -33,10 +33,10 @@
 #define DISPFACTOR 1.0e9
 
 /* Number of parameters */
-#define NUM_PARAMS 4
+#define NUM_PARAMS 9
 
 /* Refine debug */
-#define REFINE_DEBUG 1
+#define REFINE_DEBUG 0
 
 /* A simplex is an array of ten of these */
 typedef struct {
@@ -408,8 +408,6 @@ double refine_do_cell(ControlContext *ctx) {
 			d[f].dy = dly;
 			d[f].dz = dlz;
 			
-			printf("Feature %3i: %8f %8f %8f\n", f, dlx/DISPFACTOR, dly/DISPFACTOR, dlz/DISPFACTOR);
-			
 			f++;
 		
 		}
@@ -434,7 +432,7 @@ double refine_do_cell(ControlContext *ctx) {
 	memcpy(&s[9], &s[0], sizeof(SimplexVertex));	s[9].dcz = delta;
 	
 	/* Iterate */
-	maxiter = 5000;
+	maxiter = 500;
 	for ( it=0; it<maxiter; it++ ) {
 	
 		double conv;
diff --git a/src/refinetest2d.c b/src/refinetest2d.c
index 03b65d6..cdc0c8d 100644
--- a/src/refinetest2d.c
+++ b/src/refinetest2d.c
@@ -30,17 +30,17 @@ static int check_cell(Basis *cell, Basis *cell_real) {
 	
 	int fail;
 	
-//	printf("     Calculated    Actual\n");
-//	printf("---------------------------\n");
-//	printf("ax  %+8f     %+8f\n", cell->a.x/1e9, cell_real->a.x/1e9);
-//	printf("ay  %+8f     %+8f\n", cell->a.y/1e9, cell_real->a.y/1e9);
-//	printf("az  %+8f     %+8f\n", cell->a.z/1e9, cell_real->a.z/1e9);
-//	printf("bx  %+8f     %+8f\n", cell->b.x/1e9, cell_real->b.x/1e9);
-//	printf("by  %+8f     %+8f\n", cell->b.y/1e9, cell_real->b.y/1e9);
-//	printf("bz  %+8f     %+8f\n", cell->b.z/1e9, cell_real->b.z/1e9);
-//	printf("cx  %+8f     %+8f\n", cell->c.x/1e9, cell_real->c.x/1e9);
-//	printf("cy  %+8f     %+8f\n", cell->c.y/1e9, cell_real->c.y/1e9);
-//	printf("cz  %+8f     %+8f\n", cell->c.z/1e9, cell_real->c.z/1e9);
+	printf("     Refinement    Actual\n");
+	printf("---------------------------\n");
+	printf("ax  %+8f     %+8f\n", cell->a.x/1e9, cell_real->a.x/1e9);
+	printf("ay  %+8f     %+8f\n", cell->a.y/1e9, cell_real->a.y/1e9);
+	printf("az  %+8f     %+8f\n", cell->a.z/1e9, cell_real->a.z/1e9);
+	printf("bx  %+8f     %+8f\n", cell->b.x/1e9, cell_real->b.x/1e9);
+	printf("by  %+8f     %+8f\n", cell->b.y/1e9, cell_real->b.y/1e9);
+	printf("bz  %+8f     %+8f\n", cell->b.z/1e9, cell_real->b.z/1e9);
+	printf("cx  %+8f     %+8f\n", cell->c.x/1e9, cell_real->c.x/1e9);
+	printf("cy  %+8f     %+8f\n", cell->c.y/1e9, cell_real->c.y/1e9);
+	printf("cz  %+8f     %+8f\n", cell->c.z/1e9, cell_real->c.z/1e9);
 	
 	fail = 0;
 	if ( fabs(cell->a.x - cell_real->a.x) > 0.01e9 ) {
@@ -107,7 +107,7 @@ int main(int argc, char *argv[]) {
 	ctx->x_centre = 256;
 	ctx->y_centre = 256;
 	ctx->pixel_size = 5e7;
-	image_add(ctx->images, NULL, 512, 512, deg2rad(90.0), ctx);
+	image_add(ctx->images, NULL, 512, 512, deg2rad(0.0), ctx);
 	ctx->images->images[0].features = image_feature_list_new();
 	
 	/* Fudge to avoid horrifying pointer-related death */
@@ -122,7 +122,7 @@ int main(int argc, char *argv[]) {
 	/* The "real" reflections */
 	reflections_real = reflection_list_from_cell(cell_real);
 	ctx->images->images[0].features = reproject_get_reflections(&ctx->images->images[0], reflections_real);
-	printf("RT: %i test features generated.\n", ctx->images->images[0].features->n_features);
+//	printf("RT: %i test features generated.\n", ctx->images->images[0].features->n_features);
 	
 	/* The "model" cell to be refined */
 	ctx->cell = malloc(sizeof(Basis));
@@ -145,6 +145,8 @@ int main(int argc, char *argv[]) {
 	
 	if ( fail ) return 1;
 	
+	printf("\n2D refinement test OK.\n");
+	
 	return 0;
 	
 }
diff --git a/src/reproject.c b/src/reproject.c
index 439c2f7..e1ce13e 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -85,7 +85,7 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 	
 		double xl, yl, zl;
 		double a, b, c;
-		double A1, A2, s1, s2, s, temp;
+		double A1, A2, s1, s2, s;
 		
 		/* Get the coordinates of the reciprocal lattice point */
 		xl = reflection->x;
@@ -95,11 +95,9 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 		/* Next, solve the relrod equation to calculate the excitation error */
 		a = 1.0;
 		b = -2.0*(xl*nx + yl*ny + zl*nz);
-		c = xl*xl + yl*yl + zl*zl - 1.0/(image->lambda*image->lambda); /* FIXME: Don't think this is stable */
-		/* Solve the quadratic equation */
-		temp = -0.5 * (b + sign(b)*sqrt(b*b - 4*a*c));
-		A1 = temp / a;
-		A2 = c / temp;
+		c = xl*xl + yl*yl + zl*zl - 1.0/(image->lambda*image->lambda);
+		A1 = (-b + sqrt(b*b-4.0*a*c))/(2.0*a);
+		A2 = (-b - sqrt(b*b-4.0*a*c))/(2.0*a);
 		s1 = 1.0/image->lambda - A1;
 		s2 = 1.0/image->lambda - A2;
 		if ( fabs(s1) < fabs(s2) ) s = s1; else s = s2;