Rename m to P, to match documentation (including ITA)

2 files changed, 47 insertions, 37 deletions

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; h<m->rows; h++ ) vec[h] = ivec[h];
+	for ( h=0; h<cm->rows; 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; i<m->rows; i++ ) {
+		for ( i=h; i<cm->rows; 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; i<m->cols; i++ ) {
+		for ( i=0; i<cm->cols; 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; i<m->rows; i++ ) {
+		for ( i=h+1; i<cm->rows; 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; j<m->cols; j++ ) {
+			for ( j=0; j<cm->cols; 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; j<m->cols; j++ ) {
+		for ( j=i+1; j<cm->cols; 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; i<m->rows; i++ ) {
+	for ( i=0; i<P->rows; i++ ) {
 		ans[i] = rtnl_zero();
-		for ( j=0; j<m->cols; j++ ) {
+		for ( j=0; j<P->cols; 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);
 		}
 	}