/*
 * basis.c
 *
 * Handle basis structures
 *
 * (c) 2007 Thomas White <taw27@cam.ac.uk>
 *
 *  dtr - Diffraction Tomography Reconstruction
 *
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#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"
#include "utils.h"

double basis_efom(ReflectionList *reflectionlist, Basis *basis) {

	int n_indexed, n_counted;
	Reflection *cur;
	
	cur = reflectionlist->reflections;
	n_indexed = 0;
	n_counted = 0;
	while ( cur ) {

		if ( cur->type == REFLECTION_NORMAL ) {
			
			/* Can this basis "approximately" account for this reflection? */
			double det;
			double a11, a12, a13, a21, a22, a23, a31, a32, a33;
			double h, k, l;
			
			/* Set up the coordinate transform from hkl to xyz */
			a11 = basis->a.x;  a12 = basis->a.y;  a13 = basis->a.z;
			a21 = basis->b.x;  a22 = basis->b.y;  a23 = basis->b.z;
			a31 = basis->c.x;  a32 = basis->c.y;  a33 = basis->c.z;
			
			/* Invert the matrix to get hkl from xyz */
			det = a11*(a22*a33 - a23*a32) - a12*(a21*a33 - a23*a31) + a13*(a21*a32 - a22*a31);
			h = ((a22*a33-a23*a32)*cur->x + (a23*a31-a21*a33)*cur->y + (a21*a32-a22*a31)*cur->z) / det;
			k = ((a13*a32-a12*a33)*cur->x + (a11*a33-a13*a31)*cur->y + (a12*a31-a11*a32)*cur->z) / det;
			l = ((a12*a23-a13*a22)*cur->x + (a13*a21-a11*a23)*cur->y + (a11*a22-a12*a21)*cur->z) / det;
			
			/* Calculate the deviations in terms of |a|, |b| and |c| */
			h = fabs(h);  k = fabs(k);  l = fabs(l);
			h -= floor(h);  k -= floor(k);  l -= floor(l);
			if ( h == 1.0 ) h = 0.0;
			if ( k == 1.0 ) k = 0.0;
			if ( l == 1.0 ) l = 0.0;
			
			/* Define "approximately" here.  Circle in basis space becomes an ellipsoid in reciprocal space */
			if ( h*h + k*k + l*l <= 0.1*0.1*0.1 ) n_indexed++;
			
			n_counted++;
		
		}
		
		cur = cur->next;

	}
	
	return (double)n_indexed / n_counted;

}

Basis basis_add(Basis u, Basis v) {

	Basis ans;
	
	ans.a.x = u.a.x + v.a.x;	ans.a.y = u.a.y + v.a.y;	ans.a.z = u.a.z + v.a.z;
	ans.b.x = u.b.x + v.b.x;	ans.b.y = u.b.y + v.b.y;	ans.b.z = u.b.z + v.b.z;
	ans.c.x = u.c.x + v.c.x;	ans.c.y = u.c.y + v.c.y;	ans.c.z = u.c.z + v.c.z;
	
	return ans;

}

static void basis_print(Basis *cell) {

	printf("%12.8f %12.8f %12.8f\n", cell->a.x/1e9, cell->a.y/1e9, cell->a.z/1e9);
	printf("%12.8f %12.8f %12.8f\n", cell->b.x/1e9, cell->b.y/1e9, cell->b.z/1e9);
	printf("%12.8f %12.8f %12.8f\n", cell->c.x/1e9, cell->c.y/1e9, cell->c.z/1e9);	

}

static void cell_print(UnitCell *cell) {

	printf("%12.8f %12.8f %12.8f nm\n", cell->a*1e9, cell->b*1e9, cell->c*1e9);
	printf("%12.8f %12.8f %12.8f deg\n", rad2deg(cell->alpha), rad2deg(cell->beta), rad2deg(cell->gamma));

}

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;
	
//	basis->a.x = 0.5;  basis->a.y = 0.0;  basis->a.z = 0.0;
//	basis->b.x = 0.0;  basis->b.y = 0.2;  basis->b.z = 0.0;
//	basis->c.x = 0.0;  basis->c.y = 0.0;  basis->c.z = 1.0;
	printf("Reciprocal-space cell (nm^-1):\n");
	basis_print(basis);
	
	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);
	
	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);
	
	gsl_matrix_transpose(inv);
	
	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);
	
	printf("Real-space cell (nm):\n");
	printf("%12.8f %12.8f %12.8f\n", ax*1e9, ay*1e9, az*1e9);
	printf("%12.8f %12.8f %12.8f\n", bx*1e9, by*1e9, bz*1e9);
	printf("%12.8f %12.8f %12.8f\n", cx*1e9, cy*1e9, cz*1e9);
	
	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);
	
	printf("Cell parameters:\n");	
	cell_print(&cell);
	
	return cell;
	
}