Index: include/bn.h
===================================================================
--- include/bn.h	(revision 49726)
+++ include/bn.h	(working copy)
@@ -59,11 +59,12 @@
 
 
 /** @} */
+
 /*----------------------------------------------------------------------*/
+
 /** @addtogroup tol */
 /** @{ */
-/**
- * B N _ T O L
+/** B N _ T O L
  *
  * @brief Support for uniform tolerances
  *
@@ -101,6 +102,7 @@
  * tol->dist_sq = tol->dist * tol->dist;
  *@n tol->para = 1 - tol->perp;
  */
+
 struct bn_tol {
     uint32_t magic;
     double dist;		/**< @brief >= 0 */
@@ -153,7 +155,11 @@
 
 
 /** @} */
+
+
 /*----------------------------------------------------------------------*/
+
+
 /* anim.c */
 /** @addtogroup anim */
 /** @{ */
@@ -199,6 +205,12 @@
  */
 BN_EXPORT extern void anim_tran(mat_t m);
 
+
+/***************************************
+ *ANIM_MAT2* - Conversions from matrices
+ ***************************************/
+
+
 /**
  * @brief
  * Convert the rotational part of a 4x4 transformation matrix to zyx
@@ -233,6 +245,12 @@
 BN_EXPORT extern int anim_mat2quat(quat_t quat,
 				   const mat_t viewrot);
 
+
+/***************************************
+ *ANIM_*2MAT - Conversions to matrices
+ ***************************************/
+
+
 /**
  * @brief Create a premultiplication rotation matrix to turn the front
  * of an object (its x-axis) to the given yaw, pitch, and roll, which
@@ -379,6 +397,12 @@
 				    vect_t point,
 				    int end);
 
+
+/***************************************
+ * Other animation routines
+ ***************************************/
+
+
 /**
  * @brief Add pre- and post- translation to a rotation matrix.  The
  * resulting matrix has the effect of performing the first
@@ -420,7 +444,19 @@
 
 
 /*----------------------------------------------------------------------*/
-/* bn_tcl.c */
+
+
+/* tcl.c */
+
+/** @addtogroup bntcl */
+/** @{
+ *
+ * @brief
+ * Tcl interfaces to all the LIBBN math routines.
+ *
+ */
+
+
 BN_EXPORT extern int bn_decode_mat(mat_t m,
 				   const char *str);
 BN_EXPORT extern int bn_decode_quat(quat_t q,
@@ -439,15 +475,41 @@
 				      const hvect_t v);
 
 /* The presence of Tcl_Interp as an arg prevents giving arg list */
+
+/**
+ * B N _ T C L _ S E T U P
+ *@brief
+ * Add all the supported Tcl interfaces to LIBBN routines to
+ * the list of commands known by the given interpreter.
+ */
 BN_EXPORT extern void bn_tcl_setup();
+
+/**
+ * B N _ I N I T
+ *@brief
+ * Allows LIBBN to be dynamically loade to a vanilla tclsh/wish with
+ * "load /usr/brlcad/lib/libbn.so"
+ *
+ * The name of this function is specified by TCL.
+ */
 BN_EXPORT extern int Bn_Init();
+
+/**
+ * B N _ M A T _ P R I N T
+ */
 BN_EXPORT extern void bn_tcl_mat_print();
 
 /** @} */
+
+
 /*----------------------------------------------------------------------*/
+
+
 /* complex.c */
+
 /** @addtogroup complex */
 /** @{ */
+
 /*
  * Complex numbers
  */
@@ -508,22 +570,63 @@
 BN_EXPORT extern void bn_cx_sqrt(bn_complex_t *op,
 				 const bn_complex_t *ip);
 
+
 /*----------------------------------------------------------------------*/
+
+
 /* mat.c */
 /*
  * 4x4 Matrix math
  */
+
+/** @addtogroup mat */
+/** @{
+ *
+ * @brief
+ * 4 x 4 Matrix manipulation functions...
+ *
+ * Matrix array elements have the following positions in the matrix:
+ @code
+ *			|  0  1  2  3 |		| 0 |
+ *	[ 0 1 2 3 ]	|  4  5  6  7 |		| 1 |
+ *			|  8  9 10 11 |		| 2 |
+ *			| 12 13 14 15 |		| 3 |
+ *
+ * preVector (vect_t) Matrix (mat_t) postVector (vect_t)
+ @endcode
+ *
+ * TODO: need a better way to control tolerancing, either via
+ * additional tolerance parameters or perhaps providing a global
+ * tolerance state interface.
+ */
+
+
 BN_EXPORT extern const mat_t bn_mat_identity;
 
+/**
+ * B N _ M A T _ P R I N T
+ */
 BN_EXPORT extern void bn_mat_print(const char *title,
 				   const mat_t m);
+				   
 BN_EXPORT extern void bn_mat_print_guts(const char *title,
 					const mat_t m,
 					char *buf,
 					int buflen);
+					
+/**
+ * B N _ M A T _ P R I N T _ V L S
+ */
 BN_EXPORT extern void bn_mat_print_vls(const char *title,
 				       const mat_t m,
 				       struct bu_vls *vls);
+				       
+/**
+ * B N _ A T A N 2
+ *
+ * A wrapper for the system atan2().  On the Silicon Graphics, and
+ * perhaps on others, x==0 incorrectly returns infinity.
+ */
 BN_EXPORT extern double bn_atan2(double x, double y);
 
 #define bn_mat_zero(_m) { \
@@ -570,62 +673,221 @@
   #define bn_mat_copy(_d, _s)	(void)memcpy((void *)_d, (const void *)(_s), sizeof(mat_t))
 */
 
-
+/**
+ * B N _ M A T _ M U L
+ *
+ * Multiply matrix "a" by "b" and store the result in "o".
+ *
+ * This is different from multiplying "b" by "a" (most of the time!)
+ * Also, "o" must not be the same as either of the inputs.
+ */
 BN_EXPORT extern void bn_mat_mul(mat_t o,
 				 const mat_t a,
 				 const mat_t b);
+				 
+/**
+ * B N _ M A T _ M U L 2
+ *
+ * o = i * o
+ *
+ * A convenience wrapper for bn_mat_mul() to update a matrix in place.
+ * The argument ordering is confusing either way.
+ */
 BN_EXPORT extern void bn_mat_mul2(const mat_t i,
 				  mat_t o);
+				  
+/**
+ * B N _ M A T _ M U L 3
+ *
+ * o = a * b * c
+ *
+ * The output matrix may be the same as 'b' or 'c', but may not be
+ * 'a'.
+ */
 BN_EXPORT extern void bn_mat_mul3(mat_t o,
 				  const mat_t a,
 				  const mat_t b,
 				  const mat_t c);
+				  
+/**
+ * B N _ M A T _ M U L 4
+ *
+ * o = a * b * c * d
+ *
+ * The output matrix may be the same as any input matrix.
+ */
 BN_EXPORT extern void bn_mat_mul4(mat_t o,
 				  const mat_t a,
 				  const mat_t b,
 				  const mat_t c,
 				  const mat_t d);
+				  
+/**
+ * B N _ M A T X V E C
+ *
+ * Multiply the matrix "im" by the vector "iv" and store the result in
+ * the vector "ov".  Note this is post-multiply, and operates on
+ * 4-tuples.  Use MAT4X3VEC() to operate on 3-tuples.
+ */
 BN_EXPORT extern void bn_matXvec(hvect_t ov,
 				 const mat_t im,
 				 const hvect_t iv);
+
+/**
+ * B N _ M A T _ I N V
+ *
+ * The matrix pointed at by "input" is inverted and stored in the area
+ * pointed at by "output".
+ *
+ * Calls bu_bomb if matrix is singular.
+ */				 
 BN_EXPORT extern void bn_mat_inv(mat_t output,
 				 const mat_t input);
+
+/**
+ * B N _ M A T _ I N V E R S E
+ *
+ * The matrix pointed at by "input" is inverted and stored in the area
+ * pointed at by "output".
+ *
+ * Invert a 4-by-4 matrix using Algorithm 120 from ACM.  This is a
+ * modified Gauss-Jordan alogorithm Note: Inversion is done in place,
+ * with 3 work vectors
+ *
+ * @return 1 if OK.
+ * @return 0 if matrix is singular.
+ */				 
 BN_EXPORT extern int bn_mat_inverse(mat_t output,
 				    const mat_t input);
+				    
+/**
+ * B N _ V T O H _ M O V E
+ *
+ * Takes a pointer to a [x, y, z] vector, and a pointer to space for a
+ * homogeneous vector [x, y, z, w], and builds [x, y, z, 1].
+ */
 BN_EXPORT extern void bn_vtoh_move(vect_t h,
 				   const vect_t v);
+				   
+/**
+ * B N _ H T O V _ M O V E
+ *
+ * Takes a pointer to [x, y, z, w], and converts it to an ordinary
+ * vector [x/w, y/w, z/w].  Optimization for the case of w==1 is
+ * performed.
+ *
+ * FIXME: make tolerance configurable
+ */
 BN_EXPORT extern void bn_htov_move(vect_t v,
 				   const vect_t h);
+				   
+/**
+ * B N _ M A T _ T R N
+ */
 BN_EXPORT extern void bn_mat_trn(mat_t om,
 				 const mat_t im);
+				 
+/**
+ * B N _ M A T _ A E
+ *
+ * Compute a 4x4 rotation matrix given Azimuth and Elevation.
+ *
+ * Azimuth is +X, Elevation is +Z, both in degrees.
+ *
+ * Formula due to Doug Gwyn, BRL.
+ */
 BN_EXPORT extern void bn_mat_ae(mat_t m,
 				double azimuth,
 				double elev);
+				
+/**
+ * B N _ A E _ V E C
+ *
+ * Find the azimuth and elevation angles that correspond to the
+ * direction (not including twist) given by a direction vector.
+ */
 BN_EXPORT extern void bn_ae_vec(fastf_t *azp,
 				fastf_t *elp,
 				const vect_t v);
+				
+/**
+ * B N _ A E T _ V E C
+ *
+ * Find the azimuth, elevation, and twist from two vectors.  Vec_ae is
+ * in the direction of view (+z in mged view) and vec_twist points to
+ * the viewers right (+x in mged view).  Accuracy (degrees) is used to
+ * stabilze flutter between equivalent extremes of atan2(), and to
+ * snap twist to zero when elevation is near +/- 90
+ */
 BN_EXPORT extern void bn_aet_vec(fastf_t *az,
 				 fastf_t *el,
 				 fastf_t *twist,
 				 vect_t vec_ae,
 				 vect_t vec_twist,
 				 fastf_t accuracy);
+				 
+/**
+ * B N _ V E C _ A E
+ *
+ * Find a unit vector from the origin given azimuth and elevation.
+ */
 BN_EXPORT extern void bn_vec_ae(vect_t vec,
 				fastf_t az,
 				fastf_t el);
+				
+/**
+ * B N _ V E C _ A E D
+ *
+ * Find a vector from the origin given azimuth, elevation, and distance.
+ */
 BN_EXPORT extern void bn_vec_aed(vect_t vec,
 				 fastf_t az,
 				 fastf_t el,
 				 fastf_t dist);
 
+/**
+ * B N _ M A T _ A N G L E S
+ *
+ * This routine builds a Homogeneous rotation matrix, given alpha,
+ * beta, and gamma as angles of rotation, in degrees.
+ *
+ * Alpha is angle of rotation about the X axis, and is done third.
+ * Beta is angle of rotation about the Y axis, and is done second.
+ * Gamma is angle of rotation about Z axis, and is done first.
+ *
+ * FIXME: make tolerance configurable
+ */
 BN_EXPORT extern void bn_mat_angles(mat_t mat,
 				    double alpha,
 				    double beta, double ggamma);
+				    
+/**
+ * B N _ M A T _ A N G L E S _ R A D
+ *
+ * This routine builds a Homogeneous rotation matrix, given alpha,
+ * beta, and gamma as angles of rotation, in radians.
+ *
+ * Alpha is angle of rotation about the X axis, and is done third.
+ * Beta is angle of rotation about the Y axis, and is done second.
+ * Gamma is angle of rotation about Z axis, and is done first.
+ *
+ * FIXME: make tolerance configurable
+ */
 BN_EXPORT extern void bn_mat_angles_rad(mat_t mat,
 					double alpha,
 					double beta,
 					double ggamma);
 
+/**
+ * B N _ E I G E N 2 X 2
+ *
+ * Find the eigenvalues and eigenvectors of a symmetric 2x2 matrix.
+ * (a b)
+ * (b c)
+ *
+ * The eigenvalue with the smallest absolute value is returned in
+ * val1, with its eigenvector in vec1.
+ */
 BN_EXPORT extern void bn_eigen2x2(fastf_t *val1,
 				  fastf_t *val2,
 				  vect_t vec1,
@@ -634,47 +896,216 @@
 				  fastf_t b,
 				  fastf_t c);
 
+/**
+ * B N _ V E C _ P E R P
+ *
+ * Given a vector, create another vector which is perpendicular to it.
+ * The output vector will have unit length only if the input vector
+ * did.
+ *
+ * FIXME: make tolerance configurable
+ */
 BN_EXPORT extern void bn_vec_perp(vect_t new_vec,
 				  const vect_t old_vec);
+				  
+/**
+ * B N _ M A T _ F R O M T O
+ *
+ * Given two vectors, compute a rotation matrix that will transform
+ * space by the angle between the two.  There are many candidate
+ * matricies.
+ *
+ * The input 'from' and 'to' vectors need not be unit length.
+ * MAT4X3VEC(to, m, from) is the identity that is created.
+ *
+ * FIXME: make tolerance configurable
+ */
 BN_EXPORT extern void bn_mat_fromto(mat_t m,
 				    const vect_t from,
 				    const vect_t to);
+				    
+/**
+ * B N _ M A T _ X R O T
+ *
+ * Given the sin and cos of an X rotation angle, produce the rotation
+ * matrix.
+ */
 BN_EXPORT extern void bn_mat_xrot(mat_t m,
 				  double sinx,
 				  double cosx);
+				  
+/**
+ * B N _ M A T _ Y R O T
+ *
+ * Given the sin and cos of a Y rotation angle, produce the rotation
+ * matrix.
+ */
 BN_EXPORT extern void bn_mat_yrot(mat_t m,
 				  double siny,
 				  double cosy);
+				  
+/**
+ * B N _ M A T _ Z R O T
+ *
+ * Given the sin and cos of a Z rotation angle, produce the rotation
+ * matrix.
+ */
 BN_EXPORT extern void bn_mat_zrot(mat_t m,
 				  double sinz,
 				  double cosz);
+				  
+/**
+ * B N _ M A T _ L O O K A T
+ *
+ * Given a direction vector D of unit length, product a matrix which
+ * rotates that vector D onto the -Z axis.  This matrix will be
+ * suitable for use as a "model2view" matrix.
+ *
+ * XXX This routine will fail if the vector is already more or less
+ * aligned with the Z axis.
+ *
+ * This is done in several steps.
+ *
+ @code
+ 1) Rotate D about Z to match +X axis.  Azimuth adjustment.
+ 2) Rotate D about Y to match -Y axis.  Elevation adjustment.
+ 3) Rotate D about Z to make projection of X axis again point
+ in the +X direction.  Twist adjustment.
+ 4) Optionally, flip sign on Y axis if original Z becomes inverted.
+ This can be nice for static frames, but is astonishing when
+ used in animation.
+ @endcode
+*/
 BN_EXPORT extern void bn_mat_lookat(mat_t rot,
 				    const vect_t dir,
 				    int yflip);
+				    
+/**
+ * B N _ V E C _ O R T H O
+ *
+ * Given a vector, create another vector which is perpendicular to it,
+ * and with unit length.  This algorithm taken from Gift's arvec.f; a
+ * faster algorithm may be possible.
+ *
+ * FIXME: make tolerance configurable
+ */
 BN_EXPORT extern void bn_vec_ortho(vect_t out,
 				   const vect_t in);
+				   
+/**
+ * B N _ M A T _ S C A L E _ A B O U T _ P T
+ *
+ * Build a matrix to scale uniformly around a given point.
+ *
+ * @return -1 if scale is too small.
+ * @return  0 if OK.
+ *
+ * FIXME: make tolerance configurable
+ */
 BN_EXPORT extern int bn_mat_scale_about_pt(mat_t mat,
 					   const point_t pt,
 					   const double scale);
+					   
+/**
+ * B N _ M A T _ X F O R M _ A B O U T _ P T
+ *
+ * Build a matrix to apply arbitary 4x4 transformation around a given
+ * point.
+ */
 BN_EXPORT extern void bn_mat_xform_about_pt(mat_t mat,
 					    const mat_t xform,
 					    const point_t pt);
+					    
+/**
+ * B N _ M A T _ I S _ E Q U A L
+ *
+ * @return 0 When matrices are not equal
+ * @return 1 When matricies are equal
+ */
 BN_EXPORT extern int bn_mat_is_equal(const mat_t a,
 				     const mat_t b,
 				     const struct bn_tol *tol);
+				     
+/**
+ * B N _ M A T _ I S _ I D E N T I T Y
+ *
+ * This routine is intended for detecting identity matricies read in
+ * from ascii or binary files, where the numbers are pure ones or
+ * zeros.  This routine is *not* intended for tolerance-based
+ * "near-zero" comparisons; as such, it shouldn't be used on matrices
+ * which are the result of calculation.
+ *
+ *
+ * @return 0 non-identity
+ * @return 1 a perfect identity matrix
+ */
 BN_EXPORT extern int bn_mat_is_identity(const mat_t m);
+
+/**
+ * B N _ M A T _ A R B _ R O T
+ *
+ * Construct a transformation matrix for rotation about an arbitrary
+ * axis.  The axis is defined by a point (pt) and a unit direction
+ * vector (dir).  The angle of rotation is "ang"
+ *
+ * FIXME: make tolerance configurable
+ */
 BN_EXPORT extern void bn_mat_arb_rot(mat_t m,
 				     const point_t pt,
 				     const vect_t dir,
 				     const fastf_t ang);
+				     
+/**
+ * B N _ M A T _ D U P
+ *
+ * Return a pointer to a copy of the matrix in dynamically allocated
+ * memory.
+ */
 BN_EXPORT extern matp_t bn_mat_dup(const mat_t in);
+
+/**
+ * B N _ M A T _ C K
+ *
+ * Check to ensure that a rotation matrix preserves axis
+ * perpendicularily.  Note that not all matricies are rotation
+ * matricies.
+ *
+ *
+ * @return -1 FAIL
+ * @return  0 OK
+ */
 BN_EXPORT extern int bn_mat_ck(const char *title,
 			       const mat_t m);
+			       
+/**
+ * B N _ M A T _ D E T 3
+ *
+ * Calculates the determinant of the 3X3 "rotation" part of the passed
+ * matrix.
+ */
 BN_EXPORT extern fastf_t bn_mat_det3(const mat_t m);
+
+/**
+ * B N _ M A T _ D E T E R M I N A N T
+ *
+ * Calculates the determinant of the 4X4 matrix
+ */
 BN_EXPORT extern fastf_t bn_mat_determinant(const mat_t m);
 
+/**
+ * B N _ M A T _ I S _ N O N _ U N I F
+ *
+ * FIXME: make tolerance configurable
+ */
 BN_EXPORT extern int bn_mat_is_non_unif(const mat_t m);
 
+/**
+ * B N _ W R T _ P O I N T _ D I R E C
+ *
+ * Given a model-space transformation matrix "change",
+ * return a matrix which applies the change with-respect-to
+ * given "point" and "direc".
+ */
 BN_EXPORT extern void bn_wrt_point_direc(mat_t out,
 					 const mat_t change,
 					 const mat_t in,
@@ -682,10 +1113,25 @@
 					 const vect_t direc);
 
 /** @} */
+
+
 /*----------------------------------------------------------------------*/
+
+
 /* msr.c */
 /** @addtogroup msr */
-/** @{ */
+/** @{
+ *
+ * @brief
+ * Minimal Standard RANdom number generator
+ *
+ * @par From:
+ *	Stephen K. Park and Keith W. Miller
+ * @n	"Random number generators: good ones are hard to find"
+ * @n	CACM vol 31 no 10, Oct 88
+ *
+ */
+
 /*
  * Define data structures and constants for the "MSR" random number
  * package.
@@ -721,14 +1167,163 @@
     double *msr_gausses;
 };
 
+
+/*	bn_unif_init	Initialize a random number structure.
+ *
+ * @par Entry
+ *	setseed	seed to use
+ *	method  method to use to generate numbers;
+ *
+ * @par Exit
+ *	returns	a pointer to a bn_unif structure.
+ *	returns 0 on error.
+ *
+ * @par Uses
+ *	None.
+ *
+ * @par Calls
+ *	bu_malloc
+ *
+ * @par Method @code
+ *	malloc up a structure with pointers to the numbers
+ *	get space for the integer table.
+ *	get space for the floating table.
+ *	set pointer counters
+ *	set seed if one was given and setseed != 1
+ @endcode
+ *
+ */
 BN_EXPORT extern struct bn_unif *bn_unif_init(long setseed,
 					      int method);
+      
+/*	bn_unif_free	free random number table
+ *
+ */
 BN_EXPORT extern void bn_unif_free(struct bn_unif *p);
+
+/*	bn_unif_long_fill	fill a random number table.
+ *
+ * Use the msrad algorithm to fill a random number table
+ * with values from 1 to 2^31-1.  These numbers can (and are) extracted from
+ * the random number table via high speed macros and bn_unif_long_fill called
+ * when the table is exauseted.
+ *
+ * @par Entry
+ *	p	pointer to a bn_unif structure.
+ *
+ * @par Exit
+ *	if (!p) returns 1 else returns a value between 1 and 2^31-1
+ *
+ * @par Calls
+ *	None.  msran is inlined for speed reasons.
+ *
+ * @par Uses
+ *	None.
+ *
+ * @par Method @code
+ *	if (!p) return 1;
+ *	if p->msr_longs != NULL
+ *		msr_longs is reloaded with random numbers;
+ *		msr_long_ptr is set to BN_MSR_MAXTBL
+ *	endif
+ *	msr_seed is updated.
+ @endcode
+*/
 BN_EXPORT extern long bn_unif_long_fill(struct bn_unif *p);
+
+/*	bn_unif_double_fill	fill a random number table.
+ *
+ * Use the msrad algorithm to fill a random number table
+ * with values from -0.5 to 0.5.  These numbers can (and are) extracted from
+ * the random number table via high speed macros and bn_unif_double_fill
+ * called when the table is exauseted.
+ *
+ * @par Entry
+ *	p	pointer to a bn_unif structure.
+ *
+ * @par Exit
+ *	if (!p) returns 0.0 else returns a value between -0.5 and 0.5
+ *
+ * @par Calls
+ *	None.  msran is inlined for speed reasons.
+ *
+ * @par Uses
+ *	None.
+ *
+ * @par Method @code
+ *	if (!p) return (0.0)
+ *	if p->msr_longs != NULL
+ *		msr_longs is reloaded with random numbers;
+ *		msr_long_ptr is set to BN_MSR_MAXTBL
+ *	endif
+ *	msr_seed is updated.
+ @endcode
+*/
 BN_EXPORT extern double bn_unif_double_fill(struct bn_unif *p);
+
+/*	bn_gauss_init	Initialize a random number struct for gaussian
+ *	numbers.
+ *
+ * @par Entry
+ *	setseed		Seed to use.
+ *	method		method to use to generate numbers (not used)
+ *
+ * @par Exit
+ *	Returns a pointer toa bn_msr_guass structure.
+ *	returns 0 on error.
+ *
+ * @par Calls
+ *	bu_malloc
+ *
+ * @par Uses
+ *	None.
+ *
+ * @par Method @code
+ *	malloc up a structure
+ *	get table space
+ *	set seed and pointer.
+ *	if setseed != 0 then seed = setseed
+ @endcode
+*/
 BN_EXPORT extern struct bn_gauss *bn_gauss_init(long setseed,
 						int method);
+						
+/*	bn_gauss_free	free random number table
+ *
+ */
 BN_EXPORT extern void bn_gauss_free(struct bn_gauss *p);
+
+/*	bn_gauss_fill	fill a random number table.
+ *
+ * Use the msrad algorithm to fill a random number table.
+ * hese numbers can (and are) extracted from
+ * the random number table via high speed macros and bn_msr_guass_fill
+ * called when the table is exauseted.
+ *
+ * @par Entry
+ *	p	pointer to a bn_msr_guass structure.
+ *
+ * @par Exit
+ *	if (!p) returns 0.0 else returns a value with a mean of 0 and
+ *	    a variance of 1.0.
+ *
+ * @par Calls
+ *	BN_UNIF_CIRCLE to get to uniform random number whos radius is
+ *	<= 1.0. I.e. sqrt(v1*v1 + v2*v2) <= 1.0
+ *	BN_UNIF_CIRCLE is a macro which can call bn_unif_double_fill.
+ *
+ * @par Uses
+ *	None.
+ *
+ * @par Method @code
+ *	if (!p) return (0.0)
+ *	if p->msr_longs != NULL
+ *		msr_longs is reloaded with random numbers;
+ *		msr_long_ptr is set to BN_MSR_MAXTBL
+ *	endif
+ *	msr_seed is updated.
+ @endcode
+*/
 BN_EXPORT extern double bn_gauss_fill(struct bn_gauss *p);
 
 #define BN_UNIF_LONG(_p)	\
@@ -761,45 +1356,180 @@
      bn_gauss_fill(_p))
 
 /** @} */
