From 74fe33b767d9c94b5c22e274a7518937e8a7adf4 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 13 Mar 2019 15:41:38 +0100 Subject: Rename m to P, to match documentation (including ITA) --- libcrystfel/src/cell.c | 40 ++++++++++++++++++++-------------------- libcrystfel/src/rational.c | 44 +++++++++++++++++++++++++++----------------- 2 files changed, 47 insertions(+), 37 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index 63b84fc3..1dcb6bfb 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -729,14 +729,14 @@ static void maybe_eliminate(CenteringMask c, CenteringMask *cmask, Rational *nc, /* Check if the point x,y,z in the original cell matches any lattice point * in the transformed cell */ -static void check_point_fwd(RationalMatrix *m, CenteringMask *cmask, +static void check_point_fwd(RationalMatrix *P, CenteringMask *cmask, Rational x, Rational y, Rational z) { Rational c[3] = {x, y, z}; Rational nc[3]; /* Transform the lattice point */ - transform_fractional_coords_rtnl(m, c, nc); + transform_fractional_coords_rtnl(P, c, nc); /* Eliminate any centerings which don't include the transformed point */ maybe_eliminate(CMASK_P, cmask, nc, 'P'); @@ -752,14 +752,14 @@ static void check_point_fwd(RationalMatrix *m, CenteringMask *cmask, /* Check if the point x,y,z in the transformed cell matches any lattice point * in the original cell. If not, eliminate "exclude" from "*mask". */ -static void check_point_bwd(RationalMatrix *m, CenteringMask *mask, +static void check_point_bwd(RationalMatrix *P, CenteringMask *mask, char cen, CenteringMask exclude, Rational x, Rational y, Rational z) { Rational nc[3]; Rational c[3] = {x, y, z}; - transform_fractional_coords_rtnl_inverse(m, c, nc); + transform_fractional_coords_rtnl_inverse(P, c, nc); if ( !centering_has_point(cen, nc) ) { *mask |= exclude; @@ -788,19 +788,19 @@ static char cmask_decode(CenteringMask mask) } -static char determine_centering(RationalMatrix *m, char cen) +static char determine_centering(RationalMatrix *P, char cen) { CenteringMask cmask = CMASK_ALL; /* Check whether the current centering can provide all the lattice * points for the transformed cell. Eliminate any centerings for which * it can't. */ - check_point_bwd(m, &cmask, cen, CMASK_A | CMASK_F, rtnl_zero(), rtnl(1,2), rtnl(1,2)); - check_point_bwd(m, &cmask, cen, CMASK_B | CMASK_F, rtnl(1,2), rtnl_zero(), rtnl(1,2)); - check_point_bwd(m, &cmask, cen, CMASK_C | CMASK_F, rtnl(1,2), rtnl(1,2), rtnl_zero()); - check_point_bwd(m, &cmask, cen, CMASK_I, rtnl(1,2), rtnl(1,2), rtnl(1,2)); - check_point_bwd(m, &cmask, cen, CMASK_H, rtnl(2,3), rtnl(1,3), rtnl(1,3)); - check_point_bwd(m, &cmask, cen, CMASK_H, rtnl(1,3), rtnl(2,3), rtnl(2,3)); + check_point_bwd(P, &cmask, cen, CMASK_A | CMASK_F, rtnl_zero(), rtnl(1,2), rtnl(1,2)); + check_point_bwd(P, &cmask, cen, CMASK_B | CMASK_F, rtnl(1,2), rtnl_zero(), rtnl(1,2)); + check_point_bwd(P, &cmask, cen, CMASK_C | CMASK_F, rtnl(1,2), rtnl(1,2), rtnl_zero()); + check_point_bwd(P, &cmask, cen, CMASK_I, rtnl(1,2), rtnl(1,2), rtnl(1,2)); + check_point_bwd(P, &cmask, cen, CMASK_H, rtnl(2,3), rtnl(1,3), rtnl(1,3)); + check_point_bwd(P, &cmask, cen, CMASK_H, rtnl(1,3), rtnl(2,3), rtnl(2,3)); /* Check whether the current centering's lattice points will all * coincide with lattice points in the new centering. Eliminate any @@ -812,30 +812,30 @@ static char determine_centering(RationalMatrix *m, char cen) break; case 'A' : - check_point_fwd(m, &cmask, rtnl_zero(), rtnl(1,2), rtnl(1,2)); + check_point_fwd(P, &cmask, rtnl_zero(), rtnl(1,2), rtnl(1,2)); break; case 'B' : - check_point_fwd(m, &cmask, rtnl(1,2), rtnl_zero(), rtnl(1,2)); + check_point_fwd(P, &cmask, rtnl(1,2), rtnl_zero(), rtnl(1,2)); break; case 'C' : - check_point_fwd(m, &cmask, rtnl(1,2), rtnl(1,2), rtnl_zero()); + check_point_fwd(P, &cmask, rtnl(1,2), rtnl(1,2), rtnl_zero()); break; case 'I' : - check_point_fwd(m, &cmask, rtnl(1,2), rtnl(1,2), rtnl(1,2)); + check_point_fwd(P, &cmask, rtnl(1,2), rtnl(1,2), rtnl(1,2)); break; case 'F' : - check_point_fwd(m, &cmask, rtnl_zero(), rtnl(1,2), rtnl(1,2)); - check_point_fwd(m, &cmask, rtnl(1,2), rtnl_zero(), rtnl(1,2)); - check_point_fwd(m, &cmask, rtnl(1,2), rtnl(1,2), rtnl_zero()); + check_point_fwd(P, &cmask, rtnl_zero(), rtnl(1,2), rtnl(1,2)); + check_point_fwd(P, &cmask, rtnl(1,2), rtnl_zero(), rtnl(1,2)); + check_point_fwd(P, &cmask, rtnl(1,2), rtnl(1,2), rtnl_zero()); break; case 'H' : - check_point_fwd(m, &cmask, rtnl(2,3), rtnl(1,3), rtnl(1,3)); - check_point_fwd(m, &cmask, rtnl(1,3), rtnl(2,3), rtnl(2,3)); + check_point_fwd(P, &cmask, rtnl(2,3), rtnl(1,3), rtnl(1,3)); + check_point_fwd(P, &cmask, rtnl(1,3), rtnl(2,3), rtnl(2,3)); break; } diff --git a/libcrystfel/src/rational.c b/libcrystfel/src/rational.c index c1512b8d..65b209ff 100644 --- a/libcrystfel/src/rational.c +++ b/libcrystfel/src/rational.c @@ -409,7 +409,7 @@ void rtnl_mtx_free(RationalMatrix *mtx) /* rtnl_mtx_solve: - * @m: A %RationalMatrix + * @P: A %RationalMatrix * @vec: An array of %Rational * @ans: An array of %Rational in which to store the result * @@ -422,9 +422,13 @@ void rtnl_mtx_free(RationalMatrix *mtx) * of rows in @m must equal the length of @vec, but note that there is no way * for this function to check that this is the case. * + * Given that P is transformation of the unit cell axes (see ITA chapter 5.1), + * this function calculates the fractional coordinates of a point in the + * transformed axes, given its fractional coordinates in the original axes. + * * Returns: non-zero on error **/ -int transform_fractional_coords_rtnl(const RationalMatrix *m, +int transform_fractional_coords_rtnl(const RationalMatrix *P, const Rational *ivec, Rational *ans) { RationalMatrix *cm; @@ -432,28 +436,28 @@ int transform_fractional_coords_rtnl(const RationalMatrix *m, int h, k; int i; - if ( m->rows != m->cols ) return 1; + if ( P->rows != P->cols ) return 1; /* Copy the matrix and vector because the calculation will * be done in-place */ - cm = rtnl_mtx_copy(m); + cm = rtnl_mtx_copy(P); if ( cm == NULL ) return 1; - vec = malloc(m->rows*sizeof(Rational)); + vec = malloc(cm->rows*sizeof(Rational)); if ( vec == NULL ) return 1; - for ( h=0; hrows; h++ ) vec[h] = ivec[h]; + for ( h=0; hrows; h++ ) vec[h] = ivec[h]; /* Gaussian elimination with partial pivoting */ h = 0; k = 0; - while ( h<=m->rows && k<=m->cols ) { + while ( h<=cm->rows && k<=cm->cols ) { int prow = 0; Rational pval = rtnl_zero(); Rational t; /* Find the row with the largest value in column k */ - for ( i=h; irows; i++ ) { + for ( i=h; irows; i++ ) { if ( rtnl_cmp(rtnl_abs(rtnl_mtx_get(cm, i, k)), pval) > 0 ) { pval = rtnl_abs(rtnl_mtx_get(cm, i, k)); prow = i; @@ -466,7 +470,7 @@ int transform_fractional_coords_rtnl(const RationalMatrix *m, } /* Swap 'prow' with row h */ - for ( i=0; icols; i++ ) { + for ( i=0; icols; i++ ) { t = rtnl_mtx_get(cm, h, i); rtnl_mtx_set(cm, h, i, rtnl_mtx_get(cm, prow, i)); rtnl_mtx_set(cm, prow, i, t); @@ -476,7 +480,7 @@ int transform_fractional_coords_rtnl(const RationalMatrix *m, vec[h] = t; /* Divide and subtract rows below */ - for ( i=h+1; irows; i++ ) { + for ( i=h+1; irows; i++ ) { int j; Rational dval; @@ -484,7 +488,7 @@ int transform_fractional_coords_rtnl(const RationalMatrix *m, dval = rtnl_div(rtnl_mtx_get(cm, i, k), rtnl_mtx_get(cm, h, k)); - for ( j=0; jcols; j++ ) { + for ( j=0; jcols; j++ ) { Rational t = rtnl_mtx_get(cm, i, j); Rational p = rtnl_mul(dval, rtnl_mtx_get(cm, h, j)); t = rtnl_sub(t, p); @@ -504,10 +508,10 @@ int transform_fractional_coords_rtnl(const RationalMatrix *m, } /* Back-substitution */ - for ( i=m->rows-1; i>=0; i-- ) { + for ( i=cm->rows-1; i>=0; i-- ) { int j; Rational sum = rtnl_zero(); - for ( j=i+1; jcols; j++ ) { + for ( j=i+1; jcols; j++ ) { Rational av; av = rtnl_mul(rtnl_mtx_get(cm, i, j), ans[j]); sum = rtnl_add(sum, av); @@ -572,17 +576,23 @@ void rtnl_mtx_mtxmult(const RationalMatrix *A, const RationalMatrix *B, } -void transform_fractional_coords_rtnl_inverse(const RationalMatrix *m, +/* 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. + * To transform in the opposite direction, we would multiply by Q = P^-1. + * Therefore, this direction is the "easy way" and needs just a matrix + * multiplication. */ +void transform_fractional_coords_rtnl_inverse(const RationalMatrix *P, const Rational *vec, Rational *ans) { int i, j; - for ( i=0; irows; i++ ) { + for ( i=0; irows; i++ ) { ans[i] = rtnl_zero(); - for ( j=0; jcols; j++ ) { + for ( j=0; jcols; j++ ) { Rational add; - add = rtnl_mul(rtnl_mtx_get(m, i, j), vec[j]); + add = rtnl_mul(rtnl_mtx_get(P, i, j), vec[j]); ans[i] = rtnl_add(ans[i], add); } } -- cgit v1.2.3