AnalyticGeometryTool.hpp

Go to the documentation of this file.

//-------------------------------------------------------------------------
// COPYRIGHT 1998  CATERPILLAR INC.  ALL RIGHTS RESERVED
//
// Filename      : AnalyticGeometryTool.hpp
//
// Purpose       : This class performs calculations on analytic geometry
//                 (points, lines, arcs, planes, polygons).  Capabilities 
//                 include vector and point math, matrix operations, 
//                 measurements, intersections and comparison/containment checks.
//
// Special Notes : As most of these functions were taken from Caterpillar's
//                 Cmath library with minimal modification for this class,
//                 the vectors, planes, etc. are represented in non-CUBIT
//                 format.  These could be converted with some effort.
// 
//                 Other Notes:
//                 ------------
//                 1. Points and vectors are represented as double[3].
//
//                 2. Lines (unless denoted as 'segments') are defined as 
//                    a point and a vector (point-direction form).  The parametric
//                    equations of a line are:
//
//                         x = xo + t*a
//                         y = yo + t*a
//                         z = zo + t*a
//
//                 3. Planes are also defined as a point and a vector
//                    (point-normal form).  The equation of a plane is:
//
//                         a*(x-xo) + b*(y-yo) + c*(z-zo) = 0
//                      
//                     For example, the equation of a plane passing through the 
//                     point (3,-1,7) and perpendicular to the vector n = (4,2,-5) 
//                     is: 4(x-3) + 2(y+1) -5(z-7) = 0
//
//                     This can be reduced to:
//
//                         A*x + B*y + C*z + D = 0
//  
//                     For the example, the equation is: 4x + 2y -5z + 25 = 0
//
//                 4. Most functions which can accept a 3x3 matrix have
//                    an analogous 4x4 function.  This is for convenience
//                    only (the extra rows/columns are ignored within the function).
//                 
//
// Creator       : Steve Storm
//
// Creation Date : 10/16/98
//-------------------------------------------------------------------------

#ifndef ANALYTIC_GEOMETRY_TOOL_HPP
#define ANALYTIC_GEOMETRY_TOOL_HPP

#include <math.h>
#include "CubitDefines.h"
#include "CGMGeomConfigure.h"

#include <iostream>

class CubitBox;
class CubitPlane;
class CubitVector;
template <class X> class DLIList;

  // Multiplication constants
  const double AGT_E =             2.7182818284590452354;   // e
  const double AGT_LOG_2E =        1.4426950408889634074;   // log2(e)
  const double AGT_LOG_10E =       0.43429448190325182765; // log10(e)
  const double AGT_LN_2 =          0.69314718055994530942; // ln(2)
  const double AGT_LN_10 =         2.30258509299404568402; // ln(10)
  const double AGT_PI  =           3.14159265358979323846; // PI
  const double AGT_PI_DIV_2 =      1.57079632679489661923; // PI/2
  const double AGT_PI_DIV_4 =      0.78539816339744830962; // PI/4
  const double AGT_1_DIV_PI =      0.31830988618379067154; // 1/PI
  const double AGT_2_DIV_PI =      0.63661977236758134308; // 2/PI
  const double AGT_2_DIV_SQRT_PI = 1.12837916709551257390; // 2/(sqrt(PI))
  const double AGT_SQRT_2 =        1.41421356237309504880; // sqrt(2)
  const double AGT_SQRT_1_DIV_2 =  0.70710678118654752440; // sqrt(1/2)
  const double AGT_PI_DIV_180 =    0.017453292519943295;   // PI/180
  const double AGT_180_DIV_PI =    57.295779513082323;     // 180/PI 

  const double DEGtoRAD = AGT_PI_DIV_180;
  const double RADtoDEG = AGT_180_DIV_PI;

  typedef struct {
     double Vec1[3];    /* First vector that defines plane of arc */
     double Vec2[3];    /* Second vector that defines plane of arc */
     double Origin[3];  /* Center of arc (on plane) */
     double StartAngle; /* Starting angle (in radians) of arc */
     double EndAngle;   /* End angle (in radians) of arc */
     double Radius; /* Radius of arc */
  } AGT_3D_Arc;

  // Return values */
  enum AgtEquality {
     AGT_EQUAL       = 0, // define for comparisons
     AGT_LESSTHAN    = 1, // define for comparisons
     AGT_GREATERTHAN = 2  // define for comparisons
  };

  enum AgtSide {
     AGT_ON_PLANE =  0,  // define for which side of plane
     AGT_POS_SIDE =  1,  // define for which side of plane
     AGT_NEG_SIDE =  2,  // define for which side of plane
     AGT_INT_PLANE = 3,  // define for intersects plane
     AGT_CROSS = 4       // define for line crossing plane
  };

  enum AgtOrientation {
    AGT_OPP_DIR = 0,  // define for vector direction comparision
    AGT_SAME_DIR = 1  // define for vector direction comparision
  };

  enum AgtDistanceMethod {
    AGT_FRACTION = 0,  // define for distance along a curve
    AGT_ABSOLUTE = 1   // define for distance along a curve
  };

