/*
(c) Oleg V. Dolomanov 2002-3
This source code is provided without warranty of any kind.
The author is not responsible for losses of information or any other 
damage caused by this code. You use it on your own risk. You can distribute, 
modify and incorporate this code into any non-commercial packages without 
author's permission. Acknowledgements are not necessary, but would be 
appreciated.
*/

bool _fastcall TCell::InRange(TCell &C, double Dev)
{
	double 	min = 1-Dev,
			max = 1+Dev;

	if( (a >= C.a*min) && (a <= C.a*max) &&
		(b >= C.b*min) && (b <= C.b*max) &&
		(c >= C.c*min) && (c <= C.c*max) &&
		(aa >= C.aa*min) && (aa <= C.aa*max) &&
		(ab >= C.ab*min) && (ab <= C.ab*max) &&
		(ac >= C.ac*min) && (ac <= C.ac*max)	)
		return true;
	if( (b >= C.a*min) && (b <= C.a*max) &&
		(c >= C.b*min) && (c <= C.b*max) &&
		(a >= C.c*min) && (a <= C.c*max) &&
		(ab >= C.aa*min) && (aa <= C.aa*max) &&
		(ac >= C.ab*min) && (ac <= C.ab*max) &&
		(aa >= C.ac*min) && (aa <= C.ac*max)	)
		return true;
	if( (c >= C.a*min) && (c <= C.a*max) &&
		(a >= C.b*min) && (a <= C.b*max) &&
		(b >= C.c*min) && (b <= C.c*max) &&
		(ac >= C.aa*min) && (ac <= C.aa*max) &&
		(aa >= C.ab*min) && (aa <= C.ab*max) &&
		(ab >= C.ac*min) && (ab <= C.ac*max)	)
		return true;

	// compare niggli form

	if( fabs(((float)na  - (float)C.na)/(float)na) > (float)Dev )
		return false;
	if( fabs(((float)nb  - (float)C.nb)/(float)nb) > (float)Dev )
		return false;
	if( fabs(((float)nc  - (float)C.nc)/(float)nc) > (float)Dev )
		return false;
	if( fabs(((float)naa  - (float)C.naa)/(float)naa) > (float)Dev )
		return false;
	if( fabs(((float)nab  - (float)C.nab)/(float)nab) > (float)Dev )
		return false;
	if( fabs(((float)nac  - (float)C.nac)/(float)nac) > (float)Dev )
		return false;

	if( (na  >= C.na*min)  &&  (na  <= C.na*max) &&
		(nb  >= C.nb*min)  &&  (nb  <= C.nb*max) &&
		(nc  >= C.nc*min)  &&  (nc  <= C.nc*max) &&
		(naa >= C.naa*min) && (naa <= C.naa*max) &&
		(nab >= C.nab*min) && (nab <= C.nab*max) &&
		(nac >= C.nac*min) && (nac <= C.nac*max)	)
		return true;
	return false;
}

