/*----------------------------------------------------------------------------
Copyright (C) : 2002, G.J.L. Kemp, Chalmers University of Technology
----------------------------------------------------------------------------*/

/*----------------------------------------------------------------------------
File	:	mini_geom.c
Author	:	Graham Kemp

This is a cut-down version of a larger set of geometry routines written
while I was working at the University of Aberdeen in 1988, and extended
slightly in later years.  A main() function has been added with some
example calls to the geometry routines.
----------------------------------------------------------------------------*/


/*----------------------------------------------------------------------------
Include files
----------------------------------------------------------------------------*/
#include <math.h>


/*----------------------------------------------------------------------------
Definitions
----------------------------------------------------------------------------*/
#define FALSE	0
#define TRUE	1
#define PI	M_PI

#define degrees_to_radians(X)	(X * PI / 180)
#define radians_to_degrees(X)	(X * 180 / PI)

#define dot_product(V1, V2)	(V1.x * V2.x + V1.y * V2.y + V1.z * V2.z)


/*----------------------------------------------------------------------------
Type definitions
----------------------------------------------------------------------------*/
typedef double	Distance;
typedef double	Angle;


typedef struct {
	double	x, y, z;
} Point, Vector;


/*----------------------------------------------------------------------------
D_PTPT()

Return the distance between two points.
----------------------------------------------------------------------------*/
Distance d_ptpt(pt1, pt2)
	Point	pt1, pt2;
{
	double	x, y, z;

	x = pt1.x - pt2.x;
	y = pt1.y - pt2.y;
	z = pt1.z - pt2.z;
	return (sqrt(x*x + y*y + z*z));
}


/*----------------------------------------------------------------------------
MAG_VECTOR()

Return the magnitude of vector v.
----------------------------------------------------------------------------*/
double mag_vector(v)
	Vector	v;
{
	return(sqrt(v.x * v.x + v.y * v.y + v.z * v.z));
}


/*----------------------------------------------------------------------------
V_PTPT()

Return the vector from pt1 to pt2.
----------------------------------------------------------------------------*/
Vector v_ptpt(pt1, pt2)
	Point	pt1, pt2;
{
	Vector	v;

	v.x = pt2.x - pt1.x;
	v.y = pt2.y - pt1.y;
	v.z = pt2.z - pt1.z;
	return(v);
}


/*----------------------------------------------------------------------------
CROSS_PRODUCT()
----------------------------------------------------------------------------*/
Vector cross_product(v1, v2)
	Vector	v1, v2;
{
	Vector	v3;

	v3.x = v1.y * v2.z - v1.z * v2.y;
	v3.y = v1.z * v2.x - v1.x * v2.z;
	v3.z = v1.x * v2.y - v1.y * v2.x;
	return(v3);
}


/*----------------------------------------------------------------------------
A_VV()

Return the angle between the given vectors.
----------------------------------------------------------------------------*/
Angle a_vv(v1, v2)
	Vector	v1, v2;
{
	return(acos(dot_product(v1, v2) / (mag_vector(v1) * mag_vector(v2))));
}


/*----------------------------------------------------------------------------
A_PTPTPT()
----------------------------------------------------------------------------*/
Angle a_ptptpt(a, b, c)
	Point	a, b, c;
{
	return(a_vv(v_ptpt(b, a), v_ptpt(b, c)));
}


/*----------------------------------------------------------------------------
A_PTPTPTPT()

Return the torsion angle defined by the four given points.

The idea is to find the angle between the vector perpendicular to the plane
containing the first three points and the vector perpendicular to the plane
containing the last three points.
----------------------------------------------------------------------------*/
Angle a_ptptptpt(a, b, c, d)
	Point	a, b, c, d;
{

	Vector	v, vp1, vp2;
	Angle	angle;

	/* Let vp1 and vp2 be vectors perpendicular to the two planes. */
	v = v_ptpt(b, c);
	vp1 = cross_product(v_ptpt(a, b), v);
	vp2 = cross_product(v, v_ptpt(c, d));

	angle = a_vv(vp1, vp2);

	/* Check the sign of the angle. */
	if (dot_product(cross_product(vp1, vp2), v) < 0)
		angle = -angle;

	return(angle);
}


/*----------------------------------------------------------------------------
MAIN()
----------------------------------------------------------------------------*/
main()
{
	Point	thr1n;
	Point	thr1ca;
	Point	thr1c;
	Point	thr1o;
	Point	thr2n;
	Point	thr2ca;
	Point	thr2c;
	Point	thr2o;

	/*
	 * Coordinates of the main chain atoms of the first two residues of
	 * Protein Data Bank Entry 1CRN.
	 */
	thr1n.x  = 17.047;  thr1n.y  = 14.099;  thr1n.z  = 3.625;
	thr1ca.x = 16.967;  thr1ca.y = 12.784;  thr1ca.z = 4.338;
	thr1c.x  = 15.685;  thr1c.y  = 12.755;  thr1c.z  = 5.133;
	thr1o.x  = 15.268;  thr1o.y  = 13.825;  thr1o.z  = 5.594;
	thr2n.x  = 15.115;  thr2n.y  = 11.555;  thr2n.z  = 5.265;
	thr2ca.x = 13.856;  thr2ca.y = 11.469;  thr2ca.z = 6.066;
	thr2c.x  = 14.164;  thr2c.y  = 10.785;  thr2c.z  = 7.379;
	thr2o.x  = 14.993;  thr2o.y  =  9.862;  thr2o.z  = 7.443;

	printf("Distance N-CA is %f\n", d_ptpt(thr1n,thr1ca));
	printf("Distance CA-C is %f\n", d_ptpt(thr1ca,thr1c));
	printf("Distance C=O  is %f\n", d_ptpt(thr1c,thr1o));

	printf("Angle N-CA-C is %f\n",
		radians_to_degrees(a_ptptpt(thr1n,thr1ca,thr1c)));
	printf("Angle CA-C=O is %f\n",
		radians_to_degrees(a_ptptpt(thr1ca,thr1c,thr1o)));

	printf("Angle psi   of residue 1 is %f\n",
		radians_to_degrees(a_ptptptpt(thr1n,thr1ca,thr1c,thr2n)));
	printf("Angle omega of residue 1 is %f\n",
		radians_to_degrees(a_ptptptpt(thr1ca,thr1c,thr2n,thr2ca)));
	printf("Angle phi   of residue 2 is %f\n",
		radians_to_degrees(a_ptptptpt(thr1c,thr2n,thr2ca,thr2c)));
}