// Macros
#define number_sign(number) (number<0 ? -1 : 1)
#define copy_vec(vec1,vec2) copy_pnt(vec1,vec2)
#define set_vec(pnt,x,y,z) set_pnt(pnt,x,y,z)

//                               MAGIC SOFTWARE                              //
// See note later in file regarding MAGIC softare.                           //
// Minimal 2D bounding box - parameters
const int maxPartition = 32;
const int invInterp = 32;
const double angleMin = -0.25*AGT_PI;
const double angleMax = +0.25*AGT_PI;
const double angleRange = 0.5*AGT_PI;

typedef struct
{
    double x, y;
}
Point2;

typedef struct
{
    double x, y, z;
}
Point3;

typedef struct
{
    Point2 center;
    Point2 axis[2];
    double extent[2];
}
OBBox2;

typedef struct
{
    Point3 center;
    Point3 axis[3];
    double extent[3];
}
OBBox3;

// threshold on determinant to conclude lines/rays/segments are parallel
const double par_tolerance = 1e-06;
#define DIST(x) (Sqrt(Abs(x)))

// Line in 2D
typedef struct
{
    // line is Dot(N,X) = c, N not necessarily unit length
    Point2 N;
    double c;
}
AgtLine;

// Triangle in 2D
typedef struct
{
    Point2 v[3];
}
Triangle;

// Linked list of triangle
typedef struct AgtTriList
{
    Triangle* tri;
    AgtTriList* next;
}
AgtTriList;

//---------------------------------------------------------------------------
// Lines in 3D
//---------------------------------------------------------------------------
typedef struct
{
    // Line is L(t) = b+t*m for any real-valued t
    // Ray has constraint t >= 0, b is the origin of the ray
    // Line segment has constraint 0 <= t <= 1, b and b+m are end points

    Point3 b, m;
}
Line3;
//---------------------------------------------------------------------------

//---------------------------------------------------------------------------
// Circles in 3D
//---------------------------------------------------------------------------
typedef struct
{
    // Plane containing circle is Dot(N,X-C) = 0 where X is any point in the
    // plane.  Vectors U, V, and N form an orthonormal right-handed set
    // (matrix [U V N] is orthonormal and has determinant 1).  Circle within
    // the plane is parameterized by X = C + R*(cos(A)*U + sin(A)*V) where
    // A is an angle in [0,2*pi).

    Point3 U, V, N;  // coordinate system of plane containing circle
    Point3 C;  // center
    double R;  // radius > 0
}
Circle3;
//---------------------------------------------------------------------------

//---------------------------------------------------------------------------
// Triangles in 3D
//---------------------------------------------------------------------------
typedef struct
{
    // Triangle points are tri(s,t) = b+s*e0+t*e1 where 0 <= s <= 1,
    // 0 <= t <= 1, and 0 <= s+t <= 1.

    // If the vertices of the triangle are v0, v1, and v2, then b = v0,
    // e0 = v1 - v0, and e1 = v2 - v0.

    Point3 b, e0, e1;
}
Triangle3;
//---------------------------------------------------------------------------

