/*
 * crystal.c
 *
 * Crystallographic stuff
 *
 * (c) 2009 Thomas White <taw27@cam.ac.uk>
 *
 * Triclinator - solve nasty triclinic unit cells
 *
 */


#include <math.h>

#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) );
}


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);
}


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);

	/* Bodge to put the axes in the right place */
	return M_PI - 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);
	}
}