+
+
 /*----------------------------------------------------------------------*/
+
+
 /* noise.c */
+
 /* @addtogroup noise */
-/** @{ */
+/** @{
+ *
+ * These noise functions provide mostly random noise at the integer
+ * lattice points.  The functions should be evaluated at non-integer
+ * locations for their nature to be realized.
+ *
+ * Conatins contributed code from:
+ * F. Kenton Musgrave
+ * Robert Skinner
+ *
+ */
+
 /*
  * fractal noise support
  */
 
+
 BN_EXPORT extern void bn_noise_init();
+
+/**
+ *@brief
+ * Robert Skinner's Perlin-style "Noise" function
+ *
+ * Results are in the range [-0.5 .. 0.5].  Unlike many
+ * implementations, this function provides random noise at the integer
+ * lattice values.  However this produces much poorer quality and
+ * should be avoided if possible.
+ *
+ * The power distribution of the result has no particular shape,
+ * though it isn't as flat as the literature would have one believe.
+ */
 BN_EXPORT extern double bn_noise_perlin(point_t pt);
-/* FIXME: Why isn't the result listed first? */
+
+/**
+ * Vector-valued "Noise"
+ */
 BN_EXPORT extern void bn_noise_vec(point_t point,
 				   point_t result);
+
+/**
+ * @brief
+ * Procedural fBm evaluated at "point"; returns value stored in
+ * "value".
+ *
+ * @param point          location to sample noise
+ * @param ``h_val''      fractal increment parameter
+ * @param ``lacunarity'' gap between successive frequencies
+ * @param ``octaves''  	 number of frequencies in the fBm
+ *
+ * The spectral properties of the result are in the APPROXIMATE range
+ * [-1..1] Depending upon the number of octaves computed, this range
+ * may be exceeded.  Applications should clamp or scale the result to
+ * their needs.  The results have a more-or-less gaussian
+ * distribution.  Typical results for 1M samples include:
+ *
+ * @li Min -1.15246
+ * @li Max  1.23146
+ * @li Mean -0.0138744
+ * @li s.d.  0.306642
+ * @li Var 0.0940295
+ *
+ * The function call pow() is relatively expensive.  Therfore, this
+ * function pre-computes and saves the spectral weights in a table for
+ * re-use in successive invocations.
+ */
 BN_EXPORT extern double bn_noise_fbm(point_t point,
 				     double h_val,
 				     double lacunarity,
 				     double octaves);
+
+/**
+ * @brief
+ * Procedural turbulence evaluated at "point";
+ *
+ * @return turbulence value for point
+ *
+ * @param point          location to sample noise at
+ * @param ``h_val''      fractal increment parameter
+ * @param ``lacunarity'' gap between successive frequencies
+ * @param ``octaves''    number of frequencies in the fBm
+ *
+ * The result is characterized by sharp, narrow trenches in low values
+ * and a more fbm-like quality in the mid-high values.  Values are in
+ * the APPROXIMATE range [0 .. 1] depending upon the number of octaves
+ * evaluated.  Typical results:
+ @code
+ * Min  0.00857137
+ * Max  1.26712
+ * Mean 0.395122
+ * s.d. 0.174796
+ * Var  0.0305536
+ @endcode
+ * The function call pow() is relatively expensive.  Therfore, this
+ * function pre-computes and saves the spectral weights in a table for
+ * re-use in successive invocations.
+ */
 BN_EXPORT extern double bn_noise_turb(point_t point,
 				      double h_val,
 				      double lacunarity,
 				      double octaves);
+
+/**
+ *
+ * From "Texturing and Modeling, A Procedural Approach" 2nd ed
+ */
 BN_EXPORT extern double bn_noise_mf(point_t point,
 				    double h_val,
 				    double lacunarity,
 				    double octaves,
 				    double offset);
+
+/**
+ *@brief
+ * A ridged noise pattern
+ *
+ * From "Texturing and Modeling, A Procedural Approach" 2nd ed p338
+ */
 BN_EXPORT extern double bn_noise_ridged(point_t point,
 					double h_val,
 					double lacunarity,
 					double octaves,
 					double offset);
 
+
 /*----------------------------------------------------------------------*/
+
+
 /* plane.c */
 /*
  * Plane/line/point calculations
  */
 
+/** @addtogroup plane */
+/** @{
+ *
+ * @brief
+ * Some useful routines for dealing with planes and lines.
+ *
+ */
 
+
+
+/** B N _ D I S T S Q _ L I N E 3 _ L I N E 3
+ *@brief
+ * Calculate the square of the distance of closest approach for two
+ * lines.
+ *
+ * The lines are specifed as a point and a vector each.  The vectors
+ * need not be unit length.  P and d define one line; Q and e define
+ * the other.
+ *
+ * @return 0 - normal return
+ * @return 1 - lines are parallel, dist[0] is set to 0.0
+ *
+ * Output values:
+ * dist[0] is the parametric distance along the first line P + dist[0] * d of the PCA
+ * dist[1] is the parametric distance along the second line Q + dist[1] * e of the PCA
+ * dist[3] is the square of the distance between the points of closest approach
+ * pt1 is the point of closest approach on the first line
+ * pt2 is the point of closest approach on the second line
+ *
+ * This algoritm is based on expressing the distance sqaured, taking
+ * partials with respect to the two unknown parameters (dist[0] and
+ * dist[1]), seeting the two partails equal to 0, and solving the two
+ * simutaneous equations
+ */
 BN_EXPORT extern int bn_distsq_line3_line3(fastf_t dist[3],
 					   point_t P,
 					   vect_t d,
@@ -808,6 +1538,31 @@
 					   point_t pt1,
 					   point_t pt2);
 
+/**
+ * Find the distance from a point P to a line described by the
+ * endpoint A and direction dir, and the point of closest approach
+ * (PCA).
+ *
+ @code
+ 			P
+ 		       *
+ 		      /.
+ 		     / .
+ 		    /  .
+ 		   /   . (dist)
+ 		  /    .
+ 		 /     .
+ 		*------*-------->
+ 		A      PCA    dir
+ @endcode
+ * There are three distinct cases, with these return codes -
+ *   0 => P is within tolerance of point A.  *dist = 0, pca=A.
+ *   1 => P is within tolerance of line.  *dist = 0, pca=computed.
+ *   2 => P is "above/below" line.  *dist=|PCA-P|, pca=computed.
+ *
+ * TODO: For efficiency, a version of this routine that provides the
+ * distance squared would be faster.
+ */
 BN_EXPORT extern int bn_dist_pt3_line3(fastf_t *dist,
 				       point_t pca,
 				       const point_t a,
@@ -815,6 +1570,27 @@
 				       const vect_t dir,
 				       const struct bn_tol *tol);
 
+/**
+ * calculate intersection or closest approach of a line and a line
+ * segement.
+ *
+ * returns:
+ *	-2 -> line and line segment are parallel and collinear.
+ *	-1 -> line and line segment are parallel, not collinear.
+ *	 0 -> intersection between points a and b.
+ *	 1 -> intersection outside a and b.
+ *	 2 -> closest approach is between a and b.
+ *	 3 -> closest approach is outside a and b.
+ *
+ * dist[0] is actual distance from p in d direction to
+ * closest portion of segment.
+ * dist[1] is ratio of distance from a to b (0.0 at a, and 1.0 at b),
+ * dist[1] may be less than 0 or greater than 1.
+ * For return values less than 0, closest approach is defined as
+ * closest to point p.
+ * Direction vector, d, must be unit length.
+ *
+ */
 BN_EXPORT extern int bn_dist_line3_lseg3(fastf_t *dist,
 					 const fastf_t *p,
 					 const fastf_t *d,
@@ -822,6 +1598,25 @@
 					 const fastf_t *b,
 					 const struct bn_tol *tol);
 
+/**
+ * Calculate closest approach of two lines
+ *
+ * returns:
+ *	-2 -> lines are parallel and do not intersect
+ *	-1 -> lines are parallel and collinear
+ *	 0 -> lines intersect
+ *	 1 -> lines do not intersect
+ *
+ * For return values less than zero, dist is not set.  For return
+ * valuse of 0 or 1, dist[0] is the distance from p1 in the d1
+ * direction to the point of closest approach for that line.  Similar
+ * for the second line.  d1 and d2 must be unit direction vectors.
+ *
+ * XXX How is this different from bn_isect_line3_line3() ?
+ * XXX Why are the calling sequences just slightly different?
+ * XXX Can we pick the better one, and get rid of the other one?
+ * XXX If not, can we document how they differ?
+ */
 BN_EXPORT extern int bn_dist_line3_line3(fastf_t dist[2],
 					 const point_t p1,
 					 const point_t p2,
@@ -829,93 +1624,561 @@
 					 const vect_t d2,
 					 const struct bn_tol *tol);
 
+
+/**
+ * B N _ D I S T _ P T 3 _ L S E G 3
+ *@brief
+ * Find the distance from a point P to a line segment described by the
+ * two endpoints A and B, and the point of closest approach (PCA).
+ @verbatim
+ *
+ *			P
+ *		       *
+ *		      /.
+ *		     / .
+ *		    /  .
+ *		   /   . (dist)
+ *		  /    .
+ *		 /     .
+ *		*------*--------*
+ *		A      PCA	B
+ @endverbatim
+ *
+ * @return 0	P is within tolerance of lseg AB.  *dist isn't 0: (SPECIAL!!!)
+ *		  *dist = parametric dist = |PCA-A| / |B-A|.  pca=computed.
+ * @return 1	P is within tolerance of point A.  *dist = 0, pca=A.
+ * @return 2	P is within tolerance of point B.  *dist = 0, pca=B.
+ * @return 3	P is to the "left" of point A.  *dist=|P-A|, pca=A.
+ * @return 4	P is to the "right" of point B.  *dist=|P-B|, pca=B.
+ * @return 5	P is "above/below" lseg AB.  *dist=|PCA-P|, pca=computed.
+ *
+ * This routine was formerly called bn_dist_pt_lseg().
+ *
+ * XXX For efficiency, a version of this routine that provides the
+ * XXX distance squared would be faster.
+ */
 BN_EXPORT extern int bn_dist_pt3_lseg3(fastf_t *dist,
 				       point_t pca,
 				       const point_t a,
 				       const point_t b,
 				       const point_t p,
 				       const struct bn_tol *tol);
+				
+/**
+ * B N _ 3 P T S _ C O L L I N E A R
+ * @brief
+ * Check to see if three points are collinear.
+ *
+ * The algorithm is designed to work properly regardless of the order
+ * in which the points are provided.
+ *
+ * @return 1	If 3 points are collinear
+ * @return 0	If they are not
+ */
 BN_EXPORT extern int bn_3pts_collinear(point_t a,
 				       point_t b,
 				       point_t c,
 				       const struct bn_tol *tol);
+				       
+/**
+ * B N _ P T 3 _ P T 3 _ E Q U A L
+ *
+ * @return 1	if the two points are equal, within the tolerance
+ * @return 0	if the two points are not "the same"
+ */
 BN_EXPORT extern int bn_pt3_pt3_equal(const point_t a,
 				      const point_t b,
 				      const struct bn_tol *tol);
+
+/**
+ * B N _ D I S T _ P T 2 _ L S E G 2
+ *@brief
+ * Find the distance from a point P to a line segment described by the
+ * two endpoints A and B, and the point of closest approach (PCA).
+ @verbatim
+ *			P
+ *		       *
+ *		      /.
+ *		     / .
+ *		    /  .
+ *		   /   . (dist)
+ *		  /    .
+ *		 /     .
+ *		*------*--------*
+ *		A      PCA	B
+ @endverbatim
+ * There are six distinct cases, with these return codes -
+ * @return 0	P is within tolerance of lseg AB.  *dist isn't 0: (SPECIAL!!!)
+ *		  *dist = parametric dist = |PCA-A| / |B-A|.  pca=computed.
+ * @return 1	P is within tolerance of point A.  *dist = 0, pca=A.
+ * @return 2	P is within tolerance of point B.  *dist = 0, pca=B.
+ * @return 3	P is to the "left" of point A.  *dist=|P-A|**2, pca=A.
+ * @return 4	P is to the "right" of point B.  *dist=|P-B|**2, pca=B.
+ * @return 5	P is "above/below" lseg AB.  *dist=|PCA-P|**2, pca=computed.
+ *
+ *
+ * Patterned after bn_dist_pt3_lseg3().
+ */
 BN_EXPORT extern int bn_dist_pt2_lseg2(fastf_t *dist_sq,
 				       fastf_t pca[2],
 				       const point_t a,
 				       const point_t b,
 				       const point_t p,
 				       const struct bn_tol *tol);
+
+/**
+ * B N _ I S E C T _ L S E G 3 _ L S E G 3
+ *@brief
+ * Intersect two 3D line segments, defined by two points and two
+ * vectors.  The vectors are unlikely to be unit length.
+ *
+ *
+ * @return -3	missed
+ * @return -2	missed (line segments are parallel)
+ * @return -1	missed (colinear and non-overlapping)
+ * @return 0	hit (line segments colinear and overlapping)
+ * @return 1	hit (normal intersection)
+ *
+ * @param[out] dist
+ *	The value at dist[] is set to the parametric distance of the
+ *	intercept dist[0] is parameter along p, range 0 to 1, to
+ *	intercept.  dist[1] is parameter along q, range 0 to 1, to
+ *	intercept.  If within distance tolerance of the endpoints,
+ *	these will be exactly 0.0 or 1.0, to ease the job of caller.
+ *
+ *      CLARIFICATION: This function 'bn_isect_lseg3_lseg3'
+ *      returns distance values scaled where an intersect at the start
+ *      point of the line segement (within tol->dist) results in 0.0
+ *      and when the intersect is at the end point of the line
+ *      segement (within tol->dist), the result is 1.0.  Intersects
+ *      before the start point return a negative distance.  Intersects
+ *      after the end point result in a return value > 1.0.
+ *
+ * Special note: when return code is "0" for co-linearity, dist[1] has
+ * an alternate interpretation: it's the parameter along p (not q)
+ * which takes you from point p to the point (q + qdir), i.e., it's
+ * the endpoint of the q linesegment, since in this case there may be
+ * *two* intersections, if q is contained within span p to (p + pdir).
+ *
+ * @param p	point 1
+ * @param pdir	direction-1
+ * @param q	point 2
+ * @param qdir	direction-2
+ * @param tol	tolerance values
+ */
 BN_EXPORT extern int bn_isect_lseg3_lseg3(fastf_t *dist,
 					      const point_t p, const vect_t pdir,
 					      const point_t q, const vect_t qdir,
 					      const struct bn_tol *tol);
+
+/**
+ * B N _ I S E C T _ L I N E 3 _ L I N E 3
+ *
+ * Intersect two line segments, each in given in parametric form:
+ *
+ * X = p0 + pdist * pdir_i   (i.e. line p0->p1)
+ * and
+ * X = q0 + qdist * qdir_i   (i.e. line q0->q1)
+ *
+ * The input vectors 'pdir_i' and 'qdir_i' must NOT be unit vectors
+ * for this function to work correctly. The magnitude of the direction
+ * vectors indicates line segment length.
+ *
+ * The 'pdist' and 'qdist' values returned from this function are the
+ * actual distance to the intersect (i.e. not scaled). Distances may
+ * be negative, see below.
+ *
+ * @return -2	no intersection, lines are parallel.
+ * @return -1	no intersection
+ * @return 0	lines are co-linear (pdist and qdist returned) (see below)
+ * @return 1	intersection found  (pdist and qdist returned) (see below)
+ *
+ * @param	p0	point 1
+ * @param	pdir_i 	direction 1
+ * @param	q0	point 2
+ * @param	qdir_i 	direction 2
+ * @param tol	tolerance values
+ * @param[out]	pdist, qdist (distances to intersection) (see below)
+ *
+ *		When return = 1, pdist is the distance along line p0->p1 to the
+ *		intersect with line q0->q1. If the intersect is along p0->p1 but
+ *		in the opposite direction of vector pdir_i (i.e. occuring before
+ *		p0 on line p0->p1) then the distance will be negative. The value
+ *		if qdist is the same as pdist except it is the distance along line
+ *		q0->q1 to the intersect with line p0->p1.
+ *
+ *		When return code = 0 for co-linearity, pdist and qdist have a
+ *		different meaning. pdist is the distance from point p0 to point q0
+ *		and qdist is the distance from point p0 to point q1. If point q0
+ *		occurs before point p0 on line segment p0->p1 then pdist will be
+ *		negative. The same occurs for the distance to point q1.
+ */
 BN_EXPORT extern int bn_isect_line3_line3(fastf_t *s, fastf_t *t,
 					      const point_t p0,
 					      const vect_t u,
 					      const point_t q0,
 					      const vect_t v,
 					      const struct bn_tol *tol);
+
+/**
+ * B N _ 2 L I N E 3 _ C O L I N E A R
+ * @brief
+ * Returns non-zero if the 3 lines are colinear to within tol->dist
+ * over the given distance range.
+ *
+ * Range should be at least one model diameter for most applications.
+ * 1e5 might be OK for a default for "vehicle sized" models.
+ *
+ * The direction vectors do not need to be unit length.
+ */
 BN_EXPORT extern int bn_2line3_colinear(const point_t p1,
 					const vect_t d1,
 					const point_t p2,
 					const vect_t d2,
 					double range,
 					const struct bn_tol *tol);
+
+/**
+ * B N _ I S E C T _ P T 2 _ L S E G 2
+ * @brief
+ * Intersect a point P with the line segment defined by two distinct
+ * points A and B.
+ *
+ * @return -2	P on line AB but outside range of AB,
+ *			dist = distance from A to P on line.
+ * @return -1	P not on line of AB within tolerance
+ * @return 1	P is at A
+ * @return 2	P is at B
+ * @return 3	P is on AB, dist = distance from A to P on line.
+ @verbatim
+ B *
+ |
+ P'*-tol-*P
+ |    /  _
+ dist  /   /|
+ |  /   /
+ | /   / AtoP
+ |/   /
+ A *   /
+
+ tol = distance limit from line to pt P;
+ dist = distance from A to P'
+ @endverbatim
+*/					
 BN_EXPORT extern int bn_isect_pt2_lseg2(fastf_t *dist,
 					const point_t a,
 					const point_t b,
 					const point_t p,
 					const struct bn_tol *tol);
+
+/**
+ * B N _ I S E C T _ L I N E 2 _ L S E G 2
+ *@brief
+ * Intersect a line in parametric form:
+ *
+ * X = P + t * D
+ *
+ * with a line segment defined by two distinct points A and B=(A+C).
+ *
+ * XXX probably should take point B, not vector C.  Sigh.
+ *
+ * @return -4	A and B are not distinct points
+ * @return -3	Lines do not intersect
+ * @return -2	Intersection exists, but outside segemnt, < A
+ * @return -1	Intersection exists, but outside segment, > B
+ * @return 0	Lines are co-linear (special meaning of dist[1])
+ * @return 1	Intersection at vertex A
+ * @return 2	Intersection at vertex B (A+C)
+ * @return 3	Intersection between A and B
+ *
+ * Implicit Returns -
+ * @param dist	When explicit return >= 0, t is the parameter that describes
+ *		the intersection of the line and the line segment.
+ *		The actual intersection coordinates can be found by
+ *		solving P + t * D.  However, note that for return codes
+ *		1 and 2 (intersection exactly at a vertex), it is
+ *		strongly recommended that the original values passed in
+ *		A or B are used instead of solving P + t * D, to prevent
+ *		numeric error from creeping into the position of
+ *		the endpoints.
+ *
+ * @param p	point of first line
+ * @param d	direction of first line
+ * @param a	point of second line
+ * @param c	direction of second line
+ * @param tol	tolerance values
+ */
 BN_EXPORT extern int bn_isect_line2_lseg2(fastf_t *dist,
 					  const point_t p,
 					  const vect_t d,
 					  const point_t a,
 					  const vect_t c,
 					  const struct bn_tol *tol);
+
+/**
+ * B N _ I S E C T _ L S E G 2 _ L S E G 2
+ *@brief
+ * Intersect two 2D line segments, defined by two points and two
+ * vectors.  The vectors are unlikely to be unit length.
+ *
+ * @return -2	missed (line segments are parallel)
+ * @return -1	missed (colinear and non-overlapping)
+ * @return 0	hit (line segments colinear and overlapping)
+ * @return 1	hit (normal intersection)
+ *
+ * @param dist  The value at dist[] is set to the parametric distance of the
+ *		intercept.
+ *@n	dist[0] is parameter along p, range 0 to 1, to intercept.
+ *@n	dist[1] is parameter along q, range 0 to 1, to intercept.
+ *@n	If within distance tolerance of the endpoints, these will be
+ *	exactly 0.0 or 1.0, to ease the job of caller.
+ *
+ * Special note: when return code is "0" for co-linearity, dist[1] has
+ * an alternate interpretation: it's the parameter along p (not q)
+ * which takes you from point p to the point (q + qdir), i.e., its
+ * the endpoint of the q linesegment, since in this case there may be
+ * *two* intersections, if q is contained within span p to (p + pdir).
+ * And either may be -10 if the point is outside the span.
+ *
+ * @param p	point 1
+ * @param pdir	direction1
+ * @param q	point 2
+ * @param qdir	direction2
+ * @param tol	tolerance values
+ */
 BN_EXPORT extern int bn_isect_lseg2_lseg2(fastf_t *dist,
 					  const point_t p,
 					  const vect_t pdir,
 					  const point_t q,
 					  const vect_t qdir,
 					  const struct bn_tol *tol);
+					  
+/**
+ * B N _ I S E C T _ L I N E 2 _ L I N E 2
+ *
+ * Intersect two lines, each in given in parametric form:
+ @verbatim
+
+ X = P + t * D
+ and
+ X = A + u * C
+
+ @endverbatim
+ *
+ * While the parametric form is usually used to denote a ray (ie,
+ * positive values of the parameter only), in this case the full line
+ * is considered.
+ *
+ * The direction vectors C and D need not have unit length.
+ *
+ * @return -1	no intersection, lines are parallel.
+ * @return 0	lines are co-linear
+ *@n			dist[0] gives distance from P to A,
+ *@n			dist[1] gives distance from P to (A+C) [not same as below]
+ * @return 1	intersection found (t and u returned)
+ *@n			dist[0] gives distance from P to isect,
+ *@n			dist[1] gives distance from A to isect.
+ *
+ * @param dist When explicit return > 0, dist[0] and dist[1] are the
+ * line parameters of the intersection point on the 2 rays.  The
+ * actual intersection coordinates can be found by substituting either
+ * of these into the original ray equations.
+ *
+ * @param p	point of first line
+ * @param d	direction of first line
+ * @param a	point of second line
+ * @param c	direction of second line
+ * @param tol	tolerance values
+ *
+ * Note that for lines which are very nearly parallel, but not quite
+ * parallel enough to have the determinant go to "zero", the
+ * intersection can turn up in surprising places.  (e.g. when
+ * det=1e-15 and det1=5.5e-17, t=0.5)
+ */
 BN_EXPORT extern int bn_isect_line2_line2(fastf_t *dist,
 					  const point_t p,
 					  const vect_t d,
 					  const point_t a,
 					  const vect_t c,
 					  const struct bn_tol *tol);
+					  
+/**
+ * B N _ D I S T _ P T 3 _ P T 3
+ * @brief
+ * Returns distance between two points.
+ */
 BN_EXPORT extern double bn_dist_pt3_pt3(const point_t a,
 					const point_t b);
+					
+/**
+ * B N _ 3 P T S _ D I S T I N C T
+ *
+ * Check to see if three points are all distinct, i.e., ensure that
+ * there is at least sqrt(dist_tol_sq) distance between every pair of
+ * points.
+ *
+ * @return 1 If all three points are distinct
+ * @return 0 If two or more points are closer together than dist_tol_sq
+ */
 BN_EXPORT extern int bn_3pts_distinct(const point_t a,
 				      const point_t b,
 				      const point_t c,
 				      const struct bn_tol *tol);
+				      
+/**
+ * B N _ N P T S _ D I S T I N C T
+ *
+ * Check to see if the points are all distinct, i.e., ensure that
+ * there is at least sqrt(dist_tol_sq) distance between every pair of
+ * points.
+ *
+ * @return 1 If all the points are distinct
+ * @return 0 If two or more points are closer together than dist_tol_sq
+ */
 BN_EXPORT extern int bn_npts_distinct(const int npts,
 				      const point_t *pts,
 				      const struct bn_tol *tol);
