/////////////////////////////////////////////////////////////////////////////////////////////////////
//
//	FileName:	FracOrth.cpp
//	Author	:	Michael Y. Polyakov
//	email	:	myp@andrew.cmu.edu	or  mikepolyakov@hotmail.com
//	Website	:	www.angelfire.com/linux/myp
//	Date	:	7/30/2002	
//
//	class FracOrth: converts between fractional and orthonormal coordinates. This is very small and basic
//		class. For a full Crystallographic C++ library written by Vincent Favre-Nicolin see this link:
//				http://objcryst.sourceforge.net/
//////////////////////////////////////////////////////////////////////////////////////////////////////

#include "FracOrth.h"
#include <math.h>

//Class FracOrth

//constructors

//a,b,c - orthonormal unit vectors
FracOrth::FracOrth(float a, float b, float c)
{  
	Init(a, b, c, M_PI/2.0, M_PI/2.0, M_PI/2.0);
}

//all unit cell dimensions are specified (angles are in degrees)
FracOrth::FracOrth(float a, float b, float c, float alpha, float beta, float gamma)
{
	Init(a, b, c, alpha*M_PI/180.0, beta*M_PI/180.0, gamma*M_PI/180.0);
}

void FracOrth::SetDim(float a, float b, float c)
{
	Init(a, b, c, ang[0], ang[1], ang[2]);
}

void FracOrth::SetAngles(float alpha, float beta, float gamma)
{
	Init(dim[0], dim[1], dim[2], alpha*M_PI/180.0, beta*M_PI/180.0, gamma*M_PI/180.0);
}

//Converts from fractional to orthonormal coordinates
void FracOrth::FracToOrth(float &x, float &y, float &z)
{
	float pp[3] = {x,y,z};
	x = M[0][0]*pp[0] + M[0][1]*pp[1] + M[0][2]*pp[2];
	y = M[1][1]*pp[1] + M[1][2]*pp[2];
	z = M[2][2]*pp[2];
}

//Converts from orthonormal to fractional coordinates
void FracOrth::OrthToFrac(float &x, float &y, float &z)
{
	float p[3] = {x,y,z};
	x = Min[0][0]*p[0] + Min[0][1]*p[1] + Min[0][2]*p[2];
	y = Min[1][1]*p[1] + Min[1][2]*p[2];
	z = Min[2][2]*p[2];
}

void FracOrth::GetDim(float &a, float &b, float &c)
{
	a = dim[0]; b = dim[1]; c = dim[2];
}

void FracOrth::GetAngles(float &alpha, float &beta, float &gamma)
{
	alpha = ang[0]; beta = ang[1]; gamma = ang[2];
}

float FracOrth::GetVolume()
{
	return V;
}

float* FracOrth::GetFracToOrthM()
{
	return (float*)M;
}

float* FracOrth::GetOrthToFracM()
{
	return (float*)Min;
}

//private function - initializes data
void FracOrth::Init(float a, float b, float c, float alpha, float beta, float gamma)
{
	dim[0] = a; dim[1] = b; dim[2] = c;
	ang[0] = alpha; ang[1] = beta; ang[2] = gamma;
	//volume of a unit cell
	V = dim[0]*dim[1]*dim[2]*sqrt(1-cos(ang[0])*cos(ang[0]) - cos(ang[1])*cos(ang[1])
		- cos(ang[2])*cos(ang[2]) + 2*dim[1]*dim[2]*cos(ang[0])*cos(ang[1])*cos(ang[2]));
	//Transformation matrix from p' into p (fractinal to orthonormal)
	M[0][0] = dim[0];
	M[0][1] = dim[1]*cos(ang[2]);
	M[0][2] = dim[2]*cos(ang[1]);
	M[1][0] = 0.0;
	M[1][1] = dim[1]*sin(ang[2]);
	M[1][2] = dim[2]*(cos(ang[0])-cos(ang[1])*cos(ang[2]))/sin(ang[2]);
	M[2][0] = 0.0;
	M[2][1] = 0.0;
	M[2][2] = V / (dim[0]*dim[1]*sin(ang[2]));
	//Tranformation matrix from p into p' (orthonormal to fractional)
	Min[0][0] = 1.0 / dim[0];
	Min[0][1] = -cos(ang[2])/(dim[0]*sin(ang[2]));
	Min[0][2] = 1.0 / V * ( (dim[1]*dim[2]*cos(ang[2])*(cos(ang[0])-cos(ang[1])*cos(ang[2]))) / sin(ang[2]) - 
							dim[1]*dim[2]*cos(ang[1])*sin(ang[2]) );
	Min[1][0] = 0.0;
	Min[1][1] = 1.0 / (dim[1]*sin(ang[2]));
	Min[1][2] = -dim[0]*dim[2]*(cos(ang[0])-cos(ang[1])*cos(ang[2])) / (V * sin(ang[2]));
	Min[2][0] = 0.0;
	Min[2][1] = 0.0;
	Min[2][2] = dim[0]*dim[1]*sin(ang[2]) / V;
}