From 244303f588f9c4797836e062d0576d85a027ab2a Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 4 Mar 2009 17:05:42 +0000 Subject: Initial import --- .gitignore | 2 + Makefile | 16 +++ crystal.c | 110 +++++++++++++++++++ crystal.h | 33 ++++++ main.c | 334 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ triclinator.run | 32 ++++++ util.c | 39 +++++++ util.h | 32 ++++++ 8 files changed, 598 insertions(+) create mode 100644 .gitignore create mode 100644 Makefile create mode 100644 crystal.c create mode 100644 crystal.h create mode 100644 main.c create mode 100644 triclinator.run create mode 100644 util.c create mode 100644 util.h 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 + * + * Triclinator - solve nasty triclinic unit cells + * + */ + + +#include + +#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 + * + * 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 + * + * Triclinator - solve nasty triclinic unit cells + * + */ + + +#include +#include +#include +#include +#include +#include +#include +#include + +#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; in; + 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; ix, 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; idx, 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 + * + * Triclinator - solve nasty triclinic unit cells + * + */ + +#include + +#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 + * + * 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 */ -- cgit v1.2.3