*** transform.c	Mon Oct 11 15:00:58 1993
--- transform.c.frac	Mon Oct 11 14:47:18 1993
***************
*** 305,310 ****
--- 305,418 ----
  	MatrixCopy(&ttmp, inverse);
  }
  
+ /* find eigenvalues of positive definite matrices, return FALSE is singular. */
+ int MatrixEigenvalues(trans, lambda)
+ RSMatrix *trans;
+ Float lambda[3];
+ {
+         Float a[3], maxval, a11, a12, a21, a22;
+         Float *A[3];
+ 
+         a[0] = fabs(trans->matrix[0][0]); maxval = a[0];
+         a[1] = fabs(trans->matrix[1][0]); if (a[1] > maxval) maxval = a[1];
+         a[2] = fabs(trans->matrix[2][0]); if (a[2] > maxval) maxval = a[2];
+   
+         if (maxval == a[0]) {
+                 A[0] = trans->matrix[0];
+                 A[1] = trans->matrix[1];
+                 A[2] = trans->matrix[2];
+         } else if (maxval == a[1]) {
+                 A[0] = trans->matrix[1];
+                 A[1] = trans->matrix[0];
+                 A[2] = trans->matrix[2];
+         } else {
+                 A[0] = trans->matrix[2];
+                 A[1] = trans->matrix[1];
+                 A[2] = trans->matrix[0];
+         }
+ 
+         lambda[0] = A[0][0];
+         if (lambda[0] == 0.) return FALSE;
+ 
+         a[1] = A[1][0] / lambda[0];
+         a[2] = A[2][0] / lambda[0];
+         a11 = A[1][1] - a[1] * A[0][1];        
+         a21 = A[2][1] - a[2] * A[0][1];        
+         a12 = A[1][2] - a[1] * A[0][2];        
+         a22 = A[2][2] - a[2] * A[0][2];        
+ 
+         if (fabs(a11) > fabs(a21)) {
+                 lambda[1] = a11;
+                 if (lambda[1] == 0.) return FALSE;
+                 lambda[2] = a22 - a21/a11 * a12;
+         } else {
+                 lambda[1] = a21;
+                 if (lambda[1] == 0.) return FALSE;
+                 lambda[2] = a12 - a11/a21 * a22;
+         }
+ 
+         if (lambda[2] == 0.) return FALSE;
+         return TRUE;
+ }
+ 
+ /* Lipschitz number for euclidian norm: largest scaling factor */
+ Float Lipschitz(lambda)
+ Float lambda[3];
+ {
+         Float maxlambda;
+ 
+         maxlambda = lambda[0];
+         if (lambda[1] > maxlambda) maxlambda = lambda[1];
+         if (lambda[2] > maxlambda) maxlambda = lambda[2];
+         return maxlambda;
+ }
+ 
+ /* returns TRUE if the three scaling factors are equal */
+ int Uniform(lambda)
+ Float lambda[3];
+ {
+         return (equal(lambda[0], lambda[1]) && equal(lambda[1], lambda[2]));
+ }
+ 
+ /* scaling factors */
+ int MatrixScaling(trans, lambda)
+ RSMatrix *trans;
+ Float lambda[3];
+ {
+         RSMatrix Q2;
+         int i;
+ 
+ /* Q-R decomposition: we are only interested in Q**2 */
+         Q2.matrix[0][0] = trans->matrix[0][0] * trans->matrix[0][0] +
+                           trans->matrix[0][1] * trans->matrix[0][1] +
+                           trans->matrix[0][2] * trans->matrix[0][2] ;
+         Q2.matrix[0][1] = trans->matrix[0][0] * trans->matrix[1][0] +
+                           trans->matrix[0][1] * trans->matrix[1][1] +
+                           trans->matrix[0][2] * trans->matrix[1][2] ;
+         Q2.matrix[0][2] = trans->matrix[0][0] * trans->matrix[2][0] +
+                           trans->matrix[0][1] * trans->matrix[2][1] +
+                           trans->matrix[0][2] * trans->matrix[2][2] ;
+         Q2.matrix[1][1] = trans->matrix[1][0] * trans->matrix[1][0] +
+                           trans->matrix[1][1] * trans->matrix[1][1] +
+                           trans->matrix[1][2] * trans->matrix[1][2] ;
+         Q2.matrix[1][2] = trans->matrix[1][0] * trans->matrix[2][0] +
+                           trans->matrix[1][1] * trans->matrix[2][1] +
+                           trans->matrix[1][2] * trans->matrix[2][2] ;
+         Q2.matrix[2][2] = trans->matrix[2][0] * trans->matrix[2][0] +
+                           trans->matrix[2][1] * trans->matrix[2][1] +
+                           trans->matrix[2][2] * trans->matrix[2][2] ;
+         Q2.matrix[1][0] = Q2.matrix[0][1];
+         Q2.matrix[2][0] = Q2.matrix[0][2];
+         Q2.matrix[2][1] = Q2.matrix[1][2];
+ 
+ /* scaling factors are square root of eigenvalues of Q**2 */
+         if (!MatrixEigenvalues(&Q2, lambda)) return FALSE;
+         lambda[0] = sqrt(lambda[0]);
+         lambda[1] = sqrt(lambda[1]);
+         lambda[2] = sqrt(lambda[2]);
+         return TRUE;
+ }
+ 
  /*
   * Apply a transformation to a point (translation affects the point).
   */
***************
*** 399,407 ****
  Ray *ray;
  RSMatrix *trans;
  {
  	PointTransform(&ray->pos, trans);
  	VecTransform(&ray->dir, trans);
! 	return VecNormalize(&ray->dir);
  }
  
  void
--- 507,519 ----
  Ray *ray;
  RSMatrix *trans;
  {
+         Float stretch;
  	PointTransform(&ray->pos, trans);
  	VecTransform(&ray->dir, trans);
! 
! 	stretch = VecNormalize(&ray->dir);
!         ray->stretch *= stretch;
!         return stretch;
  }
  
  void
***************
*** 433,437 ****
--- 545,552 ----
  {
  	TransCopy(list, result);
  	for (list = list->next; list; list = list->next)
+ /* it was: 
  		TransCompose(list, result, result);
+  * PhB made it: */
+ 		TransCompose(result, list, result);
  }