#ifndef lint
static char *sccsid = "@(#)compute_abch.c	1.1 (UKC) 9/12/87";
#endif  lint
/*
 *	Given a list of triangles with A,B,C set, fill in the
 *	a, b, c, h fields which define the plane containing the triangle.
 *
 *	This is derived from page 92 of Ammeraal's book.
 */

#include <stdio.h>
#include "types.h"
extern double sqrt();

compute_abch()
{
	double xA, yA, zA;	/* copies of the coordinates of the triangle */
	double xB, yB, zB;
	double xC, yC, zC;
	double a, b, c, h, r;
	trilist_t *tp;

	for (tp=trilist; tp != NULLTRI; tp = tp->link) {
		xA = tp->A.x; yA = tp->A.y; zA = tp->A.z;
		xB = tp->B.x; yB = tp->B.y; zB = tp->B.z;
		xC = tp->C.x; yC = tp->C.y; zC = tp->C.z;
		
		a =   yA*(zB-zC) - yB*(zA-zC) + yC*(zA-zB);
		b = -(xA*(zB-zC) - xB*(zA-zC) + xC*(zA-zB));
		c =   xA*(yB-yC) - xB*(yA-yC) + xC*(yA-yB);
		h = xA * (yB*zC - yC*zB) -
		    xB * (yA*zC - yC*zA) +
		    xC * (yA*zB - yB*zA);
		
		if (h != 0.0) {
			r = sqrt(a*a + b*b + c*c);
			a = a/r; b = b/r; c = c/r; h = h/r;
		} else {
			/* the triangle is edge-on and so will be eliminated */
		}

		tp->a = a; tp->b = b; tp->c = c; tp->h = h;
	}
}