aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--CMakeLists.txt10
-rw-r--r--doc/man/cell_tool.1147
-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
-rw-r--r--src/cell_tool.c576
-rw-r--r--src/get_hkl.c4
-rw-r--r--src/whirligig.c20
-rw-r--r--tests/CMakeLists.txt4
-rw-r--r--tests/centering_check.c22
-rw-r--r--tests/rational_check.c194
-rw-r--r--tests/transformation_check.c227
25 files changed, 3257 insertions, 783 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 8888be0b..38b1723b 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -360,6 +360,15 @@ target_link_libraries(make_pixelmap ${COMMON_LIBRARIES} ${HDF5_C_LIBRARIES})
list(APPEND CRYSTFEL_EXECUTABLES make_pixelmap)
# ----------------------------------------------------------------------
+# cell_tool
+
+set(CELL_TOOL_SOURCES src/cell_tool.c)
+add_executable(cell_tool ${CELL_TOOL_SOURCES})
+target_include_directories(cell_tool PRIVATE ${COMMON_INCLUDES})
+target_link_libraries(cell_tool ${COMMON_LIBRARIES})
+list(APPEND CRYSTFEL_EXECUTABLES cell_tool)
+
+# ----------------------------------------------------------------------
# Install targets
set_target_properties(${CRYSTFEL_EXECUTABLES}
@@ -391,6 +400,7 @@ install(FILES
doc/man/render_hkl.1
doc/man/whirligig.1
doc/man/make_pixelmap.1
+ doc/man/cell_tool.1
DESTINATION ${CMAKE_INSTALL_MANDIR}/man1
)
diff --git a/doc/man/cell_tool.1 b/doc/man/cell_tool.1
new file mode 100644
index 00000000..576943d0
--- /dev/null
+++ b/doc/man/cell_tool.1
@@ -0,0 +1,147 @@
+.\"
+.\" cell_tool man page
+.\"
+.\" Copyright © 2015-2019 Deutsches Elektronen-Synchrotron DESY,
+.\" a research centre of the Helmholtz Association.
+.\"
+.\" Part of CrystFEL - crystallography with a FEL
+.\"
+
+.TH CELL_TOOL 1
+.SH NAME
+cell_tool \- manipulate unit cells
+.SH SYNOPSIS
+.PP
+\fBcell_tool --find-ambi \fImy_structure.cell \fR[\fB-y \fImypointgroup\fR] [\fB--tolerance=\fItols\fR]
+.PP
+\fBcell_tool --uncenter \fImy_structure.cell \fR[\fB-o \fIoutput.cell\fR]
+.PP
+\fBcell_tool --rings \fImy_structure.cell \fR[\fB--highres=\fIangstroms\fR]
+.PP
+\fBcell_tool --compare-cell \fIreference.cell \fImy_structure.cell \fR[\fB--tolerance=\fItols\fR]
+.PP
+\fBcell_tool --transform=\fIop\fR \fImy_structure.cell
+.PP
+\fBcell_tool --help\fI
+
+.SH DESCRIPTION
+\fBcell_tool\fR performs various manipulations on unit cells, including generating power ring positions, comparing one unit cell to another, calculating a primitive unit cell from a centered one and searching for indexing ambiguities.
+.PP
+The unit cell can be given as a CrystFEL unit cell file, or alternatively as a PDB file.
+
+.SH CALCULATING POWDER RING POSITIONS
+.PP
+\fBcell_tool --rings \fImy_structure.cell \fR[--highres=\fIangstroms\fR] [-y \fIpg\fR]
+.PP
+This will generate a list of d-spacings and hkl values for the powder rings given by the unit cell file. Note that screw axis and glide plane absences will not be taken into account, so some rings may be absent depending on the space group.
+.PP
+If you additionally specify the point group using \fB-y\fR (see 'man crystfel' for how to specify point groups), symmetrically equivalent rings will be combined and multiplicities calculated.
+
+.SH GENERATING A PRIMITIVE UNIT CELL
+.PP
+\fBcell_tool --uncenter \fImy_structure.cell \fR[-o \fIoutput.cell\fR]
+.PP
+This will generate a primitive unit cell representing the same lattice as the input. Add the \fB-o\fR option to write the result to a new unit cell file.
+.PP
+There are an infinite number of primitive unit cell for any lattice. This program generates only one of them.
+
+.SH COMPARING UNIT CELLS
+.PP
+\fBcell_tool --compare-cell \fIreference.cell my_structure.cell \fR[\fB--tolerance=\fItols\fR]
+.PP
+The program will compare the two cells, and report if \fImy_structure.cell\fR can be made to look similar to \fIreference.cell\fR applying any transformation matrix.
+.PP
+The tolerance \fItols\fR is given as lengthtol,angtol, in percent and degrees respectively, which will be applied to the real-space unit cell axis lengths and angles.
+
+.SH TRANSFORMING A UNIT CELL
+.PP
+\fBcell_tool --transform=\fIop\fR \fImy_structure.cell
+.PP
+The program will transform the unit cell according to \fIop\fR. Example: \fB--transform=b,c,a\fR means to permute the axes such that the new \fIa\fR axis matches the old \fIb\fR axis, and so on.
+
+.SH FINDING INDEXING AMBIGUITIES
+.PP
+\fBcell_tool --find-ambi \fImy_structure.cell \fR[-y \fIpg\fR] [\fB--tolerance=\fItols\fR]
+.PP
+The program will report all transformation matrices which produce a similar unit cell, to within the specified tolerance. The tolerance \fItols\fR is given as lengthtol,angtol, in percent and degrees respectively, which will be applied to the real-space unit cell axis lengths and angles.
+.PP
+If you additionally give the true symmetry using \fB-y\fR, the program will calculate the ambiguity operators, i.e. the operations which are not symmetry operators of the structure, but which nevertheless leave the lattice looking the same.
+.PP
+\fBExample 1: Merohedral indexing ambiguity in photosystem I\fR
+
+The space group of photosystem I crystals as described by PDB code 1JB0 is P63,
+so the point group is '6':
+
+$ cell_tool --find-ambi 1JB0.pdb -y 6
+.nf
+[...]
+Observed symmetry operations:
+ 1 : hkl -h-k,k,-l -h-k,h,l -h,-k,l -h,h+k,-l
+ -k,-h,-l -k,h+k,l k,-h-k,l k,h,-l h,-h-k,-l
+ h+k,-h,l h+k,-k,-l
+Ambiguity operations:
+ 1 -> 6 : hkl -h-k,k,-l
+.fi
+
+There are 12 reflections which cannot be distinguished between by the lattice alone, but only 6 of those are true symmetry equivalents according to the structure. The transformation describing the indexing ambiguity as follows: "A reflection hkl will be confused with one with indices -h-k,k,-l". Had the point group of the crystals been '622', there would have been no indexing ambiguity (try it!).
+
+.PP
+\fBExample 2: No indexing ambiguity in lysozyme\fR
+
+The space group of lysozyme crystals as described by PDB code 1VDS is P 43 21 2, so the point group is '422':
+
+.nf
+$ cell_tool --find-ambi 1VDS.pdb -y 422
+[...]
+Observed symmetry operations:
+ 1 : hkl -h,-k,l -h,k,-l -k,-h,-l -k,h,l
+ k,-h,l k,h,-l h,-k,-l
+Ambiguity operations:
+ 1 -> 422 : hkl
+.fi
+
+All of the apparently equivalent reflections are true symmetry equivalents according to point group 422, so there is no indexing ambiguity. The only operation produced by left coset decomposition is the identity operation.
+
+.PP
+\fBExample 3: "Accidental" ambiguity in myoglobin\fR
+
+The space group of myoglobin crystals as described by PDB code 3VAU is P2, so the point group is '2'. Note that the "unique axis b" convention has been used for 3VAU (see "man crystfel" for information about specifying point groups):
+
+.nf
+$ cell_tool --find-ambi 3VAU.pdb -y 2_uab
+[...]
+ a b c alpha beta gamma
+ 3.51 2.84 6.29 nm 90.00 105.50 90.00 deg
+[...]
+ a b c alpha beta gamma
+ 3.51 2.84 6.33 nm 90.00 106.84 90.00 deg
+[...]
+Observed symmetry operations:
+ 1 : hkl -h,-k,h+l -h,k,-l h,-k,-h-l
+Ambiguity operations:
+ 1 -> 2 : hkl -h,-k,h+l
+.fi
+
+The transformations '-h,-k,h+l' and 'h,-k,-h-l', which correspond to indexing "diagonally", produce cells which look very similar to the original cell - a difference of only 0.4A and 1.34 degrees. These two transformations are themselves related by a twofold rotation, which is a true symmetry of this crystal structure. There is therefore only one ambiguity transformation. The transformation is strange because it isn't one of the symmetries displayed by a monoclinic lattice in general. This ambiguity has arisen because of of the particular unit cell parameters for this structure.
+
+.SH AUTHOR
+This page was written by Thomas White.
+
+.SH REPORTING BUGS
+Report bugs to <taw@physics.org>, or visit <http://www.desy.de/~twhite/crystfel>.
+
+.SH COPYRIGHT AND DISCLAIMER
+Copyright © 2015-2019 Deutsches Elektronen-Synchrotron DESY, a research centre of the Helmholtz Association.
+.P
+cell-tool, and this manual, are part of CrystFEL.
+.P
+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.
+.P
+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.
+.P
+You should have received a copy of the GNU General Public License along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+
+.SH SEE ALSO
+.BR crystfel (7),
+.BR indexamajig (1),
+.BR get_hkl (1)
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);
}
diff --git a/src/cell_tool.c b/src/cell_tool.c
new file mode 100644
index 00000000..652df885
--- /dev/null
+++ b/src/cell_tool.c
@@ -0,0 +1,576 @@
+/*
+ * cell_tool.c
+ *
+ * Unit cell tool
+ *
+ * Copyright © 2018 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2012-2018 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 <stdarg.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <unistd.h>
+#include <getopt.h>
+#include <assert.h>
+
+#include "cell.h"
+#include "cell-utils.h"
+#include "reflist-utils.h"
+#include "reflist.h"
+
+
+static void show_help(const char *s)
+{
+ printf("Syntax: %s [options]\n\n", s);
+ printf(
+"Unit cell manipulation tool.\n"
+"\n"
+" -h, --help Display this help message.\n"
+" -p, --pdb=<file> Get unit cell from <file> (PDB or CrystFEL format).\n"
+" -o <file> Output unit cell file.\n"
+"\n"
+" Actions:\n"
+" --find-ambi Find indexing ambiguities for the cell.\n"
+" --uncenter Calculate a primitive cell.\n"
+" --rings Calculate powder ring positions.\n"
+" --compare-cell <file> Compare unit cell with cell from <file>.\n"
+" --cell-choices Calculate all three cell choices for monoclinic C cell.\n"
+" --transform=<op> Transform unit cell.\n"
+"\n"
+" -y <pointgroup> Real point group of the structure.\n"
+" --tolerance=<tol> Set the tolerances for cell comparison.\n"
+" Default: 5,1.5 (axis percentage, angle deg).\n"
+" --highres=n Resolution limit (Angstroms) for --rings\n"
+);
+}
+
+
+static int comparecells(UnitCell *cell, const char *comparecell,
+ double ltl, double atl)
+{
+ UnitCell *cell2;
+ RationalMatrix *m;
+
+ cell2 = load_cell_from_file(comparecell);
+ if ( cell2 == NULL ) {
+ ERROR("Failed to load unit cell from '%s'\n", comparecell);
+ return 1;
+ }
+ if ( validate_cell(cell2) > 1 ) {
+ ERROR("Comparison cell is invalid.\n");
+ return 1;
+ }
+ STATUS("------------------> The reference unit cell:\n");
+ cell_print(cell2);
+
+ STATUS("------------------> The comparison results:\n");
+ if ( !compare_reindexed_cell_parameters(cell, cell2, ltl, atl, &m) ) {
+ STATUS("No relationship found between lattices.\n");
+ return 0;
+ } else {
+ UnitCell *trans;
+ STATUS("Relationship found. To become similar to the reference"
+ " cell, the input cell should be transformed by:\n");
+ rtnl_mtx_print(m);
+ STATUS("Transformed version of input unit cell:\n");
+ trans = cell_transform_rational(cell, m);
+ cell_print(trans);
+ cell_free(trans);
+ STATUS("NB transformed cell might not really be triclinic, "
+ "it's just that I don't (yet) know how to work out what "
+ "it is.\n");
+
+ }
+
+ rtnl_mtx_free(m);
+
+ return 0;
+}
+
+
+struct sortmerefl {
+ signed int h;
+ signed int k;
+ signed int l;
+ double resolution;
+ int multi;
+};
+
+
+static int cmpres(const void *av, const void *bv)
+{
+ const struct sortmerefl *a = av;
+ const struct sortmerefl *b = bv;
+ return a->resolution > b->resolution;
+}
+
+
+static int all_rings(UnitCell *cell, SymOpList *sym, double mres)
+{
+ double ax, ay, az;
+ double bx, by, bz;
+ double cx, cy, cz;
+ int hmax, kmax, lmax;
+ signed int h, k, l;
+ RefList *list;
+ int i, n;
+ RefListIterator *iter;
+ Reflection *refl;
+ struct sortmerefl *sortus;
+
+ cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+ hmax = mres * modulus(ax, ay, az);
+ kmax = mres * modulus(bx, by, bz);
+ lmax = mres * modulus(cx, cy, cz);
+ list = reflist_new();
+ for ( h=-hmax; h<=hmax; h++ ) {
+ for ( k=-kmax; k<=kmax; k++ ) {
+ for ( l=-lmax; l<=lmax; l++ ) {
+
+ signed int ha, ka, la;
+
+ if ( forbidden_reflection(cell, h, k, l) ) continue;
+ if ( 2.0*resolution(cell, h, k, l) > mres ) continue;
+
+ if ( sym != NULL ) {
+
+ Reflection *refl;
+
+ get_asymm(sym, h, k, l, &ha, &ka, &la);
+ refl = find_refl(list, ha, ka, la);
+ if ( refl == NULL ) {
+ refl = add_refl(list, ha, ka, la);
+ set_redundancy(refl, 1);
+ } else {
+ set_redundancy(refl, get_redundancy(refl)+1);
+ }
+
+ } else {
+ Reflection *refl;
+ refl = add_refl(list, h, k, l);
+ set_redundancy(refl, 1);
+ }
+
+ }
+ }
+ }
+
+ n = num_reflections(list);
+ sortus = malloc(n*sizeof(struct sortmerefl));
+
+ i = 0;
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) )
+ {
+ signed int h, k, l;
+
+ get_indices(refl, &h, &k, &l);
+
+ sortus[i].h = h;
+ sortus[i].k = k;
+ sortus[i].l = l;
+ sortus[i].resolution = 2.0*resolution(cell, h, k, l); /* one over d */
+ sortus[i].multi = get_redundancy(refl);
+ i++;
+
+ }
+
+ qsort(sortus, n, sizeof(struct sortmerefl), cmpres);
+
+ STATUS("\nAll powder rings up to %f Ångstrøms.\n", 1e+10/mres);
+ STATUS("Note that screw axis or glide plane absences are not "
+ "omitted from this list.\n");
+ STATUS("\n d (Å) 1/d (m^-1) h k l multiplicity\n");
+ STATUS("------------------------------------------------------\n");
+ for ( i=0; i<n; i++ ) {
+ printf("%10.3f %10.3e %4i %4i %4i m = %i\n",
+ 1e10/sortus[i].resolution, sortus[i].resolution,
+ sortus[i].h, sortus[i].k, sortus[i].l,
+ sortus[i].multi);
+ }
+
+ return 0;
+}
+
+
+static int find_ambi(UnitCell *cell, SymOpList *sym, double ltl, double atl)
+{
+ SymOpList *amb;
+ SymOpList *ops;
+ signed int i[9];
+ const int maxorder = 3;
+
+ ops = get_pointgroup("1");
+ if ( ops == NULL ) return 1;
+ set_symmetry_name(ops, "Observed");
+
+ if ( sym == NULL ) {
+ sym = get_pointgroup("1");
+ }
+
+ STATUS("Looking for ambiguities up to %ix each lattice length.\n", maxorder);
+ STATUS("This will take about 30 seconds. Please wait...\n");
+
+ for ( i[0]=-maxorder; i[0]<=+maxorder; i[0]++ ) {
+ for ( i[1]=-maxorder; i[1]<=+maxorder; i[1]++ ) {
+ for ( i[2]=-maxorder; i[2]<=+maxorder; i[2]++ ) {
+ for ( i[3]=-maxorder; i[3]<=+maxorder; i[3]++ ) {
+ for ( i[4]=-maxorder; i[4]<=+maxorder; i[4]++ ) {
+ for ( i[5]=-maxorder; i[5]<=+maxorder; i[5]++ ) {
+ for ( i[6]=-maxorder; i[6]<=+maxorder; i[6]++ ) {
+ for ( i[7]=-maxorder; i[7]<=+maxorder; i[7]++ ) {
+ for ( i[8]=-maxorder; i[8]<=+maxorder; i[8]++ ) {
+
+ UnitCell *nc;
+ IntegerMatrix *m;
+ int j, k;
+ int l = 0;
+
+ m = intmat_new(3, 3);
+ for ( j=0; j<3; j++ ) {
+ for ( k=0; k<3; k++ ) {
+ intmat_set(m, j, k, i[l++]);
+ }
+ }
+
+ if ( intmat_det(m) != +1 ) continue;
+
+ nc = cell_transform_intmat(cell, m);
+
+ if ( compare_cell_parameters(cell, nc, ltl, atl) ) {
+ if ( !intmat_is_identity(m) ) add_symop(ops, m);
+ STATUS("-----------------------------------------------"
+ "-------------------------------------------\n");
+ cell_print(nc);
+ intmat_print(m);
+ } else {
+ intmat_free(m);
+ }
+
+ cell_free(nc);
+
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ STATUS("Observed symmetry operations:\n");
+ describe_symmetry(ops);
+
+ amb = get_ambiguities(ops, sym);
+ if ( amb == NULL ) {
+ STATUS("No ambiguities (or error calculating them)\n");
+ } else {
+ STATUS("Ambiguity operations:\n");
+ describe_symmetry(amb);
+ free_symoplist(amb);
+ }
+
+ free_symoplist(ops);
+
+ return 0;
+}
+
+
+static int uncenter(UnitCell *cell, const char *out_file)
+{
+ UnitCell *cnew;
+ IntegerMatrix *C;
+ RationalMatrix *Ci;
+
+ cnew = uncenter_cell(cell, &C, &Ci);
+
+ STATUS("------------------> The primitive unit cell:\n");
+ cell_print(cnew);
+
+ STATUS("------------------> The centering transformation:\n");
+ intmat_print(C);
+
+ STATUS("------------------> The un-centering transformation:\n");
+ rtnl_mtx_print(Ci);
+
+ if ( out_file != NULL ) {
+ FILE *fh = fopen(out_file, "w");
+ if ( fh == NULL ) {
+ ERROR("Failed to open '%s'\n", out_file);
+ return 1;
+ }
+ write_cell(cnew, fh);
+ fclose(fh);
+ }
+
+ return 0;
+}
+
+
+static int transform(UnitCell *cell, const char *trans_str,
+ const char *out_file)
+{
+ RationalMatrix *trans;
+ Rational det;
+ UnitCell *nc;
+
+ trans = parse_cell_transformation(trans_str);
+ if ( trans == NULL ) {
+ ERROR("Invalid cell transformation '%s'\n", trans_str);
+ return 1;
+ }
+
+ nc = cell_transform_rational(cell, trans);
+
+ STATUS("------------------> The transformation matrix:\n");
+ rtnl_mtx_print(trans);
+ det = rtnl_mtx_det(trans);
+ STATUS("Determinant = %s\n", rtnl_format(det));
+ if ( rtnl_cmp(det, rtnl_zero()) == 0 ) {
+ ERROR("Singular transformation matrix - cannot transform.\n");
+ return 1;
+ }
+
+ STATUS("------------------> The transformed unit cell:\n");
+ cell_print(nc);
+ STATUS("NB transformed cell might not really be triclinic, "
+ "it's just that I don't (yet) know how to work out what "
+ "it is.\n");
+
+ if ( out_file != NULL ) {
+ FILE *fh = fopen(out_file, "w");
+ if ( fh == NULL ) {
+ ERROR("Failed to open '%s'\n", out_file);
+ return 1;
+ }
+ write_cell(nc, fh);
+ fclose(fh);
+ }
+
+ return 0;
+}
+
+
+static int cell_choices(UnitCell *cell)
+{
+ if ( cell_get_lattice_type(cell) != L_MONOCLINIC ) {
+ ERROR("Cell must be monoclinic to use --cell-choices\n");
+ return 1;
+ }
+
+ if ( cell_get_unique_axis(cell) == 'b' ) {
+ transform(cell, "-a-c,b,a", NULL);
+ transform(cell, "c,b,-a-c", NULL);
+ } else {
+ ERROR("Sorry, --cell-choices only supports unique axis b.\n");
+ return 1;
+ }
+
+ return 0;
+}
+
+
+enum {
+ CT_NOTHING,
+ CT_FINDAMBI,
+ CT_UNCENTER,
+ CT_RINGS,
+ CT_COMPARE,
+ CT_CHOICES,
+ CT_TRANSFORM,
+};
+
+
+int main(int argc, char *argv[])
+{
+ int c;
+ char *cell_file = NULL;
+ UnitCell *cell;
+ char *toler = NULL;
+ float ltl = 0.05; /* fraction */
+ float atl = deg2rad(1.5); /* radians */
+ char *sym_str = NULL;
+ SymOpList *sym = NULL;
+ int mode = CT_NOTHING;
+ char *comparecell = NULL;
+ char *out_file = NULL;
+ float highres;
+ double rmax = 1/(2.0e-10);
+ char *trans_str = NULL;
+
+ /* Long options */
+ const struct option longopts[] = {
+ {"help", 0, NULL, 'h'},
+ {"pdb", 1, NULL, 'p'},
+ {"tolerance", 1, NULL, 2},
+ {"output", 1, NULL, 'o'},
+
+ /* Modes of operation */
+ {"find-ambi", 0, &mode, CT_FINDAMBI},
+ {"uncenter", 0, &mode, CT_UNCENTER},
+ {"uncentre", 0, &mode, CT_UNCENTER},
+ {"rings", 0, &mode, CT_RINGS},
+ {"compare-cell", 1, NULL, 3},
+ {"cell-choices", 0, &mode, CT_CHOICES},
+ {"transform", 1, NULL, 4},
+ {"highres", 1, NULL, 5},
+
+ {0, 0, NULL, 0}
+ };
+
+ /* Short options */
+ while ((c = getopt_long(argc, argv, "hp:y:o:",
+ longopts, NULL)) != -1) {
+
+ switch (c) {
+
+ case 'h' :
+ show_help(argv[0]);
+ return 0;
+
+ case 'p' :
+ cell_file = strdup(optarg);
+ break;
+
+ case 'o' :
+ out_file = strdup(optarg);
+ break;
+
+ case 'y' :
+ sym_str = strdup(optarg);
+ break;
+
+ case 2 :
+ toler = strdup(optarg);
+ break;
+
+ case 3 :
+ comparecell = strdup(optarg);
+ mode = CT_COMPARE;
+ break;
+
+ case 4 :
+ trans_str = strdup(optarg);
+ mode = CT_TRANSFORM;
+ break;
+
+ case 5 :
+ if ( sscanf(optarg, "%e", &highres) != 1 ) {
+ ERROR("Invalid value for --highres\n");
+ return 1;
+ }
+ rmax = 1.0 / (highres/1e10);
+ break;
+
+ case 0 :
+ break;
+
+ default :
+ return 1;
+
+ }
+
+ }
+
+ /* If there's a parameter left over, we assume it's the unit cell */
+ if ( (argc > optind) && (cell_file == NULL) ) {
+ cell_file = strdup(argv[optind++]);
+ }
+
+ /* If there's STILL a parameter left over, complain*/
+ if ( argc > optind ) {
+ ERROR("Excess command-line arguments:\n");
+ do {
+ ERROR("'%s'\n", argv[optind++]);
+ } while ( argc > optind );
+ return 1;
+ }
+
+ if ( cell_file == NULL ) {
+ ERROR("You must give a filename for the unit cell PDB file.\n");
+ return 1;
+ }
+ STATUS("Input unit cell: %s\n", cell_file);
+ cell = load_cell_from_file(cell_file);
+ if ( cell == NULL ) {
+ ERROR("Failed to load cell from '%s'\n", cell_file);
+ return 1;
+ }
+ free(cell_file);
+
+ if ( toler != NULL ) {
+ int i;
+ int ncomma = 0;
+ size_t l = strlen(toler);
+ for ( i=0; i<l; i++ ) if ( toler[i] == ',' ) ncomma++;
+ if ( ncomma != 1 ) {
+ ERROR("Invalid parameters for --tolerance. "
+ "Should be: --tolerance=lengthtol,angtol "
+ "(percent,degrees)\n");
+ return 1;
+ }
+ if ( sscanf(toler, "%f,%f", &ltl, &atl) != 2 ) {
+ ERROR("Invalid parameters for --tolerance\n");
+ return 1;
+ }
+ ltl /= 100.0; /* percent to fraction */
+ atl = deg2rad(atl);
+ free(toler);
+ }
+
+ STATUS("------------------> The input unit cell:\n");
+ cell_print(cell);
+
+ if ( validate_cell(cell) > 1 ) {
+ ERROR("Cell is invalid.\n");
+ return 1;
+ }
+
+ if ( sym_str != NULL ) {
+ sym = get_pointgroup(sym_str);
+ if ( sym == NULL ) return 1;
+ free(sym_str);
+ }
+
+ if ( mode == CT_NOTHING ) {
+ ERROR("Please specify mode of operation (see --help)\n");
+ return 1;
+ }
+
+ if ( mode == CT_FINDAMBI ) return find_ambi(cell, sym, ltl, atl);
+ if ( mode == CT_UNCENTER ) return uncenter(cell, out_file);
+ if ( mode == CT_RINGS ) return all_rings(cell, sym, rmax);
+ if ( mode == CT_COMPARE ) return comparecells(cell, comparecell, ltl, atl);
+ if ( mode == CT_TRANSFORM ) return transform(cell, trans_str, out_file);
+ if ( mode == CT_CHOICES ) return cell_choices(cell);
+
+ return 1;
+}
diff --git a/src/get_hkl.c b/src/get_hkl.c
index efd06c94..2a17aab5 100644
--- a/src/get_hkl.c
+++ b/src/get_hkl.c
@@ -580,6 +580,10 @@ int main(int argc, char *argv[])
if ( reindex_str != NULL ) {
reindex = parse_symmetry_operations(reindex_str);
if ( reindex == NULL ) return 1;
+ if ( num_equivs(reindex, NULL) != 1 ) {
+ ERROR("Please provide only ONE reindexing operation\n");
+ return 1;
+ }
set_symmetry_name(reindex, "Reindex");
}
diff --git a/src/whirligig.c b/src/whirligig.c
index f5a435c2..cc2b531e 100644
--- a/src/whirligig.c
+++ b/src/whirligig.c
@@ -84,7 +84,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];
@@ -309,9 +309,9 @@ static IntegerMatrix *try_all(struct window *win, int n1, int n2,
for ( i=0; i<i1->n_crystals; i++ ) {
for ( j=0; j<i2->n_crystals; j++ ) {
- if ( compare_cells(crystal_get_cell(i1->crystals[i]),
- crystal_get_cell(i2->crystals[j]),
- 0.1, deg2rad(5.0), &m) )
+ if ( compare_reindexed_cell_parameters_and_orientation(crystal_get_cell(i1->crystals[i]),
+ crystal_get_cell(i2->crystals[j]),
+ 0.1, deg2rad(5.0), &m) )
{
if ( !crystal_used(win, n1, i)
&& !crystal_used(win, n2, j) )
@@ -364,22 +364,20 @@ static int try_join(struct window *win, int sn)
int j;
Crystal *cr;
UnitCell *ref;
- UnitCellTransformation *tfn;
const int sp = win->join_ptr - 1;
/* Get the appropriately transformed cell from the last crystal in this
* series */
- tfn = tfn_from_intmat(win->mat[sn][sp]);
cr = win->img[sp].crystals[win->ser[sn][sp]];
- ref = cell_transform(crystal_get_cell(cr), tfn);
- tfn_free(tfn);
+ ref = cell_transform_intmat(crystal_get_cell(cr), win->mat[sn][sp]);
for ( j=0; j<win->img[win->join_ptr].n_crystals; j++ ) {
Crystal *cr2;
cr2 = win->img[win->join_ptr].crystals[j];
- if ( compare_cells(ref, crystal_get_cell(cr2),
- 0.1, deg2rad(5.0),
- &win->mat[sn][win->join_ptr]) ) {
+ if ( compare_reindexed_cell_parameters_and_orientation(ref, crystal_get_cell(cr2),
+ 0.1, deg2rad(5.0),
+ &win->mat[sn][win->join_ptr]) )
+ {
win->ser[sn][win->join_ptr] = j;
cell_free(ref);
return 1;
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 02f3ca02..beb141c3 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -83,3 +83,7 @@ if (HAVE_OPENCL)
add_test(gpu_sim_check gpu_sim_check)
endif (HAVE_OPENCL)
+add_executable(rational_check rational_check.c)
+target_include_directories(rational_check PRIVATE ${COMMON_INCLUDES})
+target_link_libraries(rational_check ${COMMON_LIBRARIES})
+add_test(rational_check rational_check)
diff --git a/tests/centering_check.c b/tests/centering_check.c
index 887311dd..cc553dd0 100644
--- a/tests/centering_check.c
+++ b/tests/centering_check.c
@@ -61,7 +61,8 @@ static int check_centering(double a, double b, double c,
{
UnitCell *cell, *cref;
UnitCell *n;
- UnitCellTransformation *t;
+ IntegerMatrix *C;
+ RationalMatrix *Ci;
int fail = 0;
int i;
double asx, asy, asz;
@@ -86,14 +87,19 @@ static int check_centering(double a, double b, double c,
cell_free(cref);
check_cell(cell, "Input");
- n = uncenter_cell(cell, &t);
+ n = uncenter_cell(cell, &C, &Ci);
if ( n != NULL ) {
- STATUS("Transformation was:\n");
- tfn_print(t);
+ STATUS("The back transformation is:\n");
+ intmat_print(C);
if ( check_cell(n, "Output") ) fail = 1;
- if ( !fail ) cell_print(n);
+ if ( intmat_is_identity(C) && ((cen=='P') || (cen=='R')) ) {
+ STATUS("Identity, as expected.\n");
+ } else {
+ cell_print(n);
+ }
} else {
- fail = 1;
+ STATUS("************* Could not uncenter.\n");
+ return 1;
}
cell_get_reciprocal(cell, &asx, &asy, &asz,
@@ -299,10 +305,6 @@ int main(int argc, char *argv[])
fail += check_centering(20e-10, 20e-10, 40e-10, 90.0, 90.0, 120.0,
L_HEXAGONAL, 'H', 'c', rng);
- /* Cubic P */
- fail += check_centering(30e-10, 30e-10, 30e-10, 90.0, 90.0, 90.0,
- L_CUBIC, 'P', '*', rng);
-
/* Cubic I */
fail += check_centering(30e-10, 30e-10, 30e-10, 90.0, 90.0, 90.0,
L_CUBIC, 'I', '*', rng);
diff --git a/tests/rational_check.c b/tests/rational_check.c
new file mode 100644
index 00000000..484264d9
--- /dev/null
+++ b/tests/rational_check.c
@@ -0,0 +1,194 @@
+/*
+ * rational_check.c
+ *
+ * Check that rational numbers work
+ *
+ * Copyright © 2012-2019 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2012-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 <stdarg.h>
+#include <gsl/gsl_rng.h>
+
+#include <rational.h>
+#include <utils.h>
+
+
+static Rational gen_rtnl(gsl_rng *rng, double *dptr, int sz)
+{
+ Rational r;
+ signed int num, den;
+ do {
+ num = gsl_rng_uniform_int(rng, 2*sz) - sz;
+ den = gsl_rng_uniform_int(rng, 2*sz) - sz;
+ } while ( den == 0 );
+ r = rtnl(num, den);
+ *dptr = (double)num/den;
+ return r;
+}
+
+
+static int test_rational(gsl_rng *rng)
+{
+ Rational r1, r2, r3;
+ double rd1, rd2, rd3;
+ char op;
+
+ r1 = gen_rtnl(rng, &rd1, 100);
+ r2 = gen_rtnl(rng, &rd2, 100);
+
+ switch ( gsl_rng_uniform_int(rng, 3) ) {
+
+ case 0:
+ r3 = rtnl_mul(r1, r2);
+ rd3 = rd1*rd2;
+ op = '*';
+ break;
+
+ case 1:
+ r3 = rtnl_div(r1, r2);
+ rd3 = rd1/rd2;
+ op = '/';
+ break;
+
+ case 2:
+ r3 = rtnl_add(r1, r2);
+ rd3 = rd1+rd2;
+ op = '+';
+ break;
+
+ case 3:
+ r3 = rtnl_sub(r1, r2);
+ rd3 = rd1-rd2;
+ op = '-';
+ break;
+
+ default:
+ abort();
+
+ }
+
+ if ( isinf(rd3) && isinf(rtnl_as_double(r3)) ) return 0;
+
+ if ( fabs(rtnl_as_double(r3) - rd3) > 0.001 ) {
+ ERROR("%10s %c %10s = %10s", rtnl_format(r1), op,
+ rtnl_format(r2), rtnl_format(r3));
+ ERROR(" ---> %10f %c %10f = %10f\n", rd1, op, rd2, rd3);
+ return 1;
+ }
+
+ return 0;
+}
+
+
+static int test_rational_matrix(gsl_rng *rng)
+{
+ const int size = 3;
+ RationalMatrix *rm;
+ Rational *rvec;
+ Rational *rans;
+ gsl_matrix *gm;
+ gsl_vector *gvec;
+ gsl_vector *gans;
+ int i, j;
+ int filt;
+ int fail = 0;
+
+ rm = rtnl_mtx_new(size, size);
+ gm = gsl_matrix_alloc(size, size);
+ rvec = malloc(size*sizeof(Rational));
+ gvec = gsl_vector_alloc(size);
+ rans = malloc(size*sizeof(Rational));
+ if ( (rm==NULL) || (gm==NULL) || (rvec==NULL)
+ || (gvec==NULL) || (rans==NULL) ) return 1;
+
+ /* Find a matrix equation which actually works */
+ do {
+ for ( i=0; i<size; i++ ) {
+ double d;
+ Rational r;
+ for ( j=0; j<size; j++ ) {
+ r = gen_rtnl(rng, &d, 5);
+ rtnl_mtx_set(rm, i, j, r);
+ gsl_matrix_set(gm ,i, j, d);
+ }
+ r = gen_rtnl(rng, &d, 5);
+ rvec[i] = r;
+ gsl_vector_set(gvec ,i, d);
+ }
+
+ gans = solve_svd(gvec, gm, &filt, 0);
+ } while ( (gans==NULL) || (filt!=0) || (isnan(gsl_vector_get(gans,0))) );
+
+ transform_fractional_coords_rtnl(rm, rvec, rans);
+
+ for ( i=0; i<size; i++ ) {
+ if ( fabs(rtnl_as_double(rans[i]) - gsl_vector_get(gans, i) > 0.001) )
+ {
+ STATUS("%25s %10f %10f\n", rtnl_format(rans[i]),
+ rtnl_as_double(rans[i]),
+ gsl_vector_get(gans, i));
+ fail = 1;
+ }
+ }
+
+ gsl_matrix_free(gm);
+ rtnl_mtx_free(rm);
+ free(rans);
+ free(rvec);
+ gsl_vector_free(gvec);
+ gsl_vector_free(gans);
+ return fail;
+}
+
+
+int main(int argc, char *argv[])
+{
+ int fail = 0;
+ gsl_rng *rng;
+ int i;
+
+ rng = gsl_rng_alloc(gsl_rng_mt19937);
+ gsl_set_error_handler_off();
+
+ STATUS("Arithmetic test...\n");
+ for ( i=0; i<1000; i++ ) {
+ fail += test_rational(rng);
+ }
+
+ STATUS("Matrix test...\n");
+ for ( i=0; i<1000; i++ ) {
+ fail += test_rational_matrix(rng);
+ }
+
+ gsl_rng_free(rng);
+
+ if ( fail != 0 ) return 1;
+ return 0;
+}
diff --git a/tests/transformation_check.c b/tests/transformation_check.c
index 2e7e7378..2790d6f1 100644
--- a/tests/transformation_check.c
+++ b/tests/transformation_check.c
@@ -37,6 +37,7 @@
#include <cell.h>
#include <cell-utils.h>
+#include <symmetry.h>
#define MAX_REFLS (10*1024)
@@ -134,52 +135,61 @@ static int compare_rvecs(struct rvec *a, int na, struct rvec *b, int nb)
}
-static int check_transformation(UnitCell *cell, UnitCellTransformation *tfn,
- int pred_test, UnitCell *ct)
+static int check_same_reflections(UnitCell *cell, UnitCell *cnew)
{
- UnitCell *cnew, *cback;
- UnitCellTransformation *inv;
- double a[9], b[9];
- int i;
- int fail = 0;
struct rvec *vecs;
struct rvec *tvecs;
int na, nb;
- STATUS("-----------------------\n");
- if ( ct == NULL ) {
- cnew = cell_transform(cell, tfn);
+ /* Check that the two cells predict the same reflections */
+ vecs = all_refls(cell, 1e9, &na);
+ tvecs = all_refls(cnew, 1e9, &nb);
+ if ( compare_rvecs(vecs, na, tvecs, nb)
+ || compare_rvecs(tvecs, nb, vecs, na) )
+ {
+ ERROR("********************************************** ");
+ ERROR("Transformed cell didn't predict the same reflections\n");
+ //printf("---\n");
+ //for ( i=0; i<na; i++ ) {
+ // printf("%e %e %e\n", vecs[i].u, vecs[i].v, vecs[i].w);
+ //}
+ //printf("---\n");
+ //for ( i=0; i<nb; i++ ) {
+ // printf("%e %e %e\n", tvecs[i].u, tvecs[i].v, tvecs[i].w);
+ //}
+ return 1;
} else {
- cnew = ct;
+ STATUS("The cells predict the same reflections.\n");
}
- cback = cell_transform_inverse(cnew, tfn);
+ free(vecs);
+ free(tvecs);
+ return 0;
+}
+
+
+static int check_transformation(UnitCell *cell, IntegerMatrix *tfn,
+ int pred_test)
+{
+ UnitCell *cnew, *cback;
+ double a[9], b[9];
+ int i;
+ int fail = 0;
+
+ STATUS("-----------------------\n");
+ cnew = cell_transform_intmat(cell, tfn);
+ cback = cell_transform_intmat_inverse(cnew, tfn);
+ STATUS("----> Before transformation:\n");
cell_print(cell);
- tfn_print(tfn);
+ STATUS("----> The transformation matrix:\n");
+ intmat_print(tfn);
+ STATUS("----> After transformation:\n");
cell_print(cnew);
+ STATUS("----> After back transformation:\n");
+ cell_print(cback);
if ( pred_test ) {
- /* Check that the two cells predict the same reflections */
- vecs = all_refls(cell, 1e9, &na);
- tvecs = all_refls(cnew, 1e9, &nb);
- if ( compare_rvecs(vecs, na, tvecs, nb)
- || compare_rvecs(tvecs, nb, vecs, na) )
- {
- ERROR("Transformed cell didn't predict the same reflections\n");
- //printf("---\n");
- //for ( i=0; i<na; i++ ) {
- // printf("%e %e %e\n", vecs[i].u, vecs[i].v, vecs[i].w);
- //}
- //printf("---\n");
- //for ( i=0; i<nb; i++ ) {
- // printf("%e %e %e\n", tvecs[i].u, tvecs[i].v, tvecs[i].w);
- //}
- return 1;
- } else {
- STATUS("The cells predict the same reflections.\n");
- }
- free(vecs);
- free(tvecs);
+ check_same_reflections(cell, cnew);
} else {
STATUS("Cells not expected to predict the same reflections.\n");
}
@@ -193,18 +203,18 @@ static int check_transformation(UnitCell *cell, UnitCellTransformation *tfn,
&b[6], &b[7], &b[8]);
for ( i=0; i<9; i++ ) {
if ( !tolerance(a[i], b[i]) ) {
+ //STATUS("%e %e\n", a[i], b[i]);
fail = 1;
- STATUS("%e %e\n", a[i], b[i]);
}
}
if ( fail ) {
+ ERROR("********************************************** ");
ERROR("Original cell not recovered after transformation:\n");
- cell_print(cell);
- tfn_print(tfn);
- inv = tfn_inverse(tfn);
- tfn_print(inv);
+ STATUS("----> After transformation and transformation back:\n");
cell_print(cback);
+ } else {
+ STATUS("The original cell was recovered after inverse transform.\n");
}
return fail;
@@ -214,21 +224,73 @@ static int check_transformation(UnitCell *cell, UnitCellTransformation *tfn,
static int check_uncentering(UnitCell *cell)
{
UnitCell *ct;
- UnitCellTransformation *tr;
+ IntegerMatrix *C;
+ RationalMatrix *Ci;
+ UnitCell *cback;
+ double a[9], b[9];
+ int i;
+ int fail = 0;
+
+ STATUS("-----------------------\n");
+
+ STATUS("----> Before transformation:\n");
+ cell_print_full(cell);
+
+ ct = uncenter_cell(cell, &C, &Ci);
+ if ( ct == NULL ) return 1;
+
+ STATUS("----> The primitive unit cell:\n");
+ cell_print(ct);
- ct = uncenter_cell(cell, &tr);
- return check_transformation(cell, tr, 1, ct);
+ STATUS("----> The matrix to put the centering back:\n");
+ intmat_print(C);
+
+ STATUS("----> The recovered centered cell:\n");
+ cback = cell_transform_intmat(ct, C);
+ cell_print(cback);
+
+ cell_get_cartesian(cell, &a[0], &a[1], &a[2],
+ &a[3], &a[4], &a[5],
+ &a[6], &a[7], &a[8]);
+ cell_get_cartesian(cback, &b[0], &b[1], &b[2],
+ &b[3], &b[4], &b[5],
+ &b[6], &b[7], &b[8]);
+ for ( i=0; i<9; i++ ) {
+ if ( fabs(a[i] - b[i]) > 1e-12 ) {
+ fail = 1;
+ }
+ }
+
+ if ( fail ) {
+ ERROR("********************************************** ");
+ ERROR("Original cell not recovered after back transformation\n");
+ }
+
+ fail += check_same_reflections(cell, ct);
+ cell_free(ct);
+ cell_free(cback);
+
+ return fail;
}
-static int check_identity(UnitCell *cell, UnitCellTransformation *tfn)
+static int check_identity(UnitCell *cell, IntegerMatrix *tfn)
{
UnitCell *cnew;
double a[9], b[9];
int i;
int fail = 0;
- cnew = cell_transform(cell, tfn);
+ STATUS("-----------------------\n");
+
+ cnew = cell_transform_intmat(cell, tfn);
+
+ STATUS("----> Before identity transformation:\n");
+ cell_print(cell);
+ STATUS("----> The identity transformation matrix:\n");
+ intmat_print(tfn);
+ STATUS("----> After identity transformation:\n");
+ cell_print(cnew);
cell_get_cartesian(cell, &a[0], &a[1], &a[2],
&a[3], &a[4], &a[5],
@@ -239,14 +301,15 @@ static int check_identity(UnitCell *cell, UnitCellTransformation *tfn)
for ( i=0; i<9; i++ ) {
if ( !within_tolerance(a[i], b[i], 0.1) ) {
fail = 1;
- STATUS("%e %e\n", a[i], b[i]);
+ //STATUS("%e %e\n", a[i], b[i]);
}
}
if ( fail ) {
- ERROR("Original cell not recovered after transformation:\n");
+ ERROR("********************************************** ");
+ ERROR("Original cell not recovered after identity transformation:\n");
cell_print(cell);
- tfn_print(tfn);
+ intmat_print(tfn);
cell_print(cnew);
}
@@ -258,7 +321,8 @@ int main(int argc, char *argv[])
{
int fail = 0;
UnitCell *cell, *cref;
- UnitCellTransformation *tfn;
+ IntegerMatrix *tfn;
+ IntegerMatrix *part1, *part2;
gsl_rng *rng;
rng = gsl_rng_alloc(gsl_rng_mt19937);
@@ -273,53 +337,52 @@ int main(int argc, char *argv[])
if ( cell == NULL ) return 1;
cell_free(cref);
+ tfn = intmat_identity(3);
+
/* Permutation of axes */
- tfn = tfn_identity();
if ( tfn == NULL ) return 1;
- tfn_combine(tfn, tfn_vector(0,1,0),
- tfn_vector(0,0,1),
- tfn_vector(1,0,0));
- fail += check_transformation(cell, tfn, 1, NULL);
- tfn_free(tfn);
+ intmat_set_all_3x3(tfn, 0,0,1,
+ 1,0,0,
+ 0,1,0);
+ fail += check_transformation(cell, tfn, 1);
/* Doubling of cell in one direction */
- tfn = tfn_identity();
if ( tfn == NULL ) return 1;
- tfn_combine(tfn, tfn_vector(2,0,0),
- tfn_vector(0,1,0),
- tfn_vector(0,0,1));
- fail += check_transformation(cell, tfn, 0, NULL);
- tfn_free(tfn);
-
- /* Diagonal supercell */
- tfn = tfn_identity();
+ intmat_set_all_3x3(tfn, 2,0,0,
+ 0,1,0,
+ 0,0,1);
+ fail += check_transformation(cell, tfn, 0);
+
+ /* Shearing */
if ( tfn == NULL ) return 1;
- tfn_combine(tfn, tfn_vector(1,1,0),
- tfn_vector(0,1,0),
- tfn_vector(0,0,1));
- fail += check_transformation(cell, tfn, 1, NULL);
- tfn_free(tfn);
+ intmat_set_all_3x3(tfn, 1,0,0,
+ 1,1,0,
+ 0,0,1);
+ fail += check_transformation(cell, tfn, 1);
/* Crazy */
- tfn = tfn_identity();
if ( tfn == NULL ) return 1;
- tfn_combine(tfn, tfn_vector(1,1,0),
- tfn_vector(0,1,0),
- tfn_vector(0,1,1));
- fail += check_transformation(cell, tfn, 0, NULL);
- tfn_free(tfn);
+ intmat_set_all_3x3(tfn, 1,0,0,
+ 1,1,1,
+ 0,0,1);
+ fail += check_transformation(cell, tfn, 0);
/* Identity in two parts */
- tfn = tfn_identity();
+ part1 = intmat_identity(3);
+ part2 = intmat_identity(3);
if ( tfn == NULL ) return 1;
- tfn_combine(tfn, tfn_vector(0,0,1),
- tfn_vector(0,1,0),
- tfn_vector(-1,0,0));
- tfn_combine(tfn, tfn_vector(0,0,-1),
- tfn_vector(0,1,0),
- tfn_vector(1,0,0));
+ intmat_set_all_3x3(part1, 0,0,-1,
+ 0,1,0,
+ 1,0,0);
+ intmat_set_all_3x3(part2, 0,0,1,
+ 0,1,0,
+ -1,0,0);
+ tfn = intmat_intmat_mult(part1, part2);
fail += check_identity(cell, tfn);
- tfn_free(tfn);
+ intmat_free(part1);
+ intmat_free(part2);
+
+ intmat_free(tfn);
/* Check some uncentering transformations */
cref = cell_new_from_parameters(50e-10, 50e-10, 50e-10,
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-30 13:23:42 +0000