+				      
+/**
+ * B N _ M K _ P L A N E _ 3 P T S
+ *
+ * Find the equation of a plane that contains three points.  Note that
+ * normal vector created is expected to point out (see vmath.h), so
+ * the vector from A to C had better be counter-clockwise (about the
+ * point A) from the vector from A to B.  This follows the BRL-CAD
+ * outward-pointing normal convention, and the right-hand rule for
+ * cross products.
+ *
+ @verbatim
+ *
+ *			C
+ *	                *
+ *	                |\
+ *	                | \
+ *	   ^     N      |  \
+ *	   |      \     |   \
+ *	   |       \    |    \
+ *	   |C-A     \   |     \
+ *	   |         \  |      \
+ *	   |          \ |       \
+ *	               \|        \
+ *	                *---------*
+ *	                A         B
+ *			   ----->
+ *		            B-A
+ @endverbatim
+ *
+ * If the points are given in the order A B C (e.g.,
+ * *counter*-clockwise), then the outward pointing surface normal:
+ *
+ * N = (B-A) x (C-A).
+ *
+ *  @return 0	OK
+ *  @return -1	Failure.  At least two of the points were not distinct,
+ *		or all three were colinear.
+ *
+ * @param[out]	plane	The plane equation is stored here.
+ * @param[in]	a	point 1
+ * @param[in]	b	point 2
+ * @param[in]	c	point 3
+ * @param[in]	tol	Tolerance values for doing calcualtion
+ */
 BN_EXPORT extern int bn_mk_plane_3pts(plane_t plane,
 				      const point_t a,
 				      const point_t b,
 				      const point_t c,
 				      const struct bn_tol *tol);
+
+/**
+ * B N _ M K P O I N T _ 3 P L A N E S
+ *@brief
+ * Given the description of three planes, compute the point of intersection, if
+ * any. The direction vectors of the planes need not be of unit length.
+ *
+ * Find the solution to a system of three equations in three unknowns:
+ @verbatim
+ * Px * Ax + Py * Ay + Pz * Az = -A3;
+ * Px * Bx + Py * By + Pz * Bz = -B3;
+ * Px * Cx + Py * Cy + Pz * Cz = -C3;
+ *
+ * OR
+ *
+ * [ Ax  Ay  Az ]   [ Px ]   [ -A3 ]
+ * [ Bx  By  Bz ] * [ Py ] = [ -B3 ]
+ * [ Cx  Cy  Cz ]   [ Pz ]   [ -C3 ]
+ *
+ @endverbatim
+ *
+ * @return 0	OK
+ * @return -1	Failure.  Intersection is a line or plane.
+ *
+ * @param[out] pt	The point of intersection is stored here.
+ * @param	a	plane 1
+ * @param	b	plane 2
+ * @param	c	plane 3
+ */	      
 BN_EXPORT extern int bn_mkpoint_3planes(point_t pt,
 					const plane_t a,
 					const plane_t b,
 					const plane_t c);
+	
+/**
+ * B N _ I S E C T _ L I N E 3 _ P L A N E
+ *
+ * Intersect an infinite line (specified in point and direction vector
+ * form) with a plane that has an outward pointing normal.  The
+ * direction vector need not have unit length.  The first three
+ * elements of the plane equation must form a unit lengh vector.
+ *
+ * @return -2	missed (ray is outside halfspace)
+ * @return -1	missed (ray is inside)
+ * @return 0	line lies on plane
+ * @return 1	hit (ray is entering halfspace)
+ * @return 2	hit (ray is leaving)
+ *
+ * @param[out]	dist	set to the parametric distance of the intercept
+ * @param[in]	pt	origin of ray
+ * @param[in]	dir	direction of ray
+ * @param[in]	plane	equation of plane
+ * @param[in]	tol	tolerance values
+ */				
 BN_EXPORT extern int bn_isect_line3_plane(fastf_t *dist,
 					  const point_t pt,
 					  const vect_t dir,
 					  const plane_t plane,
 					  const struct bn_tol *tol);
+					  
+/**
+ * B N _ I S E C T _ 2 P L A N E S
+ *@brief
+ * Given two planes, find the line of intersection between them, if
+ * one exists.  The line of intersection is returned in parametric
+ * line (point & direction vector) form.
+ *
+ * In order that all the geometry under consideration be in "front" of
+ * the ray, it is necessary to pass the minimum point of the model
+ * RPP.  If this convention is unnecessary, just pass (0, 0, 0) as
+ * rpp_min.
+ *
+ * @return 0	success, line of intersection stored in 'pt' and 'dir'
+ * @return -1	planes are coplanar
+ * @return -2	planes are parallel but not coplanar
+ * @return -3	error, should be intersection but unable to find
+ *
+ * @param[out]	pt	Starting point of line of intersection
+ * @param[out]	dir	Direction vector of line of intersection (unit length)
+ * @param[in]	a	plane 1 (normal must be unit length)
+ * @param[in]	b	plane 2 (normal must be unit length)
+ * @param[in]	rpp_min	minimum poit of model RPP
+ * @param[in]	tol	tolerance values
+ */
 BN_EXPORT extern int bn_isect_2planes(point_t pt,
 				      vect_t dir,
 				      const plane_t a,
 				      const plane_t b,
 				      const vect_t rpp_min,
 				      const struct bn_tol *tol);
+				      
 BN_EXPORT extern int bn_isect_2lines(fastf_t *t,
 				     fastf_t *u,
 				     const point_t p,
@@ -923,68 +2186,309 @@
 				     const point_t a,
 				     const vect_t c,
 				     const struct bn_tol *tol);
+				     
+/**
+ * B N _ I S E C T _ L I N E _ L S E G
+ *@brief
+ * Intersect a line in parametric form:
+ *
+ * X = P + t * D
+ *
+ * with a line segment defined by two distinct points A and B.
+ *
+ *
+ * @return -4	A and B are not distinct points
+ * @return -3	Intersection exists, < A (t is returned)
+ * @return -2	Intersection exists, > B (t is returned)
+ * @return -1	Lines do not intersect
+ * @return 0	Lines are co-linear (t for A is returned)
+ * @return 1	Intersection at vertex A
+ * @return 2	Intersection at vertex B
+ * @return 3	Intersection between A and B
+ *
+ * @par Implicit Returns -
+ *
+ * t When explicit return >= 0, t is the parameter that describes the
+ * intersection of the line and the line segment.  The actual
+ * intersection coordinates can be found by solving P + t * D.
+ * However, note that for return codes 1 and 2 (intersection exactly
+ * at a vertex), it is strongly recommended that the original values
+ * passed in A or B are used instead of solving P + t * D, to prevent
+ * numeric error from creeping into the position of the endpoints.
+ *
+ * XXX should probably be called bn_isect_line3_lseg3()
+ * XXX should probably be changed to return dist[2]
+ */
 BN_EXPORT extern int bn_isect_line_lseg(fastf_t *t, const point_t p,
 					const vect_t d,
 					const point_t a,
 					const point_t b,
 					const struct bn_tol *tol);
+					
+/**
+ * B N _ D I S T _ L I N E 3_ P T 3
+ *@brief
+ * Given a parametric line defined by PT + t * DIR and a point A,
+ * return the closest distance between the line and the point.
+ *
+ *  'dir' need not have unit length.
+ *
+ * Find parameter for PCA along line with unitized DIR:
+ * d = VDOT(f, dir) / MAGNITUDE(dir);
+ * Find distance g from PCA to A using Pythagoras:
+ * g = sqrt(MAGSQ(f) - d**2)
+ *
+ * Return -
+ * Distance
+ */
 BN_EXPORT extern double bn_dist_line3_pt3(const point_t pt,
 					  const vect_t dir,
 					  const point_t a);
+					  
+/**
+ * B N _ D I S T S Q _ L I N E 3 _ P T 3
+ *
+ * Given a parametric line defined by PT + t * DIR and a point A,
+ * return the square of the closest distance between the line and the
+ * point.
+ *
+ * 'dir' need not have unit length.
+ *
+ * Return -
+ * Distance squared
+ */
 BN_EXPORT extern double bn_distsq_line3_pt3(const point_t pt,
 					    const vect_t dir,
 					    const point_t a);
+					    
+/**
+ * B N _ D I S T _ L I N E _ O R I G I N
+ *@brief
+ * Given a parametric line defined by PT + t * DIR, return the closest
+ * distance between the line and the origin.
+ *
+ * 'dir' need not have unit length.
+ *
+ * @return Distance
+ */
 BN_EXPORT extern double bn_dist_line_origin(const point_t pt,
 					    const vect_t dir);
+					    
+/**
+ * B N _ D I S T _ L I N E 2 _ P O I N T 2
+ *@brief
+ * Given a parametric line defined by PT + t * DIR and a point A,
+ * return the closest distance between the line and the point.
+ *
+ * 'dir' need not have unit length.
+ *
+ * @return Distance
+ */
 BN_EXPORT extern double bn_dist_line2_point2(const point_t pt,
 					     const vect_t dir,
 					     const point_t a);
+					     
+/**
+ * B N _ D I S T S Q _ L I N E 2 _ P O I N T 2
+ *@brief
+ * Given a parametric line defined by PT + t * DIR and a point A,
+ * return the closest distance between the line and the point,
+ * squared.
+ *
+ * 'dir' need not have unit length.
+ *
+ * @return
+ * Distance squared
+ */
 BN_EXPORT extern double bn_distsq_line2_point2(const point_t pt,
 					       const vect_t dir,
 					       const point_t a);
+					       
+/**
+ * B N _ A R E A _ O F _ T R I A N G L E
+ *@brief
+ * Returns the area of a triangle. Algorithm by Jon Leech 3/24/89.
+ */
 BN_EXPORT extern double bn_area_of_triangle(const point_t a,
 					    const point_t b,
 					    const point_t c);
+					    
+/**
+ * B N _ I S E C T _ P T _ L S E G
+ *@brief
+ * Intersect a point P with the line segment defined by two distinct
+ * points A and B.
+ *
+ * @return -2	P on line AB but outside range of AB,
+ * 			dist = distance from A to P on line.
+ * @return -1	P not on line of AB within tolerance
+ * @return 1	P is at A
+ * @return 2	P is at B
+ * @return 3	P is on AB, dist = distance from A to P on line.
+ @verbatim
+ B *
+ |
+ P'*-tol-*P
+ |    /   _
+ dist  /   /|
+ |  /   /
+ | /   / AtoP
+ |/   /
+ A *   /
+
+ tol = distance limit from line to pt P;
+ dist = parametric distance from A to P' (in terms of A to B)
+ @endverbatim
+ *
+ * @param p	point
+ * @param a	start of lseg
+ * @param b	end of lseg
+ * @param tol	tolerance values
+ * @param[out] dist	parametric distance from A to P' (in terms of A to B)
+ */
 BN_EXPORT extern int bn_isect_pt_lseg(fastf_t *dist,
 				      const point_t a,
 				      const point_t b,
 				      const point_t p,
 				      const struct bn_tol *tol);
+				      
+
 BN_EXPORT extern double bn_dist_pt_lseg(point_t pca,
 					const point_t a,
 					const point_t b,
 					const point_t p,
 					const struct bn_tol *tol);
+					
+/**
+ * B N _ R O T A T E _ B B O X
+ *@brief
+ * Transform a bounding box (RPP) by the given 4x4 matrix.  There are
+ * 8 corners to the bounding RPP.  Each one needs to be transformed
+ * and min/max'ed.  This is not minimal, but does fully contain any
+ * internal object, using an axis-aligned RPP.
+ */
 BN_EXPORT extern void bn_rotate_bbox(point_t omin,
 				     point_t omax,
 				     const mat_t mat,
 				     const point_t imin,
 				     const point_t imax);
+				     
+/**
+ * B N _ R O T A T E _ P L A N E
+ *@brief
+ * Transform a plane equation by the given 4x4 matrix.
+ */
 BN_EXPORT extern void bn_rotate_plane(plane_t oplane,
 				      const mat_t mat,
 				      const plane_t iplane);
+				      
+/**
+ * B N _ C O P L A N A R
+ *@brief
+ * Test if two planes are identical.  If so, their dot products will
+ * be either +1 or -1, with the distance from the origin equal in
+ * magnitude.
+ *
+ * @return -1	not coplanar, parallel but distinct
+ * @return 0	not coplanar, not parallel.  Planes intersect.
+ * @return 1	coplanar, same normal direction
+ * @return 2	coplanar, opposite normal direction
+ */
 BN_EXPORT extern int bn_coplanar(const plane_t a,
 				 const plane_t b,
 				 const struct bn_tol *tol);
+				 
+/**
+ * B N _ A N G L E _ M E A S U R E
+ *
+ * Using two perpendicular vectors (x_dir and y_dir) which lie in the
+ * same plane as 'vec', return the angle (in radians) of 'vec' from
+ * x_dir, going CCW around the perpendicular x_dir CROSS y_dir.
+ *
+ * Trig note -
+ *
+ * theta = atan2(x, y) returns an angle in the range -pi to +pi.
+ * Here, we need an angle in the range of 0 to 2pi.  This could be
+ * implemented by adding 2pi to theta when theta is negative, but this
+ * could have nasty numeric ambiguity right in the vicinity of theta =
+ * +pi, which is a very critical angle for the applications using this
+ * routine.
+ *
+ * So, an alternative formulation is to compute gamma = atan2(-x, -y),
+ * and then theta = gamma + pi.  Now, any error will occur in the
+ * vicinity of theta = 0, which can be handled much more readily.
+ *
+ * If theta is negative, or greater than two pi, wrap it around.
+ * These conditions only occur if there are problems in atan2().
+ *
+ * @return vec == x_dir returns 0,
+ * @return vec == y_dir returns pi/2,
+ * @return vec == -x_dir returns pi,
+ * @return vec == -y_dir returns 3*pi/2.
+ *
+ * In all cases, the returned value is between 0 and bn_twopi.
+ */
 BN_EXPORT extern double bn_angle_measure(vect_t vec,
 					 const vect_t x_dir,
 					 const vect_t y_dir);
+					 
+/**
+ * B N _ D I S T _ P T 3 _ A L O N G _ L I N E 3
+ *@brief
+ * Return the parametric distance t of a point X along a line defined
+ * as a ray, i.e. solve X = P + t * D.  If the point X does not lie on
+ * the line, then t is the distance of the perpendicular projection of
+ * point X onto the line.
+ */
 BN_EXPORT extern double bn_dist_pt3_along_line3(const point_t p,
 						const vect_t d,
 						const point_t x);
+						
+/**
+ * B N _ D I S T _ P T 2 _ A L O N G _ L I N E 2
+ *@brief
+ * Return the parametric distance t of a point X along a line defined
+ * as a ray, i.e. solve X = P + t * D.  If the point X does not lie on
+ * the line, then t is the distance of the perpendicular projection of
+ * point X onto the line.
+ */
 BN_EXPORT extern double bn_dist_pt2_along_line2(const point_t p,
 						const vect_t d,
 						const point_t x);
+						
+/**
+ *
+ * @return 1	if left <= mid <= right
+ * @return 0	if mid is not in the range.
+ */
 BN_EXPORT extern int bn_between(double left,
 				double mid,
 				double right,
 				const struct bn_tol *tol);
+				
+/**
+ * B N _ D O E S _ R A Y _ I S E C T _ T R I
+ *
+ * @return 0	No intersection
+ * @return 1	Intersection, 'inter' has intersect point.
+ */
 BN_EXPORT extern int bn_does_ray_isect_tri(const point_t pt,
 					   const vect_t dir,
 					   const point_t V,
 					   const point_t A,
 					   const point_t B,
 					   point_t inter);
+
+/**
+ * B N _ H L F _ C L A S S
+ *@brief
+ * Classify a halfspace, specified by its plane equation, against a
+ * bounding RPP.
+ *
+ * @return BN_CLASSIFY_INSIDE
+ * @return BN_CLASSIFY_OVERLAPPING
+ * @return BN_CLASSIFY_OUTSIDE
+ */
 BN_EXPORT extern int bn_hlf_class(const plane_t half_eqn,
 				  const vect_t min, const vect_t max,
 				  const struct bn_tol *tol);
@@ -994,18 +2498,60 @@
 #define BN_CLASSIFY_OVERLAPPING   0x0002
 #define BN_CLASSIFY_OUTSIDE       0x0003
 
+
+/**
+ * B N _ I S E C T _ P L A N E S
+ *@brief
+ * Calculates the point that is the minimum distance from all the
+ * planes in the "planes" array.  If the planes intersect at a single
+ * point, that point is the solution.
+ *
+ * The method used here is based on:
+
+ * An expression for the distance from a point to a plane is:
+ * VDOT(pt, plane)-plane[H].
+ * Square that distance and sum for all planes to get the "total"
+ * distance.
+ * For minimum total distance, the partial derivatives of this
+ * expression (with respect to x, y, and z) must all be zero.
+ * This produces a set of three equations in three unknowns (x, y, z).
+
+ * This routine sets up the three equations as [matrix][pt] = [hpq]
+ * and solves by inverting "matrix" into "inverse" and
+ * [pt] = [inverse][hpq].
+ *
+ * There is likely a more economical solution rather than matrix
+ * inversion, but bn_mat_inv was handy at the time.
+ *
+ * Checks if these planes form a singular matrix and returns.
+ *
+ * @return 0 - all is well
+ * @return 1 - planes form a singular matrix (no solution)
+ */
 BN_EXPORT extern int bn_isect_planes(point_t pt,
 				     const plane_t planes[],
 				     const size_t pl_count);
 
 /** @} */
+
+
 /*----------------------------------------------------------------------*/
+
+
 /* poly.c */
+
 /** @addtogroup poly */
-/** @{ */
+/** @{
+ *
+ * Library for dealing with polynomials.
+ *
+ */
 
+
 /* This could be larger, or even dynamic... */
 #define BN_MAX_POLY_DEGREE 4	/* Maximum Poly Order */
+
+
 /**
  * Polynomial data type
  * bn_poly->cf[n] corresponds to X^n
@@ -1015,28 +2561,114 @@
     size_t dgr;
     double cf[BN_MAX_POLY_DEGREE+1];
 }  bn_poly_t;
+
 #define BN_CK_POLY(_p) BU_CKMAG(_p, BN_POLY_MAGIC, "struct bn_poly")
 #define BN_POLY_NULL   ((struct bn_poly *)NULL)
 
+
+
+/**
+ * bn_poly_mul
+ *
+ * @brief 
+ * multiply two polynomials
+ */
 BN_EXPORT extern struct bn_poly *bn_poly_mul(struct bn_poly *product,
 					     const struct bn_poly *m1,
 					     const struct bn_poly *m2);
+					     
+/**
+ * bn_poly_scale
+ * 
+ * @brief
+ * scale a polynomial
+ */
 BN_EXPORT extern struct bn_poly *bn_poly_scale(struct bn_poly *eqn,
 					       double factor);
+					       
+/**
+ * bn_poly_add
+ * @brief
+ * add two polynomials
+ */
 BN_EXPORT extern struct bn_poly *bn_poly_add(struct bn_poly *sum,
 					     const struct bn_poly *poly1,
 					     const struct bn_poly *poly2);
+					     
+/**
+ * bn_poly_sub
+ * @brief
+ * subtract two polynomials
+ */
 BN_EXPORT extern struct bn_poly *bn_poly_sub(struct bn_poly *diff,
 					     const struct bn_poly *poly1,
 					     const struct bn_poly *poly2);
+					     
+/**
+ * s y n D i v
+ * @brief
+ * Divides any polynomial into any other polynomial using synthetic
+ * division.  Both polynomials must have real coefficients.
+ */
 BN_EXPORT extern void bn_poly_synthetic_division(struct bn_poly *quo,
 						 struct bn_poly *rem,
 						 const struct bn_poly *dvdend,
 						 const struct bn_poly *dvsor);
+						 
+/**
+ * b n _ p o l y _ q u a d r a t i c _ r o o t s
+ *@brief
+ * Uses the quadratic formula to find the roots (in `complex' form) of
+ * any quadratic equation with real coefficients.
+ *
+ *	@return 1 for success
+ *	@return 0 for fail.
+ */
 BN_EXPORT extern int bn_poly_quadratic_roots(struct bn_complex roots[],
 					     const struct bn_poly *quadrat);
+					     
+/**
+ * b n _ p o l y _ c u b i c _ r o o t s
+ *@brief
+ * Uses the cubic formula to find the roots (in `complex' form)
+ * of any cubic equation with real coefficients.
+ *
+ * to solve a polynomial of the form:
+ *
+ * X**3 + c1*X**2 + c2*X + c3 = 0,
+ *
+ * first reduce it to the form:
+ *
+ * Y**3 + a*Y + b = 0,
+ *
+ * where
+ * Y = X + c1/3,
+ * and
+ * a = c2 - c1**2/3,
+ * b = (2*c1**3 - 9*c1*c2 + 27*c3)/27.
+ *
+ * Then we define the value delta,   D = b**2/4 + a**3/27.
+ *
+ * If D > 0, there will be one real root and two conjugate
+ * complex roots.
+ * If D = 0, there will be three real roots at least two of
+ * which are equal.
+ * If D < 0, there will be three unequal real roots.
+ *
+ * Returns 1 for success, 0 for fail.
+ */
 BN_EXPORT extern int bn_poly_cubic_roots(struct bn_complex roots[],
 					 const struct bn_poly *eqn);
+					 
+/**
+ * b n _ p o l y _ q u a r t i c _ r o o t s
+ *@brief
+ * Uses the quartic formula to find the roots (in `complex' form)
+ * of any quartic equation with real coefficients.
+ *
+ *	@return 1 for success
+ *	@return 0 for fail.
+ */
 BN_EXPORT extern int bn_poly_quartic_roots(struct bn_complex roots[],
 					   const struct bn_poly *eqn);
 BN_EXPORT extern int bn_poly_findroot(bn_poly_t *eqn, 
@@ -1047,26 +2679,51 @@
 						  bn_complex_t *b,
 						  bn_complex_t *c,
 						  bn_complex_t *d);
+						  
 BN_EXPORT extern int bn_poly_checkroots(bn_poly_t *eqn,
 					bn_complex_t *roots,
 					int nroots);
+					
 BN_EXPORT extern void bn_poly_deflate(bn_poly_t *oldP,
 				      bn_complex_t *root);
+				      
 BN_EXPORT extern int bn_poly_roots(bn_poly_t *eqn,
 				   bn_complex_t roots[],
 				   const char *name);
+				   
+/**
+ * b n _ p r _ p o l y
+ *
+ * Print out the polynomial.
+ */
 BN_EXPORT extern void bn_pr_poly(const char *title,
 				 const struct bn_poly *eqn);
+				 
+/**
+ * b n _ p r _ r o o t s
+ *
+ * Print out the roots of a given polynomial (complex numbers)
+ */
 BN_EXPORT extern void bn_pr_roots(const char *title,
 				  const struct bn_complex roots[],
 				  int n);
 
 /** @} */
+
+
 /*----------------------------------------------------------------------*/
+
+
 /* multipoly.c */
+
 /** @addtogroup multipoly */
-/** @{ */
+/** @{
+ *
+ *	Library for dealing with bivariate polynomials.
+ *
+ */
 
+
 /**
  * Polynomial data type
  */
@@ -1078,47 +2735,178 @@
 }  bn_multipoly_t;
 
 /** @} */
+
+
 /*----------------------------------------------------------------------*/
+
+
 /* qmath.c */
-/** @addtogroup mat */
-/** @{ */
+
 /*
  * Quaternion support
  */
 
+/** @addtogroup mat */
+/** @{
+ *
+ * @brief
+ *  Quaternion math routines.
+ *
+ *  Unit Quaternions:
+ *	Q = [ r, a ]	where r = cos(theta/2) = rotation amount
+ *			    |a| = sin(theta/2) = rotation axis
+ *
+ *      If a = 0 we have the reals; if one coord is zero we have
+ *	 complex numbers (2D rotations).
+ *
+ *  [r, a][s, b] = [rs - a.b, rb + sa + axb]
+ *
+ *       -1
+ *  [r, a]   = (r - a) / (r^2 + a.a)
+ *
+ *  Powers of quaternions yield incremental rotations,
+ *   e.g. Q^3 is rotated three times as far as Q.
+ *
+ *  Some operations on quaternions:
+ *            -1
+ *   [0, P'] = Q  [0, P]Q		Rotate a point P by quaternion Q
+ *                     -1  a
+ *   slerp(Q, R, a) = Q(Q  R)	Spherical linear interp: 0 < a < 1
+ *
+ *   bisect(P, Q) = (P + Q) / |P + Q|	Great circle bisector
+ *
+ *  Additions inspired by "Quaternion Calculus For Animation" by Ken Shoemake,
+ *  SIGGRAPH '89 course notes for "Math for SIGGRAPH", May 1989.
+ *
+ */
+
+
+
+/**
+ *			Q U A T _ M A T 2 Q U A T
+ *@brief
+ *  Convert Matrix to Quaternion.
+ */
 BN_EXPORT extern void quat_mat2quat(quat_t quat,
 				    const mat_t mat);
+				    
+/**
+ *			Q U A T _ Q U A T 2 M A T
+ *@brief
+ *  Convert Quaternion to Matrix.
+ *
+ * NB: This only works for UNIT quaternions.  We may get imaginary results
+ *   otherwise.  We should normalize first (still yields same rotation).
+ */
 BN_EXPORT extern void quat_quat2mat(mat_t mat,
 				    const quat_t quat);
+				    
+/**
+ *			Q U A T _ D I S T A N C E
+ *@brief
+ * Gives the euclidean distance between two quaternions.
+ */
 BN_EXPORT extern double quat_distance(const quat_t q1,
 				      const quat_t q2);
+
+/**
+ *			Q U A T _ D O U B L E
+ *@brief
+ * Gives the quaternion point representing twice the rotation
+ *   from q1 to q2.
+ *   Needed for patching Bezier curves together.
+ *   A rather poor name admittedly.
+ */				      
 BN_EXPORT extern void quat_double(quat_t qout,
 				  const quat_t q1,
 				  const quat_t q2);
