diff options
author | Thomas White <taw@physics.org> | 2019-03-09 11:42:55 +0100 |
---|---|---|
committer | Thomas White <taw@physics.org> | 2019-03-11 16:49:37 +0100 |
commit | cb51ba0e53425ae663db313d730c1df644e4f76b (patch) | |
tree | d380d61eca0f6729b3af97a5065c16c8e0c5d19c | |
parent | 7d5453e30c18a59b27b486298c3589046adbc308 (diff) |
Change matrix notation to match ITA chapter 5.1
-rw-r--r-- | libcrystfel/src/cell-utils.c | 72 | ||||
-rw-r--r-- | libcrystfel/src/cell.c | 60 | ||||
-rw-r--r-- | libcrystfel/src/integer_matrix.c | 15 | ||||
-rw-r--r-- | libcrystfel/src/rational.c | 4 | ||||
-rw-r--r-- | libcrystfel/src/symmetry.c | 6 | ||||
-rw-r--r-- | libcrystfel/src/symop.y | 12 | ||||
-rw-r--r-- | tests/transformation_check.c | 28 |
7 files changed, 102 insertions, 95 deletions
diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index cc17f108..605aa449 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -405,6 +405,10 @@ static IntegerMatrix *centering_transformation(UnitCell *in, ua = cell_get_unique_axis(in); cen = cell_get_centering(in); + /* Write the matrices exactly as they appear in ITA Table 5.1.3.1. + * C is "P", and Ci is "Q=P^-1". Vice-versa if the transformation + * should go the opposite way to what's written in the first column. */ + if ( (cen=='P') || (cen=='R') ) { *new_centering = cen; *new_latt = lt; @@ -451,12 +455,12 @@ static IntegerMatrix *centering_transformation(UnitCell *in, if ( (lt == L_HEXAGONAL) && (cen == 'H') && (ua == 'c') ) { /* Obverse setting */ - C = intmat_create_3x3(1, -1, 0, - 0, 1, -1, - 1, 1, 1); - Ci = create_rtnl_mtx( 2,3, 1,3, 1,3, - -1,3, 1,3, 1,3, - -1,3, -2,3, 1,3); + C = intmat_create_3x3( 1, 0, 1, + -1, 1, 1, + 0, -1, 1); + Ci = create_rtnl_mtx( 2,3, -1,3, -1,3, + 1,3, 1,3, -2,3, + 1,3, 1,3, 1,3); assert(lt == L_HEXAGONAL); assert(ua == 'c'); *new_latt = L_RHOMBOHEDRAL; @@ -465,12 +469,12 @@ static IntegerMatrix *centering_transformation(UnitCell *in, } if ( cen == 'A' ) { - C = intmat_create_3x3(0, 0, -1, - 1, 1, 0, - 1, -1, 0); - Ci = create_rtnl_mtx( 0,1, 1,2, 1,2, + C = intmat_create_3x3( 1, 0, 0, + 0, 1, 1, + 0, -1, 1); + Ci = create_rtnl_mtx( 1,1, 0,1, 0,1, 0,1, 1,2, -1,2, - -1,1, 0,1, 0,1); + 0,1, 1,2, 1,2); if ( lt == L_ORTHORHOMBIC ) { *new_latt = L_MONOCLINIC; *new_centering = 'P'; @@ -483,12 +487,12 @@ static IntegerMatrix *centering_transformation(UnitCell *in, } if ( cen == 'B' ) { - C = intmat_create_3x3(1, 1, 0, - 0, 0, 1, - 1, -1, 0); - Ci = create_rtnl_mtx( 1,2, 0,1, 1,2, - 1,2, 0,1, -1,2, - 0,1, 1,1, 0,1); + C = intmat_create_3x3( 1, 0, 1, + 0, 1, 0, + -1, 0, 1); + Ci = create_rtnl_mtx( 1,2, 0,1, -1,2, + 0,1, 1,1, 0,1, + 1,2, 0,1, 1,2); if ( lt == L_ORTHORHOMBIC ) { *new_latt = L_MONOCLINIC; *new_centering = 'P'; @@ -501,11 +505,11 @@ static IntegerMatrix *centering_transformation(UnitCell *in, } if ( cen == 'C' ) { - C = intmat_create_3x3(1, -1, 0, - 1, 1, 0, - 0, 0, 1); - Ci = create_rtnl_mtx( 1,2, 1,2, 0,1, - -1,2, 1,2, 0,1, + C = intmat_create_3x3( 1, 1, 0, + -1, 1, 0, + 0, 0, 1); + Ci = create_rtnl_mtx( 1,2, -1,2, 0,1, + 1,2, 1,2, 0,1, 0,1, 0,1, 1,1); if ( lt == L_ORTHORHOMBIC ) { *new_latt = L_MONOCLINIC; @@ -2047,13 +2051,13 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, /* Form the matrix using the first candidate for c */ rtnl_mtx_set(m, 0, 0, cand_a[3*ia+0]); - rtnl_mtx_set(m, 0, 1, cand_a[3*ia+1]); - rtnl_mtx_set(m, 0, 2, cand_a[3*ia+2]); - rtnl_mtx_set(m, 1, 0, cand_b[3*ib+0]); + rtnl_mtx_set(m, 1, 0, cand_a[3*ia+1]); + rtnl_mtx_set(m, 2, 0, cand_a[3*ia+2]); + rtnl_mtx_set(m, 0, 1, cand_b[3*ib+0]); rtnl_mtx_set(m, 1, 1, cand_b[3*ib+1]); - rtnl_mtx_set(m, 1, 2, cand_b[3*ib+2]); - rtnl_mtx_set(m, 2, 0, cand_c[3*ic+0]); - rtnl_mtx_set(m, 2, 1, cand_c[3*ic+1]); + rtnl_mtx_set(m, 2, 1, cand_b[3*ib+2]); + rtnl_mtx_set(m, 0, 2, cand_c[3*ic+0]); + rtnl_mtx_set(m, 1, 2, cand_c[3*ic+1]); rtnl_mtx_set(m, 2, 2, cand_c[3*ic+2]); /* Check angle between a and b */ @@ -2066,13 +2070,13 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in, for ( ic=0; ic<ncand_c; ic++ ) { rtnl_mtx_set(m, 0, 0, cand_a[3*ia+0]); - rtnl_mtx_set(m, 0, 1, cand_a[3*ia+1]); - rtnl_mtx_set(m, 0, 2, cand_a[3*ia+2]); - rtnl_mtx_set(m, 1, 0, cand_b[3*ib+0]); + rtnl_mtx_set(m, 1, 0, cand_a[3*ia+1]); + rtnl_mtx_set(m, 2, 0, cand_a[3*ia+2]); + rtnl_mtx_set(m, 0, 1, cand_b[3*ib+0]); rtnl_mtx_set(m, 1, 1, cand_b[3*ib+1]); - rtnl_mtx_set(m, 1, 2, cand_b[3*ib+2]); - rtnl_mtx_set(m, 2, 0, cand_c[3*ic+0]); - rtnl_mtx_set(m, 2, 1, cand_c[3*ic+1]); + rtnl_mtx_set(m, 2, 1, cand_b[3*ib+2]); + rtnl_mtx_set(m, 0, 2, cand_c[3*ic+0]); + rtnl_mtx_set(m, 1, 2, cand_c[3*ic+1]); rtnl_mtx_set(m, 2, 2, cand_c[3*ic+2]); if ( rtnl_cmp(rtnl_mtx_det(m),rtnl_zero()) == 0 ) continue; diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index 130a24b1..5e00ac9c 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -368,13 +368,13 @@ static int cell_invert(double ax, double ay, double az, return 1; } gsl_matrix_set(m, 0, 0, ax); - gsl_matrix_set(m, 0, 1, bx); - gsl_matrix_set(m, 0, 2, cx); gsl_matrix_set(m, 1, 0, ay); - gsl_matrix_set(m, 1, 1, by); - gsl_matrix_set(m, 1, 2, cy); gsl_matrix_set(m, 2, 0, az); + gsl_matrix_set(m, 0, 1, bx); + gsl_matrix_set(m, 1, 1, by); gsl_matrix_set(m, 2, 1, bz); + gsl_matrix_set(m, 0, 2, cx); + gsl_matrix_set(m, 1, 2, cy); gsl_matrix_set(m, 2, 2, cz); /* Invert */ @@ -410,13 +410,13 @@ static int cell_invert(double ax, double ay, double az, gsl_matrix_transpose(inv); *asx = gsl_matrix_get(inv, 0, 0); - *bsx = gsl_matrix_get(inv, 0, 1); - *csx = gsl_matrix_get(inv, 0, 2); *asy = gsl_matrix_get(inv, 1, 0); - *bsy = gsl_matrix_get(inv, 1, 1); - *csy = gsl_matrix_get(inv, 1, 2); *asz = gsl_matrix_get(inv, 2, 0); + *bsx = gsl_matrix_get(inv, 0, 1); + *bsy = gsl_matrix_get(inv, 1, 1); *bsz = gsl_matrix_get(inv, 2, 1); + *csx = gsl_matrix_get(inv, 0, 2); + *csy = gsl_matrix_get(inv, 1, 2); *csz = gsl_matrix_get(inv, 2, 2); gsl_matrix_free(inv); @@ -630,27 +630,27 @@ UnitCell *cell_transform_gsl_direct(UnitCell *in, gsl_matrix *m) c = gsl_matrix_alloc(3, 3); gsl_matrix_set(c, 0, 0, asx); - gsl_matrix_set(c, 0, 1, asy); - gsl_matrix_set(c, 0, 2, asz); - gsl_matrix_set(c, 1, 0, bsx); + gsl_matrix_set(c, 1, 0, asy); + gsl_matrix_set(c, 2, 0, asz); + gsl_matrix_set(c, 0, 1, bsx); gsl_matrix_set(c, 1, 1, bsy); - gsl_matrix_set(c, 1, 2, bsz); - gsl_matrix_set(c, 2, 0, csx); - gsl_matrix_set(c, 2, 1, csy); + gsl_matrix_set(c, 2, 1, bsz); + gsl_matrix_set(c, 0, 2, csx); + gsl_matrix_set(c, 1, 2, csy); gsl_matrix_set(c, 2, 2, csz); res = gsl_matrix_calloc(3, 3); - gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m, c, 0.0, res); + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, c, m, 0.0, res); out = cell_new_from_cell(in); cell_set_cartesian(out, gsl_matrix_get(res, 0, 0), - gsl_matrix_get(res, 0, 1), - gsl_matrix_get(res, 0, 2), gsl_matrix_get(res, 1, 0), - gsl_matrix_get(res, 1, 1), - gsl_matrix_get(res, 1, 2), gsl_matrix_get(res, 2, 0), + gsl_matrix_get(res, 0, 1), + gsl_matrix_get(res, 1, 1), gsl_matrix_get(res, 2, 1), + gsl_matrix_get(res, 0, 2), + gsl_matrix_get(res, 1, 2), gsl_matrix_get(res, 2, 2)); gsl_matrix_free(res); @@ -673,27 +673,27 @@ UnitCell *cell_transform_gsl_reciprocal(UnitCell *in, gsl_matrix *m) c = gsl_matrix_alloc(3, 3); gsl_matrix_set(c, 0, 0, asx); - gsl_matrix_set(c, 0, 1, asy); - gsl_matrix_set(c, 0, 2, asz); - gsl_matrix_set(c, 1, 0, bsx); + gsl_matrix_set(c, 1, 0, asy); + gsl_matrix_set(c, 2, 0, asz); + gsl_matrix_set(c, 0, 1, bsx); gsl_matrix_set(c, 1, 1, bsy); - gsl_matrix_set(c, 1, 2, bsz); - gsl_matrix_set(c, 2, 0, csx); - gsl_matrix_set(c, 2, 1, csy); + gsl_matrix_set(c, 2, 1, bsz); + gsl_matrix_set(c, 0, 2, csx); + gsl_matrix_set(c, 1, 2, csy); gsl_matrix_set(c, 2, 2, csz); res = gsl_matrix_calloc(3, 3); - gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m, c, 0.0, res); + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, c, m, 0.0, res); out = cell_new_from_cell(in); cell_set_reciprocal(out, gsl_matrix_get(res, 0, 0), - gsl_matrix_get(res, 0, 1), - gsl_matrix_get(res, 0, 2), gsl_matrix_get(res, 1, 0), - gsl_matrix_get(res, 1, 1), - gsl_matrix_get(res, 1, 2), gsl_matrix_get(res, 2, 0), + gsl_matrix_get(res, 0, 1), + gsl_matrix_get(res, 1, 1), gsl_matrix_get(res, 2, 1), + gsl_matrix_get(res, 0, 2), + gsl_matrix_get(res, 1, 2), gsl_matrix_get(res, 2, 2)); gsl_matrix_free(res); diff --git a/libcrystfel/src/integer_matrix.c b/libcrystfel/src/integer_matrix.c index 6e2ade6a..1b448771 100644 --- a/libcrystfel/src/integer_matrix.c +++ b/libcrystfel/src/integer_matrix.c @@ -276,23 +276,26 @@ int intmat_solve_rational(const IntegerMatrix *m, const Rational *vec, * number of columns in @m, and the size of the result equals the number of rows * in @m. * + * The multiplication looks like this: + * (a1, a2, a3) = (vec1, vec2, vec3) m + * Therefore matching the notation in ITA chapter 5.1. + * * Returns: a newly allocated array of signed integers containing the answer, * or NULL on error. **/ signed int *intmat_intvec_mult(const IntegerMatrix *m, const signed int *vec) { signed int *ans; - unsigned int i; + unsigned int j; ans = malloc(m->rows * sizeof(signed int)); if ( ans == NULL ) return NULL; - for ( i=0; i<m->rows; i++ ) { - - unsigned int j; + for ( j=0; j<m->cols; j++ ) { - ans[i] = 0; - for ( j=0; j<m->cols; j++ ) { + unsigned int i; + ans[j] = 0; + for ( i=0; i<m->rows; i++ ) { ans[i] += intmat_get(m, i, j) * vec[j]; } diff --git a/libcrystfel/src/rational.c b/libcrystfel/src/rational.c index 20a73bc8..ab1eff91 100644 --- a/libcrystfel/src/rational.c +++ b/libcrystfel/src/rational.c @@ -551,8 +551,8 @@ void rtnl_mtx_mtxmult(const RationalMatrix *A, const RationalMatrix *B, { int i, j; - assert(ans->cols == A->cols); - assert(ans->rows == B->rows); + assert(ans->cols == B->cols); + assert(ans->rows == A->rows); assert(A->cols == B->rows); for ( i=0; i<ans->rows; i++ ) { diff --git a/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c index 9ee80a16..05ee7c6b 100644 --- a/libcrystfel/src/symmetry.c +++ b/libcrystfel/src/symmetry.c @@ -197,9 +197,9 @@ static void add_symop_v(SymOpList *ops, m = intmat_new(3, 3); assert(m != NULL); - for ( i=0; i<3; i++ ) intmat_set(m, 0, i, h[i]); - for ( i=0; i<3; i++ ) intmat_set(m, 1, i, k[i]); - for ( i=0; i<3; i++ ) intmat_set(m, 2, i, l[i]); + for ( i=0; i<3; i++ ) intmat_set(m, i, 0, h[i]); + for ( i=0; i<3; i++ ) intmat_set(m, i, 1, k[i]); + for ( i=0; i<3; i++ ) intmat_set(m, i, 2, l[i]); free(h); free(k); diff --git a/libcrystfel/src/symop.y b/libcrystfel/src/symop.y index bebf823f..be31c21d 100644 --- a/libcrystfel/src/symop.y +++ b/libcrystfel/src/symop.y @@ -101,13 +101,13 @@ symoplist: symop: axexpr COMMA axexpr COMMA axexpr { rtnl_mtx_set(m, 0, 0, $1[0]); - rtnl_mtx_set(m, 0, 1, $1[1]); - rtnl_mtx_set(m, 0, 2, $1[2]); - rtnl_mtx_set(m, 1, 0, $3[0]); + rtnl_mtx_set(m, 1, 0, $1[1]); + rtnl_mtx_set(m, 2, 0, $1[2]); + rtnl_mtx_set(m, 0, 1, $3[0]); rtnl_mtx_set(m, 1, 1, $3[1]); - rtnl_mtx_set(m, 1, 2, $3[2]); - rtnl_mtx_set(m, 2, 0, $5[0]); - rtnl_mtx_set(m, 2, 1, $5[1]); + rtnl_mtx_set(m, 2, 1, $3[2]); + rtnl_mtx_set(m, 0, 2, $5[0]); + rtnl_mtx_set(m, 1, 2, $5[1]); rtnl_mtx_set(m, 2, 2, $5[2]); } ; diff --git a/tests/transformation_check.c b/tests/transformation_check.c index ba189e4e..2790d6f1 100644 --- a/tests/transformation_check.c +++ b/tests/transformation_check.c @@ -341,9 +341,9 @@ int main(int argc, char *argv[]) /* Permutation of axes */ if ( tfn == NULL ) return 1; - intmat_set_all_3x3(tfn, 0,1,0, - 0,0,1, - 1,0,0); + intmat_set_all_3x3(tfn, 0,0,1, + 1,0,0, + 0,1,0); fail += check_transformation(cell, tfn, 1); /* Doubling of cell in one direction */ @@ -355,28 +355,28 @@ int main(int argc, char *argv[]) /* Shearing */ if ( tfn == NULL ) return 1; - intmat_set_all_3x3(tfn, 1,1,0, - 0,1,0, + intmat_set_all_3x3(tfn, 1,0,0, + 1,1,0, 0,0,1); fail += check_transformation(cell, tfn, 1); /* Crazy */ if ( tfn == NULL ) return 1; - intmat_set_all_3x3(tfn, 1,1,0, - 0,1,0, - 0,1,1); + intmat_set_all_3x3(tfn, 1,0,0, + 1,1,1, + 0,0,1); fail += check_transformation(cell, tfn, 0); /* Identity in two parts */ part1 = intmat_identity(3); part2 = intmat_identity(3); if ( tfn == NULL ) return 1; - intmat_set_all_3x3(part1, 0,0,1, - 0,1,0, - -1,0,0); - intmat_set_all_3x3(part2, 0,0,-1, - 0,1,0, - 1,0,0); + intmat_set_all_3x3(part1, 0,0,-1, + 0,1,0, + 1,0,0); + intmat_set_all_3x3(part2, 0,0,1, + 0,1,0, + -1,0,0); tfn = intmat_intmat_mult(part1, part2); fail += check_identity(cell, tfn); intmat_free(part1); |