aboutsummaryrefslogtreecommitdiff
path: root/src/cell.c
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2011-11-15 12:05:55 +0100
committerThomas White <taw@physics.org>2012-02-22 15:27:40 +0100
commit38089071300b8e04ed42236dd08d9055094fb3b8 (patch)
tree91e1487ac820eb549e7652750867cd4fec039097 /src/cell.c
parent404c612223dbfa0210902ebc5c9226927335aa65 (diff)
Introduce "libcrystfel"
Diffstat (limited to 'src/cell.c')
-rw-r--r--src/cell.c1164
1 files changed, 0 insertions, 1164 deletions
diff --git a/src/cell.c b/src/cell.c
deleted file mode 100644
index c31697bc..00000000
--- a/src/cell.c
+++ /dev/null
@@ -1,1164 +0,0 @@
-/*
- * cell.c
- *
- * Unit Cell Calculations
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-#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 "utils.h"
-#include "image.h"
-
-
-/**
- * SECTION:unitcell
- * @short_description: Unit cell
- * @title: UnitCell
- * @section_id:
- * @see_also:
- * @include: "cell.h"
- * @Image:
- *
- * This structure represents a unit cell.
- */
-
-
-
-/* Weighting factor of lengths relative to angles */
-#define LWEIGHT (10.0e-9)
-
-typedef enum {
- CELL_REP_CRYST,
- CELL_REP_CART,
- CELL_REP_RECIP
-} CellRepresentation;
-
-struct _unitcell {
-
- CellRepresentation rep;
-
- /* Crystallographic representation */
- double a; /* m */
- double b; /* m */
- double c; /* m */
- double alpha; /* Radians */
- double beta; /* Radians */
- double gamma; /* Radians */
-
- /* Cartesian representation */
- double ax; double bx; double cx;
- double ay; double by; double cy;
- double az; double bz; double cz;
-
- /* Cartesian representation of reciprocal axes */
- double axs; double bxs; double cxs;
- double ays; double bys; double cys;
- double azs; double bzs; double czs;
-
- char *pointgroup;
- char *spacegroup;
-};
-
-
-/************************** Setters and Constructors **************************/
-
-
-/**
- * cell_new:
- *
- * Create a new %UnitCell.
- *
- * Returns: the new unit cell, or NULL on failure.
- *
- */
-UnitCell *cell_new()
-{
- UnitCell *cell;
-
- cell = malloc(sizeof(UnitCell));
- if ( cell == NULL ) return NULL;
-
- cell->a = 1.0;
- cell->b = 1.0;
- cell->c = 1.0;
- cell->alpha = M_PI_2;
- cell->beta = M_PI_2;
- cell->gamma = M_PI_2;
-
- cell->rep = CELL_REP_CRYST;
-
- cell->pointgroup = strdup("1");
- cell->spacegroup = strdup("P 1");
-
- return cell;
-}
-
-
-/**
- * cell_free:
- * @cell: A %UnitCell to free.
- *
- * Frees a %UnitCell, and all internal resources concerning that cell.
- *
- */
-void cell_free(UnitCell *cell)
-{
- if ( cell == NULL ) return;
- free(cell->pointgroup);
- free(cell->spacegroup);
- free(cell);
-}
-
-
-void cell_set_parameters(UnitCell *cell, double a, double b, double c,
- double alpha, double beta, double gamma)
-{
- if ( cell == NULL ) return;
-
- cell->a = a;
- cell->b = b;
- cell->c = c;
- cell->alpha = alpha;
- cell->beta = beta;
- cell->gamma = gamma;
-
- cell->rep = CELL_REP_CRYST;
-}
-
-
-void cell_set_cartesian(UnitCell *cell,
- double ax, double ay, double az,
- double bx, double by, double bz,
- double cx, double cy, double cz)
-{
- if ( cell == NULL ) return;
-
- cell->ax = ax; cell->ay = ay; cell->az = az;
- cell->bx = bx; cell->by = by; cell->bz = bz;
- cell->cx = cx; cell->cy = cy; cell->cz = cz;
-
- cell->rep = CELL_REP_CART;
-}
-
-
-void cell_set_cartesian_a(UnitCell *cell, double ax, double ay, double az)
-{
- if ( cell == NULL ) return;
- cell->ax = ax; cell->ay = ay; cell->az = az;
- cell->rep = CELL_REP_CART;
-}
-
-
-void cell_set_cartesian_b(UnitCell *cell, double bx, double by, double bz)
-{
- if ( cell == NULL ) return;
- cell->bx = bx; cell->by = by; cell->bz = bz;
- cell->rep = CELL_REP_CART;
-}
-
-
-void cell_set_cartesian_c(UnitCell *cell, double cx, double cy, double cz)
-{
- if ( cell == NULL ) return;
- cell->cx = cx; cell->cy = cy; cell->cz = cz;
- cell->rep = CELL_REP_CART;
-}
-
-
-UnitCell *cell_new_from_parameters(double a, double b, double c,
- double alpha, double beta, double gamma)
-{
- UnitCell *cell;
-
- cell = cell_new();
- if ( cell == NULL ) return NULL;
-
- cell_set_parameters(cell, a, b, c, alpha, beta, gamma);
-
- return cell;
-}
-
-
-UnitCell *cell_new_from_reciprocal_axes(struct rvec as, struct rvec bs,
- struct rvec cs)
-{
- UnitCell *cell;
-
- cell = cell_new();
- if ( cell == NULL ) return NULL;
-
- cell->axs = as.u; cell->ays = as.v; cell->azs = as.w;
- cell->bxs = bs.u; cell->bys = bs.v; cell->bzs = bs.w;
- cell->cxs = cs.u; cell->cys = cs.v; cell->czs = cs.w;
-
- cell->rep = CELL_REP_RECIP;
-
- return cell;
-}
-
-
-UnitCell *cell_new_from_direct_axes(struct rvec a, struct rvec b, struct rvec c)
-{
- UnitCell *cell;
-
- cell = cell_new();
- if ( cell == NULL ) return NULL;
-
- cell->ax = a.u; cell->ay = a.v; cell->az = a.w;
- cell->bx = b.u; cell->by = b.v; cell->bz = b.w;
- cell->cx = c.u; cell->cy = c.v; cell->cz = c.w;
-
- cell->rep = CELL_REP_CART;
-
- return cell;
-}
-
-
-UnitCell *cell_new_from_cell(UnitCell *orig)
-{
- UnitCell *new;
- double ax, ay, az, bx, by, bz, cx, cy, cz;
-
- new = cell_new();
-
- cell_get_cartesian(orig, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
- cell_set_cartesian(new, ax, ay, az, bx, by, bz, cx, cy, cz);
- cell_set_pointgroup(new, orig->pointgroup);
- cell_set_spacegroup(new, orig->spacegroup);
-
- return new;
-}
-
-
-void cell_set_reciprocal(UnitCell *cell,
- double asx, double asy, double asz,
- double bsx, double bsy, double bsz,
- double csx, double csy, double csz)
-{
- if ( cell == NULL ) return;
-
- cell->axs = asx; cell->ays = asy; cell->azs = asz;
- cell->bxs = bsx; cell->bys = bsy; cell->bzs = bsz;
- cell->cxs = csx; cell->cys = csy; cell->czs = csz;
-
- cell->rep = CELL_REP_RECIP;
-}
-
-
-void cell_set_spacegroup(UnitCell *cell, const char *sym)
-{
- free(cell->spacegroup);
- cell->spacegroup = strdup(sym);
-}
-
-
-void cell_set_pointgroup(UnitCell *cell, const char *sym)
-{
- free(cell->pointgroup);
- cell->pointgroup = strdup(sym);
-}
-
-
-/************************* Getter helper functions ****************************/
-
-static int cell_crystallographic_to_cartesian(UnitCell *cell,
- double *ax, double *ay, double *az,
- double *bx, double *by, double *bz,
- double *cx, double *cy, double *cz)
-{
- double tmp, V, cosalphastar, cstar;
-
- /* Firstly: Get a in terms of x, y and z
- * +a (cryst) is defined to lie along +x (cart) */
- *ax = cell->a;
- *ay = 0.0;
- *az = 0.0;
-
- /* b in terms of x, y and z
- * b (cryst) is defined to lie in the xy (cart) plane */
- *bx = cell->b*cos(cell->gamma);
- *by = cell->b*sin(cell->gamma);
- *bz = 0.0;
-
- tmp = cos(cell->alpha)*cos(cell->alpha)
- + cos(cell->beta)*cos(cell->beta)
- + cos(cell->gamma)*cos(cell->gamma)
- - 2.0*cos(cell->alpha)*cos(cell->beta)*cos(cell->gamma);
- V = cell->a * cell->b * cell->c * sqrt(1.0 - tmp);
-
- cosalphastar = cos(cell->beta)*cos(cell->gamma) - cos(cell->alpha);
- cosalphastar /= sin(cell->beta)*sin(cell->gamma);
-
- cstar = (cell->a * cell->b * sin(cell->gamma))/V;
-
- /* c in terms of x, y and z */
- *cx = cell->c*cos(cell->beta);
- *cy = -cell->c*sin(cell->beta)*cosalphastar;
- *cz = 1.0/cstar;
-
- return 0;
-}
-
-
-/* Why yes, I do enjoy long argument lists...! */
-static int cell_invert(double ax, double ay, double az,
- double bx, double by, double bz,
- double cx, double cy, double cz,
- double *asx, double *asy, double *asz,
- double *bsx, double *bsy, double *bsz,
- double *csx, double *csy, double *csz)
-{
- int s;
- gsl_matrix *m;
- gsl_matrix *inv;
- gsl_permutation *perm;
-
- m = gsl_matrix_alloc(3, 3);
- if ( m == NULL ) {
- ERROR("Couldn't allocate memory for matrix\n");
- 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, 2, 1, bz);
- gsl_matrix_set(m, 2, 2, cz);
-
- /* Invert */
- perm = gsl_permutation_alloc(m->size1);
- if ( perm == NULL ) {
- ERROR("Couldn't allocate permutation\n");
- gsl_matrix_free(m);
- return 1;
- }
- inv = gsl_matrix_alloc(m->size1, m->size2);
- if ( inv == NULL ) {
- ERROR("Couldn't allocate inverse\n");
- gsl_matrix_free(m);
- gsl_permutation_free(perm);
- return 1;
- }
- if ( gsl_linalg_LU_decomp(m, perm, &s) ) {
- ERROR("Couldn't decompose matrix\n");
- gsl_matrix_free(m);
- gsl_permutation_free(perm);
- return 1;
- }
- if ( gsl_linalg_LU_invert(m, perm, inv) ) {
- ERROR("Couldn't invert matrix\n");
- gsl_matrix_free(m);
- gsl_permutation_free(perm);
- return 1;
- }
- gsl_permutation_free(perm);
- gsl_matrix_free(m);
-
- /* Transpose */
- 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);
- *bsz = gsl_matrix_get(inv, 2, 1);
- *csz = gsl_matrix_get(inv, 2, 2);
-
- gsl_matrix_free(inv);
-
- return 0;
-}
-
-
-/********************************** Getters ***********************************/
-
-int cell_get_parameters(UnitCell *cell, double *a, double *b, double *c,
- double *alpha, double *beta, double *gamma)
-{
- double ax, ay, az, bx, by, bz, cx, cy, cz;
-
- if ( cell == NULL ) return 1;
-
- switch ( cell->rep ) {
-
- case CELL_REP_CRYST:
- /* Direct response */
- *a = cell->a;
- *b = cell->b;
- *c = cell->c;
- *alpha = cell->alpha;
- *beta = cell->beta;
- *gamma = cell->gamma;
- return 0;
-
- case CELL_REP_CART:
- /* Convert cartesian -> crystallographic */
- *a = modulus(cell->ax, cell->ay, cell->az);
- *b = modulus(cell->bx, cell->by, cell->bz);
- *c = modulus(cell->cx, cell->cy, cell->cz);
-
- *alpha = angle_between(cell->bx, cell->by, cell->bz,
- cell->cx, cell->cy, cell->cz);
- *beta = angle_between(cell->ax, cell->ay, cell->az,
- cell->cx, cell->cy, cell->cz);
- *gamma = angle_between(cell->ax, cell->ay, cell->az,
- cell->bx, cell->by, cell->bz);
- return 0;
-
- case CELL_REP_RECIP:
- /* Convert reciprocal -> crystallographic.
- * Start by converting reciprocal -> cartesian */
- cell_invert(cell->axs, cell->ays, cell->azs,
- cell->bxs, cell->bys, cell->bzs,
- cell->cxs, cell->cys, cell->czs,
- &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
-
- /* Now convert cartesian -> crystallographic */
- *a = modulus(ax, ay, az);
- *b = modulus(bx, by, bz);
- *c = modulus(cx, cy, cz);
-
- *alpha = angle_between(bx, by, bz, cx, cy, cz);
- *beta = angle_between(ax, ay, az, cx, cy, cz);
- *gamma = angle_between(ax, ay, az, bx, by, bz);
- return 0;
- }
-
- return 1;
-}
-
-
-int cell_get_cartesian(UnitCell *cell,
- double *ax, double *ay, double *az,
- double *bx, double *by, double *bz,
- double *cx, double *cy, double *cz)
-{
- if ( cell == NULL ) return 1;
-
- switch ( cell->rep ) {
-
- case CELL_REP_CRYST:
- /* Convert crystallographic -> cartesian. */
- return cell_crystallographic_to_cartesian(cell,
- ax, ay, az,
- bx, by, bz,
- cx, cy, cz);
-
- case CELL_REP_CART:
- /* Direct response */
- *ax = cell->ax; *ay = cell->ay; *az = cell->az;
- *bx = cell->bx; *by = cell->by; *bz = cell->bz;
- *cx = cell->cx; *cy = cell->cy; *cz = cell->cz;
- return 0;
-
- case CELL_REP_RECIP:
- /* Convert reciprocal -> cartesian */
- return cell_invert(cell->axs, cell->ays, cell->azs,
- cell->bxs, cell->bys, cell->bzs,
- cell->cxs, cell->cys, cell->czs,
- ax, ay, az, bx, by, bz, cx, cy, cz);
-
- }
-
- return 1;
-}
-
-
-int cell_get_reciprocal(UnitCell *cell,
- double *asx, double *asy, double *asz,
- double *bsx, double *bsy, double *bsz,
- double *csx, double *csy, double *csz)
-{
- int r;
- double ax, ay, az, bx, by, bz, cx, cy, cz;
- if ( cell == NULL ) return 1;
-
- switch ( cell->rep ) {
-
- case CELL_REP_CRYST:
- /* Convert crystallographic -> reciprocal */
- r = cell_crystallographic_to_cartesian(cell,
- &ax, &ay, &az,
- &bx, &by, &bz,
- &cx, &cy, &cz);
- if ( r ) return r;
- return cell_invert(ax, ay, az,bx, by, bz, cx, cy, cz,
- asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
-
- case CELL_REP_CART:
- /* Convert cartesian -> reciprocal */
- cell_invert(cell->ax, cell->ay, cell->az,
- cell->bx, cell->by, cell->bz,
- cell->cx, cell->cy, cell->cz,
- asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
- return 0;
-
- case CELL_REP_RECIP:
- /* Direct response */
- *asx = cell->axs; *asy = cell->ays; *asz = cell->azs;
- *bsx = cell->bxs; *bsy = cell->bys; *bsz = cell->bzs;
- *csx = cell->cxs; *csy = cell->cys; *csz = cell->czs;
- return 0;
-
- }
-
- return 1;
-}
-
-
-const char *cell_get_pointgroup(UnitCell *cell)
-{
- return cell->pointgroup;
-}
-
-
-const char *cell_get_spacegroup(UnitCell *cell)
-{
- return cell->spacegroup;
-}
-
-
-
-
-
-/********************************* Utilities **********************************/
-
-static const char *cell_rep(UnitCell *cell)
-{
- switch ( cell->rep ) {
- case CELL_REP_CRYST:
- return "crystallographic, direct space";
- case CELL_REP_CART:
- return "cartesian, direct space";
- case CELL_REP_RECIP:
- return "cartesian, reciprocal space";
- }
-
- 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;
-}
-
-
-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;
-
- 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,
- 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 ltl = 5.0; /* percent */
- float angtol = deg2rad(1.5);
-
- 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, ltl) )
- continue;
-
- 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 = fabs(lengths[i] - tlen);
- if ( ncand[i] == MAX_CAND ) {
- ERROR("Too many candidates\n");
- } else {
- 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, clen;
- 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);
- clen = modulus(cx, cy, cz);
-
- /* 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 cell_set_pointgroup_from_pdb(UnitCell *cell, const char *sym)
-{
- char *new = NULL;
-
- if ( strcmp(sym, "P 1") == 0 ) new = "1";
- if ( strcmp(sym, "P 63") == 0 ) new = "6";
- if ( strcmp(sym, "P 21 21 21") == 0 ) new = "222";
- if ( strcmp(sym, "P 2 2 2") == 0 ) new = "222";
- if ( strcmp(sym, "P 43 21 2") == 0 ) new = "422";
-
- if ( new != NULL ) {
- cell_set_pointgroup(cell, new);
- } else {
- ERROR("Can't determine point group for '%s'\n", sym);
- }
-}
-
-
-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];
- char *sym;
- 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';
-
- 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));
-
- if ( strlen(line) > 65 ) {
- sym = strndup(line+55, 10);
- notrail(sym);
- cell_set_pointgroup_from_pdb(cell, sym);
- cell_set_spacegroup(cell, sym);
- free(sym);
- } else {
- cell_set_pointgroup_from_pdb(cell, "P 1");
- cell_set_spacegroup(cell, "P 1");
- ERROR("CRYST1 line without space group.\n");
- }
-
- break; /* Done */
- }
-
- } while ( rval != NULL );
-
- fclose(fh);
-
- return cell;
-}
-
-
-#ifdef GSL_FUDGE
-/* Force the linker to bring in CBLAS to make GSL happy */
-void cell_fudge_gslcblas()
-{
- STATUS("%p\n", cblas_sgemm);
-}
-#endif
-
-
-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;
-}