+				  
+/**
+ *			Q U A T _ B I S E C T
+ *@brief
+ * Gives the bisector of quaternions q1 and q2.
+ * (Could be done with quat_slerp and factor 0.5)
+ * [I believe they must be unit quaternions this to work]
+ */
 BN_EXPORT extern void quat_bisect(quat_t qout,
 				  const quat_t q1,
 				  const quat_t q2);
+				  
+/**
+ *			Q U A T _ S L E R P
+ *@brief
+ * Do Spherical Linear Interpolation between two unit quaternions
+ *  by the given factor.
+ *
+ * As f goes from 0 to 1, qout goes from q1 to q2.
+ * Code based on code by Ken Shoemake
+ */
 BN_EXPORT extern void quat_slerp(quat_t qout,
 				 const quat_t q1,
 				 const quat_t q2,
 				 double f);
+				 
+/**
+ *			Q U A T _ S B E R P
+ *@brief
+ * Spherical Bezier Interpolate between four quaternions by amount f.
+ * These are intended to be used as start and stop quaternions along
+ *   with two control quaternions chosen to match spline segments with
+ *   first order continuity.
+ *
+ *  Uses the method of successive bisection.
+ */
 BN_EXPORT extern void quat_sberp(quat_t qout,
 				 const quat_t q1,
 				 const quat_t qa,
 				 const quat_t qb,
 				 const quat_t q2,
 				 double f);
+				 
+/**
+ *			Q U A T _ M A K E _ N E A R E S T
+ *@brief
+ *  Set the quaternion q1 to the quaternion which yields the
+ *   smallest rotation from q2 (of the two versions of q1 which
+ *   produce the same orientation).
+ *
+ * Note that smallest euclidian distance implies smallest great
+ *   circle distance as well (since surface is convex).
+ */
 BN_EXPORT extern void quat_make_nearest(quat_t q1,
 					const quat_t q2);
+					
+/**
+ *			Q U A T _ P R I N T
+ */
 BN_EXPORT extern void quat_print(const char *title,
 				 const quat_t quat);
+				 
+/**
+ *			Q U A T _ E X P
+ *@brief
+ *  Exponentiate a quaternion, assuming that the scalar part is 0.
+ *  Code by Ken Shoemake.
+ */
 BN_EXPORT extern void quat_exp(quat_t out,
 			       const quat_t in);
+			       
+/**
+ *			Q U A T _ L O G
+ *@brief
+ *  Take the natural logarithm of a unit quaternion.
+ *  Code by Ken Shoemake.
+ */
 BN_EXPORT extern void quat_log(quat_t out,
 			       const quat_t in);
 /** @} */
+
+
 /*----------------------------------------------------------------------*/
+
+
 /* rand.c */
+
 /** @addtogroup rnd */
 /** @{ */
 /**
@@ -1143,9 +2931,20 @@
  * numbers from bn_rand_table do not change, because the procedural
  * noise would cease to exist.
  */
+
+
 #define BN_RAND_TABSIZE 4096
 #define BN_RAND_TABMASK 0xfff
 #define BN_RANDSEED(_i, _seed)  _i = ((unsigned)_seed) % BN_RAND_TABSIZE
+
+/**
+ * This is our table of random numbers.  Rather than calling drand48()
+ * or random() or rand() we just pick numbers out of this table.  This
+ * table has 4096 unique entries with floating point values ranging
+ * from the open interval (i.e. exclusive) 0.0 to 1.0 range.
+ *
+ * There are convenience macros for access in the bn.h header.
+ */
 BN_EXPORT extern const float bn_rand_table[BN_RAND_TABSIZE];
 
 /** BN_RANDOM always gives numbers between the open interval 0.0 to 1.0 */
@@ -1188,13 +2987,29 @@
 #define bn_tab_sin(_a)	(((_a) > 0) ? \
 			 (bn_sin_table[(int)((0.5+ (_a)*bn_sin_scale))&(BN_SINTABSIZE-1)]) :\
 			 (-bn_sin_table[(int)((0.5- (_a)*bn_sin_scale))&(BN_SINTABSIZE-1)]))
+			 
+/**
+ * table of floating point sine values in the closed (i.e. inclusive)
+ * interval -1.0 to 1.0 range.
+ */
 BN_EXPORT extern const float bn_sin_table[BN_SINTABSIZE];
 
+/**
+ *			M A T H T A B _ C O N S T A N T
+ *@brief
+ *  For benchmarking purposes, make the random number table predictable.
+ *  Setting to all zeros keeps dithered values at their original values.
+ */
 BN_EXPORT extern void bn_mathtab_constant();
 
 /** @} */
+
+
 /*----------------------------------------------------------------------*/
+
+
 /* randmt.c */
+
 /** @addtogroup rnd */
 /** @{ */
 /**
@@ -1221,8 +3036,133 @@
 
 
 /*----------------------------------------------------------------------*/
+
+
 /* wavelet.c */
 
+/** @addtogroup wavelet */
+/** @{
+ *
+ * @brief
+ *  This is a standard wavelet library that takes a given data buffer of some data
+ *  type and then performs a wavelet transform on that data.
+ *
+ * The transform
+ *  operations available are to either decompose or reconstruct a signal into its
+ *  corresponding wavelet form based on the haar wavelet.
+ *
+ *  Wavelet decompose/reconstruct operations
+ *
+ *	- bn_wlt_haar_1d_double_decompose(tbuffer, buffer, dimen, channels, limit)
+ *	- bn_wlt_haar_1d_float_decompose (tbuffer, buffer, dimen, channels, limit)
+ *	- bn_wlt_haar_1d_char_decompose  (tbuffer, buffer, dimen, channels, limit)
+ *	- bn_wlt_haar_1d_short_decompose (tbuffer, buffer, dimen, channels, limit)
+ *	- bn_wlt_haar_1d_int_decompose   (tbuffer, buffer, dimen, channels, limit)
+ *	- bn_wlt_haar_1d_long_decompose  (tbuffer, buffer, dimen, channels, limit)
+ *
+ *	- bn_wlt_haar_1d_double_reconstruct(tbuffer, buffer, dimen, channels, sub_sz, limit)
+ *	- bn_wlt_haar_1d_float_reconstruct (tbuffer, buffer, dimen, channels, sub_sz, limit)
+ *	- bn_wlt_haar_1d_char_reconstruct  (tbuffer, buffer, dimen, channels, sub_sz, limit)
+ *	- bn_wlt_haar_1d_short_reconstruct (tbuffer, buffer, dimen, channels, sub_sz, limit)
+ *	- bn_wlt_haar_1d_int_reconstruct   (tbuffer, buffer, dimen, channels, sub_sz, limit)
+ *	- bn_wlt_haar_1d_long_reconstruct  (tbuffer, buffer, dimen, channels, sub_sz, limit)
+ *
+ *	- bn_wlt_haar_2d_double_decompose(tbuffer, buffer, dimen, channels, limit)
+ *	- bn_wlt_haar_2d_float_decompose (tbuffer, buffer, dimen, channels, limit)
+ *	- bn_wlt_haar_2d_char_decompose  (tbuffer, buffer, dimen, channels, limit)
+ *	- bn_wlt_haar_2d_short_decompose (tbuffer, buffer, dimen, channels, limit)
+ *	- bn_wlt_haar_2d_int_decompose   (tbuffer, buffer, dimen, channels, limit)
+ *	- bn_wlt_haar_2d_long_decompose  (tbuffer, buffer, dimen, channels, limit)
+ *
+ *	- bn_wlt_haar_2d_double_reconstruct(tbuffer, buffer, dimen, channels, sub_sz, limit)
+ *	- bn_wlt_haar_2d_float_reconstruct (tbuffer, buffer, dimen, channels, sub_sz, limit)
+ *	- bn_wlt_haar_2d_char_reconstruct  (tbuffer, buffer, dimen, channels, sub_sz, limit)
+ *	- bn_wlt_haar_2d_short_reconstruct (tbuffer, buffer, dimen, channels, sub_sz, limit)
+ *	- bn_wlt_haar_2d_int_reconstruct   (tbuffer, buffer, dimen, channels, sub_sz, limit)
+ *	- bn_wlt_haar_2d_long_reconstruct  (tbuffer, buffer, dimen, channels, sub_sz, limit)
+ *
+ *	- bn_wlt_haar_2d_double_decompose2(tbuffer, buffer, width, height, channels, limit)
+ *	- bn_wlt_haar_2d_float_decompose2 (tbuffer, buffer, width, height, channels, limit)
+ *	- bn_wlt_haar_2d_char_decompose2  (tbuffer, buffer, width, height, channels, limit)
+ *	- bn_wlt_haar_2d_short_decompose2 (tbuffer, buffer, width, height, channels, limit)
+ *	- bn_wlt_haar_2d_int_decompose2   (tbuffer, buffer, width, height, channels, limit)
+ *	- bn_wlt_haar_2d_long_decompose2  (tbuffer, buffer, width, height, channels, limit)
+ *
+ *	- bn_wlt_haar_2d_double_reconstruct2(tbuffer, buffer, width, height, channels, sub_sz, limit)
+ *	- bn_wlt_haar_2d_float_reconstruct2 (tbuffer, buffer, width, height, channels, sub_sz, limit)
+ *	- bn_wlt_haar_2d_char_reconstruct2  (tbuffer, buffer, width, height, channels, sub_sz, limit)
+ *	- bn_wlt_haar_2d_short_reconstruct2 (tbuffer, buffer, width, height, channels, sub_sz, limit)
+ *	- bn_wlt_haar_2d_int_reconstruct2   (tbuffer, buffer, width, height, channels, sub_sz, limit)
+ *	- bn_wlt_haar_2d_long_reconstruct2  (tbuffer, buffer, width, height, channels, sub_sz, limit)
+ *
+ *
+ *  For greatest accuracy, it is preferable to convert everything to "double"
+ *  and decompose/reconstruct with that.  However, there are useful
+ *  properties to performing the decomposition and/or reconstruction in
+ *  various data types (most notably char).
+ *
+ *  Rather than define all of these routines explicitly, we define
+ *  2 macros "decompose" and "reconstruct" which embody the structure of
+ *  the function (which is common to all of them).  We then instatiate
+ *  these macros once for each of the data types.  It's ugly, but it
+ *  assures that a change to the structure of one operation type
+ *  (decompose or reconstruct) occurs for all data types.
+ *
+ *
+ *
+ *
+ *  bn_wlt_haar_1d_*_decompose(tbuffer, buffer, dimen, channels, limit)
+ *  @par Parameters:
+ *	- tbuffer     a temporary data buffer 1/2 as large as "buffer". See (1) below.
+ *	- buffer      pointer to the data to be decomposed
+ *	- dimen    the number of samples in the data buffer
+ *	- channels the number of values per sample
+ *	- limit    the extent of the decomposition
+ *
+ *  Perform a Haar wavelet decomposition on the data in buffer "buffer".  The
+ *  decomposition is done "in place" on the data, hence the values in "buffer"
+ *  are not preserved, but rather replaced by their decomposition.
+ *  The number of original samples in the buffer (parameter "dimen") and the
+ *  decomposition limit ("limit") must both be a power of 2 (e.g. 512, 1024).
+ *  The buffer is decomposed into "average" and "detail" halves until the
+ *  size of the "average" portion reaches "limit".  Simultaneous
+ *  decomposition of multi-plane (e.g. pixel) data, can be performed by
+ *  indicating the number of planes in the "channels" parameter.
+ *
+ *  (1) The process requires a temporary buffer which is 1/2 the size of the
+ *  longest span to be decomposed.  If the "tbuffer" argument is non-null then
+ *  it is a pointer to a temporary buffer.  If the pointer is NULL, then a
+ *  local temporary buffer will be allocated (and freed).
+ *
+ *  @par Examples:
+ @code
+ double dbuffer[512], cbuffer[256];
+ ...
+ bn_wlt_haar_1d_double_decompose(cbuffer, dbuffer, 512, 1, 1);
+ @endcode
+ *
+ *    performs complete decomposition on the data in array "dbuffer".
+ *
+ @code
+ double buffer[3][512];	 /_* 512 samples, 3 values/sample (e.g. RGB?)*_/
+ double tbuffer[3][256];	 /_* the temporary buffer *_/
+ ...
+ bn_wlt_haar_1d_double_decompose(tbuffer, buffer, 512, 3, 1);
+ @endcode
+ *
+ *    This will completely decompose the data in buffer.  The first sample will
+ *    be the average of all the samples.  Alternatively:
+ *
+ * bn_wlt_haar_1d_double_decompose(tbuffer, buffer, 512, 3, 64);
+ *
+ *    decomposes buffer into a 64-sample "average image" and 3 "detail" sets.
+ *
+ *  bn_wlt_haar_1d_*_reconstruct(tbuffer, buffer, dimen, channels, sub_sz, limit)
+ *
+ */
+/** @} */
+
+
 #define CK_POW_2(dimen) { \
 	register unsigned long j; \
 	register int ok; \
@@ -1462,7 +3402,10 @@
 
 
 /*----------------------------------------------------------------------*/
+
+
 /* const.c */
+
 BN_EXPORT extern const double bn_pi;
 BN_EXPORT extern const double bn_twopi;
 BN_EXPORT extern const double bn_halfpi;
@@ -1473,8 +3416,32 @@
 BN_EXPORT extern const double bn_radtodeg;
 
 /*----------------------------------------------------------------------*/
+
+
 /* tabdata.c */
 
+/** @addtogroup anim */
+/** @{ */
+
+/** @brief
+ *  Routines for processing tables (curves) of data with one independent
+ *  parameter which is common to many sets of dependent data values.
+ *
+ *  Operates on bn_table (independent var) and
+ *  bn_tabdata (dependent variable) structures.
+ *
+ *  One application is for storing spectral curves, see spectrum.c
+ *
+ *  @par Inspired by -
+ *	Roy Hall and his book "Illumination and Color in Computer
+ *@n	Generated Imagery", Springer Verlag, New York, 1989.
+ *@n	ISBN 0-387-96774-5
+ *
+ *  With thanks to Russ Moulton Jr, EOSoft Inc. for his "rad.c" module.
+ */
+/** @} */
+
+
 /**
  * T A B D A T A
  *
@@ -1513,6 +3480,7 @@
  * frequencies.
  *
  */
+
 struct bn_table {
     uint32_t magic;
     size_t nx;
@@ -1565,51 +3533,140 @@
  * Routines
  */
 
+/*
+ *			B N _ T A B L E _ F R E E
+ */
 BN_EXPORT extern void bn_table_free(struct bn_table *tabp);
+
+/*
+ *			B N _ T A B D A T A _ F R E E
+ */
 BN_EXPORT extern void bn_tabdata_free(struct bn_tabdata *data);
+
+/*
+ *			B N _ C K _ T A B L E
+ */
 BN_EXPORT extern void bn_ck_table(const struct bn_table *tabp);
+
+/*
+ *			B N _ T A B L E _ M A K E _ U N I F O R M
+ *@brief
+ *  Set up an independent "table margin" from 'first' to 'last',
+ *  inclusive, using 'num' uniformly spaced samples.  Num >= 1.
+ */
 BN_EXPORT extern struct bn_table *bn_table_make_uniform(size_t num,
 							double first,
 							double last);
+							
+/*
+ *			B N _ T A B D A T A _ A D D
+ *@brief
+ *  Sum the values from two data tables.
+ */
 BN_EXPORT extern void bn_tabdata_add(struct bn_tabdata *out,
 				     const struct bn_tabdata *in1,
 				     const struct bn_tabdata *in2);
+				     
+/*
+ *			B N _ T A B D A T A _ M U L
+ *@brief
+ *  Element-by-element multiply the values from two data tables.
+ */
 BN_EXPORT extern void bn_tabdata_mul(struct bn_tabdata *out,
 				     const struct bn_tabdata *in1,
 				     const struct bn_tabdata *in2);
+				     
+/*
+ *			B N _ T A B D A T A _ M U L 3
+ *@brief
+ *  Element-by-element multiply the values from three data tables.
+ */
 BN_EXPORT extern void bn_tabdata_mul3(struct bn_tabdata *out,
 				      const struct bn_tabdata *in1,
 				      const struct bn_tabdata *in2,
 				      const struct bn_tabdata *in3);
+				      
+/*
+ *			B N _ T A B D A T A _ I N C R _ M U L 3 _ S C A L E
+ *@brief
+ *  Element-by-element multiply the values from three data tables and a scalor.
+ *
+ *	out += in1 * in2 * in3 * scale
+ */
 BN_EXPORT extern void bn_tabdata_incr_mul3_scale(struct bn_tabdata *out,
 						 const struct bn_tabdata *in1,
 						 const struct bn_tabdata *in2,
 						 const struct bn_tabdata *in3,
 						 double scale);
+						 
+/*
+ *			B N _ T A B D A T A _ I N C R _ M U L 2 _ S C A L E
+ *@brief
+ *  Element-by-element multiply the values from two data tables and a scalor.
+ *
+ *	out += in1 * in2 * scale
+ */
 BN_EXPORT extern void bn_tabdata_incr_mul2_scale(struct bn_tabdata *out,
 						 const struct bn_tabdata *in1,
 						 const struct bn_tabdata *in2,
 						 double scale);
+						 
+/*
+ *			B N _ T A B D A T A _ S C A L E
+ *@brief
+ *  Multiply every element in a data table by a scalar value 'scale'.
+ */
 BN_EXPORT extern void bn_tabdata_scale(struct bn_tabdata *out,
 				       const struct bn_tabdata *in1,
 				       double scale);
+				       
+/*
+ *			B N _ T A B L E _ S C A L E
+ *@brief
+ *  Scale the indepentent axis of a table by 'scale'.
+ */
 BN_EXPORT extern void bn_table_scale(struct bn_table *tabp,
 				     double scale);
+				     
+/*
+ *			B N _ T A B D A T A _ J O I N 1
+ *@brief
+ *  Multiply every element in data table in2 by a scalar value 'scale',
+ *  add it to the element in in1, and store in 'out'.
+ *  'out' may overlap in1 or in2.
+ */
 BN_EXPORT extern void bn_tabdata_join1(struct bn_tabdata *out,
 				       const struct bn_tabdata *in1,
 				       double scale,
 				       const struct bn_tabdata *in2);
+				       
+/*
+ *			B N _ T A B D A T A _ J O I N 2
+ *@brief
+ *  Multiply every element in data table in2 by a scalar value 'scale2',
+ *  plus in3 * scale3, and
+ *  add it to the element in in1, and store in 'out'.
+ *  'out' may overlap in1 or in2.
+ */
 BN_EXPORT extern void bn_tabdata_join2(struct bn_tabdata *out,
 				       const struct bn_tabdata *in1,
 				       double scale2,
 				       const struct bn_tabdata *in2,
 				       double scale3,
 				       const struct bn_tabdata *in3);
+				       
+/*
+ *			B N _ T A B D A T A _ B L E N D 2
+ */
 BN_EXPORT extern void bn_tabdata_blend2(struct bn_tabdata *out,
 					double scale1,
 					const struct bn_tabdata *in1,
 					double scale2,
 					const struct bn_tabdata *in2);
+					
+/*
+ *			B N _ T A B D A T A _ B L E N D 3
+ */
 BN_EXPORT extern void bn_tabdata_blend3(struct bn_tabdata *out,
 					double scale1,
 					const struct bn_tabdata *in1,
@@ -1617,61 +3674,275 @@
 					const struct bn_tabdata *in2,
 					double scale3,
 					const struct bn_tabdata *in3);
+					
+/*
+ *			B N _ T A B D A T A _ A R E A 1
+ *@brief
+ *  Following interpretation #1, where y[j] stores the total (integral
+ *  or area) value within the interval, return the area under the whole curve.
+ *  This is simply totaling up the areas from each of the intervals.
+ */
 BN_EXPORT extern double bn_tabdata_area1(const struct bn_tabdata *in);
+
+/*
+ *			B N _ T A B D A T A _ A R E A 2
+ *@brief
+ *  Following interpretation #2, where y[j] stores the average
+ *  value for the interval, return the area under
+ *  the whole curve.  Since the iterval spacing need not be uniform,
+ *  sum the areas of the rectangles.
+ */
 BN_EXPORT extern double bn_tabdata_area2(const struct bn_tabdata *in);
+
+/*
+ *			B N _ T A B D A T A _ M U L _ A R E A 1
+ *@brief
+ *  Following interpretation #1, where y[j] stores the total (integral
+ *  or area) value within the interval, return the area under the whole curve.
+ *  This is simply totaling up the areas from each of the intervals.
+ *  The curve value is found by multiplying corresponding entries from
+ *  in1 and in2.
+ */
 BN_EXPORT extern double bn_tabdata_mul_area1(const struct bn_tabdata *in1,
 					     const struct bn_tabdata *in2);
+					     
+/*
+ *			B N _ T A B D A T A _ M U L _ A R E A 2
+ *@brief
+ *  Following interpretation #2,
+ *  return the area under the whole curve.
+ *  The curve value is found by multiplying corresponding entries from
+ *  in1 and in2.
+ */
 BN_EXPORT extern double bn_tabdata_mul_area2(const struct bn_tabdata *in1,
 					     const struct bn_tabdata *in2);
+					     
+/*
+ *			B N _ T A B L E _ L I N _ I N T E R P
+ *@brief
+ *  Return the value of the curve at independent parameter value 'wl'.
+ *  Linearly interpolate between values in the input table.
+ *  Zero is returned for values outside the sampled range.
+ */
 BN_EXPORT extern fastf_t bn_table_lin_interp(const struct bn_tabdata *samp,
 					     double wl);
+					     
+/*
+ *			B N _ T A B D A T A _ R E S A M P L E _ M A X
+ *@brief
+ *  Given a set of sampled data 'olddata', resample it for different
+ *  spacing, by linearly interpolating the values when an output span
+ *  is entirely contained within an input span, and by taking the
+ *  maximum when an output span covers more than one input span.
+ *
+ *  This assumes interpretation (2) of the data, i.e. that the values
+ *  are the average value across the interval.
+ */
 BN_EXPORT extern struct bn_tabdata *bn_tabdata_resample_max(const struct bn_table *newtable,
 							    const struct bn_tabdata *olddata);
+							    
+/*
+ *			B N _ T A B D A T A _ R E S A M P L E _ A V G
+ *@brief
+ *  Given a set of sampled data 'olddata', resample it for different
+ *  spacing, by linearly interpolating the values when an output span
+ *  is entirely contained within an input span, and by taking the
+ *  average when an output span covers more than one input span.
+ *
+ *  This assumes interpretation (2) of the data, i.e. that the values
+ *  are the average value across the interval.
+ */
 BN_EXPORT extern struct bn_tabdata *bn_tabdata_resample_avg(const struct bn_table *newtable,
 							    const struct bn_tabdata *olddata);
+							    
+/*
+ *			B N _ T A B L E _ W R I T E
+ *@brief
+ *  Write out the table structure in an ASCII file,
+ *  giving the number of values (minus 1), and the
+ *  actual values.
+ */
 BN_EXPORT extern int bn_table_write(const char *filename,
 				    const struct bn_table *tabp);
+				    
+/*
+ *			B N _ T A B L E _ R E A D
+ *@brief
+ *  Allocate and read in the independent variable values from an ASCII file,
+ *  giving the number of samples (minus 1), and the
+ *  actual values.
+ */
 BN_EXPORT extern struct bn_table *bn_table_read(const char *filename);
+
+/*
+ *			B N _ P R _ T A B L E
+ */
 BN_EXPORT extern void bn_pr_table(const char *title,
 				  const struct bn_table *tabp);
+				  
+/*
+ *			B N _ P R _ T A B D A T A
+ */
 BN_EXPORT extern void bn_pr_tabdata(const char *title,
 				    const struct bn_tabdata *data);
+				    
+/*
+ *			B N _ P R I N T _ T A B L E _ A N D _ T A B D A T A
+ *@brief
+ *  Write out a given data table into an ASCII file,
+ *  suitable for input to GNUPLOT.
+ *
+ *	(set term postscript)
+ *	(set output "|print-postscript")
+ *	(plot "filename" with lines)
+ */
 BN_EXPORT extern int bn_print_table_and_tabdata(const char *filename,
 						const struct bn_tabdata *data);
+						
+/*
+ *			B N _ R E A D _ T A B L E _ A N D _ T A B D A T A
+ *@brief
+ *  Read in a file which contains two columns of numbers, the first
+ *  column being the wavelength, the second column being the sample value
+ *  at that wavelength.
+ *
+ *  A new bn_table structure and one bn_tabdata structure
+ *  are created, a pointer to the bn_tabdata structure is returned.
+ *  The final wavelength is guessed at.
+ */
 BN_EXPORT extern struct bn_tabdata *bn_read_table_and_tabdata(const char *filename);
+
+/*
+ *			B N _ T A B D A T A _ B I N A R Y _ R E A D
+ */
 BN_EXPORT extern struct bn_tabdata *bn_tabdata_binary_read(const char *filename,
 							   size_t num,
 							   const struct bn_table *tabp);
+							   
+/*
+ *			B N _ T A B D A T A _ M A L L O C _ A R R A Y
+ *@brief
+ *  Allocate storage for, and initialize, an array of 'num' data table
+ *  structures.
+ *  This subroutine is provided because the bn_tabdata structures
+ *  are variable length.
+ */
 BN_EXPORT extern struct bn_tabdata *bn_tabdata_malloc_array(const struct bn_table *tabp,
 							    size_t num);
+							    
+/*
+ *			B N _ T A B D A T A _ C O P Y
+ */
 BN_EXPORT extern void bn_tabdata_copy(struct bn_tabdata *out,
 				      const struct bn_tabdata *in);
+				      
+/*
+ *			B N _ T A B D A T A _ D U P
+ */
 BN_EXPORT extern struct bn_tabdata *bn_tabdata_dup(const struct bn_tabdata *in);
+
+/*
+ *			B N _ T A B D A T A _ G E T _ C O N S T V A L
+ *@brief
+ *  For a given table, allocate and return a tabdata structure
+ *  with all elements initialized to 'val'.
+ */
 BN_EXPORT extern struct bn_tabdata *bn_tabdata_get_constval(double val,
 							    const struct bn_table *tabp);
+							    
+/*
+ *			B N _ T A B D A T A _ C O N S T V A L
+ *@brief
+ *  Set all the tabdata elements to 'val'
+ */
 BN_EXPORT extern void bn_tabdata_constval(struct bn_tabdata *data,
 					  double val);
+					  
+/*
+ *			B N _ T A B D A T A _ T O _ T C L
+ *@brief
+ *  Convert an bn_tabdata/bn_table pair into a Tcl compatible string
+ *  appended to a VLS.  It will have form:
+ *	x {...} y {...} nx # ymin # ymax #
+ */
 BN_EXPORT extern void bn_tabdata_to_tcl(struct bu_vls *vp,
 					const struct bn_tabdata *data);
+					
+/*
+ *			B N _ T A B D A T A _ F R O M _ A R R A Y
+ *@brief
+ *  Given an array of (x, y) pairs, build the relevant bn_table and
+ *  bn_tabdata structures.
+ *
+ *  The table is terminated by an x value <= 0.
+ *  Consistent with the interpretation of the spans,
+ *  invent a final span ending x value.
+ */
 BN_EXPORT extern struct bn_tabdata *bn_tabdata_from_array(const double *array);
+
+/*
+ *			B N _ T A B D A T A _ F R E Q _ S H I F T
+ *@brief
+ *  Shift the data by a constant offset in the independent variable
+ *  (often frequency), interpolating new sample values.
+ */
 BN_EXPORT extern void bn_tabdata_freq_shift(struct bn_tabdata *out,
 					    const struct bn_tabdata *in,
 					    double offset);
+					    
+/*
+ *			B N _ T A B L E _ I N T E R V A L _ N U M _ S A M P L E S
+ *@brief
+ *  Returns number of sample points between 'low' and 'hi', inclusive.
+ */
 BN_EXPORT extern int bn_table_interval_num_samples(const struct bn_table *tabp,
 						   double low,
 						   double hi);
+						   
+/*
+ *			B N _ T A B L E _ D E L E T E _ S A M P L E _ P T S
+ *@brief
+ *  Remove all sampling points between subscripts i and j, inclusive.
+ *  Don't bother freeing the tiny bit of storage at the end of the array.
+ *  Returns number of points removed.
+ */
 BN_EXPORT extern int bn_table_delete_sample_pts(struct bn_table *tabp,
 						unsigned int i,
 						unsigned int j);
+						
+/*
+ *			B N _ T A B L E _ M E R G E 2
+ *@brief
+ *  A new table is returned which has sample points at places from
+ *  each of the input tables.
+ */
 BN_EXPORT extern struct bn_table *bn_table_merge2(const struct bn_table *a,
 						  const struct bn_table *b);
+						  
+/*
+ *		B N _ T A B D A T A _ M K _ L I N E A R _ F I L T E R
+ *@brief
+ *  Create a filter to accept power in a given band.
+ *  The first and last filter values will be in the range 0..1,
+ *  while all the internal filter values will be 1.0,
+ *  and all samples outside the given band will be 0.0.
+ *
+ *
+ *  @return	NULL	if given band does not overlap input spectrum
+ *  @return	tabdata*
+ */
 BN_EXPORT extern struct bn_tabdata *bn_tabdata_mk_linear_filter(const struct bn_table *spectrum,
 								double lower_wavelen,
 								double upper_wavelen);
 
+
 /*----------------------------------------------------------------------*/
+
+
 /* vlist.c */
+
 #define BN_VLIST_CHUNK 35		/**< @brief 32-bit mach => just less than 1k */
+
 /**
  * B N _ V L I S T
  *
@@ -1792,10 +4063,25 @@
 
 
 /*----------------------------------------------------------------------*/
+
+
 /* vert_tree.c */
 /*
  * vertex tree support
  */
+
+/** @addtogroup vtree */
+/** @{ */
+/** @brief
+ * Routines to manage a binary search tree of vertices.
+ *
+ * The actual vertices are stored in an array
+ * for convenient use by routines such as "mk_bot".
+ * The binary search tree stores indices into the array.
+ *
+ */
+
+
 /**
  * packaging structure
  * holds all the required info for a single vertex tree
@@ -1816,14 +4102,47 @@
 
 #define BN_CK_VERT_TREE(_p) BU_CKMAG(_p, VERT_TREE_MAGIC, "vert_tree")
 
+
+/**		C R E A T E _ V E R T _ T R E E
+ *@brief
+ *	routine to create a vertex tree.
+ *
+ *	Possible refinements include specifying an initial size
+ */
 BN_EXPORT extern struct vert_root *create_vert_tree();
+
+/**		C R E A T E _ V E R T _ T R E E _ W _ N O R M S
+ *@brief
+ *	routine to create a vertex tree.
+ *
+ *	Possible refinements include specifying an initial size
+ */
 BN_EXPORT extern struct vert_root *create_vert_tree_w_norms();
+
+/**		F R E E _ V E R T_ T R E E
+ *@brief
+ *	Routine to free a vertex tree and all associated dynamic memory
+ */
 BN_EXPORT extern void free_vert_tree(struct vert_root *tree_root);
+
+/**		A D D _ V E R T
+ *@brief
+ *	Routine to add a vertex to the current list of part vertices.
+ *	The array is re-alloc'd if needed.
+ *	Returns index into the array of vertices where this vertex is stored
+ */
 BN_EXPORT extern int Add_vert(double x,
 			      double y,
 			      double z,
 			      struct vert_root *tree_root,
 			      fastf_t local_tol_sq);
+			      
+/**		A D D _ V E R T _ A N D _ N O R M
+ *@brief
+ *	Routine to add a vertex and a normal to the current list of part vertices.
+ *	The array is re-alloc'd if needed.
+ *	Returns index into the array of vertices where this vertex and normal is stored
+ */
 BN_EXPORT extern int Add_vert_and_norm(double x,
 				       double y,
 				       double z,
@@ -1832,10 +4151,25 @@
 				       double nz,
 				       struct vert_root *tree_root,
 				       fastf_t local_tol_sq);
+				       
+/**		C L E A N _ V E R T _ T R E E
+ *@brief
+ *	Routine to free the binary search tree and reset the current number of vertices.
+ *	The vertex array is left untouched, for re-use later.
+ */
 BN_EXPORT extern void clean_vert_tree(struct vert_root *tree_root);
 
+
 /*----------------------------------------------------------------------*/
+
+
 /* vectfont.c */
+
+/**
+ * @brief
+ *  Once-only setup routine
+ *  Used by libplot3/symbol.c, so it can't be static.
+ */
 BN_EXPORT extern void tp_setup();
 
 /**
@@ -1845,10 +4179,14 @@
 
 __END_DECLS
 
+
 /*----------------------------------------------------------------------*/
+
+
 /*
  * plane structures (from src/librt/plane.h)
  */
+
 /**
  * Plane structures
  * holds all the required info for geometric planes.
Index: src/libbn/tabdata.c
===================================================================
--- src/libbn/tabdata.c	(revision 49726)
+++ src/libbn/tabdata.c	(working copy)
@@ -17,6 +17,7 @@
  * License along with this file; see the file named COPYING for more
  * information.
  */
+
 /** @addtogroup anim */
 /** @{ */
 /** @file libbn/tabdata.c
@@ -49,9 +50,7 @@
 #include "bu.h"
 #include "bn.h"
 
-/*
- *			B N _ T A B L E _ F R E E
- */
+
 void
 bn_table_free(struct bn_table *tabp)
 {
@@ -62,9 +61,7 @@
     bu_free(tabp, "struct bn_table");
 }
 
-/*
- *			B N _ T A B D A T A _ F R E E
- */
+
 void
 bn_tabdata_free(struct bn_tabdata *data)
 {
@@ -78,9 +75,7 @@
     bu_free(data, "struct bn_tabdata");
 }
 
-/*
- *			B N _ C K _ T A B L E
- */
+
 void
 bn_ck_table(const struct bn_table *tabp)
 {
@@ -99,12 +94,7 @@
     }
 }
 
