From c85a8a308ec7ab50545255530da6a31b841be3ed Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 4 Feb 2010 15:35:34 +0100 Subject: Better cell reduction debugging messages --- src/cell.c | 48 +++++++++++++++++++++++------------------------- 1 file changed, 23 insertions(+), 25 deletions(-) diff --git a/src/cell.c b/src/cell.c index 815873b9..27f6e3a3 100644 --- a/src/cell.c +++ b/src/cell.c @@ -319,7 +319,6 @@ static void cell_print(UnitCell *cell) double asx, asy, asz; double bsx, bsy, bsz; double csx, csy, csz; - double lengths[3]; double angles[3]; STATUS(" a b c alpha beta gamma\n"); @@ -331,24 +330,24 @@ static void cell_print(UnitCell *cell) &bsx, &bsy, &bsz, &csx, &csy, &csz); - STATUS("astar = %10.3e %10.3e %10.3e\n", asx, asy, asz); - STATUS("bstar = %10.3e %10.3e %10.3e\n", bsx, bsy, bsz); - STATUS("cstar = %10.3e %10.3e %10.3e\n", csx, csy, csz); + STATUS("a = %10.3e %10.3e %10.3e m\n", cell->ax, cell->ay, cell->az); + STATUS("b = %10.3e %10.3e %10.3e m\n", cell->bx, cell->by, cell->bz); + STATUS("c = %10.3e %10.3e %10.3e m\n", cell->cx, cell->cy, cell->cz); - lengths[0] = modulus(asx, asy, asz); - lengths[1] = modulus(bsx, bsy, bsz); - lengths[2] = modulus(csx, csy, csz); + STATUS("astar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n", + asx, asy, asz, modulus(asx, asy, asz)); + STATUS("bstar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n", + bsx, bsy, bsz, modulus(bsx, bsy, bsz)); + STATUS("cstar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n", + csx, csy, csz, modulus(csx, csy, csz)); angles[0] = angle_between(bsx, bsy, bsz, csx, csy, csz); angles[1] = angle_between(asx, asy, asz, csx, csy, csz); angles[2] = angle_between(asx, asy, asz, bsx, bsy, bsz); - STATUS("Checking %9.3e %9.3e %9.3e m^-1\n", lengths[0], - lengths[1], - lengths[2]); - STATUS("Checking %5.2f %5.2f %5.2f deg\n", rad2deg(angles[0]), - rad2deg(angles[1]), - rad2deg(angles[2])); +// STATUS("Checking %5.2f %5.2f %5.2f deg\n", rad2deg(angles[0]), +// rad2deg(angles[1]), +// rad2deg(angles[2])); } @@ -383,10 +382,16 @@ UnitCell *match_cell(UnitCell *cell, UnitCell *template) struct cvec *cand[3]; int ncand[3] = {0,0,0}; - float ltl = 5.0; - float angtol = 5.0; + float ltl = 15.0; /* percent */ + float angtol = deg2rad(9.0); UnitCell *new_cell = NULL; + STATUS("Matching this experimental cell: --------------------------\n"); + cell_print(cell); + STATUS("With this model cell: -------------------------------------\n"); + cell_print(template); + STATUS("-----------------------------------------------------------\n"); + cell_get_reciprocal(template, &asx, &asy, &asz, &bsx, &bsy, &bsz, &csx, &csy, &csz); @@ -399,13 +404,6 @@ UnitCell *match_cell(UnitCell *cell, UnitCell *template) angles[1] = angle_between(asx, asy, asz, csx, csy, csz); angles[2] = angle_between(asx, asy, asz, bsx, bsy, bsz); - STATUS("Looking for %9.3e %9.3e %9.3e m^-1\n", lengths[0], - lengths[1], - lengths[2]); - STATUS("Looking for %5.2f %5.2f %5.2f deg\n", rad2deg(angles[0]), - rad2deg(angles[1]), - rad2deg(angles[2])); - cand[0] = malloc(MAX_CAND*sizeof(struct cvec)); cand[1] = malloc(MAX_CAND*sizeof(struct cvec)); cand[2] = malloc(MAX_CAND*sizeof(struct cvec)); @@ -476,7 +474,7 @@ UnitCell *match_cell(UnitCell *cell, UnitCell *template) cand[1][j].vec.v, cand[1][j].vec.w); /* Angle between axes 0 and 1 should be angle 2 */ - if ( !within_tolerance(ang, angles[2], angtol) ) continue; + if ( fabs(ang - angles[2]) > angtol ) continue; for ( k=0; k angtol ) continue; /* Finally, the angle between the current candidate for * axis 1 and the kth candidate for axis 2 */ @@ -496,7 +494,7 @@ UnitCell *match_cell(UnitCell *cell, UnitCell *template) cand[2][k].vec.v, cand[2][k].vec.w); /* ... it should be angle 0 ... */ - if ( !within_tolerance(ang, angles[0], angtol) ) continue; + if ( fabs(ang - angles[0]) > angtol ) continue; new_cell = cell_new_from_axes(cand[0][i].vec, cand[1][j].vec, -- cgit v1.2.3