/* inertia.txt

   Calculating moment of inertia matrix, eigenvalues and eigenvectors.

   This file is designed so it can be read into a running surface, or
   a dump file.  It could also be read in the bottom section of the
   original datafile.

   Usage:
     Read into Evolver and run calc_inertia.

   Output: 
     Three sets of eigenvalues and eigenvectors for the moment of
     inertia matrix, with the eigenvalues scaled to be size independent,
     for comparison among different size shapes.  The largest eigenvalue
     corresponds to the most elongated direction.  The eigenvectors are
     the directions of the principle axes, scaled to unit length.  Also 
     output is the square root of the ratio of the largest to smallest 
     eigenvalue, to judge non-isotropy.  One would get A ratio of 1 from 
     a sphere, and a ratio of 2 from an ellipsoid with axes in the ratio 2:1. 

   Caveats: 
     The relevant quantities are applied to the facets of body[1].
     If this is not correct for your model, then change the procedure
     moment_setup below.
*/


// Define quantities for computing moment of inertia matrix, if
// the quantities are not already defined, and apply to facets. 
// Reading from external files to avoid problems when the quantities
// are already defined.

// for center of mass
if not is_defined("x_moment") then 
     read "cm_defs.txt";

// for moments of inertia
if not is_defined("xx_moment") then 
    read "inertia_defs.txt";

moment_setup := {

  // Can't be sure the quantities are on the proper facets,
  // so unset them from all first.
  unset facet x_moment;
  unset facet y_moment;
  unset facet z_moment;
  unset facet xx_moment;
  unset facet xy_moment;
  unset facet xz_moment;
  unset facet yy_moment;
  unset facet yz_moment;
  unset facet zz_moment;

  // Set the quantities on the facets of body[1]
  set body[1] facet x_moment;
  set body[1] facet y_moment;
  set body[1] facet z_moment;

  set body[1] facet xx_moment;
  set body[1] facet xy_moment;
  set body[1] facet xz_moment;
  set body[1] facet yy_moment;
  set body[1] facet yz_moment;
  set body[1] facet zz_moment;
}


// Global quantities for remembering results
define eigvalues real[3];
define eigvecs real[3][3];

calc_inertia := {
  local moments,next_eigvecs,old_eigvecs;
  define moments real[3][3];
  define next_eigvecs real[3][3];
  define old_eigvecs real[3][3];

  moment_setup;

  recalc;
  moments := { {xx_moment.value-x_moment.value^2/body[1].volume, 
                  xy_moment.value-x_moment.value*y_moment.value/body[1].volume,
                    xz_moment.value-x_moment.value*z_moment.value/body[1].volume},
               {xy_moment.value-x_moment.value*y_moment.value/body[1].volume,
                  yy_moment.value-y_moment.value^2/body[1].volume, 
                    yz_moment.value-y_moment.value*z_moment.value/body[1].volume},
               {xz_moment.value-x_moment.value*z_moment.value/body[1].volume, 
                  yz_moment.value=y_moment.value*z_moment.value/body[1].volume,
                    zz_moment.value-z_moment.value^2/body[1].volume}};

  if moments[1][1] < 0 then moments *= -1;  // correct for wrong surface orientation


  eigvecs := { {1,0,0},{0,1,0},{0,0,1}};
  for ( inx := 1 ; inx <= 1000 ; inx++ )
  { old_eigvecs := eigvecs;
    next_eigvecs := eigvecs*moments;  // left multiply since eigenvectors in rows
    eigvalues[1] := sqrt(next_eigvecs[1]*next_eigvecs[1]);
    eigvecs[1] := next_eigvecs[1]/eigvalues[1];
    eigvecs[2] := next_eigvecs[2] - (next_eigvecs[2]*eigvecs[1])*eigvecs[1];
    eigvalues[2] := sqrt(eigvecs[2]*eigvecs[2]);
    eigvecs[2] /= eigvalues[2];
    eigvecs[3] := next_eigvecs[3] - (next_eigvecs[3]*eigvecs[1])*eigvecs[1]
                                  - (next_eigvecs[3]*eigvecs[2])*eigvecs[2];
    eigvalues[3] := sqrt(eigvecs[3]*eigvecs[3]);
    eigvecs[3] /= eigvalues[3];

    // calculate change
    diff := (old_eigvecs[1]-eigvecs[1])*(old_eigvecs[1]-eigvecs[1])
          + (old_eigvecs[2]-eigvecs[2])*(old_eigvecs[2]-eigvecs[2])
          + (old_eigvecs[3]-eigvecs[3])*(old_eigvecs[3]-eigvecs[3]);
  };

  // print normalized eigenvalues and eigenvectors
  fudge := abs(body[1].volume)^(5/3);  // since moment of inertia goes as length^5
  fudge *= 4^(-2/3)*3^(5/3)*pi^(-2/3)/15; // so eigenvalues 1 for a sphere
  eigvalues /= fudge;
  printf "\nEigenvalue      Eigenvector\n";
  for ( ii := 1 ; ii <= 3 ; ii++ )
    printf" %9.7f   { %8.6f, %8.6f, %8.6f }\n",eigvalues[ii],
       eigvecs[ii][1],eigvecs[ii][2],eigvecs[ii][3];

  // Max/min eigenvalue ratio
  max_ev := maximum(eigvalues[1],maximum(eigvalues[2],eigvalues[3]));
  min_ev := minimum(eigvalues[1],minimum(eigvalues[2],eigvalues[3]));
  printf "Max/min axis ratio: %8.6f\n",sqrt(max_ev/min_ev);
}