-/*
- *			B N _ T A B L E _ M A K E _ U N I F O R M
- *@brief
- *  Set up an independent "table margin" from 'first' to 'last',
- *  inclusive, using 'num' uniformly spaced samples.  Num >= 1.
- */
+
 struct bn_table *
 bn_table_make_uniform(size_t num, double first, double last) {
     struct bn_table	*tabp;
@@ -130,11 +120,6 @@
     return tabp;
 }
 
-/*
- *			B N _ T A B D A T A _ A D D
- *@brief
- *  Sum the values from two data tables.
- */
 void
 bn_tabdata_add(struct bn_tabdata *out, const struct bn_tabdata *in1, const struct bn_tabdata *in2)
 {
@@ -161,11 +146,6 @@
     /* VADD2N( out->y, i1->y, i2->y, in1->ny ); */
 }
 
-/*
- *			B N _ T A B D A T A _ M U L
- *@brief
- *  Element-by-element multiply the values from two data tables.
- */
 void
 bn_tabdata_mul(struct bn_tabdata *out, const struct bn_tabdata *in1, const struct bn_tabdata *in2)
 {
@@ -192,11 +172,6 @@
     /* VELMUL2N( out->y, i1->y, i2->y, in1->ny ); */
 }
 
-/*
- *			B N _ T A B D A T A _ M U L 3
- *@brief
- *  Element-by-element multiply the values from three data tables.
- */
 void
 bn_tabdata_mul3(struct bn_tabdata *out, const struct bn_tabdata *in1, const struct bn_tabdata *in2, const struct bn_tabdata *in3)
 {
@@ -225,13 +200,6 @@
     /* VELMUL3N( out->y, i1->y, i2->y, i3->y, in1->ny ); */
 }
 
-/*
- *			B N _ T A B D A T A _ I N C R _ M U L 3 _ S C A L E
- *@brief
- *  Element-by-element multiply the values from three data tables and a scalor.
- *
- *	out += in1 * in2 * in3 * scale
- */
 void
 bn_tabdata_incr_mul3_scale(struct bn_tabdata *out, const struct bn_tabdata *in1, const struct bn_tabdata *in2, const struct bn_tabdata *in3, register double scale)
 {
@@ -259,13 +227,6 @@
 	*op++ += *i1++ * *i2++ * *i3++ * scale;
 }
 
-/*
- *			B N _ T A B D A T A _ I N C R _ M U L 2 _ S C A L E
- *@brief
- *  Element-by-element multiply the values from two data tables and a scalor.
- *
- *	out += in1 * in2 * scale
- */
 void
 bn_tabdata_incr_mul2_scale(struct bn_tabdata *out, const struct bn_tabdata *in1, const struct bn_tabdata *in2, register double scale)
 {
@@ -291,11 +252,6 @@
 	*op++ += *i1++ * *i2++ * scale;
 }
 
-/*
- *			B N _ T A B D A T A _ S C A L E
- *@brief
- *  Multiply every element in a data table by a scalar value 'scale'.
- */
 void
 bn_tabdata_scale(struct bn_tabdata *out, const struct bn_tabdata *in1, register double scale)
 {
@@ -320,11 +276,6 @@
     /* VSCALEN( out->y, in->y, scale ); */
 }
 
-/*
- *			B N _ T A B L E _ S C A L E
- *@brief
- *  Scale the indepentent axis of a table by 'scale'.
- */
 void
 bn_table_scale(struct bn_table *tabp, register double scale)
 {
@@ -341,13 +292,6 @@
     /* VSCALEN( tabp->x, tabp->x, scale, tabp->nx+1 ); */
 }
 
-/*
- *			B N _ T A B D A T A _ J O I N 1
- *@brief
- *  Multiply every element in data table in2 by a scalar value 'scale',
- *  add it to the element in in1, and store in 'out'.
- *  'out' may overlap in1 or in2.
- */
 void
 bn_tabdata_join1(struct bn_tabdata *out, const struct bn_tabdata *in1, register double scale, const struct bn_tabdata *in2)
 {
@@ -376,14 +320,6 @@
     /* VJOIN1N( out->y, in1->y, scale, in2->y ); */
 }
 
-/*
- *			B N _ T A B D A T A _ J O I N 2
- *@brief
- *  Multiply every element in data table in2 by a scalar value 'scale2',
- *  plus in3 * scale3, and
- *  add it to the element in in1, and store in 'out'.
- *  'out' may overlap in1 or in2.
- */
 void
 bn_tabdata_join2(struct bn_tabdata *out, const struct bn_tabdata *in1, register double scale2, const struct bn_tabdata *in2, register double scale3, const struct bn_tabdata *in3)
 {
@@ -415,9 +351,6 @@
     /* VJOIN2N( out->y, in1->y, scale2, in2->y, scale3, in3->y ); */
 }
 
-/*
- *			B N _ T A B D A T A _ B L E N D 2
- */
 void
 bn_tabdata_blend2(struct bn_tabdata *out, register double scale1, const struct bn_tabdata *in1, register double scale2, const struct bn_tabdata *in2)
 {
@@ -446,9 +379,6 @@
     /* VBLEND2N( out->y, scale1, in1->y, scale2, in2->y ); */
 }
 
-/*
- *			B N _ T A B D A T A _ B L E N D 3
- */
 void
 bn_tabdata_blend3(struct bn_tabdata *out, register double scale1, const struct bn_tabdata *in1, register double scale2, const struct bn_tabdata *in2, register double scale3, const struct bn_tabdata *in3)
 {
@@ -481,13 +411,6 @@
     /* VBLEND3N( out->y, scale1, in1->y, scale2, in2->y, scale3, in3->y ); */
 }
 
-/*
- *			B N _ T A B D A T A _ A R E A 1
- *@brief
- *  Following interpretation #1, where y[j] stores the total (integral
- *  or area) value within the interval, return the area under the whole curve.
- *  This is simply totaling up the areas from each of the intervals.
- */
 double
 bn_tabdata_area1(const struct bn_tabdata *in)
 {
@@ -507,14 +430,6 @@
     return area;
 }
 
-/*
- *			B N _ T A B D A T A _ A R E A 2
- *@brief
- *  Following interpretation #2, where y[j] stores the average
- *  value for the interval, return the area under
- *  the whole curve.  Since the iterval spacing need not be uniform,
- *  sum the areas of the rectangles.
- */
 double
 bn_tabdata_area2(const struct bn_tabdata *in)
 {
@@ -537,15 +452,6 @@
     return area;
 }
 
-/*
- *			B N _ T A B D A T A _ M U L _ A R E A 1
- *@brief
- *  Following interpretation #1, where y[j] stores the total (integral
- *  or area) value within the interval, return the area under the whole curve.
- *  This is simply totaling up the areas from each of the intervals.
- *  The curve value is found by multiplying corresponding entries from
- *  in1 and in2.
- */
 double
 bn_tabdata_mul_area1(const struct bn_tabdata *in1, const struct bn_tabdata *in2)
 {
@@ -566,14 +472,6 @@
     return area;
 }
 
-/*
- *			B N _ T A B D A T A _ M U L _ A R E A 2
- *@brief
- *  Following interpretation #2,
- *  return the area under the whole curve.
- *  The curve value is found by multiplying corresponding entries from
- *  in1 and in2.
- */
 double
 bn_tabdata_mul_area2(const struct bn_tabdata *in1, const struct bn_tabdata *in2)
 {
@@ -631,13 +529,6 @@
     return -1;
 }
 
-/*
- *			B N _ T A B L E _ L I N _ I N T E R P
- *@brief
- *  Return the value of the curve at independent parameter value 'wl'.
- *  Linearly interpolate between values in the input table.
- *  Zero is returned for values outside the sampled range.
- */
 fastf_t
 bn_table_lin_interp(const struct bn_tabdata *samp, register double wl)
 {
@@ -681,17 +572,6 @@
     return ret;
 }
 
-/*
- *			B N _ T A B D A T A _ R E S A M P L E _ M A X
- *@brief
- *  Given a set of sampled data 'olddata', resample it for different
- *  spacing, by linearly interpolating the values when an output span
- *  is entirely contained within an input span, and by taking the
- *  maximum when an output span covers more than one input span.
- *
- *  This assumes interpretation (2) of the data, i.e. that the values
- *  are the average value across the interval.
- */
 struct bn_tabdata *
 bn_tabdata_resample_max(const struct bn_table *newtable, const struct bn_tabdata *olddata) {
     const struct bn_table	*oldtable;
@@ -761,17 +641,6 @@
     return newsamp;
 }
 
-/*
- *			B N _ T A B D A T A _ R E S A M P L E _ A V G
- *@brief
- *  Given a set of sampled data 'olddata', resample it for different
- *  spacing, by linearly interpolating the values when an output span
- *  is entirely contained within an input span, and by taking the
- *  average when an output span covers more than one input span.
- *
- *  This assumes interpretation (2) of the data, i.e. that the values
- *  are the average value across the interval.
- */
 struct bn_tabdata *
 bn_tabdata_resample_avg(const struct bn_table *newtable, const struct bn_tabdata *olddata) {
     const struct bn_table	*oldtable;
@@ -847,13 +716,6 @@
     return newsamp;
 }
 
-/*
- *			B N _ T A B L E _ W R I T E
- *@brief
- *  Write out the table structure in an ASCII file,
- *  giving the number of values (minus 1), and the
- *  actual values.
- */
 int
 bn_table_write(const char *filename, const struct bn_table *tabp)
 {
@@ -884,13 +746,6 @@
     return 0;
 }
 
-/*
- *			B N _ T A B L E _ R E A D
- *@brief
- *  Allocate and read in the independent variable values from an ASCII file,
- *  giving the number of samples (minus 1), and the
- *  actual values.
- */
 struct bn_table *
 bn_table_read(const char *filename) {
     int ret;
@@ -946,9 +801,6 @@
     return tabp;
 }
 
-/*
- *			B N _ P R _ T A B L E
- */
 void
 bn_pr_table(const char *title, const struct bn_table *tabp)
 {
@@ -962,9 +814,6 @@
     }
 }
 
-/*
- *			B N _ P R _ T A B D A T A
- */
 void
 bn_pr_tabdata(const char *title, const struct bn_tabdata *data)
 {
@@ -979,16 +828,6 @@
     bu_log("\n");
 }
 
-/*
- *			B N _ P R I N T _ T A B L E _ A N D _ T A B D A T A
- *@brief
- *  Write out a given data table into an ASCII file,
- *  suitable for input to GNUPLOT.
- *
- *	(set term postscript)
- *	(set output "|print-postscript")
- *	(plot "filename" with lines)
- */
 int
 bn_print_table_and_tabdata(const char *filename, const struct bn_tabdata *data)
 {
@@ -1022,17 +861,6 @@
     return 0;
 }
 
-/*
- *			B N _ R E A D _ T A B L E _ A N D _ T A B D A T A
- *@brief
- *  Read in a file which contains two columns of numbers, the first
- *  column being the wavelength, the second column being the sample value
- *  at that wavelength.
- *
- *  A new bn_table structure and one bn_tabdata structure
- *  are created, a pointer to the bn_tabdata structure is returned.
- *  The final wavelength is guessed at.
- */
 struct bn_tabdata *
 bn_read_table_and_tabdata(const char *filename) {
     struct bn_table	*tabp;
@@ -1095,9 +923,6 @@
     return data;
 }
 
-/*
- *			B N _ T A B D A T A _ B I N A R Y _ R E A D
- */
 struct bn_tabdata *
 bn_tabdata_binary_read(const char *filename, size_t num, const struct bn_table *tabp) {
     struct bn_tabdata	*data;
@@ -1159,14 +984,6 @@
     return data;
 }
 
-/*
- *			B N _ T A B D A T A _ M A L L O C _ A R R A Y
- *@brief
- *  Allocate storage for, and initialize, an array of 'num' data table
- *  structures.
- *  This subroutine is provided because the bn_tabdata structures
- *  are variable length.
- */
 struct bn_tabdata *
 bn_tabdata_malloc_array(const struct bn_table *tabp, size_t num) {
     struct bn_tabdata *data;
@@ -1197,9 +1014,6 @@
     return data;
 }
 
-/*
- *			B N _ T A B D A T A _ C O P Y
- */
 void
 bn_tabdata_copy(struct bn_tabdata *out, const struct bn_tabdata *in)
 {
@@ -1216,9 +1030,6 @@
     memcpy((char *)out->y, (const char *)in->y, BN_SIZEOF_TABDATA_Y(in));
 }
 
-/*
- *			B N _ T A B D A T A _ D U P
- */
 struct bn_tabdata *
 bn_tabdata_dup(const struct bn_tabdata *in) {
     struct bn_tabdata *data;
@@ -1232,12 +1043,6 @@
     return data;
 }
 
-/*
- *			B N _ T A B D A T A _ G E T _ C O N S T V A L
- *@brief
- *  For a given table, allocate and return a tabdata structure
- *  with all elements initialized to 'val'.
- */
 struct bn_tabdata *
 bn_tabdata_get_constval(double val, const struct bn_table *tabp) {
     struct bn_tabdata	*data;
@@ -1256,11 +1061,6 @@
     return data;
 }
 
-/*
- *			B N _ T A B D A T A _ C O N S T V A L
- *@brief
- *  Set all the tabdata elements to 'val'
- */
 void
 bn_tabdata_constval(struct bn_tabdata *data, double val)
 {
@@ -1276,13 +1076,6 @@
 	*op++ = val;
 }
 
-/*
- *			B N _ T A B D A T A _ T O _ T C L
- *@brief
- *  Convert an bn_tabdata/bn_table pair into a Tcl compatible string
- *  appended to a VLS.  It will have form:
- *	x {...} y {...} nx # ymin # ymax #
- */
 void
 bn_tabdata_to_tcl(struct bu_vls *vp, const struct bn_tabdata *data)
 {
@@ -1312,16 +1105,6 @@
 		  tabp->nx, minval, maxval);
 }
 
-/*
- *			B N _ T A B D A T A _ F R O M _ A R R A Y
- *@brief
- *  Given an array of (x, y) pairs, build the relevant bn_table and
- *  bn_tabdata structures.
- *
- *  The table is terminated by an x value <= 0.
- *  Consistent with the interpretation of the spans,
- *  invent a final span ending x value.
- */
 struct bn_tabdata *
 bn_tabdata_from_array(const double *array) {
     register const double *dp;
@@ -1352,12 +1135,6 @@
     return data;
 }
 
-/*
- *			B N _ T A B D A T A _ F R E Q _ S H I F T
- *@brief
- *  Shift the data by a constant offset in the independent variable
- *  (often frequency), interpolating new sample values.
- */
 void
 bn_tabdata_freq_shift(struct bn_tabdata *out, const struct bn_tabdata *in, double offset)
 {
@@ -1380,11 +1157,6 @@
     }
 }
 
-/*
- *			B N _ T A B L E _ I N T E R V A L _ N U M _ S A M P L E S
- *@brief
- *  Returns number of sample points between 'low' and 'hi', inclusive.
- */
 int
 bn_table_interval_num_samples(const struct bn_table *tabp, double low, double hi)
 {
@@ -1401,13 +1173,6 @@
     return count;
 }
 
-/*
- *			B N _ T A B L E _ D E L E T E _ S A M P L E _ P T S
- *@brief
- *  Remove all sampling points between subscripts i and j, inclusive.
- *  Don't bother freeing the tiny bit of storage at the end of the array.
- *  Returns number of points removed.
- */
 int
 bn_table_delete_sample_pts(struct bn_table *tabp, unsigned int i, unsigned int j)
 {
@@ -1433,12 +1198,6 @@
     return tokill;
 }
 
-/*
- *			B N _ T A B L E _ M E R G E 2
- *@brief
- *  A new table is returned which has sample points at places from
- *  each of the input tables.
- */
 struct bn_table *
 bn_table_merge2(const struct bn_table *a, const struct bn_table *b) {
     struct bn_table *newtab;
@@ -1481,18 +1240,6 @@
     return newtab;
 }
 
