Save lattice parameters with reflection data

git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@228 bf6ca9ba-c028-0410-8290-897cf20841d1

3 files changed, 68 insertions, 2 deletions

diff --git a/src/basis.c b/src/basis.c
index f3b9b7a..46335fa 100644
--- a/src/basis.c
+++ b/src/basis.c
@@ -14,6 +14,9 @@
 #endif
 
 #include <math.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
 
 #include "reflections.h"
 #include "basis.h"
@@ -80,3 +83,43 @@ Basis basis_add(Basis u, Basis v) {
 
 }
 
+UnitCell basis_get_cell(Basis *basis) {
+
+	UnitCell cell;
+	gsl_matrix *m;
+	gsl_matrix *inv;
+	gsl_permutation *perm;
+	double ax, ay, az, bx, by, bz, cx, cy, cz;
+	int s;
+	
+	m = gsl_matrix_alloc(3, 3);
+	gsl_matrix_set(m, 0, 0, basis->a.x);  gsl_matrix_set(m, 0, 1, basis->b.x);  gsl_matrix_set(m, 0, 2, basis->c.x);
+	gsl_matrix_set(m, 1, 0, basis->a.y);  gsl_matrix_set(m, 1, 1, basis->b.y);  gsl_matrix_set(m, 1, 2, basis->c.y);
+	gsl_matrix_set(m, 2, 0, basis->a.z);  gsl_matrix_set(m, 2, 1, basis->b.z);  gsl_matrix_set(m, 2, 2, basis->c.z);
+	
+	gsl_matrix_transpose(m);
+	
+	perm = gsl_permutation_alloc(m->size1);
+	inv = gsl_matrix_alloc(m->size1, m->size2);
+	gsl_linalg_LU_decomp(m, perm, &s);
+	gsl_linalg_LU_invert(m, perm, inv);
+	gsl_permutation_free(perm);
+	gsl_matrix_free(m);
+	
+	ax = gsl_matrix_get(inv, 0, 0);  bx = gsl_matrix_get(inv, 0, 1);  cx = gsl_matrix_get(inv, 0, 2);
+	ay = gsl_matrix_get(inv, 1, 0);  by = gsl_matrix_get(inv, 1, 1);  cy = gsl_matrix_get(inv, 1, 2);
+	az = gsl_matrix_get(inv, 2, 0);  bz = gsl_matrix_get(inv, 2, 1);  cz = gsl_matrix_get(inv, 2, 2);
+	
+	cell.a = sqrt(ax*ax + ay*ay + az*az);
+	cell.b = sqrt(bx*bx + by*by + bz*bz);
+	cell.c = sqrt(cx*cx + cy*cy + cz*cz);
+	cell.alpha = acos((bx*cx + by*cy + bz*cz)/(cell.b * cell.c));
+	cell.beta = acos((ax*cx + ay*cy + az*cz)/(cell.a * cell.c));
+	cell.gamma = acos((bx*ax + by*ay + bz*az)/(cell.b * cell.a));
+	
+	gsl_matrix_free(inv);
+	
+	return cell;
+	
+}
+
diff --git a/src/basis.h b/src/basis.h
index f8248c6..25eb0ae 100644
--- a/src/basis.h
+++ b/src/basis.h
@@ -32,8 +32,18 @@ typedef struct basis_struct {
 	Vector c;
 } Basis;
 
+typedef struct cell_struct {
+	double a;
+	double b;
+	double c;
+	double alpha;
+	double beta;
+	double gamma;
+} UnitCell;
+
 extern double basis_efom(struct reflectionlist_struct *reflectionlist, Basis *basis);
 extern Basis basis_add(Basis u, Basis v);
+extern UnitCell basis_get_cell(Basis *cell);
 
 #endif	/* BASIS_H */
 
diff --git a/src/intensities.c b/src/intensities.c
index 8c02592..d5f345d 100644
--- a/src/intensities.c
+++ b/src/intensities.c
@@ -18,6 +18,8 @@
 #include "image.h"
 #include "reproject.h"
 #include "displaywindow.h"
+#include "utils.h"
+#include "basis.h"
 
 /* Extract integrated reflection intensities by estimating the spike function
  * based on the observed intensity and the calculated excitation error from
@@ -108,12 +110,23 @@ void intensities_extract(ControlContext *ctx) {
 
 }
 
-static int intensities_do_save(ReflectionList *integrated, const char *filename) {
+static int intensities_do_save(ReflectionList *integrated, Basis *cell, const char *filename) {
 
 	FILE *fh;
 	Reflection *reflection;
+	UnitCell rcell;
 	
 	fh = fopen(filename, "w");
+	
+	rcell = basis_get_cell(cell);
+	fprintf(fh, "a %12.8f\n", rcell.a*1e9);
+	fprintf(fh, "b %12.8f\n", rcell.b*1e9);
+	fprintf(fh, "c %12.8f\n", rcell.c*1e9);
+	fprintf(fh, "alpha %12.8f\n", rad2deg(rcell.alpha));
+	fprintf(fh, "beta %12.8f\n", rad2deg(rcell.beta));
+	fprintf(fh, "gamma %12.8f\n", rad2deg(rcell.gamma));
+	
+	
 	reflection = integrated->reflections;
 	while ( reflection ) {
 		fprintf(fh, "%3i %3i %3i %12.8f\n", reflection->h, reflection->k, reflection->l, reflection->intensity);
@@ -130,7 +143,7 @@ static gint intensities_save_response(GtkWidget *widget, gint response, ControlC
 	if ( response == GTK_RESPONSE_ACCEPT ) {
 		char *filename;
 		filename = gtk_file_chooser_get_filename(GTK_FILE_CHOOSER(widget));
-		if ( intensities_do_save(ctx->integrated, filename) ) {
+		if ( intensities_do_save(ctx->integrated, ctx->cell, filename) ) {
 			displaywindow_error("Failed to save cache file.", ctx->dw);
 		}
 		g_free(filename);