//---------------------------------------------------------------------------
// Parallelograms in 3D
//---------------------------------------------------------------------------
typedef struct
{
    // Rectoid points are rect(s,t) = b+s*e0+t*e1 where 0 <= s <= 1
    // and 0 <= t <= 1.  Could have called this a parallelogram, but
    // I do not like typing long words...

    Point3 b, e0, e1;
}
Rectoid3;
//---------------------------------------------------------------------------

//---------------------------------------------------------------------------
// Planes in 3D
//---------------------------------------------------------------------------
typedef struct
{
    // plane is Dot(N,X) = d
    Point3 N;
    double d;
}
Plane3;
//---------------------------------------------------------------------------

// END OF MAGIC SOFTWARE

class CUBIT_GEOM_EXPORT AnalyticGeometryTool
{
public:

   ~AnalyticGeometryTool();
   static AnalyticGeometryTool* instance();

   static void delete_instance()
   {
     if( instance_ )
       delete instance_;
     instance_ = NULL;
   };

   //*********************************************************
   // Double numbers
   //*********************************************************

   CubitBoolean dbl_eq( double val_1, double val_2 );

   double get_epsilon();
   double set_epsilon( double new_epsilon );

   void round_near_val( double& val );

   //**************************************************************************
   // Matrices & Transforms
   //**************************************************************************
   // Note: For these functions the matrix format is defined as follows:
   //    
   //    Consider the transformation from [x,y,z] to [x',y',z']:    
   //                             _          _
   //                            | x1 y1 z1 0 |
   //    [x',y',z',1] = [x,y,z,1]| x2 y2 z2 0 |
   //                            | x3 y3 z3 0 |
   //                                | ox oy oz 1 |
   //                                 -          -
   void transform_pnt( double m[4][4], double pin[3], double pout[3] );

   void transform_vec( double m3[3][3], double vin[3], double vout[3] );

   void transform_vec( double m4[4][4], double vin[3], double vout[3] );

   void transform_line( double rot_mtx[4][4], double origin[3], double axis[3] );
   void transform_line( double rot_mtx[4][4], CubitVector &origin, CubitVector &axis );

  

   void copy_mtx( double from[3][3],double to[3][3] );
   void copy_mtx( double from[4][4], double to[4][4] );
   void copy_mtx( double from[4][4], double to[3][3] );
   void copy_mtx(double from[3][3], double to[4][4], double *origin = NULL );



   void create_rotation_mtx( double theta, double v[3], double mtx3x3[3][3] );
   void create_rotation_mtx( double theta, double v[3], double mtx4x4[4][4] );

   void add_origin_to_rotation_mtx( double rot_mtx[4][4], double origin[3] );


   void identity_mtx( double mtx3x3[3][3] );
   void identity_mtx( double mtx4x4[4][4] );









   void mtx_to_angs( double mtx3x3[3][3], double &ax, double &ay, double &az );
   void mtx_to_angs( double mtx4x4[4][4], double &ax, double &ay, double &az );

  
   // This routine rotates a 3x3 system given a transformation matrix.  The 
   // matrix can be rotated in place by sending same variable in & out, or a 
   // new matrix can be created.  In the 4x4 case, the translation portion
   // of the matrix is unchanged.
   void rotate_system( double mtx[3][3], double sys_in[3][3], 
                       double sys_out[3][3] );
   void rotate_system( double mtx[4][4], double sys_in[4][4], 
                       double sys_out[4][4] );

   double det_mtx( double m[3][3] );

   // Multiply matrices together.  If any input is NULL, the identity
   // matrix is used in its place.
   void mult_mtx( double a[3][3],double b[3][3], double d[3][3] );
   void mult_mtx( double a[4][4], double b[4][4], double d[4][4] );

   CubitStatus inv_mtx_adj( double mtx[3][3], double inv_mtx[3][3] );

   CubitStatus inv_trans_mtx( double transf[4][4], double inv_transf[4][4] );



   void vecs_to_mtx( double xvec[3], double yvec[3], double zvec[3], 
                     double matrix[3][3] );
   void vecs_to_mtx( double xvec[3], double yvec[3], double zvec[3], 
                     double origin[3], double matrix[4][4] );


   void print_mtx( double mtx[3][3] );
   void print_mtx( double mtx[4][4] );