-/*
- *		B N _ T A B D A T A _ M K _ L I N E A R _ F I L T E R
- *@brief
- *  Create a filter to accept power in a given band.
- *  The first and last filter values will be in the range 0..1,
- *  while all the internal filter values will be 1.0,
- *  and all samples outside the given band will be 0.0.
- *
- *
- *  @return	NULL	if given band does not overlap input spectrum
- *  @return	tabdata*
- */
 struct bn_tabdata *
 bn_tabdata_mk_linear_filter(const struct bn_table *spectrum, double lower_wavelen, double upper_wavelen) {
     struct bn_tabdata *filt;
Index: src/libbn/msr.c
===================================================================
--- src/libbn/msr.c	(revision 49726)
+++ src/libbn/msr.c	(working copy)
@@ -17,6 +17,7 @@
  * License along with this file; see the file named COPYING for more
  * information.
  */
+
 /** @addtogroup msr */
 /** @{ */
 /** @file libbn/msr.c
@@ -43,38 +44,13 @@
  * Note: BN_MSR_MAXTBL must be an even number, preferably a power of two.
  */
 #define	BN_MSR_MAXTBL	4096	/* Size of random number tables. */
-
-
-/*	bn_unif_init	Initialize a random number structure.
- *
- * @par Entry
- *	setseed	seed to use
- *	method  method to use to generate numbers;
- *
- * @par Exit
- *	returns	a pointer to a bn_unif structure.
- *	returns 0 on error.
- *
- * @par Uses
- *	None.
- *
- * @par Calls
- *	bu_malloc
- *
- * @par Method @code
- *	malloc up a structure with pointers to the numbers
- *	get space for the integer table.
- *	get space for the floating table.
- *	set pointer counters
- *	set seed if one was given and setseed != 1
- @endcode
- *
- */
 #define	A	16807
 #define M	2147483647
 #define DM	2147483647.0
 #define Q	127773		/* Q = M / A */
 #define R	2836		/* R = M % A */
+
+
 struct bn_unif *
 bn_unif_init(long int setseed, int method)
 {
@@ -93,35 +69,8 @@
     p->magic = BN_UNIF_MAGIC;
     return p;
 }
-
-/*	bn_unif_long_fill	fill a random number table.
- *
- * Use the msrad algorithm to fill a random number table
- * with values from 1 to 2^31-1.  These numbers can (and are) extracted from
- * the random number table via high speed macros and bn_unif_long_fill called
- * when the table is exauseted.
- *
- * @par Entry
- *	p	pointer to a bn_unif structure.
- *
- * @par Exit
- *	if (!p) returns 1 else returns a value between 1 and 2^31-1
- *
- * @par Calls
- *	None.  msran is inlined for speed reasons.
- *
- * @par Uses
- *	None.
- *
- * @par Method @code
- *	if (!p) return 1;
- *	if p->msr_longs != NULL
- *		msr_longs is reloaded with random numbers;
- *		msr_long_ptr is set to BN_MSR_MAXTBL
- *	endif
- *	msr_seed is updated.
- @endcode
-*/
+
+
 long
 bn_unif_long_fill(struct bn_unif *p)
 {
@@ -155,35 +104,8 @@
     p->msr_seed =  (test < 0) ? test+M : test;
     return p->msr_seed;
 }
-
-/*	bn_unif_double_fill	fill a random number table.
- *
- * Use the msrad algorithm to fill a random number table
- * with values from -0.5 to 0.5.  These numbers can (and are) extracted from
- * the random number table via high speed macros and bn_unif_double_fill
- * called when the table is exauseted.
- *
- * @par Entry
- *	p	pointer to a bn_unif structure.
- *
- * @par Exit
- *	if (!p) returns 0.0 else returns a value between -0.5 and 0.5
- *
- * @par Calls
- *	None.  msran is inlined for speed reasons.
- *
- * @par Uses
- *	None.
- *
- * @par Method @code
- *	if (!p) return (0.0)
- *	if p->msr_longs != NULL
- *		msr_longs is reloaded with random numbers;
- *		msr_long_ptr is set to BN_MSR_MAXTBL
- *	endif
- *	msr_seed is updated.
- @endcode
-*/
+
+
 double
 bn_unif_double_fill(struct bn_unif *p)
 {
@@ -219,9 +141,7 @@
     return((p->msr_seed - M/2) * 1.0/DM);
 }
 
-/*	bn_unif_free	free random number table
- *
- */
+
 void
 bn_unif_free(struct bn_unif *p)
 {
@@ -231,31 +151,7 @@
     bu_free(p, "bn_unif");
 }
 
-
-/*	bn_gauss_init	Initialize a random number struct for gaussian
- *	numbers.
- *
- * @par Entry
- *	setseed		Seed to use.
- *	method		method to use to generate numbers (not used)
- *
- * @par Exit
- *	Returns a pointer toa bn_msr_guass structure.
- *	returns 0 on error.
- *
- * @par Calls
- *	bu_malloc
- *
- * @par Uses
- *	None.
- *
- * @par Method @code
- *	malloc up a structure
- *	get table space
- *	set seed and pointer.
- *	if setseed != 0 then seed = setseed
- @endcode
-*/
+
 struct bn_gauss *
 bn_gauss_init(long int setseed, int method)
 {
@@ -275,38 +171,8 @@
     p->magic = BN_GAUSS_MAGIC;
     return p;
 }
-
-/*	bn_gauss_fill	fill a random number table.
- *
- * Use the msrad algorithm to fill a random number table.
- * hese numbers can (and are) extracted from
- * the random number table via high speed macros and bn_msr_guass_fill
- * called when the table is exauseted.
- *
- * @par Entry
- *	p	pointer to a bn_msr_guass structure.
- *
- * @par Exit
- *	if (!p) returns 0.0 else returns a value with a mean of 0 and
- *	    a variance of 1.0.
- *
- * @par Calls
- *	BN_UNIF_CIRCLE to get to uniform random number whos radius is
- *	<= 1.0. I.e. sqrt(v1*v1 + v2*v2) <= 1.0
- *	BN_UNIF_CIRCLE is a macro which can call bn_unif_double_fill.
- *
- * @par Uses
- *	None.
- *
- * @par Method @code
- *	if (!p) return (0.0)
- *	if p->msr_longs != NULL
- *		msr_longs is reloaded with random numbers;
- *		msr_long_ptr is set to BN_MSR_MAXTBL
- *	endif
- *	msr_seed is updated.
- @endcode
-*/
+
+
 double
 bn_gauss_fill(struct bn_gauss *p)
 {
@@ -332,9 +198,8 @@
     fac = sqrt(-2.0*log(r)/r);
     return v1*fac;
 }
-/*	bn_gauss_free	free random number table
- *
- */
+
+
 void
 bn_gauss_free(struct bn_gauss *p)
 {
Index: src/libbn/qmath.c
===================================================================
--- src/libbn/qmath.c	(revision 49726)
+++ src/libbn/qmath.c	(working copy)
@@ -61,11 +61,6 @@
 #include "bn.h"
 
 
-/**
- *			Q U A T _ M A T 2 Q U A T
- *@brief
- *  Convert Matrix to Quaternion.
- */
 void
 quat_mat2quat(register fastf_t *quat, register const fastf_t *mat)
 {
@@ -125,14 +120,7 @@
 #undef MMM
 }
 
-/**
- *			Q U A T _ Q U A T 2 M A T
- *@brief
- *  Convert Quaternion to Matrix.
- *
- * NB: This only works for UNIT quaternions.  We may get imaginary results
- *   otherwise.  We should normalize first (still yields same rotation).
- */
+
 void
 quat_quat2mat(register fastf_t *mat, register const fastf_t *quat)
 {
@@ -159,11 +147,7 @@
     mat[15] = 1.0;
 }
 
-/**
- *			Q U A T _ D I S T A N C E
- *@brief
- * Gives the euclidean distance between two quaternions.
- */
+
 double
 quat_distance(const fastf_t *q1, const fastf_t *q2)
 {
@@ -173,14 +157,7 @@
     return QMAGNITUDE( qtemp );
 }
 
-/**
- *			Q U A T _ D O U B L E
- *@brief
- * Gives the quaternion point representing twice the rotation
- *   from q1 to q2.
- *   Needed for patching Bezier curves together.
- *   A rather poor name admittedly.
- */
+
 void
 quat_double(fastf_t *qout, const fastf_t *q1, const fastf_t *q2)
 {
@@ -193,13 +170,6 @@
     QUNITIZE( qout );
 }
 
-/**
- *			Q U A T _ B I S E C T
- *@brief
- * Gives the bisector of quaternions q1 and q2.
- * (Could be done with quat_slerp and factor 0.5)
- * [I believe they must be unit quaternions this to work]
- */
 void
 quat_bisect(fastf_t *qout, const fastf_t *q1, const fastf_t *q2)
 {
@@ -207,15 +177,7 @@
     QUNITIZE( qout );
 }
 
-/**
- *			Q U A T _ S L E R P
- *@brief
- * Do Spherical Linear Interpolation between two unit quaternions
- *  by the given factor.
- *
- * As f goes from 0 to 1, qout goes from q1 to q2.
- * Code based on code by Ken Shoemake
- */
+
 void
 quat_slerp(fastf_t *qout, const fastf_t *q1, const fastf_t *q2, double f)
 {
@@ -260,16 +222,7 @@
     }
 }
 
-/**
- *			Q U A T _ S B E R P
- *@brief
- * Spherical Bezier Interpolate between four quaternions by amount f.
- * These are intended to be used as start and stop quaternions along
- *   with two control quaternions chosen to match spline segments with
- *   first order continuity.
- *
- *  Uses the method of successive bisection.
- */
+
 void
 quat_sberp(fastf_t *qout, const fastf_t *q1, const fastf_t *qa, const fastf_t *qb, const fastf_t *q2, double f)
 {
@@ -288,16 +241,7 @@
     quat_slerp( qout, p4, p5, f );
 }
 
-/**
- *			Q U A T _ M A K E _ N E A R E S T
- *@brief
- *  Set the quaternion q1 to the quaternion which yields the
- *   smallest rotation from q2 (of the two versions of q1 which
- *   produce the same orientation).
- *
- * Note that smallest euclidian distance implies smallest great
- *   circle distance as well (since surface is convex).
- */
+
 void
 quat_make_nearest(fastf_t *q1, const fastf_t *q2)
 {
@@ -314,9 +258,7 @@
     }
 }
 
-/**
- *			Q U A T _ P R I N T
- */
+
 /* DEBUG ROUTINE */
 void
 quat_print(const char *title, const fastf_t *quat)
@@ -336,12 +278,7 @@
 	     axis[X], axis[Y], axis[Z] );
 }
 
-/**
- *			Q U A T _ E X P
- *@brief
- *  Exponentiate a quaternion, assuming that the scalar part is 0.
- *  Code by Ken Shoemake.
- */
+
 void
 quat_exp(fastf_t *out, const fastf_t *in)
 {
@@ -357,12 +294,7 @@
     out[W] = cos(theta);
 }
 