bool TCell::operator == (TCell &C)
{
	if( (a == C.a) && (b == C.b) && (c == C.c) &&
		(aa == C.aa) && (ab == C.ab) && (ac == C.ac) )
		return true;
	return false;
}
void _fastcall TCell::operator = (TCell &C)
{
	a = C.a;
	b = C.b;
	c = C.c;
	aa = C.aa;
	ab = C.ab;
	ac = C.ac;

	na = C.na;
	nb = C.nb;
	nc = C.nc;
	naa = C.naa;
	nab = C.nab;
	nac = C.nac;
	Lattice = C.Lattice;
}
//---------------------------------------------------------------------------
bool _fastcall ReduceCell(TCell &C)
{
	VPoint X,Y,Z;
	int R[6];
	double Rd[6];
	double val, TX, TY, TZ, AX, AY, AZ, cosa, sina;

	bool Used;
	int Cycles;
	cosa = cos(C.ac/180.0*M_PI);
	sina = sin(C.ac/180.0*M_PI);

	TX = cos(C.ab/180.0*M_PI);
	TY = (cos(C.aa*M_PI/180.0) - cosa*cos(C.ab*M_PI/180.0))/sina;
	TZ = sqrt( 1.0 - TX*TX-TY*TY);

	switch(C.Lattice)
	{
		case 1:	// P
			X.x0 = 1;	X.y0 = 0;		X.z0 = 0;
			Y.x0 = 0;	Y.y0 = 1;		Y.z0 = 0;
			Z.x0 = 0;	Z.y0 = 0;		Z.z0 = 1;
			break;
		case 2:	// I
			X.x0 = 1;	X.y0 = 0;		X.z0 = 0;
			Y.x0 = 0;	Y.y0 = 1;		Y.z0 = 0;
			Z.x0 = 1./2;	Z.y0 = 1./2;	Z.z0 = 1./2;
			break;
		case 3:	// R
			X.x0 = 2./3;	X.y0 = 1./3;	X.z0 = 1./3;
			Y.x0 = -1./3;	Y.y0 = 1./3;	Y.z0 = 1./3;
			Z.x0 = -1./3;	Z.y0 = -2./3;	Z.z0 = 1./3;
			break;
		case 4:	// F
			X.x0 = 1./2;	X.y0 = 0;		X.z0 = 1./2;
			Y.x0 = 1./2;	Y.y0 = 1./2;	Y.z0 = 0;
			Z.x0 = 0;		Z.y0 = 1./2;	Z.z0 = 1./2;
			break;
		case 5:	// A
			X.x0 = 1;	X.y0 = 0;		X.z0 = 0;
			Y.x0 = 0;	Y.y0 = 1;		Y.z0 = 0;
			Z.x0 = 0;	Z.y0 = 1./2;	Z.z0 = 1./2;
			break;
		case 6:	// B
			X.x0 = 1./2;	X.y0 = 0;		X.z0 = 1./2;
			Y.x0 = 0;	Y.y0 = 1;		Y.z0 = 0;
			Z.x0 = 0;	Z.y0 = 0;		Z.z0 = 1;
			break;
		case 7:	// C
			X.x0 = 1;	X.y0 = 0;		X.z0 = 0;
			Y.x0 = 1./2;	Y.y0 = 1./2;	Y.z0 = 0;
			Z.x0 = 0;	Z.y0 = 0;		Z.z0 = 1;
			break;
		default:
			Application->MessageBox("Cannot reduce the cell: unknown centering...", "Error", MB_OK|MB_ICONERROR);
			return false;
	}

	AX = X.x0*C.a;	AY = X.y0*C.b;	AZ = X.z0*C.c;
	X.x0 = AX + AY*cosa+ AZ*TX;
	X.y0 = AY*sina + AZ*TY;
	X.z0 = AZ*TZ;

	AX = Y.x0*C.a;	AY = Y.y0*C.b;	AZ = Y.z0*C.c;
	Y.x0 = AX + AY*cosa + AZ*TX;
	Y.y0 = AY*sina + AZ*TY;
	Y.z0 = AZ*TZ;

	AX = Z.x0*C.a;	AY = Z.y0*C.b;	AZ = Z.z0*C.c;
	Z.x0 = AX + AY*cosa+ AZ*TX;
	Z.y0 = AY*sina + AZ*TY;
	Z.z0 = AZ*TZ;

	C.na = X.Length;
	C.nb = Y.Length;
	C.nc = Z.Length;

	C.naa = acos(Y.CAngle(Z))*180/M_PI;
	C.nab = acos(X.CAngle(Z))*180/M_PI;
	C.nac = acos(Y.CAngle(X))*180/M_PI;

	Rd[0] = qrt(C.na);
	Rd[1] = qrt(C.nb);
	Rd[2] = qrt(C.nc);

	Rd[3] = 2*C.nb*C.nc*cos(C.naa*M_PI/180);
	Rd[4] = 2*C.na*C.nc*cos(C.nab*M_PI/180);
	Rd[5] = 2*C.na*C.nb*cos(C.nac*M_PI/180);

	for( int i=0; i < 6; i++ )
		R[i] = round(Rd[i]);

	Used = true;
	Cycles = 0;
	while ( Used )
	{
		Used = false;
		Cycles ++;
		if( Cycles > 100 )
		{
			Application->MessageBox("Cannot evaluate the niggli cell...", "Error", MB_OK|MB_ICONERROR);
			return false;
		}
		if( (R[0] > R[1]) || ((R[0] == R[1])&&(fabs(R[3])> fabs(R[4]))) ) //1
		{
			val = R[0];		R[0] = R[1];			R[1] = val;
			val = R[3];		R[3] = R[4];			R[4] = val;

			val = Rd[0];		Rd[0] = Rd[1];			Rd[1] = val;
			val = Rd[3];		Rd[3] = Rd[4];			Rd[4] = val;
			Used = true;
		}
		if( (R[1] > R[2]) || ( (R[1] == R[2])&&(fabs(R[4])> fabs(R[5]))) )//2
		{
			val = R[1];		R[1] = R[2];			R[2] = val;
			val = R[4];		R[4] = R[5];			R[5] = val;

			val = Rd[1];		Rd[1] = Rd[2];			Rd[2] = val;
			val = Rd[4];		Rd[4] = Rd[5];			Rd[5] = val;
			Used = true;
			continue;
		}
		if( R[3]*R[4]*R[5] > 0 ) //3
		{
			R[3] = fabs(R[3]);		R[4] = fabs(R[4]);			R[5] = fabs(R[5]);

			Rd[3] = fabs(Rd[3]);		Rd[4] = fabs(Rd[4]); 		Rd[5] = fabs(Rd[5]);
		}
		else
		{
			R[3] = -fabs(R[3]);		R[4] = -fabs(R[4]);			R[5] = -fabs(R[5]);

			Rd[3] = -fabs(Rd[3]);		Rd[4] = -fabs(Rd[4]); 		Rd[5] = -fabs(Rd[5]);
		}
		if( (fabs(R[3]) > R[1]) || ( (R[3] == R[1])&& (2*R[4]<R[5]) ) || //5
			( (R[3]==-R[1])&&(R[5]<0)) )
		{
			R[2] = R[1]+R[2] - R[3]*sign(R[3]);
			R[4] = R[4]-R[5]*sign(R[3]);
			R[3] = R[3]-2*R[1]*sign(R[3]);

			Rd[2] = Rd[1]+Rd[2] - Rd[3]*sign(Rd[3]);
			Rd[4] = Rd[4]-Rd[5]*sign(Rd[3]);
			Rd[3] = Rd[3]-2*Rd[1]*sign(Rd[3]);
			Used = true;
			continue;
		}
		if( (fabs(R[4]) > R[0]) || ( (R[4] == R[0])&& (2*R[3]<R[5]) ) ||  //6
			( (R[4]==-R[0])&&(R[5]<0)) )
		{
			R[2] = R[0]+R[2] - R[4]*sign(R[4]);
			R[3] = R[3]-R[5]*sign(R[4]);
			R[4] = R[4]-2*R[0]*sign(R[4]);

			Rd[2] = Rd[0]+Rd[2] - Rd[4]*sign(Rd[4]);
			Rd[3] = Rd[3]-Rd[5]*sign(Rd[4]);
			Rd[4] = Rd[4]-2*Rd[0]*sign(Rd[4]);
			Used = true;
			continue;
		}
		if( (fabs(R[5]) > R[0]) || ( (R[5] == R[0])&& (2*R[3]<R[4]) ) ||  //7
			( (R[5]==-R[0])&&(R[4]<0)) )
		{
			R[1] = R[0]+R[1] - R[5]*sign(R[5]);
			R[3] = R[3]-R[4]*sign(R[5]);
			R[5] = R[5]-2*R[0]*sign(R[5]);

			Rd[1] = Rd[0]+Rd[1] - Rd[5]*sign(Rd[5]);
			Rd[3] = Rd[3]-Rd[4]*sign(Rd[5]);
			Rd[5] = Rd[5]-2*Rd[0]*sign(Rd[5]);
			Used = true;
			continue;
		}
		if( ((R[0]+R[1] + R[3]+R[4]+R[5])< 0 ) ||
			( ((R[0]+R[1] + R[3]+R[4]+R[5])==0) && (2*(R[0]+R[4])+R[5])>0 )) //8
		{
			R[2] = R[0]+R[1]+R[2]+R[3]+R[4]+R[5];
			R[3] = 2*R[1] + R[3]+R[5];
			R[4] = 2*R[0] + R[4]+R[5];

			Rd[2] = Rd[0]+Rd[1]+Rd[2]+Rd[3]+Rd[4]+Rd[5];
			Rd[3] = 2*Rd[1] + Rd[3]+Rd[5];
			Rd[4] = 2*Rd[0] + Rd[4]+Rd[5];
			Used = true;
			continue;
		}
	}

	Rd[0] = sqrt(Rd[0]);
	Rd[1] = sqrt(Rd[1]);
	Rd[2] = sqrt(Rd[2]);
	C.na = Rd[0];
	C.nb = Rd[1];
	C.nc = Rd[2];
	C.naa = acos(Rd[3]/(Rd[1]*Rd[2]*2))*180/M_PI;
	C.nab = acos(Rd[4]/(Rd[0]*Rd[2]*2))*180/M_PI;
	C.nac = acos(Rd[5]/(Rd[0]*Rd[1]*2))*180/M_PI;
	return true;
}