aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel
diff options
context:
space:
mode:
Diffstat (limited to 'libcrystfel')
-rw-r--r--libcrystfel/CMakeLists.txt12
-rw-r--r--libcrystfel/src/cell-utils.c774
-rw-r--r--libcrystfel/src/cell-utils.h31
-rw-r--r--libcrystfel/src/cell.c600
-rw-r--r--libcrystfel/src/cell.h27
-rw-r--r--libcrystfel/src/index.c6
-rw-r--r--libcrystfel/src/integer_matrix.c92
-rw-r--r--libcrystfel/src/integer_matrix.h21
-rw-r--r--libcrystfel/src/rational.c675
-rw-r--r--libcrystfel/src/rational.h107
-rw-r--r--libcrystfel/src/symmetry.c175
-rw-r--r--libcrystfel/src/symmetry.h9
-rw-r--r--libcrystfel/src/symop.l52
-rw-r--r--libcrystfel/src/symop.y140
-rw-r--r--libcrystfel/src/taketwo.c103
-rw-r--r--libcrystfel/src/xgandalf.c12
16 files changed, 2156 insertions, 680 deletions
diff --git a/libcrystfel/CMakeLists.txt b/libcrystfel/CMakeLists.txt
index 247b925e..014d2a5a 100644
--- a/libcrystfel/CMakeLists.txt
+++ b/libcrystfel/CMakeLists.txt
@@ -7,6 +7,8 @@ find_package(PINKINDEXER)
find_package(NBP)
find_package(FDIP)
find_package(ZLIB REQUIRED)
+find_package(FLEX REQUIRED)
+find_package(BISON REQUIRED)
pkg_search_module(FFTW fftw3)
set(HAVE_CURSES ${CURSES_FOUND})
@@ -30,6 +32,12 @@ endif()
configure_file(config.h.cmake.in config.h)
+bison_target(symopp src/symop.y ${CMAKE_CURRENT_BINARY_DIR}/symop-parse.c COMPILE_FLAGS --report=all)
+flex_target(symopl src/symop.l ${CMAKE_CURRENT_BINARY_DIR}/symop-lex.c
+ DEFINES_FILE ${CMAKE_CURRENT_BINARY_DIR}/symop-lex.h)
+add_flex_bison_dependency(symopl symopp)
+include_directories(${PROJECT_SOURCE_DIR}/src)
+
set(LIBCRYSTFEL_SOURCES
src/reflist.c
src/utils.c
@@ -62,6 +70,9 @@ set(LIBCRYSTFEL_SOURCES
src/peakfinder8.c
src/taketwo.c
src/xgandalf.c
+ src/rational.c
+ ${BISON_symopp_OUTPUTS}
+ ${FLEX_symopl_OUTPUTS}
)
if (HAVE_FFTW)
@@ -101,6 +112,7 @@ set(LIBCRYSTFEL_HEADERS
src/peakfinder8.h
src/taketwo.h
src/xgandalf.h
+ src/rational.h
)
add_library(${PROJECT_NAME} SHARED
diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c
index 7b1984bb..a36ebe1c 100644
--- a/libcrystfel/src/cell-utils.c
+++ b/libcrystfel/src/cell-utils.c
@@ -3,13 +3,13 @@
*
* Unit Cell utility functions
*
- * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2019 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
* Copyright © 2012 Lorenzo Galli
*
* Authors:
- * 2009-2012,2014-2017 Thomas White <taw@physics.org>
- * 2012 Lorenzo Galli
+ * 2009-2019 Thomas White <taw@physics.org>
+ * 2012 Lorenzo Galli
*
* This file is part of CrystFEL.
*
@@ -219,11 +219,6 @@ int right_handed(UnitCell *cell)
void cell_print(UnitCell *cell)
{
- double asx, asy, asz;
- double bsx, bsy, bsz;
- double csx, csy, csz;
- double a, b, c, alpha, beta, gamma;
- double ax, ay, az, bx, by, bz, cx, cy, cz;
LatticeType lt;
char cen;
@@ -251,12 +246,31 @@ void cell_print(UnitCell *cell)
}
if ( cell_has_parameters(cell) ) {
+
+ double a, b, c, alpha, beta, gamma;
cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma);
STATUS("a b c alpha beta gamma\n");
STATUS("%6.2f %6.2f %6.2f A %6.2f %6.2f %6.2f deg\n",
a*1e10, b*1e10, c*1e10,
rad2deg(alpha), rad2deg(beta), rad2deg(gamma));
+ } else {
+ STATUS("Unit cell parameters are not specified.\n");
+ }
+}
+
+
+void cell_print_full(UnitCell *cell)
+{
+
+ cell_print(cell);
+
+ if ( cell_has_parameters(cell) ) {
+
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+ double ax, ay, az, bx, by, bz, cx, cy, cz;
cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
@@ -283,8 +297,6 @@ void cell_print(UnitCell *cell)
STATUS("Cell representation is %s.\n", cell_rep(cell));
- } else {
- STATUS("Unit cell parameters are not specified.\n");
}
}
@@ -349,203 +361,218 @@ int bravais_lattice(UnitCell *cell)
}
-static UnitCellTransformation *uncentering_transformation(UnitCell *in,
- char *new_centering,
- LatticeType *new_latt)
+static RationalMatrix *create_rtnl_mtx(signed int a1, signed int a2,
+ signed int b1, signed int b2,
+ signed int c1, signed int c2,
+ signed int d1, signed int d2,
+ signed int e1, signed int e2,
+ signed int f1, signed int f2,
+ signed int g1, signed int g2,
+ signed int h1, signed int h2,
+ signed int i1, signed int i2)
+{
+ RationalMatrix *m = rtnl_mtx_new(3, 3);
+ if ( m == NULL ) return NULL;
+ rtnl_mtx_set(m, 0, 0, rtnl(a1, a2));
+ rtnl_mtx_set(m, 0, 1, rtnl(b1, b2));
+ rtnl_mtx_set(m, 0, 2, rtnl(c1, c2));
+ rtnl_mtx_set(m, 1, 0, rtnl(d1, d2));
+ rtnl_mtx_set(m, 1, 1, rtnl(e1, e2));
+ rtnl_mtx_set(m, 1, 2, rtnl(f1, f2));
+ rtnl_mtx_set(m, 2, 0, rtnl(g1, g2));
+ rtnl_mtx_set(m, 2, 1, rtnl(h1, h2));
+ rtnl_mtx_set(m, 2, 2, rtnl(i1, i2));
+ return m;
+}
+
+
+/* Given a centered cell @in, return the integer transformation matrix which
+ * turns a primitive cell into @in. Set new_centering and new_latt to the
+ * centering and lattice type of the primitive cell (usually aP, sometimes rR,
+ * rarely mP). Store the inverse matrix at pCi */
+static IntegerMatrix *centering_transformation(UnitCell *in,
+ char *new_centering,
+ LatticeType *new_latt,
+ char *new_ua,
+ RationalMatrix **pCi)
{
- UnitCellTransformation *t;
- const double OT = 1.0/3.0;
- const double TT = 2.0/3.0;
- const double H = 0.5;
LatticeType lt;
char ua, cen;
+ IntegerMatrix *C = NULL;
+ RationalMatrix *Ci = NULL;
lt = cell_get_lattice_type(in);
ua = cell_get_unique_axis(in);
cen = cell_get_centering(in);
- t = tfn_identity();
- if ( t == NULL ) return NULL;
-
- if ( ua == 'a' ) {
- tfn_combine(t, tfn_vector(0,1,0),
- tfn_vector(0,0,1),
- tfn_vector(1,0,0));
- if ( lt == L_MONOCLINIC ) {
- assert(cen != 'A');
- switch ( cen ) {
- case 'B' : cen = 'A'; break;
- case 'C' : cen = 'B'; break;
- case 'I' : cen = 'I'; break;
- }
- }
- }
-
- if ( ua == 'b' ) {
- tfn_combine(t, tfn_vector(0,0,1),
- tfn_vector(1,0,0),
- tfn_vector(0,1,0));
- if ( lt == L_MONOCLINIC ) {
- assert(cen != 'B');
- switch ( cen ) {
- case 'C' : cen = 'A'; break;
- case 'A' : cen = 'B'; break;
- case 'I' : cen = 'I'; break;
- }
- }
- }
-
- switch ( cen ) {
+ /* 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. */
- case 'P' :
+ if ( (cen=='P') || (cen=='R') ) {
+ *new_centering = cen;
*new_latt = lt;
- *new_centering = 'P';
- break;
-
- case 'R' :
- *new_latt = L_RHOMBOHEDRAL;
- *new_centering = 'R';
- break;
+ *new_ua = ua;
+ C = intmat_identity(3);
+ Ci = rtnl_mtx_identity(3);
+ }
- case 'I' :
- tfn_combine(t, tfn_vector(-H,H,H),
- tfn_vector(H,-H,H),
- tfn_vector(H,H,-H));
+ if ( cen == 'I' ) {
+ C = intmat_create_3x3(0, 1, 1,
+ 1, 0, 1,
+ 1, 1, 0);
+ Ci = create_rtnl_mtx(-1,2, 1,2, 1,2,
+ 1,2, -1,2, 1,2,
+ 1,2, 1,2, -1,2);
if ( lt == L_CUBIC ) {
*new_latt = L_RHOMBOHEDRAL;
*new_centering = 'R';
+ *new_ua = '*';
} else {
- /* Tetragonal or orthorhombic */
*new_latt = L_TRICLINIC;
*new_centering = 'P';
+ *new_ua = '*';
}
- break;
+ }
- case 'F' :
- tfn_combine(t, tfn_vector(0,H,H),
- tfn_vector(H,0,H),
- tfn_vector(H,H,0));
+ if ( cen == 'F' ) {
+ C = intmat_create_3x3(-1, 1, 1,
+ 1, -1, 1,
+ 1, 1, -1);
+ Ci = create_rtnl_mtx( 0,1, 1,2, 1,2,
+ 1,2, 0,1, 1,2,
+ 1,2, 1,2, 0,1);
if ( lt == L_CUBIC ) {
*new_latt = L_RHOMBOHEDRAL;
*new_centering = 'R';
+ *new_ua = '*';
} else {
- assert(lt == L_ORTHORHOMBIC);
*new_latt = L_TRICLINIC;
*new_centering = 'P';
+ *new_ua = '*';
}
- break;
+ }
- case 'A' :
- tfn_combine(t, tfn_vector( 1, 0, 0),
- tfn_vector( 0, H, H),
- tfn_vector( 0,-H, H));
+ if ( (lt == L_HEXAGONAL) && (cen == 'H') && (ua == 'c') ) {
+ /* Obverse setting */
+ 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;
+ *new_centering = 'R';
+ *new_ua = '*';
+ }
+
+ if ( cen == 'A' ) {
+ 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,
+ 0,1, 1,2, 1,2);
if ( lt == L_ORTHORHOMBIC ) {
*new_latt = L_MONOCLINIC;
+ *new_centering = 'P';
+ *new_ua = 'a';
} else {
*new_latt = L_TRICLINIC;
+ *new_centering = 'P';
+ *new_ua = '*';
}
- *new_centering = 'P';
- break;
+ }
- case 'B' :
- tfn_combine(t, tfn_vector( H, 0, H),
- tfn_vector( 0, 1, 0),
- tfn_vector(-H, 0, H));
+ if ( cen == 'B' ) {
+ 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';
+ *new_ua = 'b';
} else {
*new_latt = L_TRICLINIC;
+ *new_centering = 'P';
+ *new_ua = '*';
}
- *new_centering = 'P';
- break;
+ }
- case 'C' :
- tfn_combine(t, tfn_vector( H, H, 0),
- tfn_vector(-H, H, 0),
- tfn_vector( 0, 0, 1));
+ 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,
+ 0,1, 0,1, 1,1);
if ( lt == L_ORTHORHOMBIC ) {
*new_latt = L_MONOCLINIC;
+ *new_centering = 'P';
+ *new_ua = 'c';
} else {
*new_latt = L_TRICLINIC;
- }
- *new_centering = 'P';
- break;
-
- case 'H' :
- /* Obverse setting */
- tfn_combine(t, tfn_vector(TT,OT,OT),
- tfn_vector(-OT,OT,OT),
- tfn_vector(-OT,-TT,OT));
- assert(lt == L_HEXAGONAL);
- *new_latt = L_RHOMBOHEDRAL;
- *new_centering = 'R';
- break;
-
- default :
- ERROR("Invalid centering '%c'\n", cell_get_centering(in));
- return NULL;
-
- }
-
- /* Reverse the axis permutation, but only if this was not an H->R
- * transformation */
- if ( !((cen=='H') && (*new_latt == L_RHOMBOHEDRAL)) ) {
- if ( ua == 'a' ) {
- tfn_combine(t, tfn_vector(0,0,1),
- tfn_vector(1,0,0),
- tfn_vector(0,1,0));
- }
-
- if ( ua == 'b' ) {
- tfn_combine(t, tfn_vector(0,1,0),
- tfn_vector(0,0,1),
- tfn_vector(1,0,0));
+ *new_centering = 'P';
+ *new_ua = '*';
}
}
- return t;
+ *pCi = Ci;
+ return C;
}
/**
* uncenter_cell:
* @in: A %UnitCell
- * @t: Location at which to store the transformation which was used.
+ * @C: Location at which to store the centering transformation
+ * @Ci: Location at which to store the inverse centering transformation
*
- * Turns any cell into a primitive one, e.g. for comparison purposes. The
- * transformation which was used is stored at @t, which can be NULL if the
- * transformation is not required.
+ * Turns any cell into a primitive one, e.g. for comparison purposes.
*
- * Returns: a primitive version of @in in a conventional (unique axis c)
- * setting.
+ * The transformation which was used is stored at @Ci. The centering
+ * transformation, which is the transformation you should apply if you want to
+ * get back the original cell, will be stored at @C. Either or both of these
+ * can be NULL if you don't need that information.
+ *
+ * Returns: a primitive version of @in.
*
*/
-UnitCell *uncenter_cell(UnitCell *in, UnitCellTransformation **t)
+UnitCell *uncenter_cell(UnitCell *in, IntegerMatrix **pC, RationalMatrix **pCi)
{
- UnitCellTransformation *tt;
+ IntegerMatrix *C;
+ RationalMatrix *Ci;
char new_centering;
LatticeType new_latt;
+ char new_ua;
UnitCell *out;
- if ( !bravais_lattice(in) ) {
- ERROR("Cannot uncenter: not a Bravais lattice.\n");
- cell_print(in);
- return NULL;
- }
-
- tt = uncentering_transformation(in, &new_centering, &new_latt);
- if ( tt == NULL ) return NULL;
+ C = centering_transformation(in, &new_centering, &new_latt,
+ &new_ua, &Ci);
+ if ( C == NULL ) return NULL;
- out = cell_transform(in, tt);
+ out = cell_transform_rational(in, Ci);
if ( out == NULL ) return NULL;
cell_set_lattice_type(out, new_latt);
cell_set_centering(out, new_centering);
+ cell_set_unique_axis(out, new_ua);
+
+ if ( pC != NULL ) {
+ *pC = C;
+ } else {
+ intmat_free(C);
+ }
- if ( t != NULL ) {
- *t = tt;
+ if ( pCi != NULL ) {
+ *pCi = Ci;
} else {
- tfn_free(tt);
+ rtnl_mtx_free(Ci);
}
return out;
@@ -609,16 +636,16 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose,
float angtol = deg2rad(tols[3]);
UnitCell *cell;
UnitCell *template;
- UnitCellTransformation *uncentering;
+ IntegerMatrix *centering;
UnitCell *new_cell_trans;
/* "Un-center" the template unit cell to make the comparison easier */
- template = uncenter_cell(template_in, &uncentering);
+ template = uncenter_cell(template_in, &centering, NULL);
if ( template == NULL ) return NULL;
/* The candidate cell is also uncentered, because it might be centered
* if it came from (e.g.) MOSFLM */
- cell = uncenter_cell(cell_in, NULL);
+ cell = uncenter_cell(cell_in, NULL, NULL);
if ( cell == NULL ) return NULL;
if ( cell_get_reciprocal(template, &asx, &asy, &asz,
@@ -627,7 +654,7 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose,
ERROR("Couldn't get reciprocal cell for template.\n");
cell_free(template);
cell_free(cell);
- tfn_free(uncentering);
+ intmat_free(centering);
return NULL;
}
@@ -649,7 +676,7 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose,
ERROR("Couldn't get reciprocal cell.\n");
cell_free(template);
cell_free(cell);
- tfn_free(uncentering);
+ intmat_free(centering);
return NULL;
}
@@ -811,7 +838,7 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose,
/* Reverse the de-centering transformation */
if ( new_cell != NULL ) {
- new_cell_trans = cell_transform_inverse(new_cell, uncentering);
+ new_cell_trans = cell_transform_intmat(new_cell, centering);
cell_free(new_cell);
cell_set_lattice_type(new_cell_trans,
cell_get_lattice_type(template_in));
@@ -821,13 +848,13 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose,
cell_get_unique_axis(template_in));
cell_free(template);
- tfn_free(uncentering);
+ intmat_free(centering);
return new_cell_trans;
} else {
cell_free(template);
- tfn_free(uncentering);
+ intmat_free(centering);
return NULL;
}
}
@@ -850,16 +877,16 @@ UnitCell *match_cell_ab(UnitCell *cell_in, UnitCell *template_in)
int have_real_c;
UnitCell *cell;
UnitCell *template;
- UnitCellTransformation *to_given_cell;
+ IntegerMatrix *to_given_cell;
UnitCell *new_cell;
UnitCell *new_cell_trans;
/* "Un-center" the template unit cell to make the comparison easier */
- template = uncenter_cell(template_in, &to_given_cell);
+ template = uncenter_cell(template_in, &to_given_cell, NULL);
/* The candidate cell is also uncentered, because it might be centered
* if it came from (e.g.) MOSFLM */
- cell = uncenter_cell(cell_in, NULL);
+ cell = uncenter_cell(cell_in, NULL, NULL);
/* Get the lengths to match */
if ( cell_get_cartesian(template, &ax, &ay, &az,
@@ -940,7 +967,7 @@ UnitCell *match_cell_ab(UnitCell *cell_in, UnitCell *template_in)
new_cell = cell_new_from_direct_axes(real_a, real_b, real_c);
/* Reverse the de-centering transformation */
- new_cell_trans = cell_transform_inverse(new_cell, to_given_cell);
+ new_cell_trans = cell_transform_intmat_inverse(new_cell, to_given_cell);
cell_free(new_cell);
cell_set_lattice_type(new_cell, cell_get_lattice_type(template_in));
cell_set_centering(new_cell, cell_get_centering(template_in));
@@ -1443,49 +1470,6 @@ void cell_fudge_gslcblas()
}
-UnitCell *transform_cell_gsl(UnitCell *in, gsl_matrix *m)
-{
- gsl_matrix *c;
- double asx, asy, asz;
- double bsx, bsy, bsz;
- double csx, csy, csz;
- gsl_matrix *res;
- UnitCell *out;
-
- cell_get_reciprocal(in, &asx, &asy, &asz, &bsx, &bsy,
- &bsz, &csx, &csy, &csz);
-
- c = gsl_matrix_alloc(3, 3);
- gsl_matrix_set(c, 0, 0, asx);
- 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, 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);
-
- out = cell_new_from_cell(in);
- cell_set_reciprocal(out, gsl_matrix_get(res, 0, 0),
- gsl_matrix_get(res, 1, 0),
- 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);
- gsl_matrix_free(c);
- return out;
-}
-
-
/**
* rotate_cell:
* @in: A %UnitCell to rotate
@@ -1586,7 +1570,9 @@ int cell_is_sensible(UnitCell *cell)
* lattice is a conventional Bravais lattice.
* Warnings are printied if any of the checks are failed.
*
- * Returns: true if cell is invalid.
+ * Returns: zero if the cell is fine, 1 if it is unconventional but otherwise
+ * OK (e.g. left-handed or not a Bravais lattice), and 2 if there is a serious
+ * problem such as the parameters being physically impossible.
*
*/
int validate_cell(UnitCell *cell)
@@ -1596,7 +1582,7 @@ int validate_cell(UnitCell *cell)
if ( cell_has_parameters(cell) && !cell_is_sensible(cell) ) {
ERROR("WARNING: Unit cell parameters are not sensible.\n");
- err = 1;
+ err = 2;
}
if ( !bravais_lattice(cell) ) {
@@ -1620,7 +1606,7 @@ int validate_cell(UnitCell *cell)
|| ((cen == 'C') && (ua == 'c')) ) {
ERROR("WARNING: A, B or C centering matches unique"
" axis.\n");
- err = 1;
+ err = 2;
}
}
@@ -1646,7 +1632,7 @@ int forbidden_reflection(UnitCell *cell,
cen = cell_get_centering(cell);
/* Reflection conditions here must match the transformation matrices
- * in uncentering_transformation(). tests/centering_check verifies
+ * in centering_transformation(). tests/centering_check verifies
* this (amongst other things). */
if ( cen == 'P' ) return 0;
@@ -1694,6 +1680,48 @@ double cell_get_volume(UnitCell *cell)
}
+/**
+ * compare_cell_parameters:
+ * @a: A %UnitCell
+ * @b: Another %UnitCell
+ * @ltl: Maximum allowable fractional difference in axis lengths
+ * @atl: Maximum allowable difference in reciprocal angles (in radians)
+ *
+ * Compare the two unit cells. If the real space parameters match to within
+ * fractional difference @ltl, and the inter-axial angles match within @atl,
+ * and the centering matches, this function returns 1. Otherwise 0.
+ *
+ * This function considers the cell parameters and centering, but ignores the
+ * orientation of the cell. If you want to compare the orientation as well,
+ * use compare_cell_parameters_and_orientation() instead.
+ *
+ * Returns: non-zero if the cells match.
+ *
+ */
+int compare_cell_parameters(UnitCell *cell1, UnitCell *cell2,
+ float ltl, float atl)
+{
+ double a1, b1, c1, al1, be1, ga1;
+ double a2, b2, c2, al2, be2, ga2;
+
+ /* Centering must match: we don't arbitrarte primitive vs centered,
+ * different cell choices etc */
+ if ( cell_get_centering(cell1) != cell_get_centering(cell2) ) return 0;
+
+ cell_get_parameters(cell1, &a1, &b1, &c1, &al1, &be1, &ga1);
+ cell_get_parameters(cell2, &a2, &b2, &c2, &al2, &be2, &ga2);
+
+ if ( !within_tolerance(a1, a2, ltl*100.0) ) return 0;
+ if ( !within_tolerance(b1, b2, ltl*100.0) ) return 0;
+ if ( !within_tolerance(c1, c2, ltl*100.0) ) return 0;
+ if ( fabs(al1-al2) > atl ) return 0;
+ if ( fabs(be1-be2) > atl ) return 0;
+ if ( fabs(ga1-ga2) > atl ) return 0;
+
+ return 1;
+}
+
+
static double moduli_check(double ax, double ay, double az,
double bx, double by, double bz)
{
@@ -1703,28 +1731,42 @@ static double moduli_check(double ax, double ay, double az,
}
-static int cells_are_similar(UnitCell *cell1, UnitCell *cell2,
- const double ltl, const double atl)
+/**
+ * compare_cell_parameters_and_orientation:
+ * @a: A %UnitCell
+ * @b: Another %UnitCell
+ * @ltl: Maximum allowable fractional difference in reciprocal axis lengths
+ * @atl: Maximum allowable difference in reciprocal angles (in radians)
+ *
+ * Compare the two unit cells. If the axes match in length (to within
+ * fractional difference @ltl) and the axes are aligned to within @atl radians,
+ * this function returns non-zero.
+ *
+ * This function compares the orientation of the cell as well as the parameters.
+ * If you just want to see if the parameters are the same, use
+ * compare_cell_parameters() instead.
+ *
+ * The cells @a and @b must have the same centering. Otherwise, this function
+ * always returns zero.
+ *
+ * Returns: non-zero if the cells match.
+ *
+ */
+int compare_cell_parameters_and_orientation(UnitCell *cell1, UnitCell *cell2,
+ const double ltl, const double atl)
{
double asx1, asy1, asz1, bsx1, bsy1, bsz1, csx1, csy1, csz1;
double asx2, asy2, asz2, bsx2, bsy2, bsz2, csx2, csy2, csz2;
- UnitCell *pcell1, *pcell2;
-
- /* Compare primitive cells, not centered */
- pcell1 = uncenter_cell(cell1, NULL);
- pcell2 = uncenter_cell(cell2, NULL);
-
- cell_get_reciprocal(pcell1, &asx1, &asy1, &asz1,
- &bsx1, &bsy1, &bsz1,
- &csx1, &csy1, &csz1);
- cell_get_reciprocal(pcell2, &asx2, &asy2, &asz2,
- &bsx2, &bsy2, &bsz2,
- &csx2, &csy2, &csz2);
+ if ( cell_get_centering(cell1) != cell_get_centering(cell2) ) return 0;
+ cell_get_cartesian(cell1, &asx1, &asy1, &asz1,
+ &bsx1, &bsy1, &bsz1,
+ &csx1, &csy1, &csz1);
- cell_free(pcell1);
- cell_free(pcell2);
+ cell_get_cartesian(cell2, &asx2, &asy2, &asz2,
+ &bsx2, &bsy2, &bsz2,
+ &csx2, &csy2, &csz2);
if ( angle_between(asx1, asy1, asz1, asx2, asy2, asz2) > atl ) return 0;
if ( angle_between(bsx1, bsy1, bsz1, bsx2, bsy2, bsz2) > atl ) return 0;
@@ -1738,28 +1780,38 @@ static int cells_are_similar(UnitCell *cell1, UnitCell *cell2,
}
-
/**
- * compare_cells:
+ * compare_reindexed_cell_parameters_and_orientation:
* @a: A %UnitCell
* @b: Another %UnitCell
* @ltl: Maximum allowable fractional difference in reciprocal axis lengths
* @atl: Maximum allowable difference in reciprocal angles (in radians)
* @pmb: Place to store pointer to matrix
*
- * Compare the two units cells. If they agree to within @ltl and @atl, using
- * any change of axes, returns non-zero and stores the transformation to map @b
- * onto @a.
+ * Compare the two unit cells. If, using any permutation of the axes, the
+ * axes can be made to match in length (to within fractional difference @ltl)
+ * and the axes aligned to within @atl radians, this function returns non-zero
+ * and stores the transformation to map @b onto @a.
+ *
+ * Note that the orientations of the cells must match, not just the parameters.
+ * The comparison is done after reindexing using
+ * compare_cell_parameters_and_orientation().
+ *
+ * The cells @a and @b must have the same centering. Otherwise, this function
+ * always returns zero.
*
* Returns: non-zero if the cells match.
*
*/
-int compare_cells(UnitCell *a, UnitCell *b, double ltl, double atl,
- IntegerMatrix **pmb)
+int compare_reindexed_cell_parameters_and_orientation(UnitCell *a, UnitCell *b,
+ double ltl, double atl,
+ IntegerMatrix **pmb)
{
IntegerMatrix *m;
int i[9];
+ if ( cell_get_centering(a) != cell_get_centering(b) ) return 0;
+
m = intmat_new(3, 3);
for ( i[0]=-1; i[0]<=+1; i[0]++ ) {
@@ -1772,7 +1824,6 @@ int compare_cells(UnitCell *a, UnitCell *b, double ltl, double atl,
for ( i[7]=-1; i[7]<=+1; i[7]++ ) {
for ( i[8]=-1; i[8]<=+1; i[8]++ ) {
- UnitCellTransformation *tfn;
UnitCell *nc;
int j, k;
int l = 0;
@@ -1783,17 +1834,14 @@ int compare_cells(UnitCell *a, UnitCell *b, double ltl, double atl,
if ( intmat_det(m) != +1 ) continue;
- tfn = tfn_from_intmat(m);
- nc = cell_transform(b, tfn);
+ nc = cell_transform_intmat(b, m);
- if ( cells_are_similar(a, nc, ltl, atl) ) {
+ if ( compare_cell_parameters_and_orientation(a, nc, ltl, atl) ) {
if ( pmb != NULL ) *pmb = m;
- tfn_free(tfn);
cell_free(nc);
return 1;
}
- tfn_free(tfn);
cell_free(nc);
}
@@ -1809,3 +1857,275 @@ int compare_cells(UnitCell *a, UnitCell *b, double ltl, double atl,
intmat_free(m);
return 0;
}
+
+
+struct cand
+{
+ Rational abc[3];
+ double fom;
+};
+
+
+static int cmpcand(const void *av, const void *bv)
+{
+ const struct cand *a = av;
+ const struct cand *b = bv;
+ return a->fom > b->fom;
+}
+
+
+static Rational *find_candidates(double len, double *a, double *b, double *c,
+ double ltl, int *pncand)
+{
+ Rational *r;
+ struct cand *cands;
+ const int max_cand = 1024;
+ int ncand = 0;
+ Rational *rat;
+ int nrat;
+ int nrej = 0;
+ int ia, ib, ic;
+ int i;
+
+ cands = malloc(max_cand * sizeof(struct cand));
+ if ( cands == NULL ) return NULL;
+
+ rat = rtnl_list(-5, 5, 1, 4, &nrat);
+ if ( rat == NULL ) return NULL;
+
+ for ( ia=0; ia<nrat; ia++ ) {
+ for ( ib=0; ib<nrat; ib++ ) {
+ for ( ic=0; ic<nrat; ic++ ) {
+ double vec[3];
+ double abc[3];
+ double veclen;
+ abc[0] = rtnl_as_double(rat[ia]);
+ abc[1] = rtnl_as_double(rat[ib]);
+ abc[2] = rtnl_as_double(rat[ic]);
+ vec[0] = a[0]*abc[0] + b[0]*abc[1] + c[0]*abc[2];
+ vec[1] = a[1]*abc[0] + b[1]*abc[1] + c[1]*abc[2];
+ vec[2] = a[2]*abc[0] + b[2]*abc[1] + c[2]*abc[2];
+ veclen = modulus(vec[0], vec[1], vec[2]);
+ if ( within_tolerance(len, veclen, ltl*100.0) ) {
+ if ( ncand == max_cand ) {
+ nrej++;
+ } else {
+ cands[ncand].abc[0] = rat[ia];
+ cands[ncand].abc[1] = rat[ib];
+ cands[ncand].abc[2] = rat[ic];
+ cands[ncand].fom = fabs(veclen - len);
+ ncand++;
+ }
+ }
+ }
+ }
+ }
+
+ if ( nrej ) {
+ ERROR("WARNING: Too many vector candidates (%i rejected)\n", nrej);
+ }
+
+ /* Sort by difference from reference vector length */
+ qsort(cands, ncand, sizeof(struct cand), cmpcand);
+
+ r = malloc(ncand * 3 * sizeof(Rational));
+ if ( r == 0 ) return NULL;
+
+ for ( i=0; i<ncand; i++ ) {
+ r[3*i+0] = cands[i].abc[0];
+ r[3*i+1] = cands[i].abc[1];
+ r[3*i+2] = cands[i].abc[2];
+ }
+ free(cands);
+
+ *pncand = ncand;
+ return r;
+}
+
+
+static void g6_components(double *g6, double a, double b, double c,
+ double al, double be, double ga)
+{
+ g6[0] = a*a;
+ g6[1] = b*b;
+ g6[2] = c*c;
+ g6[3] = 2.0*b*c*cos(al);
+ g6[4] = 2.0*a*c*cos(be);
+ g6[5] = 2.0*a*b*cos(ga);
+}
+
+
+static double g6_distance(double a1, double b1, double c1,
+ double al1, double be1, double ga1,
+ double a2, double b2, double c2,
+ double al2, double be2, double ga2)
+{
+ double g1[6], g2[6];
+ int i;
+ double total = 0.0;
+ g6_components(g1, a1, b1, c1, al1, be1, ga1);
+ g6_components(g2, a2, b2, c2, al2, be2, ga2);
+ for ( i=0; i<6; i++ ) {
+ total += (g1[i]-g2[i])*(g1[i]-g2[i]);
+ }
+ return sqrt(total);
+}
+
+
+/**
+ * compare_reindexed_cell_parameters:
+ * @cell_in: A %UnitCell
+ * @reference_in: Another %UnitCell
+ * @ltl: Maximum allowable fractional difference in direct-space axis lengths
+ * @atl: Maximum allowable difference in direct-space angles (in radians)
+ * @pmb: Place to store pointer to matrix
+ *
+ * Compare the @cell_in with @reference_in. If @cell is a derivative lattice
+ * of @reference, within fractional axis length difference @ltl and absolute angle
+ * difference @atl (in radians), this function returns non-zero and stores the
+ * transformation which needs to be applied to @cell_in at @pmb.
+ *
+ * Note that the tolerances will be applied to the primitive unit cell. If
+ * the reference cell is centered, a primitive unit cell will first be calculated.
+ *
+ * Subject to the tolerances, this function will find the transformation which
+ * gives the best match to the reference cell, using the Euclidian norm in
+ * G6 [see e.g. Andrews and Bernstein, Acta Cryst. A44 (1988) p1009].
+ *
+ * Only the cell parameters will be compared. The relative orientations are
+ * irrelevant.
+ *
+ * Returns: non-zero if the cells match, zero for no match or error.
+ *
+ */
+int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in,
+ double ltl, double atl,
+ RationalMatrix **pmb)
+{
+ UnitCell *cell;
+ UnitCell *reference;
+ IntegerMatrix *CBint;
+ RationalMatrix *CiA;
+ RationalMatrix *CB;
+ RationalMatrix *m;
+ double a, b, c, al, be, ga;
+ double av[3], bv[3], cv[3];
+ Rational *cand_a;
+ Rational *cand_b;
+ Rational *cand_c;
+ int ncand_a, ncand_b, ncand_c;
+ int ia, ib;
+ RationalMatrix *MCiA;
+ RationalMatrix *CBMCiA;
+ double min_dist = +INFINITY;
+
+ /* Actually compare against primitive version of reference */
+ reference = uncenter_cell(reference_in, &CBint, NULL);
+ if ( reference == NULL ) return 0;
+ CB = rtnl_mtx_from_intmat(CBint);
+ intmat_free(CBint);
+
+ /* Actually compare primitive version of cell */
+ cell = uncenter_cell(cell_in, NULL, &CiA);
+ if ( cell == NULL ) return 0;
+
+ /* Get target parameters */
+ cell_get_parameters(reference, &a, &b, &c, &al, &be, &ga);
+ cell_get_cartesian(cell, &av[0], &av[1], &av[2],
+ &bv[0], &bv[1], &bv[2],
+ &cv[0], &cv[1], &cv[2]);
+
+ /* Find vectors in 'cell' with lengths close to a, b and c */
+ cand_a = find_candidates(a, av, bv, cv, ltl, &ncand_a);
+ cand_b = find_candidates(b, av, bv, cv, ltl, &ncand_b);
+ cand_c = find_candidates(c, av, bv, cv, ltl, &ncand_c);
+
+ m = rtnl_mtx_new(3, 3);
+ MCiA = rtnl_mtx_new(3, 3);
+ CBMCiA = rtnl_mtx_new(3, 3);
+ for ( ia=0; ia<ncand_a; ia++ ) {
+ for ( ib=0; ib<ncand_b; ib++ ) {
+
+ UnitCell *test;
+ double at, bt, ct, alt, bet, gat;
+ double dist;
+ int ic = 0;
+
+ /* Form the matrix using the first candidate for c */
+ rtnl_mtx_set(m, 0, 0, cand_a[3*ia+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, 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 */
+ test = cell_transform_rational(cell, m);
+ cell_get_parameters(test, &at, &bt, &ct, &alt, &bet, &gat);
+ cell_free(test);
+ if ( fabs(gat - ga) > atl ) continue;
+
+ /* Gamma OK, now look for place for c axis */
+ for ( ic=0; ic<ncand_c; ic++ ) {
+
+ rtnl_mtx_set(m, 0, 0, cand_a[3*ia+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, 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;
+
+ test = cell_transform_rational(cell, m);
+ cell_get_parameters(test, &at, &bt, &ct, &alt, &bet, &gat);
+ if ( !right_handed(test) ) {
+ cell_free(test);
+ continue;
+ }
+ if ( fabs(alt - al) > atl ) {
+ cell_free(test);
+ continue;
+ }
+ if ( fabs(bet - be) > atl ) {
+ cell_free(test);
+ continue;
+ }
+
+ dist = g6_distance(at, bt, ct, alt, bet, gat,
+ a, b, c, al, be, ga);
+ if ( dist < min_dist ) {
+ min_dist = dist;
+ rtnl_mtx_mtxmult(m, CiA, MCiA);
+ }
+
+ cell_free(test);
+
+ }
+ }
+ }
+
+ rtnl_mtx_free(m);
+ free(cand_a);
+ free(cand_b);
+ free(cand_c);
+
+ if ( isinf(min_dist) ) {
+ rtnl_mtx_free(CBMCiA);
+ rtnl_mtx_free(MCiA);
+ *pmb = NULL;
+ return 0;
+ }
+
+ /* Solution found */
+ rtnl_mtx_mtxmult(CB, MCiA, CBMCiA);
+ rtnl_mtx_free(MCiA);
+ *pmb = CBMCiA;
+ return 1;
+}
diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h
index 5e2b2825..a97f2bfb 100644
--- a/libcrystfel/src/cell-utils.h
+++ b/libcrystfel/src/cell-utils.h
@@ -3,13 +3,13 @@
*
* Unit Cell utility functions
*
- * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2019 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
* Copyright © 2012 Lorenzo Galli
*
* Authors:
- * 2009-2013,2014,2017 Thomas White <taw@physics.org>
- * 2012 Lorenzo Galli
+ * 2009-2018 Thomas White <taw@physics.org>
+ * 2012 Lorenzo Galli
*
* This file is part of CrystFEL.
*
@@ -49,9 +49,9 @@ extern double resolution(UnitCell *cell,
extern UnitCell *cell_rotate(UnitCell *in, struct quaternion quat);
extern UnitCell *rotate_cell(UnitCell *in, double omega, double phi,
double rot);
-extern UnitCell *transform_cell_gsl(UnitCell *in, gsl_matrix *m);
extern void cell_print(UnitCell *cell);
+extern void cell_print_full(UnitCell *cell);
extern UnitCell *match_cell(UnitCell *cell, UnitCell *tempcell, int verbose,
const float *ltl, int reduce);
@@ -66,7 +66,8 @@ extern int cell_is_sensible(UnitCell *cell);
extern int validate_cell(UnitCell *cell);
-extern UnitCell *uncenter_cell(UnitCell *in, UnitCellTransformation **t);
+extern UnitCell *uncenter_cell(UnitCell *in, IntegerMatrix **pC,
+ RationalMatrix **pCi);
extern int bravais_lattice(UnitCell *cell);
@@ -80,8 +81,24 @@ extern int forbidden_reflection(UnitCell *cell,
extern double cell_get_volume(UnitCell *cell);
-extern int compare_cells(UnitCell *a, UnitCell *b, double ltl, double atl,
- IntegerMatrix **pmb);
+extern int compare_cell_parameters(UnitCell *cell1, UnitCell *cell2,
+ float ltl, float atl);
+
+
+extern int compare_cell_parameters_and_orientation(UnitCell *cell1,
+ UnitCell *cell2,
+ const double ltl,
+ const double atl);
+
+extern int compare_reindexed_cell_parameters_and_orientation(UnitCell *a,
+ UnitCell *b,
+ double ltl,
+ double atl,
+ IntegerMatrix **pmb);
+
+extern int compare_reindexed_cell_parameters(UnitCell *cell, UnitCell *reference,
+ double ltl, double atl,
+ RationalMatrix **pmb);
#ifdef __cplusplus
}
diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c
index cc18b49d..ce9d37bb 100644
--- a/libcrystfel/src/cell.c
+++ b/libcrystfel/src/cell.c
@@ -45,6 +45,8 @@
#include "cell.h"
#include "utils.h"
#include "image.h"
+#include "integer_matrix.h"
+#include "rational.h"
/**
@@ -95,6 +97,20 @@ struct _unitcell {
char unique_axis;
};
+typedef enum {
+ CMASK_P = 1<<0,
+ CMASK_A = 1<<1,
+ CMASK_B = 1<<2,
+ CMASK_C = 1<<3,
+ CMASK_I = 1<<4,
+ CMASK_F = 1<<5,
+ CMASK_H = 1<<6,
+ CMASK_R = 1<<7
+} CenteringMask;
+
+#define CMASK_ALL (CMASK_P | CMASK_A | CMASK_B | CMASK_C | CMASK_I \
+ | CMASK_F | CMASK_H | CMASK_R)
+
/************************** Setters and Constructors **************************/
@@ -352,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 */
@@ -394,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);
@@ -600,333 +616,405 @@ const char *cell_rep(UnitCell *cell)
}
-struct _unitcelltransformation
+UnitCell *cell_transform_gsl_direct(UnitCell *in, gsl_matrix *m)
{
- gsl_matrix *m;
-};
-
-
-/**
- * tfn_inverse:
- * @t: A %UnitCellTransformation.
- *
- * Calculates the inverse of @t. That is, if you apply cell_transform() to a
- * %UnitCell using @t, and then apply cell_transform() to the result using
- * tfn_inverse(@t) instead of @t, you will recover the same lattice vectors
- * (but note that the lattice type, centering and unique axis information will
- * be lost).
- *
- * Returns: The inverse of @t.
- *
- */
-UnitCellTransformation *tfn_inverse(UnitCellTransformation *t)
-{
- int s;
- gsl_matrix *m;
- gsl_matrix *inv;
- gsl_permutation *perm;
- UnitCellTransformation *out;
-
- m = gsl_matrix_alloc(3, 3);
- if ( m == NULL ) return NULL;
-
- out = tfn_identity();
- if ( out == NULL ) {
- gsl_matrix_free(m);
- return NULL;
- }
-
- gsl_matrix_memcpy(m, t->m);
-
- perm = gsl_permutation_alloc(m->size1);
- if ( perm == NULL ) {
- ERROR("Couldn't allocate permutation\n");
- return NULL;
- }
- inv = gsl_matrix_alloc(m->size1, m->size2);
- if ( inv == NULL ) {
- ERROR("Couldn't allocate inverse\n");
- gsl_permutation_free(perm);
- return NULL;
- }
- if ( gsl_linalg_LU_decomp(m, perm, &s) ) {
- ERROR("Couldn't decompose matrix\n");
- gsl_permutation_free(perm);
- return NULL;
- }
- if ( gsl_linalg_LU_invert(m, perm, inv) ) {
- ERROR("Couldn't invert transformation matrix\n");
- gsl_permutation_free(perm);
- return NULL;
- }
- gsl_permutation_free(perm);
+ gsl_matrix *c;
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+ gsl_matrix *res;
+ UnitCell *out;
- gsl_matrix_free(out->m);
- gsl_matrix_free(m);
- out->m = inv;
+ cell_get_cartesian(in, &asx, &asy, &asz, &bsx, &bsy,
+ &bsz, &csx, &csy, &csz);
+
+ c = gsl_matrix_alloc(3, 3);
+ gsl_matrix_set(c, 0, 0, asx);
+ 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, 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, 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, 1, 0),
+ 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);
+ gsl_matrix_free(c);
return out;
}
-/**
- * cell_transform:
- * @cell: A %UnitCell.
- * @t: A %UnitCellTransformation.
- *
- * Applies @t to @cell. Note that the lattice type, centering and unique axis
- * information will not be preserved.
- *
- * Returns: Transformed copy of @cell.
- *
- */
-UnitCell *cell_transform(UnitCell *cell, UnitCellTransformation *t)
-{
- UnitCell *out;
- double ax, ay, az;
- double bx, by, bz;
- double cx, cy, cz;
- gsl_matrix *m;
- gsl_matrix *a;
+static int centering_has_point(char cen, Rational *p)
+{
+ /* First, put the point into the range 0..1 */
+ while ( rtnl_cmp(p[0], rtnl_zero()) < 0 ) p[0] = rtnl_add(p[0], rtnl(1, 1));
+ while ( rtnl_cmp(p[1], rtnl_zero()) < 0 ) p[1] = rtnl_add(p[1], rtnl(1, 1));
+ while ( rtnl_cmp(p[2], rtnl_zero()) < 0 ) p[2] = rtnl_add(p[2], rtnl(1, 1));
+ while ( rtnl_cmp(p[0], rtnl(1, 1)) >= 0 ) p[0] = rtnl_sub(p[0], rtnl(1, 1));
+ while ( rtnl_cmp(p[1], rtnl(1, 1)) >= 0 ) p[1] = rtnl_sub(p[1], rtnl(1, 1));
+ while ( rtnl_cmp(p[2], rtnl(1, 1)) >= 0 ) p[2] = rtnl_sub(p[2], rtnl(1, 1));
+
+ /* 0,0,0 is present in all centerings */
+ if ( (rtnl_cmp(p[0], rtnl_zero()) == 0)
+ && (rtnl_cmp(p[1], rtnl_zero()) == 0)
+ && (rtnl_cmp(p[2], rtnl_zero()) == 0) ) return 1;
+
+ /* Only I has 1/2 , 1/2, 1/2 */
+ if ( (rtnl_cmp(p[0], rtnl(1,2)) == 0)
+ && (rtnl_cmp(p[1], rtnl(1,2)) == 0)
+ && (rtnl_cmp(p[2], rtnl(1,2)) == 0)
+ && (cen == 'I') ) return 1;
+
+ /* A or F has 0 , 1/2, 1/2 */
+ if ( (rtnl_cmp(p[0], rtnl_zero()) == 0)
+ && (rtnl_cmp(p[1], rtnl(1,2)) == 0)
+ && (rtnl_cmp(p[2], rtnl(1,2)) == 0)
+ && ((cen == 'A') || (cen == 'F')) ) return 1;
+
+ /* B or F has 1/2 , 0 , 1/2 */
+ if ( (rtnl_cmp(p[0], rtnl(1,2)) == 0)
+ && (rtnl_cmp(p[1], rtnl_zero()) == 0)
+ && (rtnl_cmp(p[2], rtnl(1,2)) == 0)
+ && ((cen == 'B') || (cen == 'F')) ) return 1;
+
+ /* C or F has 1/2 , 1/2 , 0 */
+ if ( (rtnl_cmp(p[0], rtnl(1,2)) == 0)
+ && (rtnl_cmp(p[1], rtnl(1,2)) == 0)
+ && (rtnl_cmp(p[2], rtnl_zero()) == 0)
+ && ((cen == 'C') || (cen == 'F')) ) return 1;
+
+ /* H has 2/3 , 1/3 , 1/3 */
+ if ( (rtnl_cmp(p[0], rtnl(2,3)) == 0)
+ && (rtnl_cmp(p[1], rtnl(1,3)) == 0)
+ && (rtnl_cmp(p[2], rtnl(1,3)) == 0)
+ && (cen == 'H') ) return 1;
+
+ /* H has 1/3 , 2/3 , 2/3 */
+ if ( (rtnl_cmp(p[0], rtnl(1,3)) == 0)
+ && (rtnl_cmp(p[1], rtnl(2,3)) == 0)
+ && (rtnl_cmp(p[2], rtnl(2,3)) == 0)
+ && (cen == 'H') ) return 1;
- if ( t == NULL ) return NULL;
+ return 0;
+}
- out = cell_new_from_cell(cell);
- if ( out == NULL ) return NULL;
- if ( cell_get_cartesian(out, &ax, &ay, &az,
- &bx, &by, &bz,
- &cx, &cy, &cz) ) return NULL;
+static void maybe_eliminate(CenteringMask c, CenteringMask *cmask, Rational *nc,
+ char cen)
+{
+ /* Skip test if this centering isn't even a candidate */
+ if ( !(*cmask & c) ) return;
- m = gsl_matrix_alloc(3,3);
- a = gsl_matrix_calloc(3,3);
- if ( (m == NULL) || (a == NULL) ) {
- cell_free(out);
- return NULL;
+ if ( !centering_has_point(cen, nc) ) {
+ *cmask |= c;
+ *cmask ^= c;
}
+}
- gsl_matrix_set(m, 0, 0, ax);
- gsl_matrix_set(m, 0, 1, ay);
- gsl_matrix_set(m, 0, 2, az);
- gsl_matrix_set(m, 1, 0, bx);
- gsl_matrix_set(m, 1, 1, by);
- gsl_matrix_set(m, 1, 2, bz);
- gsl_matrix_set(m, 2, 0, cx);
- gsl_matrix_set(m, 2, 1, cy);
- gsl_matrix_set(m, 2, 2, cz);
- gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, t->m, m, 0.0, a);
+/* 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 *P, CenteringMask *cmask,
+ Rational x, Rational y, Rational z)
+{
+ Rational c[3] = {x, y, z};
+ Rational nc[3];
- cell_set_cartesian(out, gsl_matrix_get(a, 0, 0),
- gsl_matrix_get(a, 0, 1),
- gsl_matrix_get(a, 0, 2),
- gsl_matrix_get(a, 1, 0),
- gsl_matrix_get(a, 1, 1),
- gsl_matrix_get(a, 1, 2),
- gsl_matrix_get(a, 2, 0),
- gsl_matrix_get(a, 2, 1),
- gsl_matrix_get(a, 2, 2));
+ /* Transform the lattice point */
+ transform_fractional_coords_rtnl(P, c, nc);
- gsl_matrix_free(a);
- gsl_matrix_free(m);
-
- return out;
+ /* Eliminate any centerings which don't include the transformed point */
+ maybe_eliminate(CMASK_P, cmask, nc, 'P');
+ maybe_eliminate(CMASK_R, cmask, nc, 'R');
+ maybe_eliminate(CMASK_A, cmask, nc, 'A');
+ maybe_eliminate(CMASK_B, cmask, nc, 'B');
+ maybe_eliminate(CMASK_C, cmask, nc, 'C');
+ maybe_eliminate(CMASK_I, cmask, nc, 'I');
+ maybe_eliminate(CMASK_F, cmask, nc, 'F');
+ maybe_eliminate(CMASK_H, cmask, nc, 'H');
}
-/**
- * cell_transform_inverse:
- * @cell: A %UnitCell.
- * @t: A %UnitCellTransformation.
- *
- * Applies the inverse of @t to @cell.
- *
- * Returns: Transformed copy of @cell.
- *
- */
-UnitCell *cell_transform_inverse(UnitCell *cell, UnitCellTransformation *t)
+/* 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 *P, CenteringMask *mask,
+ char cen, CenteringMask exclude,
+ Rational x, Rational y, Rational z)
{
- UnitCellTransformation *inv;
- UnitCell *out;
+ Rational nc[3];
+ Rational c[3] = {x, y, z};
- inv = tfn_inverse(t);
- out = cell_transform(cell, inv);
- tfn_free(inv);
- return out;
+ transform_fractional_coords_rtnl_inverse(P, c, nc);
+
+ if ( !centering_has_point(cen, nc) ) {
+ *mask |= exclude;
+ *mask ^= exclude; /* Unset bits */
+ }
}
-/**
- * tfn_identity:
- *
- * Returns: A %UnitCellTransformation corresponding to an identity operation.
- *
- */
-UnitCellTransformation *tfn_identity()
+static char cmask_decode(CenteringMask mask)
{
- UnitCellTransformation *tfn;
+ char res[32];
- tfn = calloc(1, sizeof(UnitCellTransformation));
- if ( tfn == NULL ) return NULL;
+ res[0] = '\0';
- tfn->m = gsl_matrix_alloc(3, 3);
- if ( tfn->m == NULL ) {
- free(tfn);
- return NULL;
- }
-
- gsl_matrix_set_identity(tfn->m);
+ if ( mask & CMASK_H ) strcat(res, "H");
+ if ( mask & CMASK_F ) strcat(res, "F");
+ if ( mask & CMASK_I ) strcat(res, "I");
+ if ( mask & CMASK_A ) strcat(res, "A");
+ if ( mask & CMASK_B ) strcat(res, "B");
+ if ( mask & CMASK_C ) strcat(res, "C");
+ if ( mask & CMASK_P ) strcat(res, "P");
+ if ( mask & CMASK_R ) strcat(res, "R");
- return tfn;
+ if ( strlen(res) == 0 ) return '?';
+ return res[0];
}
+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(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_point_bwd(P, &cmask, cen, CMASK_ALL, rtnl(1,1), rtnl(1,1), rtnl(1,1));
+
+ /* Check whether the current centering's lattice points will all
+ * coincide with lattice points in the new centering. Eliminate any
+ * centerings for which they don't (they give "excess lattice points"). */
+ switch ( cen ) {
+
+ case 'P' :
+ case 'R' :
+ check_point_fwd(P, &cmask, rtnl(1,1), rtnl(1,1), rtnl(1,1));
+ break;
+
+ case 'A' :
+ check_point_fwd(P, &cmask, rtnl(1,1), rtnl(1,1), rtnl(1,1));
+ check_point_fwd(P, &cmask, rtnl_zero(), rtnl(1,2), rtnl(1,2));
+ break;
+
+ case 'B' :
+ check_point_fwd(P, &cmask, rtnl(1,1), rtnl(1,1), rtnl(1,1));
+ check_point_fwd(P, &cmask, rtnl(1,2), rtnl_zero(), rtnl(1,2));
+ break;
+
+ case 'C' :
+ check_point_fwd(P, &cmask, rtnl(1,1), rtnl(1,1), rtnl(1,1));
+ check_point_fwd(P, &cmask, rtnl(1,2), rtnl(1,2), rtnl_zero());
+ break;
+
+ case 'I' :
+ check_point_fwd(P, &cmask, rtnl(1,1), rtnl(1,1), rtnl(1,1));
+ check_point_fwd(P, &cmask, rtnl(1,2), rtnl(1,2), rtnl(1,2));
+ break;
+
+ case 'F' :
+ check_point_fwd(P, &cmask, rtnl(1,1), rtnl(1,1), rtnl(1,1));
+ 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(P, &cmask, rtnl(1,1), rtnl(1,1), rtnl(1,1));
+ 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;
-/**
- * tfn_from_intmat:
- * @m: An %IntegerMatrix
- *
- * Returns: A %UnitCellTransformation corresponding to @m.
- *
- */
-UnitCellTransformation *tfn_from_intmat(IntegerMatrix *m)
-{
- UnitCellTransformation *tfn;
- int i, j;
-
- tfn = tfn_identity();
- if ( tfn == NULL ) return NULL;
-
- for ( i=0; i<3; i++ ) {
- for ( j=0; j<3; j++ ) {
- gsl_matrix_set(tfn->m, i, j, intmat_get(m, i, j));
- }
}
- return tfn;
+ return cmask_decode(cmask);
}
/**
- * tfn_combine:
- * @t: A %UnitCellTransformation
- * @na: Pointer to three doubles representing naa, nab, nac
- * @nb: Pointer to three doubles representing nba, nbb, nbc
- * @nc: Pointer to three doubles representing nca, ncb, ncc
+ * cell_transform_rational:
+ * @cell: A %UnitCell.
+ * @t: A %RationalMatrix.
+ *
+ * Applies @t to @cell.
*
- * Updates @t such that it represents its previous transformation followed by
- * a new transformation, corresponding to letting a = naa*a + nab*b + nac*c.
- * Likewise, a = nba*a + nbb*b + nbc*c and c = nca*a + ncb*b + ncc*c.
+ * This function will determine the centering of the resulting unit cell,
+ * producing '?' if any lattice points cannot be accounted for. Note that if
+ * there are 'excess' lattice points in the transformed cell, the centering
+ * will still be '?' even if the lattice points for another centering are
+ * all present.
+ *
+ * The lattice type will be set to triclinic, and the unique axis to '?'.
+ *
+ * Returns: Transformed copy of @cell.
*
*/
-void tfn_combine(UnitCellTransformation *t, double *na, double *nb, double *nc)
+UnitCell *cell_transform_rational(UnitCell *cell, RationalMatrix *m)
{
- gsl_matrix *a;
- gsl_matrix *n;
+ UnitCell *out;
+ gsl_matrix *tm;
+ char ncen;
+ int i, j;
+ Rational det;
- n = gsl_matrix_alloc(3, 3);
- a = gsl_matrix_calloc(3, 3);
- if ( (n == NULL) || (a == NULL) ) {
- return;
+ if ( m == NULL ) return NULL;
+
+ det = rtnl_mtx_det(m);
+ if ( rtnl_cmp(det, rtnl_zero()) == 0 ) return NULL;
+
+ tm = gsl_matrix_alloc(3,3);
+ if ( tm == NULL ) {
+ return NULL;
}
- gsl_matrix_set(n, 0, 0, na[0]);
- gsl_matrix_set(n, 0, 1, na[1]);
- gsl_matrix_set(n, 0, 2, na[2]);
- gsl_matrix_set(n, 1, 0, nb[0]);
- gsl_matrix_set(n, 1, 1, nb[1]);
- gsl_matrix_set(n, 1, 2, nb[2]);
- gsl_matrix_set(n, 2, 0, nc[0]);
- gsl_matrix_set(n, 2, 1, nc[1]);
- gsl_matrix_set(n, 2, 2, nc[2]);
- free(na);
- free(nb);
- free(nc);
+ for ( i=0; i<3; i++ ) {
+ for ( j=0; j<3; j++ ) {
+ gsl_matrix_set(tm, i, j,
+ rtnl_as_double(rtnl_mtx_get(m, i, j)));
+ }
+ }
+
+ out = cell_transform_gsl_direct(cell, tm);
+ gsl_matrix_free(tm);
+
+ ncen = determine_centering(m, cell_get_centering(cell));
+ cell_set_centering(out, ncen);
- gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, n, t->m, 0.0, a);
+ /* FIXME: Update unique axis, lattice type */
+ cell_set_lattice_type(out, L_TRICLINIC);
+ cell_set_unique_axis(out, '?');
- gsl_matrix_free(t->m);
- t->m = a;
- gsl_matrix_free(n);
+ return out;
}
/**
- * tfn_vector:
- * @a: Amount of "a" to include in new vector
- * @b: Amount of "b" to include in new vector
- * @c: Amount of "c" to include in new vector
+ * cell_transform_intmat:
+ * @cell: A %UnitCell.
+ * @t: An %IntegerMatrix.
+ *
+ * Applies @t to @cell.
*
- * This is a convenience function to use when sending vectors to tfn_combine():
- * tfn_combine(tfn, tfn_vector(1,0,0),
- * tfn_vector(0,2,0),
- * tfn_vector(0,0,1));
+ * This is just a convenience function which turns @m into a %RationalMatrix
+ * and then calls cell_transform_rational(). See the documentation for that
+ * function for some important information.
+ *
+ * Returns: Transformed copy of @cell.
*
*/
-double *tfn_vector(double a, double b, double c)
+UnitCell *cell_transform_intmat(UnitCell *cell, IntegerMatrix *m)
{
- double *vec = malloc(3*sizeof(double));
- if ( vec == NULL ) return NULL;
- vec[0] = a; vec[1] = b; vec[2] = c;
- return vec;
+ UnitCell *ans;
+ RationalMatrix *mtx = rtnl_mtx_from_intmat(m);
+ ans = cell_transform_rational(cell, mtx);
+ rtnl_mtx_free(mtx);
+ return ans;
}
/**
- * tfn_print:
- * @t: A %UnitCellTransformation
+ * cell_transform_rational_inverse:
+ * @cell: A %UnitCell.
+ * @m: A %RationalMatrix
*
- * Prints information about @t to stderr.
+ * Applies the inverse of @m to @cell.
+ *
+ * Returns: Transformed copy of @cell.
*
*/
-void tfn_print(UnitCellTransformation *t)
+UnitCell *cell_transform_rational_inverse(UnitCell *cell, RationalMatrix *m)
{
+ UnitCell *out;
+ gsl_matrix *tm;
+ gsl_matrix *inv;
gsl_permutation *perm;
- gsl_matrix *lu;
int s;
+ int i, j;
+
+ if ( m == NULL ) return NULL;
- STATUS("New a = %+.2fa %+.2fb %+.2fc\n", gsl_matrix_get(t->m, 0, 0),
- gsl_matrix_get(t->m, 0, 1),
- gsl_matrix_get(t->m, 0, 2));
- STATUS("New b = %+.2fa %+.2fb %+.2fc\n", gsl_matrix_get(t->m, 1, 0),
- gsl_matrix_get(t->m, 1, 1),
- gsl_matrix_get(t->m, 1, 2));
- STATUS("New c = %+.2fa %+.2fb %+.2fc\n", gsl_matrix_get(t->m, 2, 0),
- gsl_matrix_get(t->m, 2, 1),
- gsl_matrix_get(t->m, 2, 2));
- lu = gsl_matrix_alloc(3, 3);
- if ( lu == NULL ) {
- ERROR("Couldn't allocate LU decomposition.\n");
- return;
+ tm = gsl_matrix_alloc(3,3);
+ if ( tm == NULL ) {
+ return NULL;
}
- gsl_matrix_memcpy(lu, t->m);
+ for ( i=0; i<3; i++ ) {
+ for ( j=0; j<3; j++ ) {
+ gsl_matrix_set(tm, i, j,
+ rtnl_as_double(rtnl_mtx_get(m, i, j)));
+ }
+ }
- perm = gsl_permutation_alloc(t->m->size1);
+ perm = gsl_permutation_alloc(3);
if ( perm == NULL ) {
- ERROR("Couldn't allocate permutation.\n");
- gsl_matrix_free(lu);
- return;
+ ERROR("Couldn't allocate permutation\n");
+ return NULL;
}
- if ( gsl_linalg_LU_decomp(lu, perm, &s) ) {
- ERROR("LU decomposition failed.\n");
+ inv = gsl_matrix_alloc(3, 3);
+ if ( inv == NULL ) {
+ ERROR("Couldn't allocate inverse\n");
+ gsl_permutation_free(perm);
+ return NULL;
+ }
+ if ( gsl_linalg_LU_decomp(tm, perm, &s) ) {
+ ERROR("Couldn't decompose matrix\n");
gsl_permutation_free(perm);
- gsl_matrix_free(lu);
- return;
+ return NULL;
+ }
+ if ( gsl_linalg_LU_invert(tm, perm, inv) ) {
+ ERROR("Couldn't invert transformation matrix\n");
+ gsl_permutation_free(perm);
+ return NULL;
}
+ gsl_permutation_free(perm);
+
+ out = cell_transform_gsl_direct(cell, inv);
+
+ /* FIXME: Update centering, unique axis, lattice type */
+
+ gsl_matrix_free(tm);
+ gsl_matrix_free(inv);
- STATUS("Transformation determinant = %.2f\n", gsl_linalg_LU_det(lu, s));
+ return out;
}
/**
- * tfn_free:
- * @t: A %UnitCellTransformation
+ * cell_transform_intmat_inverse:
+ * @cell: A %UnitCell.
+ * @m: An %IntegerMatrix
*
- * Frees all resources associated with @t.
+ * Applies the inverse of @m to @cell.
+ *
+ * Returns: Transformed copy of @cell.
*
*/
-void tfn_free(UnitCellTransformation *t)
+UnitCell *cell_transform_intmat_inverse(UnitCell *cell, IntegerMatrix *m)
{
- gsl_matrix_free(t->m);
- free(t);
+ UnitCell *ans;
+ RationalMatrix *mtx = rtnl_mtx_from_intmat(m);
+ ans = cell_transform_rational_inverse(cell, mtx);
+ rtnl_mtx_free(mtx);
+ return ans;
}
diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h
index 39e6a1ed..c868cd29 100644
--- a/libcrystfel/src/cell.h
+++ b/libcrystfel/src/cell.h
@@ -91,14 +91,6 @@ typedef enum
typedef struct _unitcell UnitCell;
-/**
- * UnitCellTransformation:
- *
- * This opaque data structure represents a tranformation of a unit cell, such
- * as a rotation or a centering operation.
- **/
-typedef struct _unitcelltransformation UnitCellTransformation;
-
#ifdef __cplusplus
extern "C" {
#endif
@@ -156,18 +148,13 @@ extern void cell_set_unique_axis(UnitCell *cell, char unique_axis);
extern const char *cell_rep(UnitCell *cell);
-extern UnitCell *cell_transform(UnitCell *cell, UnitCellTransformation *t);
-extern UnitCell *cell_transform_inverse(UnitCell *cell,
- UnitCellTransformation *t);
-
-extern UnitCellTransformation *tfn_identity(void);
-extern UnitCellTransformation *tfn_from_intmat(IntegerMatrix *m);
-extern void tfn_combine(UnitCellTransformation *t,
- double *na, double *nb, double *nc);
-extern void tfn_print(UnitCellTransformation *t);
-extern UnitCellTransformation *tfn_inverse(UnitCellTransformation *t);
-extern double *tfn_vector(double a, double b, double c);
-extern void tfn_free(UnitCellTransformation *t);
+extern UnitCell *cell_transform_gsl_direct(UnitCell *in, gsl_matrix *m);
+
+extern UnitCell *cell_transform_rational(UnitCell *cell, RationalMatrix *m);
+extern UnitCell *cell_transform_rational_inverse(UnitCell *cell, RationalMatrix *m);
+
+extern UnitCell *cell_transform_intmat(UnitCell *cell, IntegerMatrix *m);
+extern UnitCell *cell_transform_intmat_inverse(UnitCell *cell, IntegerMatrix *m);
#ifdef __cplusplus
}
diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c
index f72334c4..ad9c0de3 100644
--- a/libcrystfel/src/index.c
+++ b/libcrystfel/src/index.c
@@ -698,9 +698,9 @@ static int try_indexer(struct image *image, IndexingMethod indm,
/* Don't do similarity check against bad crystals */
if ( crystal_get_user_flag(that_cr) ) continue;
- if ( compare_cells(crystal_get_cell(cr),
- crystal_get_cell(that_cr),
- 0.1, deg2rad(0.5), NULL) )
+ if ( compare_reindexed_cell_parameters_and_orientation(crystal_get_cell(cr),
+ crystal_get_cell(that_cr),
+ 0.1, deg2rad(0.5), NULL) )
{
crystal_set_user_flag(cr, 1);
}
diff --git a/libcrystfel/src/integer_matrix.c b/libcrystfel/src/integer_matrix.c
index c31072b4..90d58b35 100644
--- a/libcrystfel/src/integer_matrix.c
+++ b/libcrystfel/src/integer_matrix.c
@@ -35,7 +35,9 @@
#include <string.h>
#include <assert.h>
+#include "rational.h"
#include "integer_matrix.h"
+#include "utils.h"
/**
@@ -95,7 +97,7 @@ IntegerMatrix *intmat_new(unsigned int rows, unsigned int cols)
*
* Returns: a newly allocated copy of @m, or NULL on error/
**/
-IntegerMatrix *intmat_copy(IntegerMatrix *m)
+IntegerMatrix *intmat_copy(const IntegerMatrix *m)
{
IntegerMatrix *p;
int i, j;
@@ -127,6 +129,19 @@ void intmat_free(IntegerMatrix *m)
}
+void intmat_size(const IntegerMatrix *m, unsigned int *rows, unsigned int *cols)
+{
+ if ( m == NULL ) {
+ *rows = 0;
+ *cols = 0;
+ return;
+ }
+
+ *rows = m->rows;
+ *cols = m->cols;
+}
+
+
/**
* intmat_set:
* @m: An %IntegerMatrix
@@ -163,32 +178,38 @@ signed int intmat_get(const IntegerMatrix *m, unsigned int i, unsigned int j)
/**
- * intmat_intvec_mult:
- * @m: An %IntegerMatrix
- * @vec: An array of signed integers
+ * transform_indices:
+ * @P: An %IntegerMatrix
+ * @hkl: An array of signed integers
+ *
+ * Apply transformation matrix P to a set of reciprocal space Miller indices.
*
- * Multiplies the matrix @m by the vector @vec. The size of @vec must equal the
- * number of columns in @m, and the size of the result equals the number of rows
- * in @m.
+ * In other words:
+ * Multiplies the matrix @P by the row vector @hkl. The size of @vec must equal
+ * the number of columns in @P, and the size of the result equals the number of
+ * rows in @P.
+ *
+ * The multiplication looks like this:
+ * (a1, a2, a3) = (hkl1, hkl2, hkl3) P
+ * 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 *transform_indices(const IntegerMatrix *P, const signed int *hkl)
{
signed int *ans;
- unsigned int i;
+ unsigned int j;
- ans = malloc(m->rows * sizeof(signed int));
+ ans = malloc(P->rows * sizeof(signed int));
if ( ans == NULL ) return NULL;
- for ( i=0; i<m->rows; i++ ) {
-
- unsigned int j;
+ for ( j=0; j<P->cols; j++ ) {
- ans[i] = 0;
- for ( j=0; j<m->cols; j++ ) {
- ans[i] += intmat_get(m, i, j) * vec[j];
+ unsigned int i;
+ ans[j] = 0;
+ for ( i=0; i<P->rows; i++ ) {
+ ans[i] += intmat_get(P, i, j) * hkl[j];
}
}
@@ -342,6 +363,9 @@ static IntegerMatrix *intmat_cofactors(const IntegerMatrix *m)
*
* Calculates the inverse of @m. Inefficiently.
*
+ * Works only if the inverse of the matrix is also an integer matrix,
+ * i.e. if the determinant of @m is +/- 1.
+ *
* Returns: the inverse of @m, or NULL on error.
**/
IntegerMatrix *intmat_inverse(const IntegerMatrix *m)
@@ -532,3 +556,39 @@ IntegerMatrix *intmat_identity(int size)
return m;
}
+
+
+/**
+ * intmat_set_all_3x3
+ * @m: An %IntegerMatrix
+ *
+ * Returns: an identity %IntegerMatrix with side length @size, or NULL on error.
+ *
+ */
+void intmat_set_all_3x3(IntegerMatrix *m,
+ signed int m11, signed int m12, signed int m13,
+ signed int m21, signed int m22, signed int m23,
+ signed int m31, signed int m32, signed int m33)
+{
+ intmat_set(m, 0, 0, m11);
+ intmat_set(m, 0, 1, m12);
+ intmat_set(m, 0, 2, m13);
+
+ intmat_set(m, 1, 0, m21);
+ intmat_set(m, 1, 1, m22);
+ intmat_set(m, 1, 2, m23);
+
+ intmat_set(m, 2, 0, m31);
+ intmat_set(m, 2, 1, m32);
+ intmat_set(m, 2, 2, m33);
+}
+
+
+IntegerMatrix *intmat_create_3x3(signed int m11, signed int m12, signed int m13,
+ signed int m21, signed int m22, signed int m23,
+ signed int m31, signed int m32, signed int m33)
+{
+ IntegerMatrix *t = intmat_new(3, 3);
+ intmat_set_all_3x3(t, m11, m12, m13, m21, m22, m23, m31, m32, m33);
+ return t;
+}
diff --git a/libcrystfel/src/integer_matrix.h b/libcrystfel/src/integer_matrix.h
index 6616b5e8..6fb3e399 100644
--- a/libcrystfel/src/integer_matrix.h
+++ b/libcrystfel/src/integer_matrix.h
@@ -40,25 +40,38 @@
**/
typedef struct _integermatrix IntegerMatrix;
+#include "rational.h"
+
#ifdef __cplusplus
extern "C" {
#endif
/* Alloc/dealloc */
extern IntegerMatrix *intmat_new(unsigned int rows, unsigned int cols);
-extern IntegerMatrix *intmat_copy(IntegerMatrix *m);
+extern IntegerMatrix *intmat_copy(const IntegerMatrix *m);
extern IntegerMatrix *intmat_identity(int size);
extern void intmat_free(IntegerMatrix *m);
/* Get/set */
+extern void intmat_size(const IntegerMatrix *m, unsigned int *rows,
+ unsigned int *cols);
+
extern void intmat_set(IntegerMatrix *m, unsigned int i, unsigned int j,
signed int v);
extern signed int intmat_get(const IntegerMatrix *m,
unsigned int i, unsigned int j);
-/* Matrix-(int)vector multiplication */
-extern signed int *intmat_intvec_mult(const IntegerMatrix *m,
- const signed int *vec);
+extern void intmat_set_all_3x3(IntegerMatrix *m,
+ signed int m11, signed int m12, signed int m13,
+ signed int m21, signed int m22, signed int m23,
+ signed int m31, signed int m32, signed int m33);
+
+extern IntegerMatrix *intmat_create_3x3(signed int m11, signed int m12, signed int m13,
+ signed int m21, signed int m22, signed int m23,
+ signed int m31, signed int m32, signed int m33);
+
+/* Matrix-vector multiplication */
+extern signed int *transform_indices(const IntegerMatrix *P, const signed int *hkl);
/* Matrix-matrix multiplication */
extern IntegerMatrix *intmat_intmat_mult(const IntegerMatrix *a,
diff --git a/libcrystfel/src/rational.c b/libcrystfel/src/rational.c
new file mode 100644
index 00000000..65b209ff
--- /dev/null
+++ b/libcrystfel/src/rational.c
@@ -0,0 +1,675 @@
+/*
+ * rational.c
+ *
+ * A small rational number library
+ *
+ * Copyright © 2019 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2019 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <assert.h>
+#include <locale.h>
+
+#include "rational.h"
+#include "integer_matrix.h"
+#include "utils.h"
+
+
+/**
+ * SECTION:rational
+ * @short_description: Rational numbers
+ * @title: Rational numbers
+ * @section_id:
+ * @see_also:
+ * @include: "rational.h"
+ * @Image:
+ *
+ * A rational number library
+ */
+
+
+/* Eucliden algorithm for finding greatest common divisor */
+static signed int gcd(signed int a, signed int b)
+{
+ while ( b != 0 ) {
+ signed int t = b;
+ b = a % b;
+ a = t;
+ }
+ return a;
+}
+
+
+static void squish(Rational *rt)
+{
+ signed int g;
+
+ if ( rt->num == 0 ) {
+ rt->den = 1;
+ return;
+ }
+
+ g = gcd(rt->num, rt->den);
+ assert(g != 0);
+ rt->num /= g;
+ rt->den /= g;
+
+ if ( rt->den < 0 ) {
+ rt->num = -rt->num;
+ rt->den = -rt->den;
+ }
+}
+
+
+Rational rtnl_zero()
+{
+ Rational r;
+ r.num = 0;
+ r.den = 1;
+ return r;
+}
+
+
+Rational rtnl(signed long long int num, signed long long int den)
+{
+ Rational r;
+ r.num = num;
+ r.den = den;
+ squish(&r);
+ return r;
+}
+
+
+double rtnl_as_double(Rational r)
+{
+ return (double)r.num/r.den;
+}
+
+
+static void overflow(long long int c, long long int a, long long int b)
+{
+ setlocale(LC_ALL, "");
+ ERROR("Overflow detected in rational number library.\n");
+ ERROR("%'lli < %'lli * %'lli\n", c, a, b);
+ abort();
+}
+
+
+static void check_overflow(long long int c, long long int a, long long int b)
+{
+ if ( (a==0) || (b==0) ) {
+ if ( c != 0 ) overflow(c,a,b);
+ } else if ( (llabs(c) < llabs(a)) || (llabs(c) < llabs(b)) ) {
+ overflow(c,a,b);
+ }
+}
+
+
+Rational rtnl_mul(Rational a, Rational b)
+{
+ Rational r;
+ r.num = a.num * b.num;
+ r.den = a.den * b.den;
+ check_overflow(r.num, a.num, b.num);
+ check_overflow(r.den, a.den, b.den);
+ squish(&r);
+ return r;
+}
+
+
+Rational rtnl_div(Rational a, Rational b)
+{
+ signed int t = b.num;
+ b.num = b.den;
+ b.den = t;
+ return rtnl_mul(a, b);
+}
+
+
+Rational rtnl_add(Rational a, Rational b)
+{
+ Rational r, trt1, trt2;
+
+ trt1.num = a.num * b.den;
+ trt2.num = b.num * a.den;
+ check_overflow(trt1.num, a.num, b.den);
+ check_overflow(trt2.num, b.num, a.den);
+
+ trt1.den = a.den * b.den;
+ trt2.den = trt1.den;
+ check_overflow(trt1.den, a.den, b.den);
+
+ r.num = trt1.num + trt2.num;
+ r.den = trt1.den;
+ squish(&r);
+ return r;
+}
+
+
+Rational rtnl_sub(Rational a, Rational b)
+{
+ b.num = -b.num;
+ return rtnl_add(a, b);
+}
+
+
+/* -1, 0 +1 respectively for a<b, a==b, a>b */
+signed int rtnl_cmp(Rational a, Rational b)
+{
+ Rational trt1, trt2;
+
+ trt1.num = a.num * b.den;
+ trt2.num = b.num * a.den;
+
+ trt1.den = a.den * b.den;
+ trt2.den = a.den * b.den;
+
+ if ( trt1.num > trt2.num ) return +1;
+ if ( trt1.num < trt2.num ) return -1;
+ return 0;
+}
+
+
+Rational rtnl_abs(Rational a)
+{
+ Rational r = a;
+ squish(&r);
+ if ( r.num < 0 ) r.num = -r.num;
+ return r;
+}
+
+
+/**
+ * rtnl_format
+ * @rt: A %Rational
+ *
+ * Formats @rt as a string
+ *
+ * Returns: a string which should be freed by the caller
+ */
+char *rtnl_format(Rational rt)
+{
+ char *v = malloc(32);
+ if ( v == NULL ) return NULL;
+ if ( rt.den == 1 ) {
+ snprintf(v, 31, "%lli", rt.num);
+ } else {
+ snprintf(v, 31, "%lli/%lli", rt.num, rt.den);
+ }
+ return v;
+}
+
+
+Rational *rtnl_list(signed int num_min, signed int num_max,
+ signed int den_min, signed int den_max,
+ int *pn)
+{
+ signed int num, den;
+ Rational *list;
+ int n = 0;
+
+ list = malloc((1+num_max-num_min)*(1+den_max-den_min)*sizeof(Rational));
+ if ( list == NULL ) return NULL;
+
+ for ( num=num_min; num<=num_max; num++ ) {
+ for ( den=den_min; den<=den_max; den++ ) {
+
+ Rational r = rtnl(num, den);
+
+ /* Denominator zero? */
+ if ( den == 0 ) continue;
+
+ /* Same as last entry? */
+ if ( (n>0) && (rtnl_cmp(list[n-1], r)==0) ) continue;
+
+ /* Can be reduced? */
+ if ( gcd(num, den) != 1 ) continue;
+
+ list[n++] = r;
+ }
+ }
+ *pn = n;
+ return list;
+}
+
+
+/**
+ * SECTION:rational_matrix
+ * @short_description: Rational matrices
+ * @title: Rational matrices
+ * @section_id:
+ * @see_also:
+ * @include: "rational.h"
+ * @Image:
+ *
+ * A rational matrix library
+ */
+
+
+struct _rationalmatrix
+{
+ unsigned int rows;
+ unsigned int cols;
+ Rational *v;
+};
+
+
+/**
+ * rtnl_mtx_new:
+ * @rows: Number of rows that the new matrix is to have
+ * @cols: Number of columns that the new matrix is to have
+ *
+ * Allocates a new %RationalMatrix with all elements set to zero.
+ *
+ * Returns: a new %RationalMatrix, or NULL on error.
+ **/
+RationalMatrix *rtnl_mtx_new(unsigned int rows, unsigned int cols)
+{
+ RationalMatrix *m;
+ int i;
+
+ m = malloc(sizeof(RationalMatrix));
+ if ( m == NULL ) return NULL;
+
+ m->v = calloc(rows*cols, sizeof(Rational));
+ if ( m->v == NULL ) {
+ free(m);
+ return NULL;
+ }
+
+ m->rows = rows;
+ m->cols = cols;
+
+ for ( i=0; i<m->rows*m->cols; i++ ) {
+ m->v[i] = rtnl_zero();
+ }
+
+ return m;
+}
+
+
+RationalMatrix *rtnl_mtx_identity(int rows)
+{
+ int i;
+ RationalMatrix *m = rtnl_mtx_new(rows, rows);
+ for ( i=0; i<rows; i++ ) {
+ rtnl_mtx_set(m, i, i, rtnl(1,1));
+ }
+ return m;
+}
+
+
+RationalMatrix *rtnl_mtx_copy(const RationalMatrix *m)
+{
+ RationalMatrix *n;
+ int i;
+
+ n = rtnl_mtx_new(m->rows, m->cols);
+ if ( n == NULL ) return NULL;
+
+ for ( i=0; i<m->rows*m->cols; i++ ) {
+ n->v[i] = m->v[i];
+ }
+
+ return n;
+}
+
+
+Rational rtnl_mtx_get(const RationalMatrix *m, int i, int j)
+{
+ assert(m != NULL);
+ return m->v[j+m->cols*i];
+}
+
+
+void rtnl_mtx_set(const RationalMatrix *m, int i, int j, Rational v)
+{
+ assert(m != NULL);
+ m->v[j+m->cols*i] = v;
+}
+
+
+RationalMatrix *rtnl_mtx_from_intmat(const IntegerMatrix *m)
+{
+ RationalMatrix *n;
+ unsigned int rows, cols;
+ int i, j;
+
+ intmat_size(m, &rows, &cols);
+ n = rtnl_mtx_new(rows, cols);
+ if ( n == NULL ) return NULL;
+
+ for ( i=0; i<rows; i++ ) {
+ for ( j=0; j<cols; j++ ) {
+ rtnl_mtx_set(n, i, j, rtnl(intmat_get(m, i, j), 1));
+ }
+ }
+
+ return n;
+}
+
+
+IntegerMatrix *intmat_from_rtnl_mtx(const RationalMatrix *m)
+{
+ IntegerMatrix *n;
+ int i, j;
+
+ n = intmat_new(m->rows, m->cols);
+ if ( n == NULL ) return NULL;
+
+ for ( i=0; i<m->rows; i++ ) {
+ for ( j=0; j<m->cols; j++ ) {
+ Rational v = rtnl_mtx_get(m, i, j);
+ squish(&v);
+ if ( v.den != 1 ) {
+ ERROR("Rational matrix can't be converted to integers\n");
+ intmat_free(n);
+ return NULL;
+ }
+ intmat_set(n, i, j, v.num);
+ }
+ }
+
+ return n;
+}
+
+
+void rtnl_mtx_free(RationalMatrix *mtx)
+{
+ if ( mtx == NULL ) return;
+ free(mtx->v);
+ free(mtx);
+}
+
+
+/* rtnl_mtx_solve:
+ * @P: A %RationalMatrix
+ * @vec: An array of %Rational
+ * @ans: An array of %Rational in which to store the result
+ *
+ * Solves the matrix equation m*ans = vec, where @ans and @vec are
+ * vectors of %Rational.
+ *
+ * In this version, @m must be square.
+ *
+ * The number of columns in @m must equal the length of @ans, and the number
+ * 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 *P,
+ const Rational *ivec, Rational *ans)
+{
+ RationalMatrix *cm;
+ Rational *vec;
+ int h, k;
+ int i;
+
+ if ( P->rows != P->cols ) return 1;
+
+ /* Copy the matrix and vector because the calculation will
+ * be done in-place */
+ cm = rtnl_mtx_copy(P);
+ if ( cm == NULL ) return 1;
+
+ vec = malloc(cm->rows*sizeof(Rational));
+ if ( vec == NULL ) return 1;
+ for ( h=0; h<cm->rows; h++ ) vec[h] = ivec[h];
+
+ /* Gaussian elimination with partial pivoting */
+ h = 0;
+ k = 0;
+ 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<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;
+ }
+ }
+
+ if ( rtnl_cmp(pval, rtnl_zero()) == 0 ) {
+ k++;
+ continue;
+ }
+
+ /* Swap 'prow' with row h */
+ 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);
+ }
+ t = vec[prow];
+ vec[prow] = vec[h];
+ vec[h] = t;
+
+ /* Divide and subtract rows below */
+ for ( i=h+1; i<cm->rows; i++ ) {
+
+ int j;
+ Rational dval;
+
+ dval = rtnl_div(rtnl_mtx_get(cm, i, k),
+ rtnl_mtx_get(cm, h, k));
+
+ 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);
+ rtnl_mtx_set(cm, i, j, t);
+ }
+
+ /* Divide the right hand side as well */
+ Rational t = vec[i];
+ Rational p = rtnl_mul(dval, vec[h]);
+ vec[i] = rtnl_sub(t, p);
+ }
+
+ h++;
+ k++;
+
+
+ }
+
+ /* Back-substitution */
+ for ( i=cm->rows-1; i>=0; i-- ) {
+ int j;
+ Rational sum = rtnl_zero();
+ 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);
+ }
+ sum = rtnl_sub(vec[i], sum);
+ ans[i] = rtnl_div(sum, rtnl_mtx_get(cm, i, i));
+ }
+
+ free(vec);
+ rtnl_mtx_free(cm);
+
+ return 0;
+}
+
+
+/**
+ * rtnl_mtx_print
+ * @m: A %RationalMatrix
+ *
+ * Prints @m to stderr.
+ *
+ */
+void rtnl_mtx_print(const RationalMatrix *m)
+{
+ unsigned int i, j;
+
+ for ( i=0; i<m->rows; i++ ) {
+
+ fprintf(stderr, "[ ");
+ for ( j=0; j<m->cols; j++ ) {
+ char *v = rtnl_format(rtnl_mtx_get(m, i, j));
+ fprintf(stderr, "%4s ", v);
+ free(v);
+ }
+ fprintf(stderr, "]\n");
+ }
+}
+
+
+void rtnl_mtx_mtxmult(const RationalMatrix *A, const RationalMatrix *B,
+ RationalMatrix *ans)
+{
+ int i, j;
+
+ assert(ans->cols == B->cols);
+ assert(ans->rows == A->rows);
+ assert(A->cols == B->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_mtx_get(B, k, j));
+ 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.
+ * 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<P->rows; i++ ) {
+ ans[i] = rtnl_zero();
+ for ( j=0; j<P->cols; j++ ) {
+ Rational add;
+ add = rtnl_mul(rtnl_mtx_get(P, i, j), vec[j]);
+ ans[i] = rtnl_add(ans[i], add);
+ }
+ }
+}
+
+
+static RationalMatrix *delete_row_and_column(const RationalMatrix *m,
+ unsigned int di, unsigned int dj)
+{
+ RationalMatrix *n;
+ unsigned int i, j;
+
+ n = rtnl_mtx_new(m->rows-1, m->cols-1);
+ if ( n == NULL ) return NULL;
+
+ for ( i=0; i<n->rows; i++ ) {
+ for ( j=0; j<n->cols; j++ ) {
+
+ Rational val;
+ unsigned int gi, gj;
+
+ gi = (i>=di) ? i+1 : i;
+ gj = (j>=dj) ? j+1 : j;
+ val = rtnl_mtx_get(m, gi, gj);
+ rtnl_mtx_set(n, i, j, val);
+
+ }
+ }
+
+ return n;
+}
+
+
+static Rational cofactor(const RationalMatrix *m,
+ unsigned int i, unsigned int j)
+{
+ RationalMatrix *n;
+ Rational t, C;
+
+ n = delete_row_and_column(m, i, j);
+ if ( n == NULL ) {
+ fprintf(stderr, "Failed to allocate matrix.\n");
+ return rtnl_zero();
+ }
+
+ /* -1 if odd, +1 if even */
+ t = (i+j) & 0x1 ? rtnl(-1, 1) : rtnl(1, 1);
+
+ C = rtnl_mul(t, rtnl_mtx_det(n));
+ rtnl_mtx_free(n);
+
+ return C;
+}
+
+
+
+Rational rtnl_mtx_det(const RationalMatrix *m)
+{
+ unsigned int i, j;
+ Rational det;
+
+ assert(m->rows == m->cols); /* Otherwise determinant doesn't exist */
+
+ if ( m->rows == 2 ) {
+ Rational a, b;
+ a = rtnl_mul(rtnl_mtx_get(m, 0, 0), rtnl_mtx_get(m, 1, 1));
+ b = rtnl_mul(rtnl_mtx_get(m, 0, 1), rtnl_mtx_get(m, 1, 0));
+ return rtnl_sub(a, b);
+ }
+
+ i = 0; /* Fixed */
+ det = rtnl_zero();
+ for ( j=0; j<m->cols; j++ ) {
+ Rational a;
+ a = rtnl_mul(rtnl_mtx_get(m, i, j), cofactor(m, i, j));
+ det = rtnl_add(det, a);
+ }
+
+ return det;
+}
diff --git a/libcrystfel/src/rational.h b/libcrystfel/src/rational.h
new file mode 100644
index 00000000..23c918cf
--- /dev/null
+++ b/libcrystfel/src/rational.h
@@ -0,0 +1,107 @@
+/*
+ * rational.h
+ *
+ * A small rational number library
+ *
+ * Copyright © 2019 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2019 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef RATIONAL_H
+#define RATIONAL_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+/**
+ * Rational
+ *
+ * The Rational is an opaque-ish data structure representing a rational number.
+ *
+ * "Opaque-ish" means that the structure isn't technically opaque, allowing you
+ * to assign and allocate them easily. But you shouldn't look at or set its
+ * contents, except by using the accessor functions.
+ **/
+typedef struct {
+ /* Private, don't modify */
+ signed long long int num;
+ signed long long int den;
+} Rational;
+
+
+/**
+ * RationalMatrix
+ *
+ * The RationalMatrix is an opaque data structure representing a matrix of
+ * rational numbers.
+ **/
+typedef struct _rationalmatrix RationalMatrix;
+
+#include "integer_matrix.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+extern Rational rtnl_zero(void);
+extern Rational rtnl(signed long long int num, signed long long int den);
+extern double rtnl_as_double(Rational r);
+
+extern Rational rtnl_mul(Rational a, Rational b);
+extern Rational rtnl_div(Rational a, Rational b);
+extern Rational rtnl_add(Rational a, Rational b);
+extern Rational rtnl_sub(Rational a, Rational b);
+
+extern signed int rtnl_cmp(Rational a, Rational b);
+extern Rational rtnl_abs(Rational a);
+
+extern char *rtnl_format(Rational rt);
+
+extern Rational *rtnl_list(signed int num_min, signed int num_max,
+ signed int den_min, signed int den_max,
+ int *pn);
+
+extern RationalMatrix *rtnl_mtx_new(unsigned int rows, unsigned int cols);
+extern RationalMatrix *rtnl_mtx_copy(const RationalMatrix *m);
+extern Rational rtnl_mtx_get(const RationalMatrix *m, int i, int j);
+extern void rtnl_mtx_set(const RationalMatrix *m, int i, int j, Rational v);
+extern RationalMatrix *rtnl_mtx_from_intmat(const IntegerMatrix *m);
+extern RationalMatrix *rtnl_mtx_identity(int rows);
+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 int transform_fractional_coords_rtnl(const RationalMatrix *P,
+ const Rational *ivec,
+ Rational *ans);
+extern void transform_fractional_coords_rtnl_inverse(const RationalMatrix *P,
+ const Rational *vec,
+ Rational *ans);
+extern void rtnl_mtx_print(const RationalMatrix *m);
+extern Rational rtnl_mtx_det(const RationalMatrix *m);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* RATIONAL_H */
diff --git a/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c
index 06cc368c..69dada16 100644
--- a/libcrystfel/src/symmetry.c
+++ b/libcrystfel/src/symmetry.c
@@ -3,12 +3,12 @@
*
* Symmetry
*
- * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2019 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2014,2016 Thomas White <taw@physics.org>
- * 2014 Kenneth Beyerlein <kenneth.beyerlein@desy.de>
+ * 2010-2019 Thomas White <taw@physics.org>
+ * 2014 Kenneth Beyerlein <kenneth.beyerlein@desy.de>
*
* This file is part of CrystFEL.
*
@@ -41,6 +41,8 @@
#include "symmetry.h"
#include "utils.h"
#include "integer_matrix.h"
+#include "symop-parse.h"
+#include "symop-lex.h"
/**
@@ -195,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);
@@ -369,13 +371,13 @@ static void expand_ops(SymOpList *s)
/* Transform all the operations in a SymOpList by a given matrix.
* The matrix must have a determinant of +/- 1 (otherwise its inverse would
* not also be an integer matrix). */
-static void transform_ops(SymOpList *s, IntegerMatrix *t)
+static void transform_ops(SymOpList *s, IntegerMatrix *P)
{
int n, i;
- IntegerMatrix *inv;
+ IntegerMatrix *Pi;
signed int det;
- det = intmat_det(t);
+ det = intmat_det(P);
if ( det == -1 ) {
ERROR("WARNING: mirrored SymOpList.\n");
} else if ( det != 1 ) {
@@ -383,8 +385,8 @@ static void transform_ops(SymOpList *s, IntegerMatrix *t)
return;
}
- inv = intmat_inverse(t);
- if ( inv == NULL ) {
+ Pi = intmat_inverse(P);
+ if ( Pi == NULL ) {
ERROR("Failed to invert matrix.\n");
return;
}
@@ -394,13 +396,13 @@ static void transform_ops(SymOpList *s, IntegerMatrix *t)
IntegerMatrix *r, *f;
- r = intmat_intmat_mult(s->ops[i], t);
+ r = intmat_intmat_mult(P, s->ops[i]);
if ( r == NULL ) {
ERROR("Matrix multiplication failed.\n");
return;
}
- f = intmat_intmat_mult(inv, r);
+ f = intmat_intmat_mult(r, Pi);
if ( f == NULL ) {
ERROR("Matrix multiplication failed.\n");
return;
@@ -413,7 +415,7 @@ static void transform_ops(SymOpList *s, IntegerMatrix *t)
}
- intmat_free(inv);
+ intmat_free(Pi);
}
@@ -1012,14 +1014,14 @@ static SymOpList *getpg_arbitrary_ua(const char *sym, size_t s)
case 'a' :
intmat_set(t, 0, 2, 1);
- intmat_set(t, 1, 1, 1);
- intmat_set(t, 2, 0, -1);
+ intmat_set(t, 1, 0, 1);
+ intmat_set(t, 2, 1, 1);
break;
case 'b' :
- intmat_set(t, 0, 0, 1);
+ intmat_set(t, 0, 1, 1);
intmat_set(t, 1, 2, 1);
- intmat_set(t, 2, 1, -1);
+ intmat_set(t, 2, 0, 1);
break;
@@ -1124,7 +1126,7 @@ static void do_op(const IntegerMatrix *op,
v[0] = h; v[1] = k; v[2] = l;
- ans = intmat_intvec_mult(op, v);
+ ans = transform_indices(op, v);
assert(ans != NULL);
*he = ans[0]; *ke = ans[1]; *le = ans[2];
@@ -1636,105 +1638,76 @@ SymOpList *get_ambiguities(const SymOpList *source, const SymOpList *target)
}
-static IntegerMatrix *parse_symmetry_operation(const char *s)
+/* Parse a single symmetry operation, e.g. 'h,-2k,(h+l)/3' */
+RationalMatrix *parse_symmetry_operation(const char *s)
{
- IntegerMatrix *m;
- char **els;
- int n, i;
+ YY_BUFFER_STATE b;
+ RationalMatrix *m;
+ int r;
+ m = rtnl_mtx_new(3, 3);
+ b = symop_scan_string(s);
+ r = symopparse(m, NULL);
+ symop_delete_buffer(b);
- n = assplode(s, ",", &els, ASSPLODE_NONE);
- if ( n != 3 ) {
- for ( i=0; i<n; i++ ) free(els[i]);
- free(els);
+ if ( r ) {
+ ERROR("Failed to parse '%s'\n", s);
+ rtnl_mtx_free(m);
return NULL;
}
- m = intmat_new(3, 3);
- if ( m == NULL ) return NULL;
-
- for ( i=0; i<n; i++ ) {
-
- int c;
- size_t cl;
- signed int nh = 0;
- signed int nk = 0;
- signed int nl = 0;
- signed int mult = 1;
- int ndigit = 0;
- signed int sign = +1;
-
- /* We have one expression something like "-2h+k" */
- cl = strlen(els[i]);
- for ( c=0; c<cl; c++ ) {
-
- if ( els[i][c] == '-' ) sign *= -1;
- if ( els[i][c] == 'h' ) {
- nh = mult*sign;
- mult = 1;
- ndigit = 0;
- sign = +1;
- }
- if ( els[i][c] == 'k' ) {
- nk = mult*sign;
- mult = 1;
- ndigit = 0;
- sign = +1;
- }
- if ( els[i][c] == 'l' ) {
- nl = mult*sign;
- mult = 1;
- ndigit = 0;
- sign = +1;
- }
- if ( isdigit(els[i][c]) ) {
- if ( ndigit > 0 ) {
- mult *= 10;
- mult += els[i][c] - '0';
- } else {
- mult *= els[i][c] - '0';
- }
- ndigit++;
- }
- }
-
- intmat_set(m, i, 0, nh);
- intmat_set(m, i, 1, nk);
- intmat_set(m, i, 2, nl);
-
- free(els[i]);
+ return m;
+}
- }
- free(els);
- return m;
+/**
+ * parse_cell_transformation
+ * @s: Textual representation of cell transformation
+ *
+ * Parses @s, for example 'a,(b+a)/2,c', and returns the corresponding
+ * %RationalMatrix.
+ *
+ * Returns: A %RationalMatrix describing the transformation, or NULL on error.
+ *
+ */
+RationalMatrix *parse_cell_transformation(const char *s)
+{
+ return parse_symmetry_operation(s);
}
+/**
+ * parse_symmetry_operations
+ * @s: Textual representation of a list of symmetry operations
+ *
+ * Parses @s, for example 'h,k,l;k,h,-l', and returns the corresponding
+ * %SymOpList
+ *
+ * Returns: A %SymOpList, or NULL on error.
+ *
+ */
SymOpList *parse_symmetry_operations(const char *s)
{
- SymOpList *sol;
- char **ops;
- int n, i;
+ YY_BUFFER_STATE b;
+ RationalMatrix *m;
+ SymOpList *list;
+ int r;
- sol = new_symoplist();
- if ( sol == NULL ) return NULL;
+ m = rtnl_mtx_new(3, 3); /* Scratch space for parser */
+ list = new_symoplist(); /* The result we want */
- n = assplode(s, ";:", &ops, ASSPLODE_NONE);
- for ( i=0; i<n; i++ ) {
- IntegerMatrix *m;
- m = parse_symmetry_operation(ops[i]);
- if ( m != NULL ) {
- add_symop(sol, m);
- } else {
- ERROR("Invalid symmetry operation '%s'\n", ops[i]);
- /* Try the next one */
- }
- free(ops[i]);
+ b = symop_scan_string(s);
+ r = symopparse(m, list);
+ symop_delete_buffer(b);
+ rtnl_mtx_free(m);
+
+ if ( r ) {
+ ERROR("Failed to parse '%s'\n", s);
+ free_symoplist(list);
+ return NULL;
}
- free(ops);
- return sol;
+ return list;
}
diff --git a/libcrystfel/src/symmetry.h b/libcrystfel/src/symmetry.h
index a61ca96f..ab2aa934 100644
--- a/libcrystfel/src/symmetry.h
+++ b/libcrystfel/src/symmetry.h
@@ -3,12 +3,12 @@
*
* Symmetry
*
- * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2019 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2014,2016 Thomas White <taw@physics.org>
- * 2014 Kenneth Beyerlein <kenneth.beyerlein@desy.de>
+ * 2010-2019 Thomas White <taw@physics.org>
+ * 2014 Kenneth Beyerlein <kenneth.beyerlein@desy.de>
*
* This file is part of CrystFEL.
*
@@ -36,6 +36,7 @@
#include "integer_matrix.h"
+#include "rational.h"
/**
* SymOpList
@@ -93,7 +94,9 @@ extern int is_centric(signed int h, signed int k, signed int l,
extern void pointgroup_warning(const char *sym);
extern void add_symop(SymOpList *ops, IntegerMatrix *m);
+extern RationalMatrix *parse_symmetry_operation(const char *s);
extern SymOpList *parse_symmetry_operations(const char *s);
+extern RationalMatrix *parse_cell_transformation(const char *s);
extern char *get_matrix_name(const IntegerMatrix *m, int row);
#ifdef __cplusplus
diff --git a/libcrystfel/src/symop.l b/libcrystfel/src/symop.l
new file mode 100644
index 00000000..ed38a65e
--- /dev/null
+++ b/libcrystfel/src/symop.l
@@ -0,0 +1,52 @@
+/*
+ * symop.l
+ *
+ * Lexical scanner for symmetry operations
+ *
+ * Copyright © 2019 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2019 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+%{
+ #define YYDEBUG 1
+ #include "rational.h"
+ #include "symop-parse.h"
+%}
+
+%option prefix="symop"
+%option noyywrap nounput noinput
+
+%%
+
+[,] { return COMMA; }
+[0-9]+ { symoplval.n = atoi(yytext); return NUMBER; }
+[/] { return DIVIDE; }
+[+] { return PLUS; }
+[-] { return MINUS; }
+[ahx] { return H; }
+[bky] { return K; }
+[clz] { return L; }
+[(] { return OPENB; }
+[)] { return CLOSEB; }
+[;] { return SEMICOLON; }
+
+%%
diff --git a/libcrystfel/src/symop.y b/libcrystfel/src/symop.y
new file mode 100644
index 00000000..be31c21d
--- /dev/null
+++ b/libcrystfel/src/symop.y
@@ -0,0 +1,140 @@
+/*
+ * symop.y
+ *
+ * Parser for symmetry operations
+ *
+ * Copyright © 2019 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2019 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+%{
+ #include <stdio.h>
+
+ #include "rational.h"
+ #include "symmetry.h"
+
+ extern int symoplex();
+ extern int symopparse(RationalMatrix *m, SymOpList *list);
+ void symoperror(RationalMatrix *m, SymOpList *list, const char *s);
+%}
+
+%define api.prefix {symop}
+
+%parse-param {RationalMatrix *m} {SymOpList *list}
+
+%code requires {
+ #include "symmetry.h"
+}
+
+%union {
+ RationalMatrix *m; /* Full rational matrix */
+ Rational rv[3]; /* Rational vector, e.g. '1/2h+3k' */
+ Rational r; /* Rational number */
+ int n; /* Just a number */
+}
+
+%token SEMICOLON
+%token COMMA
+%token NUMBER
+%token OPENB CLOSEB
+%token H K L
+
+%left PLUS MINUS
+%left DIVIDE
+%precedence MUL
+%precedence NEG
+
+%type <m> symop
+%type <rv> axexpr
+%type <rv> part
+%type <n> NUMBER
+%type <r> fraction
+
+%{
+static int try_add_symop(SymOpList *list, RationalMatrix *m, int complain)
+{
+ if ( list == NULL ) {
+ /* Only complain if this isn't the only operation provided */
+ if ( complain ) {
+ yyerror(m, list, "Must be a single symmetry operation");
+ }
+ return 1;
+ } else {
+ IntegerMatrix *im;
+ im = intmat_from_rtnl_mtx(m);
+ if ( im == NULL ) {
+ yyerror(m, list, "Symmetry operations must all be integer");
+ return 1;
+ } else {
+ add_symop(list, im);
+ }
+ }
+ return 0;
+}
+%}
+
+%%
+
+symoplist:
+ symop { try_add_symop(list, m, 0); }
+| symoplist SEMICOLON symop { if ( try_add_symop(list, m, 1) ) YYERROR; }
+;
+
+symop:
+ axexpr COMMA axexpr COMMA axexpr { rtnl_mtx_set(m, 0, 0, $1[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, 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]);
+ }
+;
+
+axexpr:
+ part { int i; for ( i=0; i<3; i++ ) $$[i] = $1[i]; }
+| axexpr PLUS axexpr { int i; for ( i=0; i<3; i++ ) $$[i] = rtnl_add($1[i], $3[i]); }
+| axexpr MINUS axexpr { int i; for ( i=0; i<3; i++ ) $$[i] = rtnl_sub($1[i], $3[i]); }
+| MINUS axexpr %prec NEG { int i; for ( i=0; i<3; i++ ) $$[i] = rtnl_sub(rtnl_zero(), $2[i]); }
+| OPENB axexpr CLOSEB { int i; for ( i=0; i<3; i++ ) $$[i] = $2[i]; }
+| axexpr DIVIDE NUMBER { int i; for ( i=0; i<3; i++ ) $$[i] = rtnl_div($1[i], rtnl($3, 1)); }
+| NUMBER axexpr %prec MUL { int i; for ( i=0; i<3; i++ ) $$[i] = rtnl_mul($2[i], rtnl($1, 1)); }
+| fraction axexpr %prec MUL { int i; for ( i=0; i<3; i++ ) $$[i] = rtnl_mul($2[i], $1); }
+;
+
+part:
+ H { $$[0] = rtnl(1, 1); $$[1] = rtnl_zero(); $$[2] = rtnl_zero(); }
+| K { $$[1] = rtnl(1, 1); $$[0] = rtnl_zero(); $$[2] = rtnl_zero(); }
+| L { $$[2] = rtnl(1, 1); $$[0] = rtnl_zero(); $$[1] = rtnl_zero(); }
+;
+
+fraction:
+ NUMBER DIVIDE NUMBER { $$ = rtnl($1, $3); }
+;
+
+%%
+
+void symoperror(RationalMatrix *m, SymOpList *list, const char *s) {
+ printf("Error: %s\n", s);
+}
diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c
index 203e5a05..933fa76a 100644
--- a/libcrystfel/src/taketwo.c
+++ b/libcrystfel/src/taketwo.c
@@ -98,6 +98,7 @@
#include <assert.h>
#include <time.h>
+#include "cell.h"
#include "cell-utils.h"
#include "index.h"
#include "taketwo.h"
@@ -355,6 +356,54 @@ static void show_rvec(struct rvec r2)
* functions called under the core functions, still specialised (Level 3)
* ------------------------------------------------------------------------*/
+/* cell_transform_gsl_direct() doesn't do quite what we want here.
+ * The matrix m should be post-multiplied by a matrix of real or reciprocal
+ * basis vectors (it doesn't matter which because it's just a rotation).
+ * M contains the basis vectors written in columns: M' = mM */
+static UnitCell *cell_post_smiley_face(UnitCell *in, gsl_matrix *m)
+{
+ gsl_matrix *c;
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+ gsl_matrix *res;
+ UnitCell *out;
+
+ cell_get_cartesian(in, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+
+ c = gsl_matrix_alloc(3, 3);
+ gsl_matrix_set(c, 0, 0, asx);
+ 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, 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);
+
+ out = cell_new_from_cell(in);
+ cell_set_cartesian(out, gsl_matrix_get(res, 0, 0),
+ gsl_matrix_get(res, 1, 0),
+ 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);
+ gsl_matrix_free(c);
+ return out;
+}
+
+
static void rotation_around_axis(struct rvec c, double th,
gsl_matrix *res)
{
@@ -463,35 +512,6 @@ static void closest_rot_mat(struct rvec vec1, struct rvec vec2,
rotation_around_axis(axis, bestAngle, twizzle);
}
-/*
-static double matrix_angle(gsl_matrix *m)
-{
- double a = gsl_matrix_get(m, 0, 0);
- double b = gsl_matrix_get(m, 1, 1);
- double c = gsl_matrix_get(m, 2, 2);
-
- double cos_t = (a + b + c - 1) / 2;
- double theta = acos(cos_t);
-
- return theta;
-}
-
-static struct rvec matrix_axis(gsl_matrix *a)
-{
- double ang = matrix_angle(a);
- double cos_t = cos(ang);
- double p = gsl_matrix_get(a, 0, 0);
- double q = gsl_matrix_get(a, 1, 1);
- double r = gsl_matrix_get(a, 2, 2);
- double x = sqrt((p - cos_t) / (1 - cos_t));
- double y = sqrt((q - cos_t) / (1 - cos_t));
- double z = sqrt((r - cos_t) / (1 - cos_t));
-
- struct rvec v = new_rvec(x, y, z);
- return v;
-}
-*/
-
static double matrix_trace(gsl_matrix *a)
{
int i;
@@ -1601,7 +1621,8 @@ static unsigned int start_seeds(gsl_matrix **rotation, struct TakeTwoCell *cell)
}
-static void set_gsl_matrix(gsl_matrix *mat, double asx, double asy, double asz,
+static void set_gsl_matrix(gsl_matrix *mat,
+ double asx, double asy, double asz,
double bsx, double bsy, double bsz,
double csx, double csy, double csz)
{
@@ -1630,15 +1651,23 @@ static int generate_rotation_sym_ops(struct TakeTwoCell *ttCell)
gsl_matrix *recip = gsl_matrix_alloc(3, 3);
gsl_matrix *cart = gsl_matrix_alloc(3, 3);
- cell_get_reciprocal(ttCell->cell, &asx, &asy, &asz, &bsx, &bsy,
- &bsz, &csx, &csy, &csz);
- set_gsl_matrix(recip, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
+ cell_get_reciprocal(ttCell->cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+
+ set_gsl_matrix(recip, asx, asy, asz,
+ asx, bsy, bsz,
+ csx, csy, csz);
+
+ cell_get_cartesian(ttCell->cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
- cell_get_cartesian(ttCell->cell, &asx, &asy, &asz, &bsx, &bsy,
- &bsz, &csx, &csy, &csz);
+ set_gsl_matrix(cart, asx, bsx, csx,
+ asy, bsy, csy,
+ asz, bsz, csz);
- set_gsl_matrix(cart, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
int i, j, k;
int numOps = num_equivs(rawList, NULL);
@@ -2058,7 +2087,7 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts,
tp->numPrevs++;
/* Prepare the solution for CrystFEL friendliness */
- result = transform_cell_gsl(cell, solution);
+ result = cell_post_smiley_face(cell, solution);
cleanup_taketwo_cell(&ttCell);
return result;
diff --git a/libcrystfel/src/xgandalf.c b/libcrystfel/src/xgandalf.c
index 5f31df33..51819956 100644
--- a/libcrystfel/src/xgandalf.c
+++ b/libcrystfel/src/xgandalf.c
@@ -46,7 +46,7 @@ struct xgandalf_private_data {
IndexingMethod indm;
UnitCell *cellTemplate;
Lattice_t sampleRealLattice_A; //same as cellTemplate
- UnitCellTransformation *uncenteringTransformation;
+ IntegerMatrix *centeringTransformation;
LatticeTransform_t latticeReductionTransform;
};
@@ -114,7 +114,7 @@ int run_xgandalf(struct image *image, void *ipriv)
if(xgandalf_private_data->cellTemplate != NULL){
restoreCell(uc, &xgandalf_private_data->latticeReductionTransform);
- UnitCell *new_cell_trans = cell_transform_inverse(uc, xgandalf_private_data->uncenteringTransformation);
+ UnitCell *new_cell_trans = cell_transform_intmat(uc, xgandalf_private_data->centeringTransformation);
cell_free(uc);
uc = new_cell_trans;
@@ -147,7 +147,7 @@ void *xgandalf_prepare(IndexingMethod *indm, UnitCell *cell,
allocReciprocalPeaks(&(xgandalf_private_data->reciprocalPeaks_1_per_A));
xgandalf_private_data->indm = *indm;
xgandalf_private_data->cellTemplate = NULL;
- xgandalf_private_data->uncenteringTransformation = NULL;
+ xgandalf_private_data->centeringTransformation = NULL;
float tolerance = xgandalf_opts->tolerance;
samplingPitch_t samplingPitch = xgandalf_opts->sampling_pitch;
@@ -157,7 +157,7 @@ void *xgandalf_prepare(IndexingMethod *indm, UnitCell *cell,
xgandalf_private_data->cellTemplate = cell;
- UnitCell* primitiveCell = uncenter_cell(cell, &xgandalf_private_data->uncenteringTransformation);
+ UnitCell* primitiveCell = uncenter_cell(cell, &xgandalf_private_data->centeringTransformation, NULL);
UnitCell *uc = cell_new_from_cell(primitiveCell);
reduceCell(primitiveCell, &xgandalf_private_data->latticeReductionTransform);
@@ -250,8 +250,8 @@ void xgandalf_cleanup(void *pp)
freeReciprocalPeaks(xgandalf_private_data->reciprocalPeaks_1_per_A);
IndexerPlain_delete(xgandalf_private_data->indexer);
- if(xgandalf_private_data->uncenteringTransformation != NULL){
- tfn_free(xgandalf_private_data->uncenteringTransformation);
+ if(xgandalf_private_data->centeringTransformation != NULL){
+ intmat_free(xgandalf_private_data->centeringTransformation);
}
free(xgandalf_private_data);
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-05 18:19:16 +0000