From 244303f588f9c4797836e062d0576d85a027ab2a Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 4 Mar 2009 17:05:42 +0000 Subject: Initial import --- main.c | 334 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 334 insertions(+) create mode 100644 main.c (limited to 'main.c') 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