Initial import

8 files changed, 598 insertions, 0 deletions

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..6003c91
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,2 @@
+*.o
+triclinator
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..3878741
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,16 @@
+triclinator: main.o crystal.o util.o
+	gcc main.o crystal.o util.o -o triclinator -lgsl -lgslcblas
+
+main.o:	main.c
+	gcc -g -c -I/usr/include main.c -o main.o
+
+crystal.o: crystal.c
+	gcc -g -c -I/usr/include crystal.c -o crystal.o
+
+util.o: util.c
+	gcc -g -c -I/usr/include util.c -o util.o
+
+clean:
+	rm -f triclinator main.o crystal.o util.o
+
+.PHONY: clean
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);
+	}
+}
diff --git a/crystal.h b/crystal.h
new file mode 100644
index 0000000..5c38da8
--- /dev/null
+++ b/crystal.h
@@ -0,0 +1,33 @@
+/*
+ * crystal.h
+ *
+ * Crystallographic stuff
+ *
+ * (c) 2009 Thomas White <taw27@cam.ac.uk>
+ *
+ * Triclinator - solve nasty triclinic unit cells
+ *
+ */
+
+#ifndef CRYSTAL_H
+#define CRYSTAL_H
+
+#define RAD2DEG(a) ((a)*180/M_PI)
+#define DEG2RAD(a) ((a)*M_PI/180)
+
+#include "util.h"
+
+typedef struct
+{
+	double a;
+	double b;
+	double c;
+	double al;
+	double be;
+	double ga;
+} Cell;
+
+/* Return what the measurement 'val' would have been if the cell were 'cell' */
+extern double crystal_calc(MVal val, Cell cell);
+
+#endif	/* CRYSTAL_H */
diff --git a/main.c b/main.c
new file mode 100644
index 0000000..f5f621b
--- /dev/null
+++ b/main.c
@@ -0,0 +1,334 @@
+/*
+ * main.c
+ *
+ * The Top Level Source File
+ *
+ * (c) 2009 Thomas White <taw27@cam.ac.uk>
+ *
+ * Triclinator - solve nasty triclinic unit cells
+ *
+ */
+
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_multifit_nlin.h>
+#include <math.h>
+
+#include "crystal.h"
+#include "util.h"
+
+#define DSTEP 0.0001
+#define ASTEP 0.00001
+
+struct data {
+	size_t n;	/* Number of measurements */
+	MVal *vals;	/* The measurements themselves */
+};
+
+
+int expb_f(const gsl_vector *x, void *data, gsl_vector *f)
+{
+	size_t n = ((struct data *)data)->n;
+	MVal *vals = ((struct data *) data)->vals;
+	Cell cell;
+	size_t i;
+
+	cell.a = gsl_vector_get(x, 0);
+	cell.b = gsl_vector_get(x, 1);
+	cell.c = gsl_vector_get(x, 2);
+	cell.al = gsl_vector_get(x, 3);
+	cell.be = gsl_vector_get(x, 4);
+	cell.ga = gsl_vector_get(x, 5);
+
+	for ( i=0; i<n; i++ ) {
+		double calc = crystal_calc(vals[i], cell);
+		gsl_vector_set(f, i, (calc - vals[i].meas)/vals[i].esd);
+	}
+
+	return GSL_SUCCESS;
+}
+
+
+int expb_df(const gsl_vector *x, void *data, gsl_matrix *J)
+{
+	size_t n = ((struct data *)data)->n;
+	MVal *vals = ((struct data *) data)->vals;
+	Cell cell;
+	size_t i;
+
+	cell.a = gsl_vector_get(x, 0);
+	cell.b = gsl_vector_get(x, 1);
+	cell.c = gsl_vector_get(x, 2);
+	cell.al = gsl_vector_get(x, 3);
+	cell.be = gsl_vector_get(x, 4);
+	cell.ga = gsl_vector_get(x, 5);
+
+	for ( i=0; i<n; i++ ) {
+
+		double calc, calc2, s;
+		Cell cell2;
+		MVal val;
+
+		val = vals[i];
+
+		calc = crystal_calc(vals[i], cell);
+		s = vals[i].esd;
+
+		cell2 = cell; cell2.a += DSTEP; calc2 = crystal_calc(val, cell2);
+		gsl_matrix_set(J, i, 0, (calc2 - calc)/(DSTEP*s)); /* df/da */
+
+		cell2 = cell; cell2.b += DSTEP; calc2 = crystal_calc(val, cell2);
+		gsl_matrix_set(J, i, 1, (calc2 - calc)/(DSTEP*s)); /* df/db */
+
+		cell2 = cell; cell2.c += DSTEP; calc2 = crystal_calc(val, cell2);
+		gsl_matrix_set(J, i, 2, (calc2 - calc)/(DSTEP*s)); /* df/dc */
+
+		cell2 = cell; cell2.al += ASTEP; calc2 = crystal_calc(val, cell2);
+		gsl_matrix_set(J, i, 3, (calc2 - calc)/(ASTEP*s)); /* df/dal */
+
+		cell2 = cell; cell2.be += ASTEP; calc2 = crystal_calc(val, cell2);
+		gsl_matrix_set(J, i, 4, (calc2 - calc)/(ASTEP*s)); /* df/dbe */
+
+		cell2 = cell; cell2.ga += ASTEP; calc2 = crystal_calc(val, cell2);
+		gsl_matrix_set(J, i, 5, (calc2 - calc)/(ASTEP*s)); /* df/dga */
+
+	}
+
+	return GSL_SUCCESS;
+}
+
+
+int expb_fdf(const gsl_vector *x, void *data, gsl_vector *f, gsl_matrix *J)
+{
+	expb_f(x, data, f);
+	expb_df(x, data, J);
+	return GSL_SUCCESS;
+}
+
+
+void print_state (size_t iter, gsl_multifit_fdfsolver *s)
+{
+	printf ("%3u %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f", iter,
+		gsl_vector_get(s->x, 0),
+		gsl_vector_get(s->x, 1),
+		gsl_vector_get(s->x, 2),
+		gsl_vector_get(s->x, 3),
+		gsl_vector_get(s->x, 4),
+		gsl_vector_get(s->x, 5));
+
+}
+
+
+int main(void)
+{
+	const gsl_multifit_fdfsolver_type *T;
+	gsl_multifit_fdfsolver *s;
+	int status;
+	unsigned int i, iter, done;
+	size_t n;
+	gsl_matrix *covar;
+	MVal *vals;
+	struct data d;
+	gsl_multifit_function_fdf f;
+	double x_init[6];
+	gsl_vector_view x = gsl_vector_view_array(x_init, 6);
+	Cell cell;
+
+	/* Get the data to be fitted */
+	done = 0; n = 0;
+	vals = NULL;
+	while ( !done ) {
+
+		signed int h, k, l;
+		char buf[64];
+		int dspacing = 0;
+
+//		printf("Enter the indicies of the first set of planes, "
+//			" in the format 'h k l'.\n");
+//		if ( n != 0 ) printf("If you have no further planes, "
+//					"just hit enter.\n");
+		printf("        First plane indicies: ");
+		if ( fgets(buf, 63, stdin) != buf ) {
+			fprintf(stderr, "Error reading from input\n");
+		}
+		if ( sscanf(buf, "%i %i %i", &h, &k, &l) != 3 ) {
+			if ( buf[0] == '\n' ) {
+				done = 1;
+				continue;
+			} else {
+				printf("Invalid input, try again.\n");
+			}
+		}
+
+		/* Allocate space for the new reading,
+			now that we know it exists */
+		vals = realloc(vals, (n+1)*sizeof(MVal));
+		vals[n].h1 = h;
+		vals[n].k1 = k;
+		vals[n].l1 = l;
+
+//		printf("Enter the indicies of the second set of planes in the "
+//			"same format.\n");
+//		if ( n != 0 ) printf("If you are trying to give a d-spacing "
+//					"only, just hit enter here.\n");
+		printf("       Second plane indicies: ");
+		if ( fgets(buf, 63, stdin) != buf ) {
+			fprintf(stderr, "Error reading from input\n");
+		}
+		if ( sscanf(buf, "%i %i %i", &h, &k, &l) != 3 ) {
+			if ( buf[0] == '\n' ) {
+				dspacing = 1;
+			} else {
+				printf("Invalid input, try again.\n");
+			}
+		}
+
+		if ( dspacing ) {
+			vals[n].h2 = 0;
+			vals[n].k2 = 0;
+			vals[n].l2 = 0;
+			vals[n].meas = read_value("                   D-spacing: ");
+		} else {
+			vals[n].h2 = h;
+			vals[n].k2 = k;
+			vals[n].l2 = l;
+			vals[n].meas = DEG2RAD(read_value("        Angle between planes: "));
+		}
+
+		vals[n].esd = read_value("Estimated standard deviation: ");
+
+		printf("--------------------------------------\n");
+		n++;
+
+	}
+
+	x_init[0] = read_value("    Initial guess for a: ");
+	x_init[1] = read_value("    Initial guess for b: ");
+	x_init[2] = read_value("    Initial guess for c: ");
+	x_init[3] = read_value("Initial guess for alpha: ");
+	x_init[4] = read_value(" Initial guess for beta: ");
+	x_init[5] = read_value("Initial guess for gamma: ");
+
+	printf("-----------------------------------------\n\n");
+
+	printf("Inputted values:\n");
+	printf("--------------------------------------------------------------\n");
+	printf(" h1  k1  l1   h2  k2  l2    Spacing   Angle          ESD\n");
+	printf("--------------------------------------------------------------\n");
+	for ( i=0; i<n; i++ ) {
+		printf("%3i %3i %3i  ", vals[i].h1, vals[i].k1, vals[i].l1);
+		if ( is_dspacing(vals[i]) ) {
+			printf("  -   -   -   %8.5f    -     ", vals[i].meas);
+		} else {
+			printf("%3i %3i %3i        -    %8.5f", vals[i].h2,
+				vals[i].k2, vals[i].l2, RAD2DEG(vals[i].meas));
+		}
+		printf(" +/- %8.5f", vals[i].esd);
+		if ( !is_dspacing(vals[i]) ) {
+			printf(" deg\n");
+		} else {
+			printf(" nm\n");
+		}
+	}
+	printf("--------------------------------------------------------------\n");
+
+	d.vals = vals;
+	d.n = n;
+
+	f.f = &expb_f;
+	f.df = &expb_df;
+	f.fdf = &expb_fdf;
+	f.n = n;
+	f.p = 6;
+	f.params = &d;
+
+	T = gsl_multifit_fdfsolver_lmsder;
+	s = gsl_multifit_fdfsolver_alloc(T, n, 6);
+	gsl_multifit_fdfsolver_set(s, &f, &x.vector);
+
+	iter = 0;
+	printf("\n");
+	printf("iter    a        b        c      alpha    beta    gamma\n");
+	printf("---------------------------------------------------------\n");
+	print_state(iter, s);
+
+	do {
+
+		iter++;
+		status = gsl_multifit_fdfsolver_iterate(s);
+
+		printf("  status = %s\n", gsl_strerror(status));
+		print_state(iter, s);
+		if ( status ) break;
+
+		status = gsl_multifit_test_delta(s->dx, s->x, 1e-6, 1e-6);
+
+	} while ( (status == GSL_CONTINUE) && (iter < 500) );
+
+	covar = gsl_matrix_alloc(6, 6);
+	gsl_multifit_covar(s->J, 0.0, covar);
+
+	#define FIT(i) gsl_vector_get(s->x, i)
+	#define ERR(i) sqrt(gsl_matrix_get(covar,i,i))
+
+	printf("\n---------------------------------------------------------\n");
+
+	{
+		double chi = gsl_blas_dnrm2(s->f);
+		double dof = n - 6;
+		double c = GSL_MAX_DBL(1, chi / sqrt(dof));
+
+		printf("chisq/dof = %g / %g = %g\n",  pow(chi, 2.0), dof,
+							pow(chi, 2.0) / dof);
+
+		printf ("a      = %8.5f +/- %8.5f\n", FIT(0), c*ERR(0));
+		printf ("b      = %8.5f +/- %8.5f\n", FIT(1), c*ERR(1));
+		printf ("c      = %8.5f +/- %8.5f\n", FIT(2), c*ERR(2));
+		printf ("alpha  = %8.5f +/- %8.5f\n", FIT(3), c*ERR(3));
+		printf (" beta  = %8.5f +/- %8.5f\n", FIT(4), c*ERR(4));
+		printf ("gamma  = %8.5f +/- %8.5f\n", FIT(5), c*ERR(5));
+	}
+
+	printf ("status = %s\n", gsl_strerror (status));
+
+	gsl_multifit_fdfsolver_free(s);
+	gsl_matrix_free(covar);
+
+	cell.a = FIT(0);
+	cell.b = FIT(1);
+	cell.c = FIT(2);
+	cell.al = FIT(3);
+	cell.be = FIT(4);
+	cell.ga = FIT(5);
+
+	printf("\nValues given by fitted cell:\n");
+	printf("--------------------------------------------------------------\n");
+	printf(" h1  k1  l1   h2  k2  l2    Spacing   Angle          ESD\n");
+	printf("--------------------------------------------------------------\n");
+	for ( i=0; i<n; i++ ) {
+		printf("%3i %3i %3i  ", vals[i].h1, vals[i].k1, vals[i].l1);
+		if ( is_dspacing(vals[i]) ) {
+			printf("  -   -   -   %8.5f    -     ",
+				crystal_calc(vals[i], cell));
+		} else {
+			printf("%3i %3i %3i        -    %8.5f", vals[i].h2,
+				vals[i].k2, vals[i].l2,
+				RAD2DEG(crystal_calc(vals[i], cell)));
+		}
+		printf(" +/- %8.5f", vals[i].esd);
+		if ( !is_dspacing(vals[i]) ) {
+			printf(" deg\n");
+		} else {
+			printf(" nm\n");
+		}
+	}
+	printf("--------------------------------------------------------------\n");
+
+
+	return 0;
+}
diff --git a/triclinator.run b/triclinator.run
new file mode 100644
index 0000000..2a6d503
--- /dev/null
+++ b/triclinator.run
@@ -0,0 +1,32 @@
+1 0 0
+
+0.799
+0.04
+0 1 0
+
+0.548
+0.03
+0 0 1
+
+0.475
+0.02
+0 1 1
+
+0.376
+0.02
+1 0 0
+0 1 0
+90.5
+0.3
+1 0 0
+0 0 1
+84.2
+0.3
+
+0.821
+0.493
+0.585
+90.0
+94.7
+90.0
+
diff --git a/util.c b/util.c
new file mode 100644
index 0000000..11b1bcf
--- /dev/null
+++ b/util.c
@@ -0,0 +1,39 @@
+/*
+ * util.c
+ *
+ * Utility stuff
+ *
+ * (c) 2009 Thomas White <taw27@cam.ac.uk>
+ *
+ * Triclinator - solve nasty triclinic unit cells
+ *
+ */
+
+#include <stdio.h>
+
+#include "util.h"
+
+double read_value(const char *text)
+{
+	while ( 1 ) {
+
+		float d;
+		char buf[64];
+
+		printf("%s", text);
+		if ( fgets(buf, 63, stdin) != buf ) {
+			fprintf(stderr, "Error reading from input\n");
+		}
+		if ( sscanf(buf, "%f", &d) != 1 ) {
+			printf("Invalid input, try again.\n");
+		} else {
+			return d;
+		}
+	}
+}
+
+int is_dspacing(MVal val)
+{
+	if ( (val.h2 == 0) && (val.k2 == 0) && (val.l2 == 0) ) return 1;
+	return 0;
+}
diff --git a/util.h b/util.h
new file mode 100644
index 0000000..20e6a37
--- /dev/null
+++ b/util.h
@@ -0,0 +1,32 @@
+/*
+ * util.h
+ *
+ * Utility stuff
+ *
+ * (c) 2009 Thomas White <taw27@cam.ac.uk>
+ *
+ * Triclinator - solve nasty triclinic unit cells
+ *
+ */
+
+#ifndef UTIL_H
+#define UTIL_H
+
+typedef struct
+{
+	signed int h1;
+	signed int k1;
+	signed int l1;
+
+	signed int h2;
+	signed int k2;
+	signed int l2;	/* Or zero */
+
+	double meas;
+	double esd;
+} MVal;
+
+extern double read_value(const char *text);
+extern int is_dspacing(MVal val);
+
+#endif	/* UTIL_H */