-/**
- *			Q U A T _ L O G
- *@brief
- *  Take the natural logarithm of a unit quaternion.
- *  Code by Ken Shoemake.
- */
+
 void
 quat_log(fastf_t *out, const fastf_t *in)
 {
Index: src/libbn/poly.c
===================================================================
--- src/libbn/poly.c	(revision 49726)
+++ src/libbn/poly.c	(working copy)
@@ -17,6 +17,7 @@
  * License along with this file; see the file named COPYING for more
  * information.
  */
+
 /** @addtogroup poly */
 /** @{ */
 /** @file libbn/poly.c
@@ -62,11 +63,6 @@
 }
 
 
-/**
- * bn_poly_mul
- *
- * @brief multiply two polynomials
- */
 struct bn_poly *
 bn_poly_mul(register struct bn_poly *product, register const struct bn_poly *m1, register const struct bn_poly *m2)
 {
@@ -116,11 +112,6 @@
 }
 
 
-/**
- * bn_poly_scale
- * @brief
- * scale a polynomial
- */
 struct bn_poly *
 bn_poly_scale(register struct bn_poly *eqn, double factor)
 {
@@ -133,11 +124,6 @@
 }
 
 
-/**
- * bn_poly_add
- * @brief
- * add two polynomials
- */
 struct bn_poly *
 bn_poly_add(register struct bn_poly *sum, register const struct bn_poly *poly1, register const struct bn_poly *poly2)
 {
@@ -167,11 +153,6 @@
 }
 
 
-/**
- * bn_poly_sub
- * @brief
- * subtract two polynomials
- */
 struct bn_poly *
 bn_poly_sub(register struct bn_poly *diff, register const struct bn_poly *poly1, register const struct bn_poly *poly2)
 {
@@ -203,12 +184,6 @@
 }
 
 
-/**
- * s y n D i v
- * @brief
- * Divides any polynomial into any other polynomial using synthetic
- * division.  Both polynomials must have real coefficients.
- */
 void
 bn_poly_synthetic_division(register struct bn_poly *quo, register struct bn_poly *rem, register const struct bn_poly *dvdend, register const struct bn_poly *dvsor)
 {
@@ -239,15 +214,6 @@
 }
 
 
-/**
- * b n _ p o l y _ q u a d r a t i c _ r o o t s
- *@brief
- * Uses the quadratic formula to find the roots (in `complex' form) of
- * any quadratic equation with real coefficients.
- *
- *	@return 1 for success
- *	@return 0 for fail.
- */
 int
 bn_poly_quadratic_roots(register struct bn_complex *roots, register const struct bn_poly *quadrat)
 {
@@ -308,36 +274,6 @@
 }
 
 
-/**
- * b n _ p o l y _ c u b i c _ r o o t s
- *@brief
- * Uses the cubic formula to find the roots (in `complex' form)
- * of any cubic equation with real coefficients.
- *
- * to solve a polynomial of the form:
- *
- * X**3 + c1*X**2 + c2*X + c3 = 0,
- *
- * first reduce it to the form:
- *
- * Y**3 + a*Y + b = 0,
- *
- * where
- * Y = X + c1/3,
- * and
- * a = c2 - c1**2/3,
- * b = (2*c1**3 - 9*c1*c2 + 27*c3)/27.
- *
- * Then we define the value delta,   D = b**2/4 + a**3/27.
- *
- * If D > 0, there will be one real root and two conjugate
- * complex roots.
- * If D = 0, there will be three real roots at least two of
- * which are equal.
- * If D < 0, there will be three unequal real roots.
- *
- * Returns 1 for success, 0 for fail.
- */
 int
 bn_poly_cubic_roots(register struct bn_complex *roots, register const struct bn_poly *eqn)
 {
@@ -436,15 +372,6 @@
 }
 
 
-/**
- * b n _ p o l y _ q u a r t i c _ r o o t s
- *@brief
- * Uses the quartic formula to find the roots (in `complex' form)
- * of any quartic equation with real coefficients.
- *
- *	@return 1 for success
- *	@return 0 for fail.
- */
 int
 bn_poly_quartic_roots(register struct bn_complex *roots, register const struct bn_poly *eqn)
 {
@@ -524,11 +451,6 @@
 }
 
 
-/**
- * b n _ p r _ p o l y
- *
- * Print out the polynomial.
- */
 void
 bn_pr_poly(const char *title, register const struct bn_poly *eqn)
 {
@@ -568,11 +490,7 @@
     bu_vls_free(&str);
 }
 
-/**
- * b n _ p r _ r o o t s
- *
- * Print out the roots of a given polynomial (complex numbers)
- */
+
 void
 bn_pr_roots(const char *title, const struct bn_complex *roots, int n)
 {
@@ -585,6 +503,7 @@
 }
 
 /** @} */
+
 /*
  * Local Variables:
  * mode: C
Index: src/libbn/vectfont.c
===================================================================
--- src/libbn/vectfont.c	(revision 49726)
+++ src/libbn/vectfont.c	(working copy)
@@ -17,6 +17,7 @@
  * License along with this file; see the file named COPYING for more
  * information.
  */
+
 /** @addtogroup plot */
 /** @{ */
 /** @file libbn/vectfont.c
@@ -56,11 +57,6 @@
 
 int *tp_cindex[256];	/* index to stroke tokens */
 
-/**
- * @brief
- *  Once-only setup routine
- *  Used by libplot3/symbol.c, so it can't be static.
- */
 void
 tp_setup(void)
 {
Index: src/libbn/multipoly.c
===================================================================
--- src/libbn/multipoly.c	(revision 49726)
+++ src/libbn/multipoly.c	(working copy)
@@ -17,6 +17,7 @@
  * License along with this file; see the file named COPYING for more
  * information.
  */
+
 /** @addtogroup multipoly */
 /** @{ */
 /** @file libbn/multipoly.c
Index: src/libbn/vert_tree.c
===================================================================
--- src/libbn/vert_tree.c	(revision 49726)
+++ src/libbn/vert_tree.c	(working copy)
@@ -17,6 +17,7 @@
  * License along with this file; see the file named COPYING for more
  * information.
  */
+
 /** @addtogroup vtree */
 /** @{ */
 /** @file libbn/vert_tree.c
@@ -73,12 +74,6 @@
 #define VERT_NODE	'n'
 
 
-/**		C R E A T E _ V E R T _ T R E E
- *@brief
- *	routine to create a vertex tree.
- *
- *	Possible refinements include specifying an initial size
- */
 struct vert_root *
 create_vert_tree()
 {
@@ -95,12 +90,6 @@
     return tree;
 }
 
-/**		C R E A T E _ V E R T _ T R E E _ W _ N O R M S
- *@brief
- *	routine to create a vertex tree.
- *
- *	Possible refinements include specifying an initial size
- */
 struct vert_root *
 create_vert_tree_w_norms()
 {
@@ -133,11 +122,6 @@
 
 }
 
-/**		C L E A N _ V E R T _ T R E E
- *@brief
- *	Routine to free the binary search tree and reset the current number of vertices.
- *	The vertex array is left untouched, for re-use later.
- */
 void
 clean_vert_tree( struct vert_root *tree_root )
 {
@@ -167,10 +151,6 @@
 
 }
 
-/**		F R E E _ V E R T_ T R E E
- *@brief
- *	Routine to free a vertex tree and all associated dynamic memory
- */
 void
 free_vert_tree( struct vert_root *vert_root )
 {
@@ -197,12 +177,6 @@
     vert_root->max_vert = 0;
 }
 
-/**		A D D _ V E R T
- *@brief
- *	Routine to add a vertex to the current list of part vertices.
- *	The array is re-alloc'd if needed.
- *	Returns index into the array of vertices where this vertex is stored
- */
 int
 Add_vert( double x, double y, double z, struct vert_root *vert_root, fastf_t local_tol_sq )
 {
@@ -322,12 +296,6 @@
     return new_leaf->vleaf.index;
 }
 
-/**		A D D _ V E R T _ A N D _ N O R M
- *@brief
- *	Routine to add a vertex and a normal to the current list of part vertices.
- *	The array is re-alloc'd if needed.
- *	Returns index into the array of vertices where this vertex and normal is stored
- */
 int
 Add_vert_and_norm( double x, double y, double z, double nx, double ny, double nz, struct vert_root *vert_root, fastf_t local_tol_sq )
 {
@@ -463,6 +431,7 @@
     /* return the index into the vertex array */
     return new_leaf->vleaf.index;
 }
+
 /** @} */
 /*
  * Local Variables:
Index: src/libbn/plane.c
===================================================================
--- src/libbn/plane.c	(revision 49726)
+++ src/libbn/plane.c	(working copy)
@@ -38,11 +38,7 @@
 
 #define UNIT_SQ_TOL 1.0e-13
 
-/**
- * B N _ D I S T _ P T 3 _ P T 3
- * @brief
- * Returns distance between two points.
- */
+
 double
 bn_dist_pt3_pt3(const fastf_t *a, const fastf_t *b)
 {
@@ -53,12 +49,6 @@
 }
 
 
-/**
- * B N _ P T 3 _ P T 3 _ E Q U A L
- *
- * @return 1	if the two points are equal, within the tolerance
- * @return 0	if the two points are not "the same"
- */
 int
 bn_pt3_pt3_equal(const fastf_t *a, const fastf_t *b, const struct bn_tol *tol)
 {
@@ -103,17 +93,6 @@
 }
 
 
-/**
- * B N _ 3 P T S _ C O L L I N E A R
- * @brief
- * Check to see if three points are collinear.
- *
- * The algorithm is designed to work properly regardless of the order
- * in which the points are provided.
- *
- * @return 1	If 3 points are collinear
- * @return 0	If they are not
- */
 int
 bn_3pts_collinear(fastf_t *a, fastf_t *b, fastf_t *c, const struct bn_tol *tol)
 {
@@ -166,16 +145,6 @@
 }
 
 
-/**
- * B N _ 3 P T S _ D I S T I N C T
- *
- * Check to see if three points are all distinct, i.e., ensure that
- * there is at least sqrt(dist_tol_sq) distance between every pair of
- * points.
- *
- * @return 1 If all three points are distinct
- * @return 0 If two or more points are closer together than dist_tol_sq
- */
 int
 bn_3pts_distinct(const fastf_t *a, const fastf_t *b, const fastf_t *c, const struct bn_tol *tol)
 {
@@ -194,16 +163,6 @@
 }
 
 
-/**
- * B N _ N P T S _ D I S T I N C T
- *
- * Check to see if the points are all distinct, i.e., ensure that
- * there is at least sqrt(dist_tol_sq) distance between every pair of
- * points.
- *
- * @return 1 If all the points are distinct
- * @return 0 If two or more points are closer together than dist_tol_sq
- */
 int
 bn_npts_distinct(const int npt, const point_t *pts, const struct bn_tol *tol)
 {
@@ -222,50 +181,6 @@
 }
 
 
-/**
- * B N _ M K _ P L A N E _ 3 P T S
- *
- * Find the equation of a plane that contains three points.  Note that
- * normal vector created is expected to point out (see vmath.h), so
- * the vector from A to C had better be counter-clockwise (about the
- * point A) from the vector from A to B.  This follows the BRL-CAD
- * outward-pointing normal convention, and the right-hand rule for
- * cross products.
- *
- @verbatim
- *
- *			C
- *	                *
- *	                |\
- *	                | \
- *	   ^     N      |  \
- *	   |      \     |   \
- *	   |       \    |    \
- *	   |C-A     \   |     \
- *	   |         \  |      \
- *	   |          \ |       \
- *	               \|        \
- *	                *---------*
- *	                A         B
- *			   ----->
- *		            B-A
- @endverbatim
- *
- * If the points are given in the order A B C (e.g.,
- * *counter*-clockwise), then the outward pointing surface normal:
- *
- * N = (B-A) x (C-A).
- *
- *  @return 0	OK
- *  @return -1	Failure.  At least two of the points were not distinct,
- *		or all three were colinear.
- *
- * @param[out]	plane	The plane equation is stored here.
- * @param[in]	a	point 1
- * @param[in]	b	point 2
- * @param[in]	c	point 3
- * @param[in]	tol	Tolerance values for doing calcualtion
- */
 int
 bn_mk_plane_3pts(fastf_t *plane,
 		 const fastf_t *a,
@@ -303,34 +218,6 @@
 }
 
 
-/**
- * B N _ M K P O I N T _ 3 P L A N E S
- *@brief
- * Given the description of three planes, compute the point of intersection, if
- * any. The direction vectors of the planes need not be of unit length.
- *
- * Find the solution to a system of three equations in three unknowns:
- @verbatim
- * Px * Ax + Py * Ay + Pz * Az = -A3;
- * Px * Bx + Py * By + Pz * Bz = -B3;
- * Px * Cx + Py * Cy + Pz * Cz = -C3;
- *
- * OR
- *
- * [ Ax  Ay  Az ]   [ Px ]   [ -A3 ]
- * [ Bx  By  Bz ] * [ Py ] = [ -B3 ]
- * [ Cx  Cy  Cz ]   [ Pz ]   [ -C3 ]
- *
- @endverbatim
- *
- * @return 0	OK
- * @return -1	Failure.  Intersection is a line or plane.
- *
- * @param[out] pt	The point of intersection is stored here.
- * @param	a	plane 1
- * @param	b	plane 2
- * @param	c	plane 3
- */
 int
 bn_mkpoint_3planes(fastf_t *pt, const fastf_t *a, const fastf_t *b, const fastf_t *c)
 {
@@ -384,17 +271,6 @@
 }
 
 
-/**
- * B N _ 2 L I N E 3 _ C O L I N E A R
- * @brief
- * Returns non-zero if the 3 lines are colinear to within tol->dist
- * over the given distance range.
- *
- * Range should be at least one model diameter for most applications.
- * 1e5 might be OK for a default for "vehicle sized" models.
- *
- * The direction vectors do not need to be unit length.
- */
 int
 bn_2line3_colinear(const fastf_t *p1,
 		   const fastf_t *d1,
@@ -442,31 +318,6 @@
 }
 
 
-/**
- * Find the distance from a point P to a line described by the
- * endpoint A and direction dir, and the point of closest approach
- * (PCA).
- *
- @code
- 			P
- 		       *
- 		      /.
- 		     / .
- 		    /  .
- 		   /   . (dist)
- 		  /    .
- 		 /     .
- 		*------*-------->
- 		A      PCA    dir
- @endcode
- * There are three distinct cases, with these return codes -
- *   0 => P is within tolerance of point A.  *dist = 0, pca=A.
- *   1 => P is within tolerance of line.  *dist = 0, pca=computed.
- *   2 => P is "above/below" line.  *dist=|PCA-P|, pca=computed.
- *
- * TODO: For efficiency, a version of this routine that provides the
- * distance squared would be faster.
- */
 int
 bn_dist_pt3_line3(fastf_t *dist, fastf_t *pca, const fastf_t *a, const fastf_t *dir, const fastf_t *p, const struct bn_tol *tol)
 {
@@ -514,25 +365,6 @@
 }
 
 
-/**
- * Calculate closest approach of two lines
- *
- * returns:
- *	-2 -> lines are parallel and do not intersect
- *	-1 -> lines are parallel and collinear
- *	 0 -> lines intersect
- *	 1 -> lines do not intersect
- *
- * For return values less than zero, dist is not set.  For return
- * valuse of 0 or 1, dist[0] is the distance from p1 in the d1
- * direction to the point of closest approach for that line.  Similar
- * for the second line.  d1 and d2 must be unit direction vectors.
- *
- * XXX How is this different from bn_isect_line3_line3() ?
- * XXX Why are the calling sequences just slightly different?
- * XXX Can we pick the better one, and get rid of the other one?
- * XXX If not, can we document how they differ?
- */
 int
 bn_dist_line3_line3(fastf_t *dist, const fastf_t *p1, const fastf_t *d1, const fastf_t *p2, const fastf_t *d2, const struct bn_tol *tol)
 {
@@ -591,27 +423,6 @@
 }
 
 
-/**
- * calculate intersection or closest approach of a line and a line
- * segement.
- *
- * returns:
- *	-2 -> line and line segment are parallel and collinear.
- *	-1 -> line and line segment are parallel, not collinear.
- *	 0 -> intersection between points a and b.
- *	 1 -> intersection outside a and b.
- *	 2 -> closest approach is between a and b.
- *	 3 -> closest approach is outside a and b.
- *
- * dist[0] is actual distance from p in d direction to
- * closest portion of segment.
- * dist[1] is ratio of distance from a to b (0.0 at a, and 1.0 at b),
- * dist[1] may be less than 0 or greater than 1.
- * For return values less than 0, closest approach is defined as
- * closest to point p.
- * Direction vector, d, must be unit length.
- *
- */
 int
 bn_dist_line3_lseg3(fastf_t *dist, const fastf_t *p, const fastf_t *d, const fastf_t *a, const fastf_t *b, const struct bn_tol *tol)
 {
@@ -665,26 +476,6 @@
 }
 
 
-/**
- * B N _ I S E C T _ L I N E 3 _ P L A N E
- *
- * Intersect an infinite line (specified in point and direction vector
- * form) with a plane that has an outward pointing normal.  The
- * direction vector need not have unit length.  The first three
- * elements of the plane equation must form a unit lengh vector.
- *
- * @return -2	missed (ray is outside halfspace)
- * @return -1	missed (ray is inside)
- * @return 0	line lies on plane
- * @return 1	hit (ray is entering halfspace)
- * @return 2	hit (ray is leaving)
- *
- * @param[out]	dist	set to the parametric distance of the intercept
- * @param[in]	pt	origin of ray
- * @param[in]	dir	direction of ray
- * @param[in]	plane	equation of plane
- * @param[in]	tol	tolerance values
- */
 int
 bn_isect_line3_plane(fastf_t *dist,
 		     const fastf_t *pt,
@@ -725,30 +516,6 @@
 }
 
 
-/**
- * B N _ I S E C T _ 2 P L A N E S
- *@brief
- * Given two planes, find the line of intersection between them, if
- * one exists.  The line of intersection is returned in parametric
- * line (point & direction vector) form.
- *
- * In order that all the geometry under consideration be in "front" of
- * the ray, it is necessary to pass the minimum point of the model
- * RPP.  If this convention is unnecessary, just pass (0, 0, 0) as
- * rpp_min.
- *
- * @return 0	success, line of intersection stored in 'pt' and 'dir'
- * @return -1	planes are coplanar
- * @return -2	planes are parallel but not coplanar
- * @return -3	error, should be intersection but unable to find
- *
- * @param[out]	pt	Starting point of line of intersection
- * @param[out]	dir	Direction vector of line of intersection (unit length)
- * @param[in]	a	plane 1 (normal must be unit length)
- * @param[in]	b	plane 2 (normal must be unit length)
- * @param[in]	rpp_min	minimum poit of model RPP
- * @param[in]	tol	tolerance values
- */
 int
 bn_isect_2planes(fastf_t *pt,
 		 fastf_t *dir,
@@ -836,48 +603,6 @@
 }
 
 
-/**
- * B N _ I S E C T _ L I N E 2 _ L I N E 2
- *
- * Intersect two lines, each in given in parametric form:
- @verbatim
-
- X = P + t * D
- and
- X = A + u * C
-
- @endverbatim
- *
- * While the parametric form is usually used to denote a ray (ie,
- * positive values of the parameter only), in this case the full line
- * is considered.
- *
- * The direction vectors C and D need not have unit length.
- *
- * @return -1	no intersection, lines are parallel.
- * @return 0	lines are co-linear
- *@n			dist[0] gives distance from P to A,
- *@n			dist[1] gives distance from P to (A+C) [not same as below]
- * @return 1	intersection found (t and u returned)
- *@n			dist[0] gives distance from P to isect,
- *@n			dist[1] gives distance from A to isect.
- *
- * @param dist When explicit return > 0, dist[0] and dist[1] are the
- * line parameters of the intersection point on the 2 rays.  The
- * actual intersection coordinates can be found by substituting either
- * of these into the original ray equations.
- *
- * @param p	point of first line
- * @param d	direction of first line
- * @param a	point of second line
- * @param c	direction of second line
- * @param tol	tolerance values
- *
- * Note that for lines which are very nearly parallel, but not quite
- * parallel enough to have the determinant go to "zero", the
- * intersection can turn up in surprising places.  (e.g. when
- * det=1e-15 and det1=5.5e-17, t=0.5)
- */
 int
 bn_isect_line2_line2(fastf_t *dist, const fastf_t *p, const fastf_t *d, const fastf_t *a, const fastf_t *c, const struct bn_tol *tol)
 /* dist[2] */
@@ -1032,43 +757,6 @@
 }
 
 
-/**
- * B N _ I S E C T _ L I N E 2 _ L S E G 2
- *@brief
- * Intersect a line in parametric form:
- *
- * X = P + t * D
- *
- * with a line segment defined by two distinct points A and B=(A+C).
- *
- * XXX probably should take point B, not vector C.  Sigh.
- *
- * @return -4	A and B are not distinct points
- * @return -3	Lines do not intersect
- * @return -2	Intersection exists, but outside segemnt, < A
- * @return -1	Intersection exists, but outside segment, > B
- * @return 0	Lines are co-linear (special meaning of dist[1])
- * @return 1	Intersection at vertex A
- * @return 2	Intersection at vertex B (A+C)
- * @return 3	Intersection between A and B
- *
- * Implicit Returns -
- * @param dist	When explicit return >= 0, t is the parameter that describes
- *		the intersection of the line and the line segment.
- *		The actual intersection coordinates can be found by
- *		solving P + t * D.  However, note that for return codes
- *		1 and 2 (intersection exactly at a vertex), it is
- *		strongly recommended that the original values passed in
- *		A or B are used instead of solving P + t * D, to prevent
- *		numeric error from creeping into the position of
- *		the endpoints.
- *
- * @param p	point of first line
- * @param d	direction of first line
- * @param a	point of second line
- * @param c	direction of second line
- * @param tol	tolerance values
- */
 int
 bn_isect_line2_lseg2(fastf_t *dist,
 		     const fastf_t *p,
@@ -1242,37 +930,6 @@
 }
 
 
-/**
- * B N _ I S E C T _ L S E G 2 _ L S E G 2
- *@brief
- * Intersect two 2D line segments, defined by two points and two
- * vectors.  The vectors are unlikely to be unit length.
- *
- * @return -2	missed (line segments are parallel)
- * @return -1	missed (colinear and non-overlapping)
- * @return 0	hit (line segments colinear and overlapping)
- * @return 1	hit (normal intersection)
- *
- * @param dist  The value at dist[] is set to the parametric distance of the
- *		intercept.
- *@n	dist[0] is parameter along p, range 0 to 1, to intercept.
- *@n	dist[1] is parameter along q, range 0 to 1, to intercept.
- *@n	If within distance tolerance of the endpoints, these will be
- *	exactly 0.0 or 1.0, to ease the job of caller.
- *
- * Special note: when return code is "0" for co-linearity, dist[1] has
- * an alternate interpretation: it's the parameter along p (not q)
- * which takes you from point p to the point (q + qdir), i.e., its
- * the endpoint of the q linesegment, since in this case there may be
- * *two* intersections, if q is contained within span p to (p + pdir).
- * And either may be -10 if the point is outside the span.
- *
- * @param p	point 1
- * @param pdir	direction1
- * @param q	point 2
- * @param qdir	direction2
- * @param tol	tolerance values
- */
 int
 bn_isect_lseg2_lseg2(fastf_t *dist,
 		     const fastf_t *p,
@@ -1344,46 +1001,6 @@
 }
 
 
-/**
- * B N _ I S E C T _ L S E G 3 _ L S E G 3
- *@brief
- * Intersect two 3D line segments, defined by two points and two
- * vectors.  The vectors are unlikely to be unit length.
- *
- *
- * @return -3	missed
- * @return -2	missed (line segments are parallel)
- * @return -1	missed (colinear and non-overlapping)
- * @return 0	hit (line segments colinear and overlapping)
- * @return 1	hit (normal intersection)
- *
- * @param[out] dist
- *	The value at dist[] is set to the parametric distance of the
- *	intercept dist[0] is parameter along p, range 0 to 1, to
- *	intercept.  dist[1] is parameter along q, range 0 to 1, to
- *	intercept.  If within distance tolerance of the endpoints,
- *	these will be exactly 0.0 or 1.0, to ease the job of caller.
- *
- *      CLARIFICATION: This function 'bn_isect_lseg3_lseg3'
- *      returns distance values scaled where an intersect at the start
- *      point of the line segement (within tol->dist) results in 0.0
- *      and when the intersect is at the end point of the line
- *      segement (within tol->dist), the result is 1.0.  Intersects
- *      before the start point return a negative distance.  Intersects
- *      after the end point result in a return value > 1.0.
- *
- * Special note: when return code is "0" for co-linearity, dist[1] has
- * an alternate interpretation: it's the parameter along p (not q)
- * which takes you from point p to the point (q + qdir), i.e., it's
- * the endpoint of the q linesegment, since in this case there may be
- * *two* intersections, if q is contained within span p to (p + pdir).
- *
- * @param p	point 1
- * @param pdir	direction-1
- * @param q	point 2
- * @param qdir	direction-2
- * @param tol	tolerance values
- */
 int
 bn_isect_lseg3_lseg3(fastf_t *dist,
 		     const fastf_t *p,
@@ -1528,48 +1145,6 @@
 }
 
 
-/**
- * B N _ I S E C T _ L I N E 3 _ L I N E 3
- *
- * Intersect two line segments, each in given in parametric form:
- *
- * X = p0 + pdist * pdir_i   (i.e. line p0->p1)
- * and
- * X = q0 + qdist * qdir_i   (i.e. line q0->q1)
- *
- * The input vectors 'pdir_i' and 'qdir_i' must NOT be unit vectors
- * for this function to work correctly. The magnitude of the direction
- * vectors indicates line segment length.
- *
- * The 'pdist' and 'qdist' values returned from this function are the
- * actual distance to the intersect (i.e. not scaled). Distances may
- * be negative, see below.
- *
- * @return -2	no intersection, lines are parallel.
- * @return -1	no intersection
- * @return 0	lines are co-linear (pdist and qdist returned) (see below)
- * @return 1	intersection found  (pdist and qdist returned) (see below)
- *
- * @param	p0	point 1
- * @param	pdir_i 	direction 1
- * @param	q0	point 2
- * @param	qdir_i 	direction 2
- * @param tol	tolerance values
- * @param[out]	pdist, qdist (distances to intersection) (see below)
- *
- *		When return = 1, pdist is the distance along line p0->p1 to the
- *		intersect with line q0->q1. If the intersect is along p0->p1 but
- *		in the opposite direction of vector pdir_i (i.e. occuring before
- *		p0 on line p0->p1) then the distance will be negative. The value
- *		if qdist is the same as pdist except it is the distance along line
- *		q0->q1 to the intersect with line p0->p1.
- *
- *		When return code = 0 for co-linearity, pdist and qdist have a
- *		different meaning. pdist is the distance from point p0 to point q0
- *		and qdist is the distance from point p0 to point q1. If point q0
- *		occurs before point p0 on line segment p0->p1 then pdist will be
- *		negative. The same occurs for the distance to point q1.
- */
 int
 bn_isect_line3_line3(fastf_t *pdist,        /* see above */
 		     fastf_t *qdist,        /* see above */
@@ -1702,38 +1277,6 @@
 }
 
 
-/**
- * B N _ I S E C T _ L I N E _ L S E G
- *@brief
- * Intersect a line in parametric form:
- *
- * X = P + t * D
- *
- * with a line segment defined by two distinct points A and B.
- *
- *
- * @return -4	A and B are not distinct points
- * @return -3	Intersection exists, < A (t is returned)
- * @return -2	Intersection exists, > B (t is returned)
- * @return -1	Lines do not intersect
- * @return 0	Lines are co-linear (t for A is returned)
- * @return 1	Intersection at vertex A
- * @return 2	Intersection at vertex B
- * @return 3	Intersection between A and B
- *
- * @par Implicit Returns -
- *
- * t When explicit return >= 0, t is the parameter that describes the
- * intersection of the line and the line segment.  The actual
- * intersection coordinates can be found by solving P + t * D.
- * However, note that for return codes 1 and 2 (intersection exactly
- * at a vertex), it is strongly recommended that the original values
- * passed in A or B are used instead of solving P + t * D, to prevent
- * numeric error from creeping into the position of the endpoints.
- *
- * XXX should probably be called bn_isect_line3_lseg3()
- * XXX should probably be changed to return dist[2]
- */
 int
 bn_isect_line_lseg(fastf_t *t, const fastf_t *p, const fastf_t *d, const fastf_t *a, const fastf_t *b, const struct bn_tol *tol)
 {
@@ -1926,22 +1469,6 @@
 }
 
 
-/**
- * B N _ D I S T _ L I N E 3_ P T 3
- *@brief
- * Given a parametric line defined by PT + t * DIR and a point A,
- * return the closest distance between the line and the point.
- *
- *  'dir' need not have unit length.
- *
- * Find parameter for PCA along line with unitized DIR:
- * d = VDOT(f, dir) / MAGNITUDE(dir);
- * Find distance g from PCA to A using Pythagoras:
- * g = sqrt(MAGSQ(f) - d**2)
- *
- * Return -
- * Distance
- */
 double
 bn_dist_line3_pt3(const fastf_t *pt, const fastf_t *dir, const fastf_t *a)
 {
@@ -1968,18 +1495,6 @@
 }
 
 
-/**
- * B N _ D I S T S Q _ L I N E 3 _ P T 3
- *
- * Given a parametric line defined by PT + t * DIR and a point A,
- * return the square of the closest distance between the line and the
- * point.
- *
- * 'dir' need not have unit length.
- *
- * Return -
- * Distance squared
- */
 double
 bn_distsq_line3_pt3(const fastf_t *pt, const fastf_t *dir, const fastf_t *a)
 {
@@ -2005,16 +1520,6 @@
 }
 
 
-/**
- * B N _ D I S T _ L I N E _ O R I G I N
- *@brief
- * Given a parametric line defined by PT + t * DIR, return the closest
- * distance between the line and the origin.
- *
- * 'dir' need not have unit length.
- *
- * @return Distance
- */
 double
 bn_dist_line_origin(const fastf_t *pt, const fastf_t *dir)
 {
@@ -2029,16 +1534,6 @@
 }
 
 
-/**
- * B N _ D I S T _ L I N E 2 _ P O I N T 2
- *@brief
- * Given a parametric line defined by PT + t * DIR and a point A,
- * return the closest distance between the line and the point.
- *
- * 'dir' need not have unit length.
- *
- * @return Distance
- */
 double
 bn_dist_line2_point2(const fastf_t *pt, const fastf_t *dir, const fastf_t *a)
 {
@@ -2055,18 +1550,6 @@
 }
 
 
-/**
- * B N _ D I S T S Q _ L I N E 2 _ P O I N T 2
- *@brief
- * Given a parametric line defined by PT + t * DIR and a point A,
- * return the closest distance between the line and the point,
- * squared.
- *
- * 'dir' need not have unit length.
- *
- * @return
- * Distance squared
- */
 double
 bn_distsq_line2_point2(const fastf_t *pt, const fastf_t *dir, const fastf_t *a)
 {
@@ -2083,11 +1566,6 @@
 }
 
 
-/**
- * B N _ A R E A _ O F _ T R I A N G L E
- *@brief
- * Returns the area of a triangle. Algorithm by Jon Leech 3/24/89.
- */
 double
 bn_area_of_triangle(register const fastf_t *a, register const fastf_t *b, register const fastf_t *c)
 {
@@ -2111,39 +1589,6 @@
 }
 
 
-/**
- * B N _ I S E C T _ P T _ L S E G
- *@brief
- * Intersect a point P with the line segment defined by two distinct
- * points A and B.
- *
- * @return -2	P on line AB but outside range of AB,
- * 			dist = distance from A to P on line.
- * @return -1	P not on line of AB within tolerance
- * @return 1	P is at A
- * @return 2	P is at B
- * @return 3	P is on AB, dist = distance from A to P on line.
- @verbatim
- B *
- |
- P'*-tol-*P
- |    /   _
- dist  /   /|
- |  /   /
- | /   / AtoP
- |/   /
- A *   /
-
- tol = distance limit from line to pt P;
- dist = parametric distance from A to P' (in terms of A to B)
- @endverbatim
- *
- * @param p	point
- * @param a	start of lseg
- * @param b	end of lseg
- * @param tol	tolerance values
- * @param[out] dist	parametric distance from A to P' (in terms of A to B)
- */
 int bn_isect_pt_lseg(fastf_t *dist,
 		     const fastf_t *a,
 		     const fastf_t *b,
@@ -2200,33 +1645,6 @@
 }
 
 
-/**
- * B N _ I S E C T _ P T 2 _ L S E G 2
- * @brief
- * Intersect a point P with the line segment defined by two distinct
- * points A and B.
- *
- * @return -2	P on line AB but outside range of AB,
- *			dist = distance from A to P on line.
- * @return -1	P not on line of AB within tolerance
- * @return 1	P is at A
- * @return 2	P is at B
- * @return 3	P is on AB, dist = distance from A to P on line.
- @verbatim
- B *
- |
- P'*-tol-*P
- |    /  _
- dist  /   /|
- |  /   /
- | /   / AtoP
- |/   /
- A *   /
-
- tol = distance limit from line to pt P;
- dist = distance from A to P'
- @endverbatim
-*/
 int
 bn_isect_pt2_lseg2(fastf_t *dist, const fastf_t *a, const fastf_t *b, const fastf_t *p, const struct bn_tol *tol)
 /* distance along line from A to P */
@@ -2291,38 +1709,6 @@
 }
 
 
-/**
- * B N _ D I S T _ P T 3 _ L S E G 3
- *@brief
- * Find the distance from a point P to a line segment described by the
- * two endpoints A and B, and the point of closest approach (PCA).
- @verbatim
- *
- *			P
- *		       *
- *		      /.
- *		     / .
- *		    /  .
- *		   /   . (dist)
- *		  /    .
- *		 /     .
- *		*------*--------*
- *		A      PCA	B
- @endverbatim
- *
- * @return 0	P is within tolerance of lseg AB.  *dist isn't 0: (SPECIAL!!!)
- *		  *dist = parametric dist = |PCA-A| / |B-A|.  pca=computed.
- * @return 1	P is within tolerance of point A.  *dist = 0, pca=A.
- * @return 2	P is within tolerance of point B.  *dist = 0, pca=B.
- * @return 3	P is to the "left" of point A.  *dist=|P-A|, pca=A.
- * @return 4	P is to the "right" of point B.  *dist=|P-B|, pca=B.
- * @return 5	P is "above/below" lseg AB.  *dist=|PCA-P|, pca=computed.
- *
- * This routine was formerly called bn_dist_pt_lseg().
- *
- * XXX For efficiency, a version of this routine that provides the
- * XXX distance squared would be faster.
- */
 int
 bn_dist_pt3_lseg3(fastf_t *dist,
 		  fastf_t *pca,
@@ -2416,35 +1802,6 @@
 }
 
 
-/**
- * B N _ D I S T _ P T 2 _ L S E G 2
- *@brief
- * Find the distance from a point P to a line segment described by the
- * two endpoints A and B, and the point of closest approach (PCA).
- @verbatim
- *			P
- *		       *
- *		      /.
- *		     / .
- *		    /  .
- *		   /   . (dist)
- *		  /    .
- *		 /     .
- *		*------*--------*
- *		A      PCA	B
- @endverbatim
- * There are six distinct cases, with these return codes -
- * @return 0	P is within tolerance of lseg AB.  *dist isn't 0: (SPECIAL!!!)
- *		  *dist = parametric dist = |PCA-A| / |B-A|.  pca=computed.
- * @return 1	P is within tolerance of point A.  *dist = 0, pca=A.
- * @return 2	P is within tolerance of point B.  *dist = 0, pca=B.
- * @return 3	P is to the "left" of point A.  *dist=|P-A|**2, pca=A.
- * @return 4	P is to the "right" of point B.  *dist=|P-B|**2, pca=B.
- * @return 5	P is "above/below" lseg AB.  *dist=|PCA-P|**2, pca=computed.
- *
- *
- * Patterned after bn_dist_pt3_lseg3().
- */
 int
 bn_dist_pt2_lseg2(fastf_t *dist_sq, fastf_t *pca, const fastf_t *a, const fastf_t *b, const fastf_t *p, const struct bn_tol *tol)
 {
@@ -2534,14 +1891,6 @@
 }
 
 
-/**
- * B N _ R O T A T E _ B B O X
- *@brief
- * Transform a bounding box (RPP) by the given 4x4 matrix.  There are
- * 8 corners to the bounding RPP.  Each one needs to be transformed
- * and min/max'ed.  This is not minimal, but does fully contain any
- * internal object, using an axis-aligned RPP.
- */
 void
 bn_rotate_bbox(fastf_t *omin, fastf_t *omax, const fastf_t *mat, const fastf_t *imin, const fastf_t *imax)
 {
@@ -2566,11 +1915,6 @@
 }
 
 
-/**
- * B N _ R O T A T E _ P L A N E
- *@brief
- * Transform a plane equation by the given 4x4 matrix.
- */
 void
 bn_rotate_plane(fastf_t *oplane, const fastf_t *mat, const fastf_t *iplane)
 {
@@ -2594,18 +1938,6 @@
 }
 
 
-/**
- * B N _ C O P L A N A R
- *@brief
- * Test if two planes are identical.  If so, their dot products will
- * be either +1 or -1, with the distance from the origin equal in
- * magnitude.
- *
- * @return -1	not coplanar, parallel but distinct
- * @return 0	not coplanar, not parallel.  Planes intersect.
- * @return 1	coplanar, same normal direction
- * @return 2	coplanar, opposite normal direction
- */
 int
 bn_coplanar(const fastf_t *a, const fastf_t *b, const struct bn_tol *tol)
 {
@@ -2644,36 +1976,6 @@
 }
 
 
-/**
- * B N _ A N G L E _ M E A S U R E
- *
- * Using two perpendicular vectors (x_dir and y_dir) which lie in the
- * same plane as 'vec', return the angle (in radians) of 'vec' from
- * x_dir, going CCW around the perpendicular x_dir CROSS y_dir.
- *
- * Trig note -
- *
- * theta = atan2(x, y) returns an angle in the range -pi to +pi.
- * Here, we need an angle in the range of 0 to 2pi.  This could be
- * implemented by adding 2pi to theta when theta is negative, but this
- * could have nasty numeric ambiguity right in the vicinity of theta =
- * +pi, which is a very critical angle for the applications using this
- * routine.
- *
- * So, an alternative formulation is to compute gamma = atan2(-x, -y),
- * and then theta = gamma + pi.  Now, any error will occur in the
- * vicinity of theta = 0, which can be handled much more readily.
- *
- * If theta is negative, or greater than two pi, wrap it around.
- * These conditions only occur if there are problems in atan2().
- *
- * @return vec == x_dir returns 0,
- * @return vec == y_dir returns pi/2,
- * @return vec == -x_dir returns pi,
- * @return vec == -y_dir returns 3*pi/2.
- *
- * In all cases, the returned value is between 0 and bn_twopi.
- */
 double
 bn_angle_measure(fastf_t *vec, const fastf_t *x_dir, const fastf_t *y_dir)
 {
@@ -2701,14 +2003,6 @@
 }
 
 
-/**
- * B N _ D I S T _ P T 3 _ A L O N G _ L I N E 3
- *@brief
- * Return the parametric distance t of a point X along a line defined
- * as a ray, i.e. solve X = P + t * D.  If the point X does not lie on
- * the line, then t is the distance of the perpendicular projection of
- * point X onto the line.
- */
 double
 bn_dist_pt3_along_line3(const fastf_t *p, const fastf_t *d, const fastf_t *x)
 {
@@ -2719,14 +2013,6 @@
 }
 
 
-/**
- * B N _ D I S T _ P T 2 _ A L O N G _ L I N E 2
- *@brief
- * Return the parametric distance t of a point X along a line defined
- * as a ray, i.e. solve X = P + t * D.  If the point X does not lie on
- * the line, then t is the distance of the perpendicular projection of
- * point X onto the line.
- */
 double
 bn_dist_pt2_along_line2(const fastf_t *p, const fastf_t *d, const fastf_t *x)
 {
@@ -2746,11 +2032,6 @@
 }
 
 
