aboutsummaryrefslogtreecommitdiff
path: root/src/cell.c
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2010-03-08 18:46:03 +0100
committerThomas White <taw@physics.org>2010-03-08 18:46:03 +0100
commit68f81a0b0950a6e269a80a57a42eca62b29f2198 (patch)
treef1ef101f95af8ecf7380139f72f1ce95b3b9dd84 /src/cell.c
parente8de19bdf531faa2e759ffd40fa54b033f2936d5 (diff)
Fix and tidy-up the cell matching
Diffstat (limited to 'src/cell.c')
-rw-r--r--src/cell.c12
1 files changed, 8 insertions, 4 deletions
diff --git a/src/cell.c b/src/cell.c
index fb621911..c564c3b5 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -401,7 +401,7 @@ UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose)
double angles[3];
struct cvec *cand[3];
UnitCell *new_cell = NULL;
- float best_fom = HUGE_VAL;
+ float best_fom = +HUGE_VALF;
int ncand[3] = {0,0,0};
float ltl = 5.0; /* percent */
float angtol = deg2rad(5.0);
@@ -457,13 +457,14 @@ UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose)
n1 *= b1; n2 *= b2; n3 *= b3;
- tx = n1*asx + n2*asy + n3*asz;
- ty = n1*bsx + n2*bsy + n3*bsz;
- tz = n1*csx + n2*csy + n3*csz;
+ tx = n1*asx + n2*bsx + n3*csx;
+ ty = n1*asy + n2*bsy + n3*csy;
+ tz = n1*asz + n2*bsz + n3*csz;
tlen = modulus(tx, ty, tz);
/* Test modulus for agreement with moduli of template */
for ( i=0; i<3; i++ ) {
+
if ( !within_tolerance(lengths[i], tlen, ltl) )
continue;
@@ -510,6 +511,7 @@ UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose)
/* Angle between axes 0 and 1 should be angle 2 */
if ( fabs(ang - angles[2]) > angtol ) continue;
+
fom1 = fabs(ang - angles[2]);
for ( k=0; k<ncand[2]; k++ ) {
@@ -526,6 +528,7 @@ UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose)
/* ... it should be angle 1 ... */
if ( fabs(ang - angles[1]) > angtol ) continue;
+
fom2 = fom1 + fabs(ang - angles[1]);
/* Finally, the angle between the current candidate for
@@ -536,6 +539,7 @@ UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose)
/* ... it should be angle 0 ... */
if ( fabs(ang - angles[0]) > angtol ) continue;
+
fom3 = fom2 + fabs(ang - angles[0]);
fom3 += LWEIGHT * (cand[0][i].fom + cand[1][j].fom
+ cand[2][k].fom);