diff options
Diffstat (limited to 'libcrystfel')
-rw-r--r-- | libcrystfel/src/cell-utils.c | 33 | ||||
-rw-r--r-- | libcrystfel/src/rational.c | 23 | ||||
-rw-r--r-- | libcrystfel/src/rational.h | 2 |
3 files changed, 46 insertions, 12 deletions
diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 6e7cebb5..dbc845bc 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -1844,6 +1844,23 @@ static void rtnl_mult_in_place(RationalMatrix *T, RationalMatrix *M) } +/* Replace the elements of T with those of T*M */ +static void rtnl_int_mult_in_place(RationalMatrix *T, IntegerMatrix *M) +{ + int i, j; + RationalMatrix *tmp; + + tmp = rtnl_mtx_new(3, 3); + rtnl_mtx_intmatmult(T, M, tmp); + for ( i=0; i<3; i++ ) { + for ( j=0; j<3; j++ ) { + rtnl_mtx_set(T, i, j, rtnl_mtx_get(tmp, i, j)); + } + } + rtnl_mtx_free(tmp); +} + + /* Cell volume from G6 components */ static double g6_volume(struct g6 g) { @@ -2243,24 +2260,16 @@ UnitCell *compare_reindexed_cell_parameters(UnitCell *cell_in, //intmat_print(P); /* Calculate combined matrix: CB.RiB.RA.CiA */ - RationalMatrix *rRA = rtnl_mtx_from_intmat(RA); - RationalMatrix *rP = rtnl_mtx_from_intmat(P); IntegerMatrix *RiB = intmat_inverse(RB); - RationalMatrix *rRiB = rtnl_mtx_from_intmat(RiB); - RationalMatrix *rCB = rtnl_mtx_from_intmat(CB); RationalMatrix *comb = rtnl_mtx_identity(3); rtnl_mult_in_place(comb, CiA); - rtnl_mult_in_place(comb, rRA); - rtnl_mult_in_place(comb, rP); - rtnl_mult_in_place(comb, rRiB); - rtnl_mult_in_place(comb, rCB); + rtnl_int_mult_in_place(comb, RA); + rtnl_int_mult_in_place(comb, P); + rtnl_int_mult_in_place(comb, RiB); + rtnl_int_mult_in_place(comb, CB); - rtnl_mtx_free(rRA); - rtnl_mtx_free(rP); intmat_free(RiB); - rtnl_mtx_free(rRiB); - rtnl_mtx_free(rCB); match = cell_transform_rational(cell_in, comb); //STATUS("Original cell transformed to look like reference:\n"); diff --git a/libcrystfel/src/rational.c b/libcrystfel/src/rational.c index ce286d7a..f8d3cb68 100644 --- a/libcrystfel/src/rational.c +++ b/libcrystfel/src/rational.c @@ -553,6 +553,29 @@ void rtnl_mtx_mtxmult(const RationalMatrix *A, const RationalMatrix *B, } +void rtnl_mtx_intmatmult(const RationalMatrix *A, const IntegerMatrix *B, + RationalMatrix *ans) +{ + int i, j; + + assert(ans->rows == A->rows); + + for ( i=0; i<ans->rows; i++ ) { + for ( j=0; j<ans->cols; j++ ) { + int k; + Rational sum = rtnl_zero(); + for ( k=0; k<A->rows; k++ ) { + Rational add; + add = rtnl_mul(rtnl_mtx_get(A, i, k), + rtnl(intmat_get(B, k, j), 1)); + sum = rtnl_add(sum, add); + } + rtnl_mtx_set(ans, i, j, sum); + } + } +} + + /* Given a "P-matrix" (see ITA chapter 5.1), calculate the fractional * coordinates of point "vec" in the original axes, given its fractional * coordinates in the transformed axes. diff --git a/libcrystfel/src/rational.h b/libcrystfel/src/rational.h index 8681bb06..012b929e 100644 --- a/libcrystfel/src/rational.h +++ b/libcrystfel/src/rational.h @@ -92,6 +92,8 @@ extern IntegerMatrix *intmat_from_rtnl_mtx(const RationalMatrix *m); extern void rtnl_mtx_free(RationalMatrix *mtx); extern void rtnl_mtx_mtxmult(const RationalMatrix *A, const RationalMatrix *B, RationalMatrix *ans); +extern void rtnl_mtx_intmatmult(const RationalMatrix *A, const IntegerMatrix *B, + RationalMatrix *ans); extern int transform_fractional_coords_rtnl(const RationalMatrix *P, const Rational *ivec, Rational *ans); |