-/**
- *
- * @return 1	if left <= mid <= right
- * @return 0	if mid is not in the range.
- */
 int
 bn_between(double left, double mid, double right, const struct bn_tol *tol)
 {
@@ -2780,12 +2061,6 @@
 }
 
 
-/**
- * B N _ D O E S _ R A Y _ I S E C T _ T R I
- *
- * @return 0	No intersection
- * @return 1	Intersection, 'inter' has intersect point.
- */
 int
 bn_does_ray_isect_tri(
     const point_t pt,
@@ -2944,16 +2219,7 @@
 }
 #endif
 
-/**
- * B N _ H L F _ C L A S S
- *@brief
- * Classify a halfspace, specified by its plane equation, against a
- * bounding RPP.
- *
- * @return BN_CLASSIFY_INSIDE
- * @return BN_CLASSIFY_OVERLAPPING
- * @return BN_CLASSIFY_OUTSIDE
- */
+
 int
 bn_hlf_class(const fastf_t *half_eqn, const fastf_t *min, const fastf_t *max, const struct bn_tol *tol)
 {
@@ -2990,30 +2256,6 @@
 }
 
 
-/** B N _ D I S T S Q _ L I N E 3 _ L I N E 3
- *@brief
- * Calculate the square of the distance of closest approach for two
- * lines.
- *
- * The lines are specifed as a point and a vector each.  The vectors
- * need not be unit length.  P and d define one line; Q and e define
- * the other.
- *
- * @return 0 - normal return
- * @return 1 - lines are parallel, dist[0] is set to 0.0
- *
- * Output values:
- * dist[0] is the parametric distance along the first line P + dist[0] * d of the PCA
- * dist[1] is the parametric distance along the second line Q + dist[1] * e of the PCA
- * dist[3] is the square of the distance between the points of closest approach
- * pt1 is the point of closest approach on the first line
- * pt2 is the point of closest approach on the second line
- *
- * This algoritm is based on expressing the distance sqaured, taking
- * partials with respect to the two unknown parameters (dist[0] and
- * dist[1]), seeting the two partails equal to 0, and solving the two
- * simutaneous equations
- */
 int
 bn_distsq_line3_line3(fastf_t *dist, fastf_t *P, fastf_t *d_in, fastf_t *Q, fastf_t *e_in, fastf_t *pt1, fastf_t *pt2)
 {
@@ -3065,35 +2307,6 @@
 }
 
 
-/**
- * B N _ I S E C T _ P L A N E S
- *@brief
- * Calculates the point that is the minimum distance from all the
- * planes in the "planes" array.  If the planes intersect at a single
- * point, that point is the solution.
- *
- * The method used here is based on:
-
- * An expression for the distance from a point to a plane is:
- * VDOT(pt, plane)-plane[H].
- * Square that distance and sum for all planes to get the "total"
- * distance.
- * For minimum total distance, the partial derivatives of this
- * expression (with respect to x, y, and z) must all be zero.
- * This produces a set of three equations in three unknowns (x, y, z).
-
- * This routine sets up the three equations as [matrix][pt] = [hpq]
- * and solves by inverting "matrix" into "inverse" and
- * [pt] = [inverse][hpq].
- *
- * There is likely a more economical solution rather than matrix
- * inversion, but bn_mat_inv was handy at the time.
- *
- * Checks if these planes form a singular matrix and returns.
- *
- * @return 0 - all is well
- * @return 1 - planes form a singular matrix (no solution)
- */
 int
 bn_isect_planes(fastf_t *pt, const fastf_t (*planes)[4], const size_t pl_count)
 {
Index: src/libbn/wavelet.c
===================================================================
--- src/libbn/wavelet.c	(revision 49726)
+++ src/libbn/wavelet.c	(working copy)
@@ -17,6 +17,7 @@
  * License along with this file; see the file named COPYING for more
  * information.
  */
+
 /** @addtogroup wavelet */
 /** @{ */
 /** @file libbn/wavelet.c
Index: src/libbn/noise.c
===================================================================
--- src/libbn/noise.c	(revision 49726)
+++ src/libbn/noise.c	(working copy)
@@ -218,18 +218,6 @@
 }
 
 
-/**
- *@brief
- * Robert Skinner's Perlin-style "Noise" function
- *
- * Results are in the range [-0.5 .. 0.5].  Unlike many
- * implementations, this function provides random noise at the integer
- * lattice values.  However this produces much poorer quality and
- * should be avoided if possible.
- *
- * The power distribution of the result has no particular shape,
- * though it isn't as flat as the literature would have one believe.
- */
 double
 bn_noise_perlin(fastf_t *point)
 {
@@ -538,32 +526,8 @@
 
     return ep;
 }
-/**
- * @brief
- * Procedural fBm evaluated at "point"; returns value stored in
- * "value".
- *
- * @param point          location to sample noise
- * @param ``h_val''      fractal increment parameter
- * @param ``lacunarity'' gap between successive frequencies
- * @param ``octaves''  	 number of frequencies in the fBm
- *
- * The spectral properties of the result are in the APPROXIMATE range
- * [-1..1] Depending upon the number of octaves computed, this range
- * may be exceeded.  Applications should clamp or scale the result to
- * their needs.  The results have a more-or-less gaussian
- * distribution.  Typical results for 1M samples include:
- *
- * @li Min -1.15246
- * @li Max  1.23146
- * @li Mean -0.0138744
- * @li s.d.  0.306642
- * @li Var 0.0940295
- *
- * The function call pow() is relatively expensive.  Therfore, this
- * function pre-computes and saves the spectral weights in a table for
- * re-use in successive invocations.
- */
+
+
 double
 bn_noise_fbm(fastf_t *point, double h_val, double lacunarity, double octaves)
 {
@@ -608,32 +572,6 @@
 } /* bn_noise_fbm() */
 
 
-/**
- * @brief
- * Procedural turbulence evaluated at "point";
- *
- * @return turbulence value for point
- *
- * @param point          location to sample noise at
- * @param ``h_val''      fractal increment parameter
- * @param ``lacunarity'' gap between successive frequencies
- * @param ``octaves''    number of frequencies in the fBm
- *
- * The result is characterized by sharp, narrow trenches in low values
- * and a more fbm-like quality in the mid-high values.  Values are in
- * the APPROXIMATE range [0 .. 1] depending upon the number of octaves
- * evaluated.  Typical results:
- @code
- * Min  0.00857137
- * Max  1.26712
- * Mean 0.395122
- * s.d. 0.174796
- * Var  0.0305536
- @endcode
- * The function call pow() is relatively expensive.  Therfore, this
- * function pre-computes and saves the spectral weights in a table for
- * re-use in successive invocations.
- */
 double
 bn_noise_turb(fastf_t *point, double h_val, double lacunarity, double octaves)
 {
@@ -703,12 +641,7 @@
 
 } /* bn_noise_turb() */
 
-/**
- *@brief
- * A ridged noise pattern
- *
- * From "Texturing and Modeling, A Procedural Approach" 2nd ed p338
- */
+
 double
 bn_noise_ridged(fastf_t *point, double h_val, double lacunarity, double octaves, double offset)
 {
@@ -763,10 +696,7 @@
     return result;
 }
 
-/**
- *
- * From "Texturing and Modeling, A Procedural Approach" 2nd ed
- */
+
 double
 bn_noise_mf(fastf_t *point, double h_val, double lacunarity, double octaves, double offset)
 {
Index: src/libbn/mat.c
===================================================================
--- src/libbn/mat.c	(revision 49726)
+++ src/libbn/mat.c	(working copy)
@@ -87,9 +87,6 @@
 }
 
 
-/**
- * B N _ M A T _ P R I N T
- */
 void
 bn_mat_print(const char *title,
 	     const mat_t m)
@@ -101,9 +98,6 @@
 }
 
 
-/**
- * B N _ M A T _ P R I N T _ V L S
- */
 void
 bn_mat_print_vls(const char *title,
 		 const mat_t m,
@@ -116,12 +110,6 @@
 }
 
 
-/**
- * B N _ A T A N 2
- *
- * A wrapper for the system atan2().  On the Silicon Graphics, and
- * perhaps on others, x==0 incorrectly returns infinity.
- */
 double
 bn_atan2(double y, double x)
 {
@@ -137,14 +125,6 @@
 }
 
 
-/**
- * B N _ M A T _ M U L
- *
- * Multiply matrix "a" by "b" and store the result in "o".
- *
- * This is different from multiplying "b" by "a" (most of the time!)
- * Also, "o" must not be the same as either of the inputs.
- */
 void
 bn_mat_mul(register mat_t o, register const mat_t a, register const mat_t b)
 {
@@ -170,14 +150,6 @@
 }
 
 
-/**
- * B N _ M A T _ M U L 2
- *
- * o = i * o
- *
- * A convenience wrapper for bn_mat_mul() to update a matrix in place.
- * The argument ordering is confusing either way.
- */
 void
 bn_mat_mul2(register const mat_t i, register mat_t o)
 {
@@ -188,14 +160,6 @@
 }
 
 
-/**
- * B N _ M A T _ M U L 3
- *
- * o = a * b * c
- *
- * The output matrix may be the same as 'b' or 'c', but may not be
- * 'a'.
- */
 void
 bn_mat_mul3(mat_t o, const mat_t a, const mat_t b, const mat_t c)
 {
@@ -206,13 +170,6 @@
 }
 
 
-/**
- * B N _ M A T _ M U L 4
- *
- * o = a * b * c * d
- *
- * The output matrix may be the same as any input matrix.
- */
 void
 bn_mat_mul4(mat_t ao, const mat_t a, const mat_t b, const mat_t c, const mat_t d)
 {
@@ -224,13 +181,6 @@
 }
 
 
-/**
- * B N _ M A T X V E C
- *
- * Multiply the matrix "im" by the vector "iv" and store the result in
- * the vector "ov".  Note this is post-multiply, and operates on
- * 4-tuples.  Use MAT4X3VEC() to operate on 3-tuples.
- */
 void
 bn_matXvec(register vect_t ov, register const mat_t im, register const vect_t iv)
 {
@@ -249,14 +199,6 @@
 }
 
 
-/**
- * B N _ M A T _ I N V
- *
- * The matrix pointed at by "input" is inverted and stored in the area
- * pointed at by "output".
- *
- * Calls bu_bomb if matrix is singular.
- */
 void
 bn_mat_inv(register mat_t output, const mat_t input)
 {
@@ -272,19 +214,6 @@
 }
 
 
-/**
- * B N _ M A T _ I N V E R S E
- *
- * The matrix pointed at by "input" is inverted and stored in the area
- * pointed at by "output".
- *
- * Invert a 4-by-4 matrix using Algorithm 120 from ACM.  This is a
- * modified Gauss-Jordan alogorithm Note: Inversion is done in place,
- * with 3 work vectors
- *
- * @return 1 if OK.
- * @return 0 if matrix is singular.
- */
 int
 bn_mat_inverse(register mat_t output, const mat_t input)
 {
@@ -370,12 +299,6 @@
 }
 
 
-/**
- * B N _ V T O H _ M O V E
- *
- * Takes a pointer to a [x, y, z] vector, and a pointer to space for a
- * homogeneous vector [x, y, z, w], and builds [x, y, z, 1].
- */
 void
 bn_vtoh_move(register vect_t h, register const vect_t v)
 {
@@ -384,15 +307,6 @@
 }
 
 
-/**
- * B N _ H T O V _ M O V E
- *
- * Takes a pointer to [x, y, z, w], and converts it to an ordinary
- * vector [x/w, y/w, z/w].  Optimization for the case of w==1 is
- * performed.
- *
- * FIXME: make tolerance configurable
- */
 void
 bn_htov_move(register vect_t v, register const vect_t h)
 {
@@ -411,9 +325,6 @@
 }
 
 
-/**
- * B N _ M A T _ T R N
- */
 void
 bn_mat_trn(mat_t om, register const mat_t im)
 {
@@ -421,15 +332,6 @@
 }
 
 
-/**
- * B N _ M A T _ A E
- *
- * Compute a 4x4 rotation matrix given Azimuth and Elevation.
- *
- * Azimuth is +X, Elevation is +Z, both in degrees.
- *
- * Formula due to Doug Gwyn, BRL.
- */
 void
 bn_mat_ae(register fastf_t *m, double azimuth, double elev)
 {
@@ -464,12 +366,6 @@
 }
 
 
-/**
- * B N _ A E _ V E C
- *
- * Find the azimuth and elevation angles that correspond to the
- * direction (not including twist) given by a direction vector.
- */
 void
 bn_ae_vec(fastf_t *azp, fastf_t *elp, const vect_t v)
 {
@@ -486,15 +382,6 @@
 }
 
 
-/**
- * B N _ A E T _ V E C
- *
- * Find the azimuth, elevation, and twist from two vectors.  Vec_ae is
- * in the direction of view (+z in mged view) and vec_twist points to
- * the viewers right (+x in mged view).  Accuracy (degrees) is used to
- * stabilze flutter between equivalent extremes of atan2(), and to
- * snap twist to zero when elevation is near +/- 90
- */
 void
 bn_aet_vec(fastf_t *az, fastf_t *el, fastf_t *twist, fastf_t *vec_ae, fastf_t *vec_twist, fastf_t accuracy)
 {
@@ -530,11 +417,6 @@
 }
 
 
-/**
- * B N _ V E C _ A E
- *
- * Find a unit vector from the origin given azimuth and elevation.
- */
 void
 bn_vec_ae(vect_t vect, fastf_t az, fastf_t el)
 {
@@ -547,11 +429,6 @@
 }
 
 
-/**
- * B N _ V E C _ A E D
- *
- * Find a vector from the origin given azimuth, elevation, and distance.
- */
 void
 bn_vec_aed(vect_t vect, fastf_t az, fastf_t el, fastf_t distance)
 {
@@ -564,18 +441,6 @@
 }
 
 
-/**
- * B N _ M A T _ A N G L E S
- *
- * This routine builds a Homogeneous rotation matrix, given alpha,
- * beta, and gamma as angles of rotation, in degrees.
- *
- * Alpha is angle of rotation about the X axis, and is done third.
- * Beta is angle of rotation about the Y axis, and is done second.
- * Gamma is angle of rotation about Z axis, and is done first.
- *
- * FIXME: make tolerance configurable
- */
 void
 bn_mat_angles(register fastf_t *mat, double alpha_in, double beta_in, double ggamma_in)
 {
@@ -634,18 +499,6 @@
 }
 
 
-/**
- * B N _ M A T _ A N G L E S _ R A D
- *
- * This routine builds a Homogeneous rotation matrix, given alpha,
- * beta, and gamma as angles of rotation, in radians.
- *
- * Alpha is angle of rotation about the X axis, and is done third.
- * Beta is angle of rotation about the Y axis, and is done second.
- * Gamma is angle of rotation about Z axis, and is done first.
- *
- * FIXME: make tolerance configurable
- */
 void
 bn_mat_angles_rad(register mat_t mat,
 		  double alpha,
@@ -687,16 +540,6 @@
 }
 
 
-/**
- * B N _ E I G E N 2 X 2
- *
- * Find the eigenvalues and eigenvectors of a symmetric 2x2 matrix.
- * (a b)
- * (b c)
- *
- * The eigenvalue with the smallest absolute value is returned in
- * val1, with its eigenvector in vec1.
- */
 void
 bn_eigen2x2(fastf_t *val1, fastf_t *val2, fastf_t *vec1, fastf_t *vec2, fastf_t a, fastf_t b, fastf_t c)
 {
@@ -741,15 +584,6 @@
 }
 
 
-/**
- * B N _ V E C _ P E R P
- *
- * Given a vector, create another vector which is perpendicular to it.
- * The output vector will have unit length only if the input vector
- * did.
- *
- * FIXME: make tolerance configurable
- */
 void
 bn_vec_perp(vect_t new, const vect_t old)
 {
@@ -769,18 +603,6 @@
 }
 
 
-/**
- * B N _ M A T _ F R O M T O
- *
- * Given two vectors, compute a rotation matrix that will transform
- * space by the angle between the two.  There are many candidate
- * matricies.
- *
- * The input 'from' and 'to' vectors need not be unit length.
- * MAT4X3VEC(to, m, from) is the identity that is created.
- *
- * FIXME: make tolerance configurable
- */
 void
 bn_mat_fromto(mat_t m, const vect_t from, const vect_t to)
 {
@@ -856,12 +678,6 @@
 }
 
 
-/**
- * B N _ M A T _ X R O T
- *
- * Given the sin and cos of an X rotation angle, produce the rotation
- * matrix.
- */
 void
 bn_mat_xrot(fastf_t *m, double sinx, double cosx)
 {
@@ -885,12 +701,6 @@
 }
 
 
-/**
- * B N _ M A T _ Y R O T
- *
- * Given the sin and cos of a Y rotation angle, produce the rotation
- * matrix.
- */
 void
 bn_mat_yrot(fastf_t *m, double siny, double cosy)
 {
@@ -914,12 +724,6 @@
 }
 
 
-/**
- * B N _ M A T _ Z R O T
- *
- * Given the sin and cos of a Z rotation angle, produce the rotation
- * matrix.
- */
 void
 bn_mat_zrot(fastf_t *m, double sinz, double cosz)
 {
@@ -943,28 +747,6 @@
 }
 
 
-/**
- * B N _ M A T _ L O O K A T
- *
- * Given a direction vector D of unit length, product a matrix which
- * rotates that vector D onto the -Z axis.  This matrix will be
- * suitable for use as a "model2view" matrix.
- *
- * XXX This routine will fail if the vector is already more or less
- * aligned with the Z axis.
- *
- * This is done in several steps.
- *
- @code
- 1) Rotate D about Z to match +X axis.  Azimuth adjustment.
- 2) Rotate D about Y to match -Y axis.  Elevation adjustment.
- 3) Rotate D about Z to make projection of X axis again point
- in the +X direction.  Twist adjustment.
- 4) Optionally, flip sign on Y axis if original Z becomes inverted.
- This can be nice for static frames, but is astonishing when
- used in animation.
- @endcode
-*/
 void
 bn_mat_lookat(mat_t rot, const vect_t dir, int yflip)
 {
@@ -1020,15 +802,6 @@
 }
 
 
-/**
- * B N _ V E C _ O R T H O
- *
- * Given a vector, create another vector which is perpendicular to it,
- * and with unit length.  This algorithm taken from Gift's arvec.f; a
- * faster algorithm may be possible.
- *
- * FIXME: make tolerance configurable
- */
 void
 bn_vec_ortho(register vect_t out, register const vect_t in)
 {
@@ -1073,16 +846,6 @@
 }
 
 
-/**
- * B N _ M A T _ S C A L E _ A B O U T _ P T
- *
- * Build a matrix to scale uniformly around a given point.
- *
- * @return -1 if scale is too small.
- * @return  0 if OK.
- *
- * FIXME: make tolerance configurable
- */
 int
 bn_mat_scale_about_pt(mat_t mat, const point_t pt, const double scale)
 {
@@ -1108,12 +871,6 @@
 }
 
 
-/**
- * B N _ M A T _ X F O R M _ A B O U T _ P T
- *
- * Build a matrix to apply arbitary 4x4 transformation around a given
- * point.
- */
 void
 bn_mat_xform_about_pt(mat_t mat, const mat_t xform, const point_t pt)
 {
@@ -1130,12 +887,6 @@
 }
 
 
-/**
- * B N _ M A T _ I S _ E Q U A L
- *
- * @return 0 When matrices are not equal
- * @return 1 When matricies are equal
- */
 int
 bn_mat_is_equal(const mat_t a, const mat_t b, const struct bn_tol *tol)
 {
@@ -1181,19 +932,6 @@
 }
 
 
-/**
- * B N _ M A T _ I S _ I D E N T I T Y
- *
- * This routine is intended for detecting identity matricies read in
- * from ascii or binary files, where the numbers are pure ones or
- * zeros.  This routine is *not* intended for tolerance-based
- * "near-zero" comparisons; as such, it shouldn't be used on matrices
- * which are the result of calculation.
- *
- *
- * @return 0 non-identity
- * @return 1 a perfect identity matrix
- */
 int
 bn_mat_is_identity(const mat_t m)
 {
@@ -1201,15 +939,6 @@
 }
 
 
-/**
- * B N _ M A T _ A R B _ R O T
- *
- * Construct a transformation matrix for rotation about an arbitrary
- * axis.  The axis is defined by a point (pt) and a unit direction
- * vector (dir).  The angle of rotation is "ang"
- *
- * FIXME: make tolerance configurable
- */
 void
 bn_mat_arb_rot(mat_t m, const point_t pt, const vect_t dir, const fastf_t ang)
 {
@@ -1265,12 +994,6 @@
 }
 
 
-/**
- * B N _ M A T _ D U P
- *
- * Return a pointer to a copy of the matrix in dynamically allocated
- * memory.
- */
 matp_t
 bn_mat_dup(const mat_t in)
 {
@@ -1282,17 +1005,6 @@
 }
 
 
-/**
- * B N _ M A T _ C K
- *
- * Check to ensure that a rotation matrix preserves axis
- * perpendicularily.  Note that not all matricies are rotation
- * matricies.
- *
- *
- * @return -1 FAIL
- * @return  0 OK
- */
 int
 bn_mat_ck(const char *title, const mat_t m)
 {
@@ -1338,12 +1050,6 @@
 }
 
 
-/**
- * B N _ M A T _ D E T 3
- *
- * Calculates the determinant of the 3X3 "rotation" part of the passed
- * matrix.
- */
 fastf_t
 bn_mat_det3(const mat_t m)
 {
@@ -1357,11 +1063,6 @@
 }
 
 
-/**
- * B N _ M A T _ D E T E R M I N A N T
- *
- * Calculates the determinant of the 4X4 matrix
- */
 fastf_t
 bn_mat_determinant(const mat_t m)
 {
@@ -1391,11 +1092,6 @@
 }
 
 
-/**
- * B N _ M A T _ I S _ N O N _ U N I F
- *
- * FIXME: make tolerance configurable
- */
 int
 bn_mat_is_non_unif (const mat_t m)
 {
@@ -1418,13 +1114,6 @@
 }
 
 
-/**
- * B N _ W R T _ P O I N T _ D I R E C
- *
- * Given a model-space transformation matrix "change",
- * return a matrix which applies the change with-respect-to
- * given "point" and "direc".
- */
 void
 bn_wrt_point_direc(mat_t out, const mat_t change, const mat_t in, const point_t point, const vect_t direc)
 {
Index: src/libbn/tcl.c
===================================================================
--- src/libbn/tcl.c	(revision 49726)
+++ src/libbn/tcl.c	(working copy)
@@ -899,9 +899,6 @@
 }
 
 
-/**
- * B N _ M A T _ P R I N T
- */
 void
 bn_tcl_mat_print(Tcl_Interp *interp,
 		 const char *title,
@@ -914,12 +911,6 @@
 }
 
 
-/**
- * B N _ T C L _ S E T U P
- *@brief
- * Add all the supported Tcl interfaces to LIBBN routines to
- * the list of commands known by the given interpreter.
- */
 void
 bn_tcl_setup(Tcl_Interp *interp)
 {
@@ -954,14 +945,6 @@
 }
 
 
-/**
- * B N _ I N I T
- *@brief
- * Allows LIBBN to be dynamically loade to a vanilla tclsh/wish with
- * "load /usr/brlcad/lib/libbn.so"
- *
- * The name of this function is specified by TCL.
- */
 int
 Bn_Init(Tcl_Interp *interp)
 {
Index: src/libbn/rand.c
===================================================================
--- src/libbn/rand.c	(revision 49726)
+++ src/libbn/rand.c	(working copy)
@@ -32,14 +32,7 @@
 #include "vmath.h"
 #include "bn.h"
 
-/**
- * This is our table of random numbers.  Rather than calling drand48()
- * or random() or rand() we just pick numbers out of this table.  This
- * table has 4096 unique entries with floating point values ranging
- * from the open interval (i.e. exclusive) 0.0 to 1.0 range.
- *
- * There are convenience macros for access in the bn.h header.
- */
+
 const float bn_rand_table[BN_RAND_TABSIZE] = {
     0.39646477f, 0.84048537f, 0.35333610f, 0.44658343f, 0.31869277f, 0.88642843f,
     0.01558285f, 0.58409022f, 0.15936863f, 0.38371587f, 0.69100437f, 0.05885891f,
@@ -729,10 +722,7 @@
 
 double bn_sin_scale = BN_SINTABSIZE / (2.0 * M_PI);
 
-/**
- * table of floating point sine values in the closed (i.e. inclusive)
- * interval -1.0 to 1.0 range.
- */
+
 const float bn_sin_table[BN_SINTABSIZE] = {
     0.0f,		0.00306796f,	0.00613588f,	0.00920375f,
     0.0122715f,	0.0153392f,	0.0184067f,	0.0214741f,
@@ -3634,12 +3624,7 @@
     9e20f, 9e20f, 9e20f, 9e20f, 9e20f, 9e20f, 9e20f, 9e20f, 9e20f, 9e20f
 };
 
-/**
- *			M A T H T A B _ C O N S T A N T
- *@brief
- *  For benchmarking purposes, make the random number table predictable.
- *  Setting to all zeros keeps dithered values at their original values.
- */
+
 void
 bn_mathtab_constant(void)
 {
