aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2021-01-11 11:25:05 +0100
committerThomas White <taw@physics.org>2021-01-11 11:25:05 +0100
commitc4325566cc37a59b726afc1034bbfb1ade6a6ca3 (patch)
tree9979ca60141da3081dd3ebb6fa55c3bea97d6da2
parent9d9e819bc535102ca52d6270b6a928e6e70539b6 (diff)
resolution(): Use reciprocal representation instead of crystallographic
This avoids a load of trigonometric functions. In combination with the new UnitCell representation caching, this gives a significant speedup for cases where resolution() is called in a loop.
-rw-r--r--libcrystfel/src/cell-utils.c30
1 files changed, 10 insertions, 20 deletions
diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c
index 6b5d7717..484f9c5c 100644
--- a/libcrystfel/src/cell-utils.c
+++ b/libcrystfel/src/cell-utils.c
@@ -617,28 +617,18 @@ UnitCell *uncenter_cell(UnitCell *in, IntegerMatrix **pC, RationalMatrix **pCi)
/* 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));
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
- 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);
+ cell_get_reciprocal(cell,
+ &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
- return oneoverd / 2;
+ return modulus(h*asx + k*bsx + l*csx,
+ h*asy + k*bsy + l*csy,
+ h*asz + k*bsz + l*csz) / 2.0;
}