int pointplane( point p1, point p2, point p3, point rn, point rc, float b )
{

//   for points p1, p2, p3, return center and plane normal vector,
//     plus the scalar b, such that the equation of the plane is
//     rn@x = b.

float anorm;

rc = (p1 + p2 + p3) * (1./3.);
// original:  rn = p1^p2 + p3^p1 + p2^p3;
rn = (p2 - p1) ^ (p3 - p1);

anorm = sqrt(rn @ rn);
if( anorm > 0.0001 ){
    anorm = 1./anorm;
} else {
    fprintf( stderr, "anorm is too small in plane\n" );  
	fprintf( stderr, "p1: %8.3f %8.3f%8.3f\n", p1.x, p1.y, p1.z );
	fprintf( stderr, "p2: %8.3f %8.3f%8.3f\n", p2.x, p2.y, p2.z );
	fprintf( stderr, "p3: %8.3f %8.3f%8.3f\n", p3.x, p3.y, p3.z );
	fprintf( stderr, "rn: %8.3f %8.3f%8.3f\n", rn.x, rn.y, rn.z );
	fprintf( stderr, "rc: %8.3f %8.3f%8.3f\n", rc.x, rc.y, rc.z );
	exit( 1 );
}
rn = rn*anorm;
b = rn @ p1;
return( 0 );

};