1 files changed, 110 insertions, 0 deletions
diff --git a/crystal.c b/crystal.c
new file mode 100644
index 0000000..065236a
--- /dev/null
+++ b/crystal.c
@@ -0,0 +1,110 @@
+/*
+ * crystal.c
+ *
+ * Crystallographic stuff
+ *
+ * (c) 2009 Thomas White <taw27@cam.ac.uk>
+ *
+ * Triclinator - solve nasty triclinic unit cells
+ *
+ */
+
+
+#include <math.h>
+
+#include "crystal.h"
+#include "util.h"
+
+
+/* Some attempt at preserving sanity */
+#define S11 s11(cell)
+#define S22 s22(cell)
+#define S33 s33(cell)
+#define S12 s12(cell)
+#define S23 s23(cell)
+#define S13 s13(cell)
+
+
+static double s11(Cell cell)
+{
+	return cell.b*cell.b*cell.c*cell.c*sin(cell.al)*sin(cell.al);
+}
+
+static double s22(Cell cell)
+{
+	return cell.a*cell.a*cell.c*cell.c*sin(cell.be)*sin(cell.be);
+}
+
+static double s33(Cell cell)
+{
+	return cell.a*cell.a*cell.b*cell.b*sin(cell.ga)*sin(cell.ga);
+}
+
+static double s12(Cell cell)
+{
+	return cell.a*cell.b*cell.c*cell.c*(cos(cell.al)*cos(cell.be) - cos(cell.ga));
+}
+
+static double s23(Cell cell)
+{
+	return cell.a*cell.a*cell.b*cell.c*(cos(cell.be)*cos(cell.ga) - cos(cell.al));
+}
+
+static double s13(Cell cell)
+{
+	return cell.a*cell.b*cell.b*cell.c*(cos(cell.ga)*cos(cell.al) - cos(cell.be));
+}
+
+
+static double volume(Cell cell)
+{
+	return cell.a*cell.b*cell.c*sqrt(1 - cos(cell.al)*cos(cell.al)
+			    - cos(cell.be)*cos(cell.be)
+			    - cos(cell.ga)*cos(cell.ga)
+			    + 2*cos(cell.al)*cos(cell.be)*cos(cell.ga) );
+}
+
+
+static double dspacing(Cell cell, double h, double k, double l)
+{
+	double sum_S_terms, one_over_Vsq, one_over_LHS;
+
+	sum_S_terms = S11*h*h + S22*k*k + S33*l*l
+			+ 2*S12*h*k + 2*S23*k*l + 2*S13*h*l;
+
+	one_over_Vsq = 1 / pow(volume(cell), 2);
+
+	one_over_LHS = 1/( one_over_Vsq*sum_S_terms );
+
+	return sqrt(one_over_LHS);
+}
+
+
+static double plane_angle(Cell cell, double h1, double k1, double l1,
+				   double h2, double k2, double l2)
+{
+	double sum_S_terms, d1, d2, dd_over_Vsq;
+
+	sum_S_terms = S11*h1*h2 + S22*k1*k2 * S33*l1*l2
+			+ (S23*(k1*l2 + k2*l1))
+			+ (S13*(l1*h2 + l2*h1))
+			+ (S12*(h1*k2 + h2*k1));
+
+	d1 = dspacing(cell, h1, k1, l1);
+	d2 = dspacing(cell, h2, k2, l2);
+	dd_over_Vsq = d1*d2 / pow(volume(cell), 2);
+
+	return acos(dd_over_Vsq * sum_S_terms);
+}
+
+
+/* Return what the measurement 'val' would have been if the cell were 'cell' */
+double crystal_calc(MVal val, Cell cell)
+{
+	if ( is_dspacing(val) ) {
+		return dspacing(cell, val.h1, val.k1, val.l1);
+	} else {
+		return plane_angle(cell, val.h1, val.k1, val.l1,
+					 val.h2, val.k2, val.l2);
+	}
+}