aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel/src
diff options
context:
space:
mode:
Diffstat (limited to 'libcrystfel/src')
-rw-r--r--libcrystfel/src/cell-utils.c72
-rw-r--r--libcrystfel/src/cell.c60
-rw-r--r--libcrystfel/src/integer_matrix.c15
-rw-r--r--libcrystfel/src/rational.c4
-rw-r--r--libcrystfel/src/symmetry.c6
-rw-r--r--libcrystfel/src/symop.y12
6 files changed, 88 insertions, 81 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]);
}
;