summaryrefslogtreecommitdiff
path: root/main.c
diff options
context:
space:
mode:
Diffstat (limited to 'main.c')
-rw-r--r--main.c334
1 files changed, 334 insertions, 0 deletions
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;
+}
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-03 05:55:59 +0000