   //**************************************************************************
   // 3D Points
   //**************************************************************************
   void copy_pnt( double pnt_in[3], double pnt_out[3] );


   void copy_pnt( double pnt_in[3], CubitVector &cubit_vec );
   void copy_pnt( CubitVector &cubit_vec, double pnt_out[3] );

   void set_pnt( double pnt[3], double x, double y, double z );

   CubitBoolean pnt_eq( double pnt1[3],double pnt2[3] );

   CubitStatus mirror_pnt( double pnt[3], double pln_orig[3], 
                           double pln_norm[3], double m_pnt[3]);

   CubitStatus next_pnt( double str_pnt[3], double vec_dir[3], double len,
                         double new_pnt[3]);

   //***************************************************************************
   // 3D Vectors
   //***************************************************************************
   CubitStatus unit_vec( double vin[3], double vout[3] );

   double dot_vec( double uval[3], double vval[3] );

   void cross_vec( double uval[3], double vval[3], double cross[3] );

   void cross_vec_unit( double uval[3], double vval[3], double cross[3] );

   void orth_vecs( double uvect[3], double vvect[3], double wvect[3] );

   double mag_vec( double vec[3] );

   CubitStatus get_vec ( double str_pnt[3], double stp_pnt[3],
                         double vector_out[3] );

   CubitStatus get_vec_unit( double str_pnt[3], double stp_pnt[3],
                             double uv_out[3] );

   void mult_vecxconst( double constant, double vec[3], double vec_out[3] );

   void reverse_vec( double vin[3],double vout[3] );

   double angle_vec_vec( double v1[3],double v2[3] );
   
   //***************************************************************************
   // Distances
   //***************************************************************************
   double dist_pnt_pnt( double pnt1[3], double pnt2[3] );

   double dist_pln_pln( double pln_1_orig[3], double pln_1_norm[3], 
                        double pln_2_orig[3], double pln_2_norm[3],
                        AgtSide *side = NULL, AgtOrientation *orien = NULL,
                        unsigned short *status = NULL );
   
   //***************************************************************************
   // Intersections
   //***************************************************************************
   CubitStatus int_ln_pln( double ln_orig[3], double ln_vec[3], double pln_orig[3],
                           double pln_norm[3], double int_pnt[3] );
   int int_ln_ln( double p1[3], double v1[3], double p2[3], double v2[3], 
                  double int_pnt1[3], double int_pnt2[3] );

   int int_pnt_ln( double pnt[3], double ln_orig[3],  double ln_vec[3],
                            double int_pnt[3] );

   int int_pnt_pln( double pnt[3], double pln_orig[3], double pln_norm[3], 
                             double pln_int[3] );

   CubitStatus int_pln_pln( double pln_1_orig[3], double pln_1_norm[3],
                            double pln_2_orig[3], double pln_2_norm[3],
                            double ln_orig[3], double ln_vec[3] );
   
   //**************************************************************************
   // Comparison/Containment Tests
   //**************************************************************************
   CubitBoolean is_vec_par( double vec_1[3], double vec_2[3] );

   CubitBoolean is_vec_perp( double vec_1[3],double vec_2[3]);

   CubitBoolean is_vecs_same_dir( double vec_1[3], double vec_2[3] );

   CubitBoolean is_pnt_on_ln( double pnt[3], double ln_orig[3], double ln_vec[3] );

   CubitBoolean is_pnt_on_ln_seg( double pnt1[3], double end1[3], double end2[3] );

   CubitBoolean is_pnt_on_pln( double pnt[3], double pln_orig[3], double pln_norm[3] );

   CubitBoolean is_ln_on_pln( double ln_orig[3], double ln_vec[3],
                              double pln_orig[3], double pln_norm[3] );


   CubitBoolean is_pln_on_pln( double pln_orig1[3], double pln_norm1[3],
                               double pln_orig2[3], double pln_norm2[3] );

   //**************************************************************************
   // Arcs/Circles
   //**************************************************************************






