aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel
diff options
context:
space:
mode:
Diffstat (limited to 'libcrystfel')
-rw-r--r--libcrystfel/Makefile.am5
-rw-r--r--libcrystfel/src/cell-utils.c863
-rw-r--r--libcrystfel/src/cell-utils.h58
-rw-r--r--libcrystfel/src/cell.c825
-rw-r--r--libcrystfel/src/cell.h22
-rw-r--r--libcrystfel/src/geometry.c1
-rw-r--r--libcrystfel/src/index.c1
-rw-r--r--libcrystfel/src/peaks.c1
-rw-r--r--libcrystfel/src/reax.c1
-rw-r--r--libcrystfel/src/reflist-utils.c1
10 files changed, 933 insertions, 845 deletions
diff --git a/libcrystfel/Makefile.am b/libcrystfel/Makefile.am
index d55f6e74..28d483ac 100644
--- a/libcrystfel/Makefile.am
+++ b/libcrystfel/Makefile.am
@@ -7,7 +7,8 @@ libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \
src/beam-parameters.c src/geometry.c src/statistics.c \
src/symmetry.c src/stream.c src/peaks.c \
src/reflist-utils.c src/filters.c \
- src/render.c src/index.c src/dirax.c src/mosflm.c
+ src/render.c src/index.c src/dirax.c src/mosflm.c \
+ src/cell-utils.c
if HAVE_FFTW
libcrystfel_la_SOURCES += src/reax.c
@@ -22,7 +23,7 @@ libcrystfel_la_include_HEADERS = src/beam-parameters.h src/hdf5-file.h \
src/geometry.h src/peaks.h src/stream.h \
src/render.h src/index.h src/image.h \
src/filters.h src/dirax.h src/mosflm.h \
- src/index-priv.h src/reax.h
+ src/index-priv.h src/reax.h src/cell-utils.h
INCLUDES = "-I$(top_srcdir)/data"
AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" -I$(top_builddir)/lib -Wall
diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c
new file mode 100644
index 00000000..c30d9a03
--- /dev/null
+++ b/libcrystfel/src/cell-utils.c
@@ -0,0 +1,863 @@
+/*
+ * cell-utils.c
+ *
+ * Unit Cell utility functions
+ *
+ * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ * Copyright © 2012 Lorenzo Galli
+ *
+ * Authors:
+ * 2009-2012 Thomas White <taw@physics.org>
+ * 2012 Lorenzo Galli
+ *
+ * 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 <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
+
+#include "cell.h"
+#include "cell-utils.h"
+#include "utils.h"
+#include "image.h"
+
+
+/* Weighting factor of lengths relative to angles */
+#define LWEIGHT (10.0e-9)
+
+
+/**
+ * cell_rotate:
+ * @in: A %UnitCell to rotate
+ * @quat: A %quaternion
+ *
+ * Rotate a %UnitCell using a %quaternion.
+ *
+ * Returns: a newly allocated rotated copy of @in.
+ *
+ */
+UnitCell *cell_rotate(UnitCell *in, struct quaternion quat)
+{
+ struct rvec a, b, c;
+ struct rvec an, bn, cn;
+ UnitCell *out = cell_new_from_cell(in);
+
+ cell_get_cartesian(in, &a.u, &a.v, &a.w,
+ &b.u, &b.v, &b.w,
+ &c.u, &c.v, &c.w);
+
+ an = quat_rot(a, quat);
+ bn = quat_rot(b, quat);
+ cn = quat_rot(c, quat);
+
+ cell_set_cartesian(out, an.u, an.v, an.w,
+ bn.u, bn.v, bn.w,
+ cn.u, cn.v, cn.w);
+
+ return out;
+}
+
+
+static const char *str_lattice(LatticeType l)
+{
+ switch ( l )
+ {
+ case L_TRICLINIC : return "triclinic";
+ case L_MONOCLINIC : return "monoclinic";
+ case L_ORTHORHOMBIC : return "orthorhombic";
+ case L_TETRAGONAL : return "tetragonal";
+ case L_RHOMBOHEDRAL : return "rhombohedral";
+ case L_HEXAGONAL : return "hexagonal";
+ case L_CUBIC : return "cubic";
+ }
+
+ return "unknown lattice";
+}
+
+
+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;
+
+ STATUS("%s %c\n", str_lattice(cell_get_lattice_type(cell)),
+ cell_get_centering(cell));
+
+ cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma);
+
+ STATUS(" a b c alpha beta gamma\n");
+ STATUS("%5.2f %5.2f %5.2f nm %6.2f %6.2f %6.2f deg\n",
+ a*1e9, b*1e9, c*1e9,
+ rad2deg(alpha), rad2deg(beta), rad2deg(gamma));
+
+ cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+
+ STATUS("a = %10.3e %10.3e %10.3e m\n", ax, ay, az);
+ STATUS("b = %10.3e %10.3e %10.3e m\n", bx, by, bz);
+ STATUS("c = %10.3e %10.3e %10.3e m\n", cx, cy, cz);
+
+ cell_get_reciprocal(cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+
+ STATUS("astar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n",
+ asx, asy, asz, modulus(asx, asy, asz));
+ STATUS("bstar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n",
+ bsx, bsy, bsz, modulus(bsx, bsy, bsz));
+ STATUS("cstar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n",
+ csx, csy, csz, modulus(csx, csy, csz));
+
+ STATUS("Point group: %s\n", cell_get_pointgroup(cell));
+ STATUS("Cell representation is %s.\n", cell_rep(cell));
+}
+
+
+#define MAX_CAND (1024)
+
+static int right_handed(struct rvec a, struct rvec b, struct rvec c)
+{
+ struct rvec aCb;
+ double aCb_dot_c;
+
+ /* "a" cross "b" */
+ aCb.u = a.v*b.w - a.w*b.v;
+ aCb.v = - (a.u*b.w - a.w*b.u);
+ aCb.w = a.u*b.v - a.v*b.u;
+
+ /* "a cross b" dot "c" */
+ aCb_dot_c = aCb.u*c.u + aCb.v*c.v + aCb.w*c.w;
+
+ if ( aCb_dot_c > 0.0 ) return 1;
+ return 0;
+}
+
+
+struct cvec {
+ struct rvec vec;
+ float na;
+ float nb;
+ float nc;
+ float fom;
+};
+
+
+static int same_vector(struct cvec a, struct cvec b)
+{
+ if ( a.na != b.na ) return 0;
+ if ( a.nb != b.nb ) return 0;
+ if ( a.nc != b.nc ) return 0;
+ return 1;
+}
+
+
+/* Attempt to make 'cell' fit into 'template' somehow */
+UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose,
+ const float *tols, int reduce)
+{
+ signed int n1l, n2l, n3l;
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+ int i, j;
+ double lengths[3];
+ double angles[3];
+ struct cvec *cand[3];
+ UnitCell *new_cell = NULL;
+ float best_fom = +999999999.9; /* Large number.. */
+ int ncand[3] = {0,0,0};
+ signed int ilow, ihigh;
+ float angtol = deg2rad(tols[3]);
+
+ if ( cell_get_reciprocal(template, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz) ) {
+ ERROR("Couldn't get reciprocal cell for template.\n");
+ return NULL;
+ }
+
+ lengths[0] = modulus(asx, asy, asz);
+ lengths[1] = modulus(bsx, bsy, bsz);
+ lengths[2] = modulus(csx, csy, csz);
+
+ angles[0] = angle_between(bsx, bsy, bsz, csx, csy, csz);
+ angles[1] = angle_between(asx, asy, asz, csx, csy, csz);
+ angles[2] = angle_between(asx, asy, asz, bsx, bsy, bsz);
+
+ cand[0] = malloc(MAX_CAND*sizeof(struct cvec));
+ cand[1] = malloc(MAX_CAND*sizeof(struct cvec));
+ cand[2] = malloc(MAX_CAND*sizeof(struct cvec));
+
+ if ( cell_get_reciprocal(cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz) ) {
+ ERROR("Couldn't get reciprocal cell.\n");
+ return NULL;
+ }
+
+ if ( reduce ) {
+ ilow = -2; ihigh = 4;
+ } else {
+ ilow = 0; ihigh = 1;
+ }
+
+ /* Negative values mean 1/n, positive means n, zero means zero */
+ for ( n1l=ilow; n1l<=ihigh; n1l++ ) {
+ for ( n2l=ilow; n2l<=ihigh; n2l++ ) {
+ for ( n3l=ilow; n3l<=ihigh; n3l++ ) {
+
+ float n1, n2, n3;
+ signed int b1, b2, b3;
+
+ n1 = (n1l>=0) ? (n1l) : (1.0/n1l);
+ n2 = (n2l>=0) ? (n2l) : (1.0/n2l);
+ n3 = (n3l>=0) ? (n3l) : (1.0/n3l);
+
+ if ( !reduce ) {
+ if ( n1l + n2l + n3l > 1 ) continue;
+ }
+
+ /* 'bit' values can be +1 or -1 */
+ for ( b1=-1; b1<=1; b1+=2 ) {
+ for ( b2=-1; b2<=1; b2+=2 ) {
+ for ( b3=-1; b3<=1; b3+=2 ) {
+
+ double tx, ty, tz;
+ double tlen;
+ int i;
+
+ n1 *= b1; n2 *= b2; n3 *= b3;
+
+ tx = n1*asx + n2*bsx + n3*csx;
+ ty = n1*asy + n2*bsy + n3*csy;
+ tz = n1*asz + n2*bsz + n3*csz;
+ tlen = modulus(tx, ty, tz);
+
+ /* Test modulus for agreement with moduli of template */
+ for ( i=0; i<3; i++ ) {
+
+ if ( !within_tolerance(lengths[i], tlen,
+ tols[i]) )
+ {
+ continue;
+ }
+
+ if ( ncand[i] == MAX_CAND ) {
+ ERROR("Too many cell candidates - ");
+ ERROR("consider tightening the unit ");
+ ERROR("cell tolerances.\n");
+ } else {
+
+ double fom;
+
+ fom = fabs(lengths[i] - tlen);
+
+ cand[i][ncand[i]].vec.u = tx;
+ cand[i][ncand[i]].vec.v = ty;
+ cand[i][ncand[i]].vec.w = tz;
+ cand[i][ncand[i]].na = n1;
+ cand[i][ncand[i]].nb = n2;
+ cand[i][ncand[i]].nc = n3;
+ cand[i][ncand[i]].fom = fom;
+
+ ncand[i]++;
+
+ }
+ }
+
+ }
+ }
+ }
+
+ }
+ }
+ }
+
+ if ( verbose ) {
+ STATUS("Candidates: %i %i %i\n", ncand[0], ncand[1], ncand[2]);
+ }
+
+ for ( i=0; i<ncand[0]; i++ ) {
+ for ( j=0; j<ncand[1]; j++ ) {
+
+ double ang;
+ int k;
+ float fom1;
+
+ if ( same_vector(cand[0][i], cand[1][j]) ) continue;
+
+ /* Measure the angle between the ith candidate for axis 0
+ * and the jth candidate for axis 1 */
+ ang = angle_between(cand[0][i].vec.u, cand[0][i].vec.v,
+ cand[0][i].vec.w, cand[1][j].vec.u,
+ cand[1][j].vec.v, cand[1][j].vec.w);
+
+ /* Angle between axes 0 and 1 should be angle 2 */
+ if ( fabs(ang - angles[2]) > angtol ) continue;
+
+ fom1 = fabs(ang - angles[2]);
+
+ for ( k=0; k<ncand[2]; k++ ) {
+
+ float fom2, fom3;
+
+ if ( same_vector(cand[1][j], cand[2][k]) ) continue;
+
+ /* Measure the angle between the current candidate for
+ * axis 0 and the kth candidate for axis 2 */
+ ang = angle_between(cand[0][i].vec.u, cand[0][i].vec.v,
+ cand[0][i].vec.w, cand[2][k].vec.u,
+ cand[2][k].vec.v, cand[2][k].vec.w);
+
+ /* ... it should be angle 1 ... */
+ if ( fabs(ang - angles[1]) > angtol ) continue;
+
+ fom2 = fom1 + fabs(ang - angles[1]);
+
+ /* Finally, the angle between the current candidate for
+ * axis 1 and the kth candidate for axis 2 */
+ ang = angle_between(cand[1][j].vec.u, cand[1][j].vec.v,
+ cand[1][j].vec.w, cand[2][k].vec.u,
+ cand[2][k].vec.v, cand[2][k].vec.w);
+
+ /* ... it should be angle 0 ... */
+ if ( fabs(ang - angles[0]) > angtol ) continue;
+
+ /* Unit cell must be right-handed */
+ if ( !right_handed(cand[0][i].vec, cand[1][j].vec,
+ cand[2][k].vec) ) continue;
+
+ fom3 = fom2 + fabs(ang - angles[0]);
+ fom3 += LWEIGHT * (cand[0][i].fom + cand[1][j].fom
+ + cand[2][k].fom);
+
+ if ( fom3 < best_fom ) {
+ if ( new_cell != NULL ) free(new_cell);
+ new_cell = cell_new_from_reciprocal_axes(
+ cand[0][i].vec, cand[1][j].vec,
+ cand[2][k].vec);
+ best_fom = fom3;
+ }
+
+ }
+
+ }
+ }
+
+ free(cand[0]);
+ free(cand[1]);
+ free(cand[2]);
+
+ return new_cell;
+}
+
+
+
+UnitCell *match_cell_ab(UnitCell *cell, UnitCell *template)
+{
+ double ax, ay, az;
+ double bx, by, bz;
+ double cx, cy, cz;
+ int i;
+ double lengths[3];
+ int used[3];
+ struct rvec real_a, real_b, real_c;
+ struct rvec params[3];
+ double alen, blen;
+ float ltl = 5.0; /* percent */
+ int have_real_a;
+ int have_real_b;
+ int have_real_c;
+
+ /* Get the lengths to match */
+ if ( cell_get_cartesian(template, &ax, &ay, &az,
+ &bx, &by, &bz,
+ &cx, &cy, &cz) )
+ {
+ ERROR("Couldn't get cell for template.\n");
+ return NULL;
+ }
+ alen = modulus(ax, ay, az);
+ blen = modulus(bx, by, bz);
+
+ /* Get the lengths from the cell and turn them into anonymous vectors */
+ if ( cell_get_cartesian(cell, &ax, &ay, &az,
+ &bx, &by, &bz,
+ &cx, &cy, &cz) )
+ {
+ ERROR("Couldn't get cell.\n");
+ return NULL;
+ }
+ lengths[0] = modulus(ax, ay, az);
+ lengths[1] = modulus(bx, by, bz);
+ lengths[2] = modulus(cx, cy, cz);
+ used[0] = 0; used[1] = 0; used[2] = 0;
+ params[0].u = ax; params[0].v = ay; params[0].w = az;
+ params[1].u = bx; params[1].v = by; params[1].w = bz;
+ params[2].u = cx; params[2].v = cy; params[2].w = cz;
+
+ real_a.u = 0.0; real_a.v = 0.0; real_a.w = 0.0;
+ real_b.u = 0.0; real_b.v = 0.0; real_b.w = 0.0;
+ real_c.u = 0.0; real_c.v = 0.0; real_c.w = 0.0;
+
+ /* Check each vector against a and b */
+ have_real_a = 0;
+ have_real_b = 0;
+ for ( i=0; i<3; i++ ) {
+ if ( within_tolerance(lengths[i], alen, ltl)
+ && !used[i] && !have_real_a )
+ {
+ used[i] = 1;
+ memcpy(&real_a, &params[i], sizeof(struct rvec));
+ have_real_a = 1;
+ }
+ if ( within_tolerance(lengths[i], blen, ltl)
+ && !used[i] && !have_real_b )
+ {
+ used[i] = 1;
+ memcpy(&real_b, &params[i], sizeof(struct rvec));
+ have_real_b = 1;
+ }
+ }
+
+ /* Have we matched both a and b? */
+ if ( !(have_real_a && have_real_b) ) return NULL;
+
+ /* "c" is "the other one" */
+ have_real_c = 0;
+ for ( i=0; i<3; i++ ) {
+ if ( !used[i] ) {
+ memcpy(&real_c, &params[i], sizeof(struct rvec));
+ have_real_c = 1;
+ }
+ }
+
+ if ( !have_real_c ) {
+ ERROR("Huh? Couldn't find the third vector.\n");
+ ERROR("Matches: %i %i %i\n", used[0], used[1], used[2]);
+ return NULL;
+ }
+
+ /* Flip c if not right-handed */
+ if ( !right_handed(real_a, real_b, real_c) ) {
+ real_c.u = -real_c.u;
+ real_c.v = -real_c.v;
+ real_c.w = -real_c.w;
+ }
+
+ return cell_new_from_direct_axes(real_a, real_b, real_c);
+}
+
+
+/* Return sin(theta)/lambda = 1/2d. Multiply by two if you want 1/d */
+double resolution(UnitCell *cell, signed int h, signed int k, signed int l)
+{
+ double a, b, c, alpha, beta, gamma;
+
+ cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma);
+
+ const double Vsq = a*a*b*b*c*c*(1 - cos(alpha)*cos(alpha)
+ - cos(beta)*cos(beta)
+ - cos(gamma)*cos(gamma)
+ + 2*cos(alpha)*cos(beta)*cos(gamma) );
+
+ const double S11 = b*b*c*c*sin(alpha)*sin(alpha);
+ const double S22 = a*a*c*c*sin(beta)*sin(beta);
+ const double S33 = a*a*b*b*sin(gamma)*sin(gamma);
+ const double S12 = a*b*c*c*(cos(alpha)*cos(beta) - cos(gamma));
+ const double S23 = a*a*b*c*(cos(beta)*cos(gamma) - cos(alpha));
+ const double S13 = a*b*b*c*(cos(gamma)*cos(alpha) - cos(beta));
+
+ const double brackets = S11*h*h + S22*k*k + S33*l*l
+ + 2*S12*h*k + 2*S23*k*l + 2*S13*h*l;
+ const double oneoverdsq = brackets / Vsq;
+ const double oneoverd = sqrt(oneoverdsq);
+
+ return oneoverd / 2;
+}
+
+
+static void determine_lattice(UnitCell *cell,
+ const char *as, const char *bs, const char *cs,
+ const char *als, const char *bes, const char *gas)
+{
+ int n_right;
+
+ /* Rhombohedral or cubic? */
+ if ( (strcmp(as, bs) == 0) && (strcmp(as, cs) == 0) ) {
+
+ if ( (strcmp(als, " 90.00") == 0)
+ && (strcmp(bes, " 90.00") == 0)
+ && (strcmp(gas, " 90.00") == 0) )
+ {
+ /* Cubic. Unique axis irrelevant. */
+ cell_set_lattice_type(cell, L_CUBIC);
+ return;
+ }
+
+ if ( (strcmp(als, bes) == 0) && (strcmp(als, gas) == 0) ) {
+ /* Rhombohedral. Unique axis irrelevant. */
+ cell_set_lattice_type(cell, L_RHOMBOHEDRAL);
+ return;
+ }
+
+ }
+
+ if ( (strcmp(als, " 90.00") == 0)
+ && (strcmp(bes, " 90.00") == 0)
+ && (strcmp(gas, " 90.00") == 0) )
+ {
+ if ( strcmp(bs, cs) == 0 ) {
+ /* Tetragonal, unique axis a */
+ cell_set_lattice_type(cell, L_TETRAGONAL);
+ cell_set_unique_axis(cell, 'a');
+ return;
+ }
+
+ if ( strcmp(as, cs) == 0 ) {
+ /* Tetragonal, unique axis b */
+ cell_set_lattice_type(cell, L_TETRAGONAL);
+ cell_set_unique_axis(cell, 'b');
+ return;
+ }
+
+ if ( strcmp(as, bs) == 0 ) {
+ /* Tetragonal, unique axis c */
+ cell_set_lattice_type(cell, L_TETRAGONAL);
+ cell_set_unique_axis(cell, 'c');
+ return;
+ }
+
+ /* Orthorhombic. Unique axis irrelevant, but point group
+ * can have different orientations. */
+ cell_set_lattice_type(cell, L_ORTHORHOMBIC);
+ return;
+ }
+
+ n_right = 0;
+ if ( strcmp(als, " 90.00") == 0 ) n_right++;
+ if ( strcmp(bes, " 90.00") == 0 ) n_right++;
+ if ( strcmp(gas, " 90.00") == 0 ) n_right++;
+
+ /* Hexgonal or monoclinic? */
+ if ( n_right == 2 ) {
+
+ if ( (strcmp(als, " 120.00") == 0)
+ && (strcmp(bs, cs) == 0) )
+ {
+ /* Hexagonal, unique axis a */
+ cell_set_lattice_type(cell, L_HEXAGONAL);
+ cell_set_unique_axis(cell, 'a');
+ return;
+ }
+
+ if ( (strcmp(bes, " 120.00") == 0)
+ && (strcmp(as, cs) == 0) )
+ {
+ /* Hexagonal, unique axis b */
+ cell_set_lattice_type(cell, L_HEXAGONAL);
+ cell_set_unique_axis(cell, 'b');
+ return;
+ }
+
+ if ( (strcmp(gas, " 120.00") == 0)
+ && (strcmp(as, bs) == 0) )
+ {
+ /* Hexagonal, unique axis c */
+ cell_set_lattice_type(cell, L_HEXAGONAL);
+ cell_set_unique_axis(cell, 'c');
+ return;
+ }
+
+ if ( strcmp(als, " 90.00") != 0 ) {
+ /* Monoclinic, unique axis a */
+ cell_set_lattice_type(cell, L_MONOCLINIC);
+ cell_set_unique_axis(cell, 'a');
+ return;
+ }
+
+ if ( strcmp(bes, " 90.00") != 0 ) {
+ /* Monoclinic, unique axis b */
+ cell_set_lattice_type(cell, L_MONOCLINIC);
+ cell_set_unique_axis(cell, 'b');
+ return;
+ }
+
+ if ( strcmp(gas, " 90.00") != 0 ) {
+ /* Monoclinic, unique axis c */
+ cell_set_lattice_type(cell, L_MONOCLINIC);
+ cell_set_unique_axis(cell, 'c');
+ return;
+ }
+ }
+
+ /* Triclinic, unique axis irrelevant. */
+ cell_set_lattice_type(cell, L_TRICLINIC);
+}
+
+
+UnitCell *load_cell_from_pdb(const char *filename)
+{
+ FILE *fh;
+ char *rval;
+ UnitCell *cell = NULL;
+
+ fh = fopen(filename, "r");
+ if ( fh == NULL ) {
+ ERROR("Couldn't open '%s'\n", filename);
+ return NULL;
+ }
+
+ do {
+
+ char line[1024];
+
+ rval = fgets(line, 1023, fh);
+
+ if ( strncmp(line, "CRYST1", 6) == 0 ) {
+
+ float a, b, c, al, be, ga;
+ char as[10], bs[10], cs[10];
+ char als[8], bes[8], gas[8];
+ int r;
+
+ memcpy(as, line+6, 9); as[9] = '\0';
+ memcpy(bs, line+15, 9); bs[9] = '\0';
+ memcpy(cs, line+24, 9); cs[9] = '\0';
+ memcpy(als, line+33, 7); als[7] = '\0';
+ memcpy(bes, line+40, 7); bes[7] = '\0';
+ memcpy(gas, line+47, 7); gas[7] = '\0';
+
+ STATUS("'%s' '%s' '%s'\n", as, bs, cs);
+ STATUS("'%s' '%s' '%s'\n", als, bes, gas);
+
+ r = sscanf(as, "%f", &a);
+ r += sscanf(bs, "%f", &b);
+ r += sscanf(cs, "%f", &c);
+ r += sscanf(als, "%f", &al);
+ r += sscanf(bes, "%f", &be);
+ r += sscanf(gas, "%f", &ga);
+
+ if ( r != 6 ) {
+ STATUS("Couldn't understand CRYST1 line.\n");
+ continue;
+ }
+
+ cell = cell_new_from_parameters(a*1e-10,
+ b*1e-10, c*1e-10,
+ deg2rad(al),
+ deg2rad(be),
+ deg2rad(ga));
+
+ determine_lattice(cell, as, bs, cs, als, bes, gas);
+
+ if ( strlen(line) > 65 ) {
+ cell_set_centering(cell, line[55]);
+ } else {
+ cell_set_pointgroup(cell, "1");
+ ERROR("CRYST1 line without centering.\n");
+ }
+
+ break; /* Done */
+ }
+
+ } while ( rval != NULL );
+
+ fclose(fh);
+
+ /* FIXME: Turn "H" centered cells into "R" cells */
+
+ validate_cell(cell);
+
+ return cell;
+}
+
+
+/* Force the linker to bring in CBLAS to make GSL happy */
+void cell_fudge_gslcblas()
+{
+ STATUS("%p\n", cblas_sgemm);
+}
+
+
+UnitCell *rotate_cell(UnitCell *in, double omega, double phi, double rot)
+{
+ UnitCell *out;
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+ double xnew, ynew, znew;
+
+ cell_get_reciprocal(in, &asx, &asy, &asz, &bsx, &bsy,
+ &bsz, &csx, &csy, &csz);
+
+ /* Rotate by "omega" about +z (parallel to c* and c unless triclinic) */
+ xnew = asx*cos(omega) + asy*sin(omega);
+ ynew = -asx*sin(omega) + asy*cos(omega);
+ znew = asz;
+ asx = xnew; asy = ynew; asz = znew;
+ xnew = bsx*cos(omega) + bsy*sin(omega);
+ ynew = -bsx*sin(omega) + bsy*cos(omega);
+ znew = bsz;
+ bsx = xnew; bsy = ynew; bsz = znew;
+ xnew = csx*cos(omega) + csy*sin(omega);
+ ynew = -csx*sin(omega) + csy*cos(omega);
+ znew = csz;
+ csx = xnew; csy = ynew; csz = znew;
+
+ /* Rotate by "phi" about +x (not parallel to anything specific) */
+ xnew = asx;
+ ynew = asy*cos(phi) + asz*sin(phi);
+ znew = -asy*sin(phi) + asz*cos(phi);
+ asx = xnew; asy = ynew; asz = znew;
+ xnew = bsx;
+ ynew = bsy*cos(phi) + bsz*sin(phi);
+ znew = -bsy*sin(phi) + bsz*cos(phi);
+ bsx = xnew; bsy = ynew; bsz = znew;
+ xnew = csx;
+ ynew = csy*cos(phi) + csz*sin(phi);
+ znew = -csy*sin(phi) + csz*cos(phi);
+ csx = xnew; csy = ynew; csz = znew;
+
+ /* Rotate by "rot" about the new +z (in-plane rotation) */
+ xnew = asx*cos(rot) + asy*sin(rot);
+ ynew = -asx*sin(rot) + asy*cos(rot);
+ znew = asz;
+ asx = xnew; asy = ynew; asz = znew;
+ xnew = bsx*cos(rot) + bsy*sin(rot);
+ ynew = -bsx*sin(rot) + bsy*cos(rot);
+ znew = bsz;
+ bsx = xnew; bsy = ynew; bsz = znew;
+ xnew = csx*cos(rot) + csy*sin(rot);
+ ynew = -csx*sin(rot) + csy*cos(rot);
+ znew = csz;
+ csx = xnew; csy = ynew; csz = znew;
+
+ out = cell_new_from_cell(in);
+ cell_set_reciprocal(out, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
+
+ return out;
+}
+
+
+int cell_is_sensible(UnitCell *cell)
+{
+ double a, b, c, al, be, ga;
+
+ cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga);
+ if ( al + be + ga >= 2.0*M_PI ) return 0;
+ if ( al + be - ga >= 2.0*M_PI ) return 0;
+ if ( al - be + ga >= 2.0*M_PI ) return 0;
+ if ( - al + be + ga >= 2.0*M_PI ) return 0;
+ if ( al + be + ga <= 0.0 ) return 0;
+ if ( al + be - ga <= 0.0 ) return 0;
+ if ( al - be + ga <= 0.0 ) return 0;
+ if ( - al + be + ga <= 0.0 ) return 0;
+ if ( isnan(al) ) return 0;
+ if ( isnan(be) ) return 0;
+ if ( isnan(ga) ) return 0;
+ return 1;
+}
+
+
+static int bravais_lattice(UnitCell *cell)
+{
+ LatticeType lattice = cell_get_lattice_type(cell);
+ char centering = cell_get_centering(cell);
+ char ua = cell_get_unique_axis(cell);
+
+ switch ( centering )
+ {
+ case 'P' :
+ return 1;
+
+ case 'A' :
+ case 'B' :
+ case 'C' :
+ if ( (lattice != L_MONOCLINIC)
+ && (lattice != L_ORTHORHOMBIC) )
+ {
+ return 0;
+ }
+ if ( (ua=='a') && (centering=='A') ) return 1;
+ if ( (ua=='b') && (centering=='B') ) return 1;
+ if ( (ua=='c') && (centering=='C') ) return 1;
+ return 0;
+
+ case 'I' :
+ if ( (lattice == L_ORTHORHOMBIC)
+ || (lattice == L_TETRAGONAL)
+ || (lattice == L_CUBIC) )
+ {
+ return 1;
+ }
+ return 0;
+
+ case 'F' :
+ if ( (lattice == L_ORTHORHOMBIC) || (lattice == L_CUBIC) ) {
+ return 1;
+ }
+ return 0;
+
+ case 'H' :
+ if ( lattice == L_HEXAGONAL ) return 1;
+ return 0;
+
+ default :
+ return 0;
+ }
+}
+
+
+/**
+ * validate_cell:
+ * @cell: A %UnitCell to validate
+ *
+ * Perform some checks for crystallographic validity @cell, such as that the
+ * lattice is a conventional Bravais lattice.
+ * Warnings are printied if any of the checks are failed.
+ *
+ */
+void validate_cell(UnitCell *cell)
+{
+ int err = 0;
+
+ if ( !cell_is_sensible(cell) ) {
+ ERROR("Warning: Unit cell parameters are not sensible.\n");
+ err = 1;
+ }
+
+ if ( !bravais_lattice(cell) ) {
+ ERROR("Warning: Unit cell is not a conventional Bravais"
+ " lattice.\n");
+ err = 1;
+ }
+
+ if ( err ) cell_print(cell);
+}
diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h
new file mode 100644
index 00000000..8d8cfc6d
--- /dev/null
+++ b/libcrystfel/src/cell-utils.h
@@ -0,0 +1,58 @@
+/*
+ * cell-utils.h
+ *
+ * Unit Cell utility functions
+ *
+ * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ * Copyright © 2012 Lorenzo Galli
+ *
+ * Authors:
+ * 2009-2012 Thomas White <taw@physics.org>
+ * 2012 Lorenzo Galli
+ *
+ * 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 CELL_UTILS_H
+#define CELL_UTILS_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+extern double resolution(UnitCell *cell,
+ signed int h, signed int k, signed int l);
+
+extern UnitCell *cell_rotate(UnitCell *in, struct quaternion quat);
+extern UnitCell *rotate_cell(UnitCell *in, double omega, double phi,
+ double rot);
+
+extern void cell_print(UnitCell *cell);
+
+extern UnitCell *match_cell(UnitCell *cell, UnitCell *tempcell, int verbose,
+ const float *ltl, int reduce);
+
+extern UnitCell *match_cell_ab(UnitCell *cell, UnitCell *tempcell);
+
+extern UnitCell *load_cell_from_pdb(const char *filename);
+
+extern int cell_is_sensible(UnitCell *cell);
+
+extern void validate_cell(UnitCell *cell);
+
+#endif /* CELL_UTILS_H */
diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c
index 8f380128..64bcd3f4 100644
--- a/libcrystfel/src/cell.c
+++ b/libcrystfel/src/cell.c
@@ -1,7 +1,7 @@
/*
* cell.c
*
- * Unit Cell Calculations
+ * A class representing a unit cell
*
* Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
@@ -60,10 +60,6 @@
*/
-
-/* Weighting factor of lengths relative to angles */
-#define LWEIGHT (10.0e-9)
-
typedef enum {
CELL_REP_CRYST,
CELL_REP_CART,
@@ -603,11 +599,7 @@ char cell_get_unique_axis(UnitCell *cell)
}
-
-
-/********************************* Utilities **********************************/
-
-static const char *cell_rep(UnitCell *cell)
+const char *cell_rep(UnitCell *cell)
{
switch ( cell->rep ) {
@@ -624,816 +616,3 @@ static const char *cell_rep(UnitCell *cell)
return "unknown";
}
-
-
-/**
- * cell_rotate:
- * @in: A %UnitCell to rotate
- * @quat: A %quaternion
- *
- * Rotate a %UnitCell using a %quaternion.
- *
- * Returns: a newly allocated rotated copy of @in.
- *
- */
-UnitCell *cell_rotate(UnitCell *in, struct quaternion quat)
-{
- struct rvec a, b, c;
- struct rvec an, bn, cn;
- UnitCell *out = cell_new_from_cell(in);
-
- cell_get_cartesian(in, &a.u, &a.v, &a.w,
- &b.u, &b.v, &b.w,
- &c.u, &c.v, &c.w);
-
- an = quat_rot(a, quat);
- bn = quat_rot(b, quat);
- cn = quat_rot(c, quat);
-
- cell_set_cartesian(out, an.u, an.v, an.w,
- bn.u, bn.v, bn.w,
- cn.u, cn.v, cn.w);
-
- return out;
-}
-
-
-static const char *str_lattice(LatticeType l)
-{
- switch ( l )
- {
- case L_TRICLINIC : return "triclinic";
- case L_MONOCLINIC : return "monoclinic";
- case L_ORTHORHOMBIC : return "orthorhombic";
- case L_TETRAGONAL : return "tetragonal";
- case L_RHOMBOHEDRAL : return "rhombohedral";
- case L_HEXAGONAL : return "hexagonal";
- case L_CUBIC : return "cubic";
- }
-
- return "unknown lattice";
-}
-
-
-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;
-
- STATUS("%s %c\n", str_lattice(cell_get_lattice_type(cell)),
- cell_get_centering(cell));
-
- cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma);
-
- STATUS(" a b c alpha beta gamma\n");
- STATUS("%5.2f %5.2f %5.2f nm %6.2f %6.2f %6.2f deg\n",
- a*1e9, b*1e9, c*1e9,
- rad2deg(alpha), rad2deg(beta), rad2deg(gamma));
-
- cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
-
- STATUS("a = %10.3e %10.3e %10.3e m\n", ax, ay, az);
- STATUS("b = %10.3e %10.3e %10.3e m\n", bx, by, bz);
- STATUS("c = %10.3e %10.3e %10.3e m\n", cx, cy, cz);
-
- cell_get_reciprocal(cell, &asx, &asy, &asz,
- &bsx, &bsy, &bsz,
- &csx, &csy, &csz);
-
- STATUS("astar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n",
- asx, asy, asz, modulus(asx, asy, asz));
- STATUS("bstar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n",
- bsx, bsy, bsz, modulus(bsx, bsy, bsz));
- STATUS("cstar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n",
- csx, csy, csz, modulus(csx, csy, csz));
-
- STATUS("Point group: %s\n", cell_get_pointgroup(cell));
- STATUS("Cell representation is %s.\n", cell_rep(cell));
-}
-
-
-#define MAX_CAND (1024)
-
-static int right_handed(struct rvec a, struct rvec b, struct rvec c)
-{
- struct rvec aCb;
- double aCb_dot_c;
-
- /* "a" cross "b" */
- aCb.u = a.v*b.w - a.w*b.v;
- aCb.v = - (a.u*b.w - a.w*b.u);
- aCb.w = a.u*b.v - a.v*b.u;
-
- /* "a cross b" dot "c" */
- aCb_dot_c = aCb.u*c.u + aCb.v*c.v + aCb.w*c.w;
-
- if ( aCb_dot_c > 0.0 ) return 1;
- return 0;
-}
-
-
-struct cvec {
- struct rvec vec;
- float na;
- float nb;
- float nc;
- float fom;
-};
-
-
-static int same_vector(struct cvec a, struct cvec b)
-{
- if ( a.na != b.na ) return 0;
- if ( a.nb != b.nb ) return 0;
- if ( a.nc != b.nc ) return 0;
- return 1;
-}
-
-
-/* Attempt to make 'cell' fit into 'template' somehow */
-UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose,
- const float *tols, int reduce)
-{
- signed int n1l, n2l, n3l;
- double asx, asy, asz;
- double bsx, bsy, bsz;
- double csx, csy, csz;
- int i, j;
- double lengths[3];
- double angles[3];
- struct cvec *cand[3];
- UnitCell *new_cell = NULL;
- float best_fom = +999999999.9; /* Large number.. */
- int ncand[3] = {0,0,0};
- signed int ilow, ihigh;
- float angtol = deg2rad(tols[3]);
-
- if ( cell_get_reciprocal(template, &asx, &asy, &asz,
- &bsx, &bsy, &bsz,
- &csx, &csy, &csz) ) {
- ERROR("Couldn't get reciprocal cell for template.\n");
- return NULL;
- }
-
- lengths[0] = modulus(asx, asy, asz);
- lengths[1] = modulus(bsx, bsy, bsz);
- lengths[2] = modulus(csx, csy, csz);
-
- angles[0] = angle_between(bsx, bsy, bsz, csx, csy, csz);
- angles[1] = angle_between(asx, asy, asz, csx, csy, csz);
- angles[2] = angle_between(asx, asy, asz, bsx, bsy, bsz);
-
- cand[0] = malloc(MAX_CAND*sizeof(struct cvec));
- cand[1] = malloc(MAX_CAND*sizeof(struct cvec));
- cand[2] = malloc(MAX_CAND*sizeof(struct cvec));
-
- if ( cell_get_reciprocal(cell, &asx, &asy, &asz,
- &bsx, &bsy, &bsz,
- &csx, &csy, &csz) ) {
- ERROR("Couldn't get reciprocal cell.\n");
- return NULL;
- }
-
- if ( reduce ) {
- ilow = -2; ihigh = 4;
- } else {
- ilow = 0; ihigh = 1;
- }
-
- /* Negative values mean 1/n, positive means n, zero means zero */
- for ( n1l=ilow; n1l<=ihigh; n1l++ ) {
- for ( n2l=ilow; n2l<=ihigh; n2l++ ) {
- for ( n3l=ilow; n3l<=ihigh; n3l++ ) {
-
- float n1, n2, n3;
- signed int b1, b2, b3;
-
- n1 = (n1l>=0) ? (n1l) : (1.0/n1l);
- n2 = (n2l>=0) ? (n2l) : (1.0/n2l);
- n3 = (n3l>=0) ? (n3l) : (1.0/n3l);
-
- if ( !reduce ) {
- if ( n1l + n2l + n3l > 1 ) continue;
- }
-
- /* 'bit' values can be +1 or -1 */
- for ( b1=-1; b1<=1; b1+=2 ) {
- for ( b2=-1; b2<=1; b2+=2 ) {
- for ( b3=-1; b3<=1; b3+=2 ) {
-
- double tx, ty, tz;
- double tlen;
- int i;
-
- n1 *= b1; n2 *= b2; n3 *= b3;
-
- tx = n1*asx + n2*bsx + n3*csx;
- ty = n1*asy + n2*bsy + n3*csy;
- tz = n1*asz + n2*bsz + n3*csz;
- tlen = modulus(tx, ty, tz);
-
- /* Test modulus for agreement with moduli of template */
- for ( i=0; i<3; i++ ) {
-
- if ( !within_tolerance(lengths[i], tlen,
- tols[i]) )
- {
- continue;
- }
-
- if ( ncand[i] == MAX_CAND ) {
- ERROR("Too many cell candidates - ");
- ERROR("consider tightening the unit ");
- ERROR("cell tolerances.\n");
- } else {
-
- double fom;
-
- fom = fabs(lengths[i] - tlen);
-
- cand[i][ncand[i]].vec.u = tx;
- cand[i][ncand[i]].vec.v = ty;
- cand[i][ncand[i]].vec.w = tz;
- cand[i][ncand[i]].na = n1;
- cand[i][ncand[i]].nb = n2;
- cand[i][ncand[i]].nc = n3;
- cand[i][ncand[i]].fom = fom;
-
- ncand[i]++;
-
- }
-
- }
-
- }
- }
- }
-
- }
- }
- }
-
- if ( verbose ) {
- STATUS("Candidates: %i %i %i\n", ncand[0], ncand[1], ncand[2]);
- }
-
- for ( i=0; i<ncand[0]; i++ ) {
- for ( j=0; j<ncand[1]; j++ ) {
-
- double ang;
- int k;
- float fom1;
-
- if ( same_vector(cand[0][i], cand[1][j]) ) continue;
-
- /* Measure the angle between the ith candidate for axis 0
- * and the jth candidate for axis 1 */
- ang = angle_between(cand[0][i].vec.u, cand[0][i].vec.v,
- cand[0][i].vec.w, cand[1][j].vec.u,
- cand[1][j].vec.v, cand[1][j].vec.w);
-
- /* Angle between axes 0 and 1 should be angle 2 */
- if ( fabs(ang - angles[2]) > angtol ) continue;
-
- fom1 = fabs(ang - angles[2]);
-
- for ( k=0; k<ncand[2]; k++ ) {
-
- float fom2, fom3;
-
- if ( same_vector(cand[1][j], cand[2][k]) ) continue;
-
- /* Measure the angle between the current candidate for
- * axis 0 and the kth candidate for axis 2 */
- ang = angle_between(cand[0][i].vec.u, cand[0][i].vec.v,
- cand[0][i].vec.w, cand[2][k].vec.u,
- cand[2][k].vec.v, cand[2][k].vec.w);
-
- /* ... it should be angle 1 ... */
- if ( fabs(ang - angles[1]) > angtol ) continue;
-
- fom2 = fom1 + fabs(ang - angles[1]);
-
- /* Finally, the angle between the current candidate for
- * axis 1 and the kth candidate for axis 2 */
- ang = angle_between(cand[1][j].vec.u, cand[1][j].vec.v,
- cand[1][j].vec.w, cand[2][k].vec.u,
- cand[2][k].vec.v, cand[2][k].vec.w);
-
- /* ... it should be angle 0 ... */
- if ( fabs(ang - angles[0]) > angtol ) continue;
-
- /* Unit cell must be right-handed */
- if ( !right_handed(cand[0][i].vec, cand[1][j].vec,
- cand[2][k].vec) ) continue;
-
- fom3 = fom2 + fabs(ang - angles[0]);
- fom3 += LWEIGHT * (cand[0][i].fom + cand[1][j].fom
- + cand[2][k].fom);
-
- if ( fom3 < best_fom ) {
- if ( new_cell != NULL ) free(new_cell);
- new_cell = cell_new_from_reciprocal_axes(
- cand[0][i].vec, cand[1][j].vec,
- cand[2][k].vec);
- best_fom = fom3;
- }
-
- }
-
- }
- }
-
- free(cand[0]);
- free(cand[1]);
- free(cand[2]);
-
- return new_cell;
-}
-
-
-
-UnitCell *match_cell_ab(UnitCell *cell, UnitCell *template)
-{
- double ax, ay, az;
- double bx, by, bz;
- double cx, cy, cz;
- int i;
- double lengths[3];
- int used[3];
- struct rvec real_a, real_b, real_c;
- struct rvec params[3];
- double alen, blen;
- float ltl = 5.0; /* percent */
- int have_real_a;
- int have_real_b;
- int have_real_c;
-
- /* Get the lengths to match */
- if ( cell_get_cartesian(template, &ax, &ay, &az,
- &bx, &by, &bz,
- &cx, &cy, &cz) )
- {
- ERROR("Couldn't get cell for template.\n");
- return NULL;
- }
- alen = modulus(ax, ay, az);
- blen = modulus(bx, by, bz);
-
- /* Get the lengths from the cell and turn them into anonymous vectors */
- if ( cell_get_cartesian(cell, &ax, &ay, &az,
- &bx, &by, &bz,
- &cx, &cy, &cz) )
- {
- ERROR("Couldn't get cell.\n");
- return NULL;
- }
- lengths[0] = modulus(ax, ay, az);
- lengths[1] = modulus(bx, by, bz);
- lengths[2] = modulus(cx, cy, cz);
- used[0] = 0; used[1] = 0; used[2] = 0;
- params[0].u = ax; params[0].v = ay; params[0].w = az;
- params[1].u = bx; params[1].v = by; params[1].w = bz;
- params[2].u = cx; params[2].v = cy; params[2].w = cz;
-
- real_a.u = 0.0; real_a.v = 0.0; real_a.w = 0.0;
- real_b.u = 0.0; real_b.v = 0.0; real_b.w = 0.0;
- real_c.u = 0.0; real_c.v = 0.0; real_c.w = 0.0;
-
- /* Check each vector against a and b */
- have_real_a = 0;
- have_real_b = 0;
- for ( i=0; i<3; i++ ) {
- if ( within_tolerance(lengths[i], alen, ltl)
- && !used[i] && !have_real_a )
- {
- used[i] = 1;
- memcpy(&real_a, &params[i], sizeof(struct rvec));
- have_real_a = 1;
- }
- if ( within_tolerance(lengths[i], blen, ltl)
- && !used[i] && !have_real_b )
- {
- used[i] = 1;
- memcpy(&real_b, &params[i], sizeof(struct rvec));
- have_real_b = 1;
- }
- }
-
- /* Have we matched both a and b? */
- if ( !(have_real_a && have_real_b) ) return NULL;
-
- /* "c" is "the other one" */
- have_real_c = 0;
- for ( i=0; i<3; i++ ) {
- if ( !used[i] ) {
- memcpy(&real_c, &params[i], sizeof(struct rvec));
- have_real_c = 1;
- }
- }
-
- if ( !have_real_c ) {
- ERROR("Huh? Couldn't find the third vector.\n");
- ERROR("Matches: %i %i %i\n", used[0], used[1], used[2]);
- return NULL;
- }
-
- /* Flip c if not right-handed */
- if ( !right_handed(real_a, real_b, real_c) ) {
- real_c.u = -real_c.u;
- real_c.v = -real_c.v;
- real_c.w = -real_c.w;
- }
-
- return cell_new_from_direct_axes(real_a, real_b, real_c);
-}
-
-
-/* Return sin(theta)/lambda = 1/2d. Multiply by two if you want 1/d */
-double resolution(UnitCell *cell, signed int h, signed int k, signed int l)
-{
- double a, b, c, alpha, beta, gamma;
-
- cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma);
-
- const double Vsq = a*a*b*b*c*c*(1 - cos(alpha)*cos(alpha)
- - cos(beta)*cos(beta)
- - cos(gamma)*cos(gamma)
- + 2*cos(alpha)*cos(beta)*cos(gamma) );
-
- const double S11 = b*b*c*c*sin(alpha)*sin(alpha);
- const double S22 = a*a*c*c*sin(beta)*sin(beta);
- const double S33 = a*a*b*b*sin(gamma)*sin(gamma);
- const double S12 = a*b*c*c*(cos(alpha)*cos(beta) - cos(gamma));
- const double S23 = a*a*b*c*(cos(beta)*cos(gamma) - cos(alpha));
- const double S13 = a*b*b*c*(cos(gamma)*cos(alpha) - cos(beta));
-
- const double brackets = S11*h*h + S22*k*k + S33*l*l
- + 2*S12*h*k + 2*S23*k*l + 2*S13*h*l;
- const double oneoverdsq = brackets / Vsq;
- const double oneoverd = sqrt(oneoverdsq);
-
- return oneoverd / 2;
-}
-
-
-static void determine_lattice(UnitCell *cell,
- const char *as, const char *bs, const char *cs,
- const char *als, const char *bes, const char *gas)
-{
- int n_right;
-
- /* Rhombohedral or cubic? */
- if ( (strcmp(as, bs) == 0) && (strcmp(as, cs) == 0) ) {
-
- if ( (strcmp(als, " 90.00") == 0)
- && (strcmp(bes, " 90.00") == 0)
- && (strcmp(gas, " 90.00") == 0) )
- {
- /* Cubic. Unique axis irrelevant. */
- cell_set_lattice_type(cell, L_CUBIC);
- return;
- }
-
- if ( (strcmp(als, bes) == 0) && (strcmp(als, gas) == 0) ) {
- /* Rhombohedral. Unique axis irrelevant. */
- cell_set_lattice_type(cell, L_RHOMBOHEDRAL);
- return;
- }
-
- }
-
- if ( (strcmp(als, " 90.00") == 0)
- && (strcmp(bes, " 90.00") == 0)
- && (strcmp(gas, " 90.00") == 0) )
- {
- if ( strcmp(bs, cs) == 0 ) {
- /* Tetragonal, unique axis a */
- cell_set_lattice_type(cell, L_TETRAGONAL);
- cell_set_unique_axis(cell, 'a');
- return;
- }
-
- if ( strcmp(as, cs) == 0 ) {
- /* Tetragonal, unique axis b */
- cell_set_lattice_type(cell, L_TETRAGONAL);
- cell_set_unique_axis(cell, 'b');
- return;
- }
-
- if ( strcmp(as, bs) == 0 ) {
- /* Tetragonal, unique axis c */
- cell_set_lattice_type(cell, L_TETRAGONAL);
- cell_set_unique_axis(cell, 'c');
- return;
- }
-
- /* Orthorhombic. Unique axis irrelevant, but point group
- * can have different orientations. */
- cell_set_lattice_type(cell, L_ORTHORHOMBIC);
- return;
- }
-
- n_right = 0;
- if ( strcmp(als, " 90.00") == 0 ) n_right++;
- if ( strcmp(bes, " 90.00") == 0 ) n_right++;
- if ( strcmp(gas, " 90.00") == 0 ) n_right++;
-
- /* Hexgonal or monoclinic? */
- if ( n_right == 2 ) {
-
- if ( (strcmp(als, " 120.00") == 0)
- && (strcmp(bs, cs) == 0) )
- {
- /* Hexagonal, unique axis a */
- cell_set_lattice_type(cell, L_HEXAGONAL);
- cell_set_unique_axis(cell, 'a');
- return;
- }
-
- if ( (strcmp(bes, " 120.00") == 0)
- && (strcmp(as, cs) == 0) )
- {
- /* Hexagonal, unique axis b */
- cell_set_lattice_type(cell, L_HEXAGONAL);
- cell_set_unique_axis(cell, 'b');
- return;
- }
-
- if ( (strcmp(gas, " 120.00") == 0)
- && (strcmp(as, bs) == 0) )
- {
- /* Hexagonal, unique axis c */
- cell_set_lattice_type(cell, L_HEXAGONAL);
- cell_set_unique_axis(cell, 'c');
- return;
- }
-
- if ( strcmp(als, " 90.00") != 0 ) {
- /* Monoclinic, unique axis a */
- cell_set_lattice_type(cell, L_MONOCLINIC);
- cell_set_unique_axis(cell, 'a');
- return;
- }
-
- if ( strcmp(bes, " 90.00") != 0 ) {
- /* Monoclinic, unique axis b */
- cell_set_lattice_type(cell, L_MONOCLINIC);
- cell_set_unique_axis(cell, 'b');
- return;
- }
-
- if ( strcmp(gas, " 90.00") != 0 ) {
- /* Monoclinic, unique axis c */
- cell_set_lattice_type(cell, L_MONOCLINIC);
- cell_set_unique_axis(cell, 'c');
- return;
- }
- }
-
- /* Triclinic, unique axis irrelevant. */
- cell_set_lattice_type(cell, L_TRICLINIC);
-}
-
-
-UnitCell *load_cell_from_pdb(const char *filename)
-{
- FILE *fh;
- char *rval;
- UnitCell *cell = NULL;
-
- fh = fopen(filename, "r");
- if ( fh == NULL ) {
- ERROR("Couldn't open '%s'\n", filename);
- return NULL;
- }
-
- do {
-
- char line[1024];
-
- rval = fgets(line, 1023, fh);
-
- if ( strncmp(line, "CRYST1", 6) == 0 ) {
-
- float a, b, c, al, be, ga;
- char as[10], bs[10], cs[10];
- char als[8], bes[8], gas[8];
- int r;
-
- memcpy(as, line+6, 9); as[9] = '\0';
- memcpy(bs, line+15, 9); bs[9] = '\0';
- memcpy(cs, line+24, 9); cs[9] = '\0';
- memcpy(als, line+33, 7); als[7] = '\0';
- memcpy(bes, line+40, 7); bes[7] = '\0';
- memcpy(gas, line+47, 7); gas[7] = '\0';
-
- STATUS("'%s' '%s' '%s'\n", as, bs, cs);
- STATUS("'%s' '%s' '%s'\n", als, bes, gas);
-
- r = sscanf(as, "%f", &a);
- r += sscanf(bs, "%f", &b);
- r += sscanf(cs, "%f", &c);
- r += sscanf(als, "%f", &al);
- r += sscanf(bes, "%f", &be);
- r += sscanf(gas, "%f", &ga);
-
- if ( r != 6 ) {
- STATUS("Couldn't understand CRYST1 line.\n");
- continue;
- }
-
- cell = cell_new_from_parameters(a*1e-10,
- b*1e-10, c*1e-10,
- deg2rad(al),
- deg2rad(be),
- deg2rad(ga));
-
- determine_lattice(cell, as, bs, cs, als, bes, gas);
-
- if ( strlen(line) > 65 ) {
- cell_set_centering(cell, line[55]);
- } else {
- cell_set_pointgroup(cell, "1");
- cell_set_spacegroup(cell, "P 1");
- ERROR("CRYST1 line without centering.\n");
- }
-
- break; /* Done */
- }
-
- } while ( rval != NULL );
-
- fclose(fh);
-
- validate_cell(cell);
-
- return cell;
-}
-
-
-/* Force the linker to bring in CBLAS to make GSL happy */
-void cell_fudge_gslcblas()
-{
- STATUS("%p\n", cblas_sgemm);
-}
-
-
-UnitCell *rotate_cell(UnitCell *in, double omega, double phi, double rot)
-{
- UnitCell *out;
- double asx, asy, asz;
- double bsx, bsy, bsz;
- double csx, csy, csz;
- double xnew, ynew, znew;
-
- cell_get_reciprocal(in, &asx, &asy, &asz, &bsx, &bsy,
- &bsz, &csx, &csy, &csz);
-
- /* Rotate by "omega" about +z (parallel to c* and c unless triclinic) */
- xnew = asx*cos(omega) + asy*sin(omega);
- ynew = -asx*sin(omega) + asy*cos(omega);
- znew = asz;
- asx = xnew; asy = ynew; asz = znew;
- xnew = bsx*cos(omega) + bsy*sin(omega);
- ynew = -bsx*sin(omega) + bsy*cos(omega);
- znew = bsz;
- bsx = xnew; bsy = ynew; bsz = znew;
- xnew = csx*cos(omega) + csy*sin(omega);
- ynew = -csx*sin(omega) + csy*cos(omega);
- znew = csz;
- csx = xnew; csy = ynew; csz = znew;
-
- /* Rotate by "phi" about +x (not parallel to anything specific) */
- xnew = asx;
- ynew = asy*cos(phi) + asz*sin(phi);
- znew = -asy*sin(phi) + asz*cos(phi);
- asx = xnew; asy = ynew; asz = znew;
- xnew = bsx;
- ynew = bsy*cos(phi) + bsz*sin(phi);
- znew = -bsy*sin(phi) + bsz*cos(phi);
- bsx = xnew; bsy = ynew; bsz = znew;
- xnew = csx;
- ynew = csy*cos(phi) + csz*sin(phi);
- znew = -csy*sin(phi) + csz*cos(phi);
- csx = xnew; csy = ynew; csz = znew;
-
- /* Rotate by "rot" about the new +z (in-plane rotation) */
- xnew = asx*cos(rot) + asy*sin(rot);
- ynew = -asx*sin(rot) + asy*cos(rot);
- znew = asz;
- asx = xnew; asy = ynew; asz = znew;
- xnew = bsx*cos(rot) + bsy*sin(rot);
- ynew = -bsx*sin(rot) + bsy*cos(rot);
- znew = bsz;
- bsx = xnew; bsy = ynew; bsz = znew;
- xnew = csx*cos(rot) + csy*sin(rot);
- ynew = -csx*sin(rot) + csy*cos(rot);
- znew = csz;
- csx = xnew; csy = ynew; csz = znew;
-
- out = cell_new_from_cell(in);
- cell_set_reciprocal(out, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
-
- return out;
-}
-
-
-int cell_is_sensible(UnitCell *cell)
-{
- double a, b, c, al, be, ga;
-
- cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga);
- if ( al + be + ga >= 2.0*M_PI ) return 0;
- if ( al + be - ga >= 2.0*M_PI ) return 0;
- if ( al - be + ga >= 2.0*M_PI ) return 0;
- if ( - al + be + ga >= 2.0*M_PI ) return 0;
- if ( al + be + ga <= 0.0 ) return 0;
- if ( al + be - ga <= 0.0 ) return 0;
- if ( al - be + ga <= 0.0 ) return 0;
- if ( - al + be + ga <= 0.0 ) return 0;
- if ( isnan(al) ) return 0;
- if ( isnan(be) ) return 0;
- if ( isnan(ga) ) return 0;
- return 1;
-}
-
-
-static int bravais_lattice(UnitCell *cell)
-{
- LatticeType lattice = cell_get_lattice_type(cell);
- char centering = cell_get_centering(cell);
- char ua = cell_get_unique_axis(cell);
-
- switch ( centering )
- {
- case 'P' :
- return 1;
-
- case 'A' :
- case 'B' :
- case 'C' :
- if ( (lattice != L_MONOCLINIC)
- && (lattice != L_ORTHORHOMBIC) )
- {
- return 0;
- }
- if ( (ua=='a') && (centering=='A') ) return 1;
- if ( (ua=='b') && (centering=='B') ) return 1;
- if ( (ua=='c') && (centering=='C') ) return 1;
- return 0;
-
- case 'I' :
- if ( (lattice == L_ORTHORHOMBIC)
- || (lattice == L_TETRAGONAL)
- || (lattice == L_CUBIC) )
- {
- return 1;
- }
- return 0;
-
- case 'F' :
- if ( (lattice == L_ORTHORHOMBIC) || (lattice == L_CUBIC) ) {
- return 1;
- }
- return 0;
-
- case 'H' :
- if ( lattice == L_HEXAGONAL ) return 1;
- return 0;
-
- default :
- return 0;
- }
-}
-
-
-/**
- * validate_cell:
- * @cell: A %UnitCell to validate
- *
- * Perform some checks for crystallographic validity @cell, such as that the
- * lattice is a conventional Bravais lattice.
- * Warnings are printied if any of the checks are failed.
- *
- */
-void validate_cell(UnitCell *cell)
-{
- int err = 0;
-
- if ( !cell_is_sensible(cell) ) {
- ERROR("Warning: Unit cell parameters are not sensible.\n");
- err = 1;
- }
-
- if ( !bravais_lattice(cell) ) {
- ERROR("Warning: Unit cell is not a conventional Bravais"
- " lattice.\n");
- err = 1;
- }
-
- if ( err ) cell_print(cell);
-}
diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h
index 851237cf..e4e9d6ce 100644
--- a/libcrystfel/src/cell.h
+++ b/libcrystfel/src/cell.h
@@ -1,7 +1,7 @@
/*
* cell.h
*
- * Unit Cell Calculations
+ * A class representing a unit cell
*
* Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
@@ -128,24 +128,6 @@ extern void cell_set_centering(UnitCell *cell, char centering);
extern char cell_get_unique_axis(UnitCell *cell);
extern void cell_set_unique_axis(UnitCell *cell, char unique_axis);
-extern double resolution(UnitCell *cell,
- signed int h, signed int k, signed int l);
-
-extern UnitCell *cell_rotate(UnitCell *in, struct quaternion quat);
-extern UnitCell *rotate_cell(UnitCell *in, double omega, double phi,
- double rot);
-
-extern void cell_print(UnitCell *cell);
-
-extern UnitCell *match_cell(UnitCell *cell, UnitCell *tempcell, int verbose,
- const float *ltl, int reduce);
-
-extern UnitCell *match_cell_ab(UnitCell *cell, UnitCell *tempcell);
-
-extern UnitCell *load_cell_from_pdb(const char *filename);
-
-extern int cell_is_sensible(UnitCell *cell);
-
-extern void validate_cell(UnitCell *cell);
+extern const char *cell_rep(UnitCell *cell);
#endif /* CELL_H */
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index 24669aef..091b4fed 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -36,6 +36,7 @@
#include "utils.h"
#include "cell.h"
+#include "cell-utils.h"
#include "image.h"
#include "peaks.h"
#include "beam-parameters.h"
diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c
index 102a5312..7d912902 100644
--- a/libcrystfel/src/index.c
+++ b/libcrystfel/src/index.c
@@ -49,6 +49,7 @@
#include "index-priv.h"
#include "reax.h"
#include "geometry.h"
+#include "cell-utils.h"
/* Base class constructor for unspecialised indexing private data */
diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
index b87eb56e..f7f6c650 100644
--- a/libcrystfel/src/peaks.c
+++ b/libcrystfel/src/peaks.c
@@ -51,6 +51,7 @@
#include "filters.h"
#include "reflist-utils.h"
#include "beam-parameters.h"
+#include "cell-utils.h"
/* Degree of polarisation of X-ray beam */
diff --git a/libcrystfel/src/reax.c b/libcrystfel/src/reax.c
index 5cfa908a..8bd9f29f 100644
--- a/libcrystfel/src/reax.c
+++ b/libcrystfel/src/reax.c
@@ -47,6 +47,7 @@
#include "utils.h"
#include "peaks.h"
#include "cell.h"
+#include "cell-utils.h"
#include "index.h"
#include "index-priv.h"
diff --git a/libcrystfel/src/reflist-utils.c b/libcrystfel/src/reflist-utils.c
index 2e6715f2..f2292929 100644
--- a/libcrystfel/src/reflist-utils.c
+++ b/libcrystfel/src/reflist-utils.c
@@ -35,6 +35,7 @@
#include "reflist.h"
#include "cell.h"
+#include "cell-utils.h"
#include "utils.h"
#include "reflist-utils.h"
#include "symmetry.h"
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-05 20:44:30 +0000