From 9775637ff55ba497168d4778e2c1d8ca4b6f481a Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Tue, 10 May 2011 16:09:07 +0200
Subject: Update debugging output

---
 src/post-refinement.c | 65 ++++++++++++++++++++++++++++++++++++++++-----------
 1 file changed, 52 insertions(+), 13 deletions(-)

(limited to 'src')

diff --git a/src/post-refinement.c b/src/post-refinement.c
index 71156026..1c98346f 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -185,7 +185,9 @@ static void apply_cell_shift(UnitCell *cell, int k, double shift)
 	cell_get_reciprocal(cell, &asx, &asy, &asz,
 	                          &bsx, &bsy, &bsz,
 	                          &csx, &csy, &csz);
-	STATUS("Shifting by %e\n", shift);
+
+	//STATUS("Shifting by %e\n", shift);
+
 	switch ( k )
 	{
 		case REF_ASX :  asx += shift;  break;
@@ -281,11 +283,11 @@ static double pr_iterate(struct image *image, const RefList *full,
 
 		/* Find the full version */
 		get_indices(refl, &ha, &ka, &la);
-		if ( (ha!=23) || (ka!=12) || (la!=3) ) continue;
 		match = find_refl(full, ha, ka, la);
-		assert(match != NULL);  /* Never happens because all scalable
-		                         * reflections had their LSQ intensities
-		                         * calculated in lsq_intensities(). */
+		if ( match == NULL ) continue;
+		/* Some reflections may have recently become scalable, but
+		 * scale_intensities() might not yet have been called, so the
+		 * full version may not have been calculated yet. */
 		I_full = get_intensity(match);
 
 		/* Actual measurement of this reflection from this pattern? */
@@ -326,9 +328,8 @@ static double pr_iterate(struct image *image, const RefList *full,
 		}
 
 	}
-	double tg = gsl_matrix_get(M, 0, 0);
-	STATUS("total gradient = %e\n", tg);
-	show_matrix_eqn(M, v, NUM_PARAMS);
+	//STATUS("total gradient = %e\n",  gsl_matrix_get(M, 0, 0));
+	//show_matrix_eqn(M, v, NUM_PARAMS);
 
 	shifts = gsl_vector_alloc(NUM_PARAMS);
 	max_shift = 0.0;
@@ -352,6 +353,28 @@ static double pr_iterate(struct image *image, const RefList *full,
 }
 
 
+static void check_scalable(RefList *check)
+{
+	Reflection *refl;
+	RefListIterator *iter;
+
+	STATUS("------ Scalability check ------\n");
+	for ( refl = first_refl(check, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) )
+	{
+
+		signed int h, k, l;
+		if ( get_scalable(refl) ) continue;
+
+		get_indices(refl, &h, &k, &l);
+		STATUS("%3i %3i %3i is not scalable.\n", h, k, l);
+
+	}
+	STATUS("------ End of scalability check ------\n");
+}
+
+
 static double mean_partial_dev(struct image *image,
                                const RefList *full, const char *sym)
 {
@@ -373,11 +396,15 @@ static double mean_partial_dev(struct image *image,
 		if ( !get_scalable(refl) ) continue;
 
 		get_indices(refl, &h, &k, &l);
-		if ( (h!=23) || (k!=12) || (l!=3) ) continue;
 		assert ((h!=0) || (k!=0) || (l!=0));
 
+		if ( !get_scalable(refl) ) continue;
+
 		full_version = find_refl(full, h, k, l);
-		assert(full_version != NULL);
+		if ( full_version == NULL ) continue;
+		/* Some reflections may have recently become scalable, but
+		 * scale_intensities() might not yet have been called, so the
+		 * full version may not have been calculated yet. */
 
 		G = image->osf;
 		p = get_partiality(refl);
@@ -431,20 +458,28 @@ void pr_refine(struct image *image, const RefList *full, const char *sym)
 {
 	double max_shift, dev;
 	int i;
+	double ax, ay, az;
+	double bx, by, bz;
+	double cx, cy, cz;
+
+	cell_get_reciprocal(image->indexed_cell, &ax, &ay, &az,
+	                    &bx, &by, &bz, &cx, &cy, &cz);
+	STATUS("At start of post refinement, ax*=%e\n", ax);
+
+	dev = mean_partial_dev(image, full, sym);
+	//STATUS("PR starting dev = %5.2f\n", dev);
 
 	/* FIXME: This is for debugging */
 	//plot_curve(image, full, sym);
 	//return;
 
-	dev = mean_partial_dev(image, full, sym);
-	STATUS("PR starting dev = %5.2f\n", dev);
-
 	i = 0;
 	do {
 
 		double dev;
 
 		max_shift = pr_iterate(image, full, sym);
+
 		update_partialities(image, sym, NULL, NULL, NULL, NULL);
 
 		dev = mean_partial_dev(image, full, sym);
@@ -454,4 +489,8 @@ void pr_refine(struct image *image, const RefList *full, const char *sym)
 		i++;
 
 	} while ( (max_shift > 0.01) && (i < MAX_CYCLES) );
+
+	cell_get_reciprocal(image->indexed_cell, &ax, &ay, &az,
+	                    &bx, &by, &bz, &cx, &cy, &cz);
+	STATUS("At end of post refinement, ax*=%e\n", ax);
 }
-- 
cgit v1.2.3