   void setup_arc( double vec_1[3], double vec_2[3], double origin[3],
                  double start_angle, double end_angle, // Angles in radians
                  double radius, AGT_3D_Arc &arc );
   void setup_arc( CubitVector& vec_1, CubitVector& vec_2, 
                   CubitVector origin, double start_angle, // Angles in radians
                   double end_angle, double radius,
                   AGT_3D_Arc &arc );
  


   void get_arc_xyz( AGT_3D_Arc &arc, double param, double pnt[3] );
   void get_arc_xyz( AGT_3D_Arc &arc, double param, CubitVector& pnt );

   int get_num_circle_tess_pnts( double radius, double len_tolerance = 1e-4 );

   //**************************************************************************
   // Miscellaneous
   //**************************************************************************

   void get_pln_orig_norm( double A, double B, double C, double D, 
                           double pln_orig[3], double pln_norm[3] = NULL );
   void get_box_corners( double box_min[3],double box_max[3],double c[8][3] );











   CubitStatus min_pln_box_int_corners( const CubitPlane& plane, 
                                        const CubitBox& box,
                                        int extension_type, double extension,
                                        CubitVector& p1, CubitVector& p2,
                                        CubitVector& p3, CubitVector& p4,
                                        CubitBoolean silent = CUBIT_FALSE );
   CubitStatus min_pln_box_int_corners( CubitVector& vec1,
                                        CubitVector& vec2,
                                        CubitVector& vec3,
                                        CubitVector& box_min, CubitVector& box_max,
                                        int extension_type, double extension,
                                        CubitVector& p1, CubitVector& p2, 
                                        CubitVector& p3, CubitVector& p4,
                                        CubitBoolean silent = CUBIT_FALSE );
   CubitStatus min_pln_box_int_corners( double plane_norm[3], double plane_coeff,
                                       double box_min[3], double box_max[3],
                                       int extension_type, double extension,
                                       double p1[3], double p2[3],   // OUT
                                       double p3[3], double p4[3],   // OUT
                                       CubitBoolean silent = CUBIT_FALSE); 



   CubitStatus get_tight_bounding_box( DLIList<CubitVector*> &point_list,                                       
                                       CubitVector& p1, CubitVector& p2,
                                       CubitVector& p3, CubitVector& p4,
                                       CubitVector& p5, CubitVector& p6,
                                       CubitVector& p7, CubitVector& p8);
   CubitStatus get_tight_bounding_box( DLIList<CubitVector*> &point_list,                              
                                       CubitVector &center,
                                       CubitVector axes[3],
                                       CubitVector &extension );

  










   CubitStatus min_cyl_box_int( double radius, CubitVector& axis, CubitVector& center,
                                CubitBox& box,
                                int extension_type, double extension,
                                CubitVector &start, CubitVector &end,
                                int num_tess_pnts = 2048 );
   CubitStatus min_cyl_box_int( double radius, double axis[3], double center[3], 
                                double box_min[3], double box_max[3], 
                                int extension_type, double extension, 
                                double start[3], double end[3],
                                int num_tess_pnts = 2048 );

   double MinTriangleTriangle (Triangle3& tri0, 
                               Triangle3& tri1, 
                              double& s, double& t, double& u, double& v);
   
   double AreaIntersection (Triangle &tri1, Triangle &tri2 );

   double Angle( Triangle3& tri1, Triangle3& tri2 );

   double MinPointTriangle (const Point3& p, const Triangle3& tri, double& s, double& t);

   void Normal( Triangle3& tri, double normal[3] );

   double ProjectedOverlap( Triangle3& tri1, Triangle3& tri2, bool draw_overlap = false );

protected:
   AnalyticGeometryTool();
   //- Class Constructor. (Not callable by user code. Class is constructed
   //- by the {instance()} member function.
   
private:
   
   static AnalyticGeometryTool* instance_;

