aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel/src/cell-utils.c
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2017-03-06 17:14:54 +0100
committerThomas White <taw@physics.org>2017-09-07 16:29:54 +0200
commit141fb2cd3af6ea14098cc65ba3c64f93f72c9cf3 (patch)
treefbfbbb0ffe94dcf6d4667b7d8f6ac7ca8f99cd37 /libcrystfel/src/cell-utils.c
parentd2d7ecbebca7b5a063e8a8487f33feb2b235e462 (diff)
Add compare_cells() (and use it in whirligig)
Diffstat (limited to 'libcrystfel/src/cell-utils.c')
-rw-r--r--libcrystfel/src/cell-utils.c121
1 files changed, 119 insertions, 2 deletions
diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c
index 6610abc3..41436bb4 100644
--- a/libcrystfel/src/cell-utils.c
+++ b/libcrystfel/src/cell-utils.c
@@ -3,12 +3,12 @@
*
* Unit Cell utility functions
*
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
* Copyright © 2012 Lorenzo Galli
*
* Authors:
- * 2009-2012,2014-2015 Thomas White <taw@physics.org>
+ * 2009-2012,2014-2017 Thomas White <taw@physics.org>
* 2012 Lorenzo Galli
*
* This file is part of CrystFEL.
@@ -1686,3 +1686,120 @@ double cell_get_volume(UnitCell *cell)
return 1/rec_volume;
}
+
+
+static double moduli_check(double ax, double ay, double az,
+ double bx, double by, double bz)
+{
+ double ma = modulus(ax, ay, az);
+ double mb = modulus(bx, by, bz);
+ return fabs(ma-mb)/ma;
+}
+
+
+static int cells_are_similar(UnitCell *cell1, UnitCell *cell2,
+ const double ltl, const double atl)
+{
+ double asx1, asy1, asz1, bsx1, bsy1, bsz1, csx1, csy1, csz1;
+ double asx2, asy2, asz2, bsx2, bsy2, bsz2, csx2, csy2, csz2;
+ UnitCell *pcell1, *pcell2;
+
+ /* Compare primitive cells, not centered */
+ pcell1 = uncenter_cell(cell1, NULL);
+ pcell2 = uncenter_cell(cell2, NULL);
+
+ cell_get_reciprocal(pcell1, &asx1, &asy1, &asz1,
+ &bsx1, &bsy1, &bsz1,
+ &csx1, &csy1, &csz1);
+
+ cell_get_reciprocal(pcell2, &asx2, &asy2, &asz2,
+ &bsx2, &bsy2, &bsz2,
+ &csx2, &csy2, &csz2);
+
+
+ cell_free(pcell1);
+ cell_free(pcell2);
+
+ if ( angle_between(asx1, asy1, asz1, asx2, asy2, asz2) > atl ) return 0;
+ if ( angle_between(bsx1, bsy1, bsz1, bsx2, bsy2, bsz2) > atl ) return 0;
+ if ( angle_between(csx1, csy1, csz1, csx2, csy2, csz2) > atl ) return 0;
+
+ if ( moduli_check(asx1, asy1, asz1, asx2, asy2, asz2) > ltl ) return 0;
+ if ( moduli_check(bsx1, bsy1, bsz1, bsx2, bsy2, bsz2) > ltl ) return 0;
+ if ( moduli_check(csx1, csy1, csz1, csx2, csy2, csz2) > ltl ) return 0;
+
+ return 1;
+}
+
+
+
+/**
+ * compare_cells:
+ * @a: A %UnitCell
+ * @b: Another %UnitCell
+ * @ltl: Maximum allowable fractional difference in reciprocal axis lengths
+ * @atl: Maximum allowable difference in reciprocal angles (in radians)
+ * @pmb: Place to store pointer to matrix
+ *
+ * Compare the two units cells. If they agree to within @ltl and @atl, using
+ * any change of axes, returns non-zero and stores the transformation to map @b
+ * onto @a.
+ *
+ * Returns: non-zero if the cells match.
+ *
+ */
+int compare_cells(UnitCell *a, UnitCell *b, double ltl, double atl,
+ IntegerMatrix **pmb)
+{
+ IntegerMatrix *m;
+ int i[9];
+
+ m = intmat_new(3, 3);
+
+ for ( i[0]=-1; i[0]<=+1; i[0]++ ) {
+ for ( i[1]=-1; i[1]<=+1; i[1]++ ) {
+ for ( i[2]=-1; i[2]<=+1; i[2]++ ) {
+ for ( i[3]=-1; i[3]<=+1; i[3]++ ) {
+ for ( i[4]=-1; i[4]<=+1; i[4]++ ) {
+ for ( i[5]=-1; i[5]<=+1; i[5]++ ) {
+ for ( i[6]=-1; i[6]<=+1; i[6]++ ) {
+ for ( i[7]=-1; i[7]<=+1; i[7]++ ) {
+ for ( i[8]=-1; i[8]<=+1; i[8]++ ) {
+
+ UnitCellTransformation *tfn;
+ UnitCell *nc;
+ int j, k;
+ int l = 0;
+
+ for ( j=0; j<3; j++ )
+ for ( k=0; k<3; k++ )
+ intmat_set(m, j, k, i[l++]);
+
+ if ( intmat_det(m) != +1 ) continue;
+
+ tfn = tfn_from_intmat(m);
+ nc = cell_transform(b, tfn);
+
+ if ( cells_are_similar(a, nc, ltl, atl) ) {
+ if ( pmb != NULL ) *pmb = m;
+ tfn_free(tfn);
+ cell_free(nc);
+ return 1;
+ }
+
+ tfn_free(tfn);
+ cell_free(nc);
+
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ intmat_free(m);
+ return 0;
+}