   double agtEpsilon;  // default = 1e-6


//                               MAGIC SOFTWARE
//This code was obtained from Dave Eberly, at www.magic-software.com.  It
//has been modified only to be placed into AnalyticGeometryTool.  This was
//done for convenience only.  An alternate solution would be to create a 
//separate library for these functions.
//
//Steve Storm, [email protected], 3-27-99
//MAGIC is an acronym for My Alternate Graphics and Image Code. The
//initial code base originated in 1991 as an attempt to answer questions that
//arise in my favorite news group, comp.graphics.algorithms, and has been
//growing ever since. Magic is intended to provide free source code for solving
//problems that commonly arise in computer graphics and image analysis. While
//the code at this web site is free, additional conditions are: 
//
// *  You may distribute the original source code to others at no charge. You
//    got it for free, so do not charge others for it. 
// *  You may modify the original source code and distribute it to others at no
//    charge. The modified code must be documented to indicate that it is not
//    part of the original package. I do not want folks to blame me for bugs
//    introduced by modifications. I do accept blame for bugs that are my
//    doing and will do my best to fix them. 
// *  You may use this code for non-commercial purposes. You may also
//    incorporate this code into commercial packages. However, you may not
//    sell any of your source code which contains my original and/or modified
//    source code (see items 1 and 2 in this list of conditions). In such a case,
//    you would need to factor out my code and freely distribute it. Send me
//    email if you use it. I am always interested in knowing where and how
//    my code is being used. 
// *  The original code comes with absolutely no warranty and no guarantee is
//    made that the code is bug-free. Caveat emptor. 
//
// Dave Eberly, [email protected], www.magic-software.com
//FILE: minbox2.h
   double Area (int N, Point2* pt, double angle);
   void MinimalBoxForAngle (int N, Point2* pt, double angle, OBBox2& box);
   OBBox2 MinimalBox2 (int N, Point2* pt);
   // Functions to find minimal area rectangle that surrounds a set of points.

#if 1
//FILE: minbox3.h
   void MatrixToAngleAxis( double** R, double& angle, double axis[3] );
   void AngleAxisToMatrix( double angle, double axis[3], double R[3][3] );
   double Volume( int N, Point3* pt, double angle[3] );
   void MinimalBoxForAngles( int N, Point3* pt, double angle[3], OBBox3& box );
   void GetInterval( double A[3], double D[3], double& tmin, double& tmax );
   void Combine( double result[3], double A[3], double t, double D[3] );
   double MinimizeOnInterval( int N, Point3* pt, double A[3], double D[3] );
   double MinimizeOnLattice( int N, Point3* pt, double A[3], int layers,
                             double thickness );
   void InitialGuess( int N, Point3* pt, double angle[3] );
   OBBox3 MinimalBox3( int N, Point3* pt );
   // Functions to find minimal volume box that surrounds a set of points.
#endif

   double Abs (double x);
   double ACos (double x);
   double ATan2 (double y, double x);
   double Cos (double x);
   double Sign (double x);
   double Sin (double x);
   double Sqrt (double x);
   double UnitRandom ();
   double Dot (const Point2& p, const Point2& q);
   double Length (const Point2& p);
   CubitBoolean Unitize (Point2& p );
   double Dot (const Point3& p, const Point3& q);
   Point3 Cross (const Point3& p, const Point3& q);
   double Length (const Point3& p);
   CubitBoolean Unitize (Point3& p );
   // Supporting code for triangle calculations

   double MinLineSegmentLineSegment (const Line3& seg0, const Line3& seg1, double& s, double& t);

   //get intersection between bounding box and plane
   int get_plane_bbox_intersections( double box_min[3],
                                     double box_max[3],
                                     double pln_orig[3],
                                     double pln_norm[3],
                                     double int_array[6][3] );
// FILE: pt3tri3.cpp


// FILE: lin3tri3.cpp
   double MinLineSegmentTriangle (const Line3& seg, const Triangle3& tri, double& r, double& s, double& t);
   // Code to support distance between triangles

//FILE: triasect.cpp
   AgtLine EdgeToLine (Point2* v0, Point2* v1);
   void TriangleLines (Triangle* T, AgtLine line[3]);
   AgtTriList* SplitAndDecompose (Triangle* T, AgtLine* line);
   AgtTriList* Intersection (Triangle* T0, Triangle* T1);
   double AreaTriangle (Triangle* T);
   double Area (AgtTriList* list);
   double Area (Triangle &tri1, Triangle &tri2 );
   // Code to support area of overlap of two triangles
};

#if 1
// FILE: eigen.h
class mgcEigen
{
public:
        mgcEigen (int _size);
        ~mgcEigen ();

        mgcEigen& Matrix (float** inmat);

private:
        int size;
        float** mat;
        float* diag;
        float* subd;

        // Householder reduction to tridiagonal form
        void Tridiagonal2 (float** pmat, float* pdiag, float* psubd);
        void Tridiagonal3 (float** pmat, float* pdiag, float* psubd);
        void Tridiagonal4 (float** pmat, float* pdiag, float* psubd);
        void TridiagonalN (int n, float** mat, float* diag, float* subd);

        // QL algorithm with implicit shifting, applies to tridiagonal matrices
        void QLAlgorithm (int n, float* pdiag, float* psubd, float** pmat);

        // sort eigenvalues from largest to smallest
        void DecreasingSort (int n, float* eigval, float** eigvec);

        // sort eigenvalues from smallest to largest
        void IncreasingSort (int n, float* eigval, float** eigvec);

// error handling
public:
        static int verbose1;
        static unsigned error;
        static void Report (std::ostream& ostr);
private:
        static const unsigned invalid_size;
        static const unsigned allocation_failed;
        static const unsigned ql_exceeded;
        static const char* message[3];
        static int Number (unsigned single_error);
        static void Report (unsigned single_error);
};

// FILE: eigen.h
class mgcEigenD
{
public:
        mgcEigenD (int _size);
        ~mgcEigenD ();

    // set the matrix for eigensolving
        double& Matrix (int row, int col) { return mat[row][col]; }
        mgcEigenD& Matrix (double** inmat);

        // get the results of eigensolving
        double Eigenvector (int row, int col) { return mat[row][col]; }
        const double** Eigenvector () { return (const double**) mat; }

        void EigenStuff3 ();  // uses TriDiagonal3

private:
        int size;
        double** mat;
        double* diag;
        double* subd;

        // Householder reduction to tridiagonal form
        void Tridiagonal2 (double** mat, double* diag, double* subd);
        void Tridiagonal3 (double** mat, double* diag, double* subd);
        void Tridiagonal4 (double** mat, double* diag, double* subd);
        void TridiagonalN (int n, double** mat, double* diag, double* subd);

        // QL algorithm with implicit shifting, applies to tridiagonal matrices
        void QLAlgorithm (int n, double* diag, double* subd, double** mat);

        // sort eigenvalues from largest to smallest
        void DecreasingSort (int n, double* eigval, double** eigvec);

        // sort eigenvalues from smallest to largest
        void IncreasingSort (int n, double* eigval, double** eigvec);

// error handling
public:
        static int verbose1;
        static unsigned error;
        static void Report (std::ostream& ostr);
private:
        static const unsigned invalid_size;
        static const unsigned allocation_failed;
        static const unsigned ql_exceeded;
        static const char* message[3];
        static int Number (unsigned single_error);
        static void Report (unsigned single_error);
};
#endif

//---------------------------------------------------------------------------
// Wrapped math functions
//---------------------------------------------------------------------------
inline double AnalyticGeometryTool::Abs (double x)
{
    return double(fabs(x));
}

inline double AnalyticGeometryTool::ATan2 (double y, double x)
{
    return double(atan2(y,x));
}

inline double AnalyticGeometryTool::Sqrt (double x)
{
    return double(sqrt(x));
}

inline double AnalyticGeometryTool::UnitRandom ()
{
    return double(rand())/double(RAND_MAX);
}

//---------------------------------------------------------------------------
// Points in 2D
//---------------------------------------------------------------------------

inline double AnalyticGeometryTool::Dot (const Point3& p, const Point3& q)
{
    return double(p.x*q.x + p.y*q.y + p.z*q.z);
}


inline double AnalyticGeometryTool::set_epsilon( double new_epsilon )
{ 
   double old_epsilon = agtEpsilon;
   agtEpsilon = new_epsilon;
   return old_epsilon;
}

inline CubitBoolean AnalyticGeometryTool::dbl_eq(double val_1,double val_2)
{ 
   CubitBoolean result;
   if (fabs(val_1 - val_2) < agtEpsilon)
      result = CUBIT_TRUE;
   else
      result = CUBIT_FALSE;
   return result;
}

#endif