#ifndef GR_GRCURVEGEOM_H
#define GR_GRCURVEGEOM_H 1

#include "CubitDefines.h"
#include "CubitVector.hpp"
#include "DLIList.hpp"
#include "GR_config.h"
#include "GR_misc.h"

class CubitBox;

class GRCurveGeom {

 public:

  typedef enum curve_type {UNKNOWN, LINE, CUBIC, ARC, BEZIER, SPLINE};

 protected:

  int n_pts;           //Number of points defining the curve
  CubitVector* points; //Points locations
  int n_dim;           //Curve geometric dimension

  //Functor used to compare parameters on a curve. Two params
  //that are almost equal (to a small tolerance) are considered equal.
  struct param_compare {
    bool operator()(const double& d1, const double& d2) const {
      if(iFuzzyComp(d1, d2) == -1) { return true; }
      else { return false; }      
    }
  };
   
 private:

  //Disabling copy construction and operator=.
  GRCurveGeom(const GRCurveGeom& C) 
    { assert(0); n_pts = 0; points = NULL; }  
  GRCurveGeom& operator=(const GRCurveGeom& C) 
    { assert(0); n_pts = 0; points = NULL; return *this; }
 
 public:
  
  GRCurveGeom() : n_pts(0), points(NULL), n_dim(0) {}
  virtual ~GRCurveGeom() {}

  //Returns the curve number of topological dimensions.
  int topo_dim() const { return 1; }

  //Returns the curve number of geometric dimensions.
  virtual int geom_dim() const { return n_dim; }

  //Returns the type of geometric entity.
  virtual curve_type get_curve_type() const = 0;

  //Returns the min and max param for this curve;
  virtual double min_param() const = 0;
  virtual double max_param() const = 0;

  //Calculates the bounding box for this curve.
  virtual CubitBox bounding_box() const = 0;

  //Returns parameter at a given coordinate.
  virtual double param_at_coord(const CubitVector& coord) const = 0;

  //Gets the parameter located arc length "length" from "param_start" 
  virtual double param_at_arc_length(const double param_start, const double length) const = 0;

  //Returns coordinate at a given parameter.
  virtual void coord_at_param(const double param, CubitVector& coord) const = 0;

  //Given a coordinate near the curve, computes neareast coord on the curve.
  virtual void closest_coord_on_curve(const CubitVector& coord,
				      CubitVector& closest) const = 0;

  //Given a coordinate near the curve, computes neareast parameter on the curve.
  virtual double closest_param_on_curve(const CubitVector& coord) const = 0;

  //Returns true if coord is directly on the curve.
  virtual bool coord_on_curve(const CubitVector& coord,
			      double* const param_X = NULL,
			      double* const param_Y = NULL,
			      double* const param_Z = NULL) const = 0; 

  //Returns the first and second derivative with respect to the parameter
  virtual void first_deriv(const double param, CubitVector& FD) const = 0;
  virtual void first_deriv(const CubitVector& coord, CubitVector& FD) const = 0;

  virtual void second_deriv(const double param, CubitVector& SD) const = 0;
  virtual void second_deriv(const CubitVector& coord, CubitVector& SD) const = 0;

  //Curvature at param and at coord (located on the curve).
  virtual double curvature(const double param) const = 0;
  virtual double curvature(const CubitVector& coord) const = 0;

  //Curvature vector
  virtual void curvature(const double param, CubitVector& curvature) const = 0;
  virtual void curvature(const CubitVector& coord, CubitVector& Curvature) const = 0;

  //Unit normal points to the left side of the curve.
  //ONLY WORKS IN 2D!!!
  virtual void unit_normal(const double param, CubitVector& normal) const = 0;
  virtual void unit_normal(const CubitVector& coord, CubitVector& normal) const = 0;
  
  //Tangent vector has the same "sense" as the curve.
  virtual void unit_tangent(const double param, CubitVector& tangent) const = 0;
  virtual void unit_tangent(const CubitVector& coord, CubitVector& tangent) const = 0;

  //Total length, length between dParam1 and dParam2, between adCoord1 and adCoord2.
  virtual double arc_length() const = 0;
  virtual double arc_length(const double param1, const double param2) const = 0;
  virtual double arc_length(const CubitVector& coord1, const CubitVector& coord2) const = 0;

  //Rerturns the center point based on arc length. Either based on total
  //arc length or between params or points.
  virtual void center_point(CubitVector& mid_point) const = 0;
  virtual void center_point(const double param1, const double param2,
			    CubitVector& mid_point) const = 0;
  virtual void center_point(const CubitVector& coord1, const CubitVector& coord2,
			    CubitVector& mid_point) const = 0;

  //Rerturns the point located in the middle of the parameterization range.
  void mid_point(CubitVector& mp) const;
  void mid_point(const double param1, const double param2,
		 CubitVector& mp) const;
  void mid_point(const CubitVector& coord1, const CubitVector& coord2,
		 CubitVector& mp) const;

  //Total variation of the tangent angle
  virtual double TVT() const = 0;
  virtual double TVT(const double param1, const double param2) const = 0;  
  virtual double TVT(const CubitVector& coord1, const CubitVector& coord2) const = 0;

  //Returns the local interior extrema for the curve.
  virtual void interior_extrema(DLIList<CubitVector*>& point_list, 
				CubitSense return_sense) const = 0;

  //Gets points[pt] in some form or another.
  virtual double X(const int pt) const = 0;
  virtual double Y(const int pt) const = 0;
  virtual double Z(const int pt) const = 0;
  virtual void get_coord(const int pt, CubitVector& coord) const = 0;
  virtual const CubitVector& get_coord(const int pt) const = 0;

  //Sets Points[iPt] to new location
  virtual void set_X(const int pt, const double X) = 0;
  virtual void set_Y(const int pt, const double Y) = 0;
  virtual void set_Z(const int pt, const double Z) = 0;
  virtual void set_coord(const int pt, const CubitVector& new_coord) = 0; 

  virtual bool qValid() const = 0;  

};



class GRLine : public GRCurveGeom {
 
 public:

  GRLine();
  GRLine(const CubitVector& coord_start, const CubitVector& coord_end);
  GRLine(const GRLine& L);
  GRLine& operator=(const GRLine& L);
  virtual ~GRLine();

  GRCurveGeom::curve_type get_curve_type() const { return GRCurveGeom::LINE; }

  double min_param() const { return 0.; }
  double max_param() const { return 1.; }
 
  CubitBox bounding_box() const;
  
  double param_at_coord(const CubitVector& coord) const;
  
  double param_at_arc_length(const double param_start, const double length) const;

  void coord_at_param(const double param, CubitVector& coord) const;
 
  void closest_coord_on_curve(const CubitVector& coord,
			      CubitVector& closest) const;
  
  double closest_param_on_curve(const CubitVector& coord) const;

  bool coord_on_curve(const CubitVector& coord,
		      double* const param_X = NULL,
		      double* const param_Y = NULL,
		      double* const param_Z = NULL) const; 
  
  void first_deriv(const double param, CubitVector& FD) const;
  void first_deriv(const CubitVector& coord, CubitVector& FD) const;

  void second_deriv(const double param, CubitVector& SD) const;
  void second_deriv(const CubitVector& coord, CubitVector& SD) const;

  double curvature(const double param) const;
  double curvature(const CubitVector& coord) const;
  void curvature(const double param, CubitVector& curvature) const;
  void curvature(const CubitVector& coord, CubitVector& Curvature) const;
  
  void unit_normal(const double param, CubitVector& normal) const;
  void unit_normal(const CubitVector& coord, CubitVector& normal) const;
    
  void unit_tangent(const double param, CubitVector& tangent) const;
  void unit_tangent(const CubitVector& coord, CubitVector& tangent) const;
  
  double arc_length() const;
  double arc_length(const double param1, const double param2) const;
  double arc_length(const CubitVector& coord1, const CubitVector& coord2) const;

  void center_point(CubitVector& mid_point) const;
  void center_point(const double param1, const double param2,
		    CubitVector& mid_point) const;
  void center_point(const CubitVector& coord1, const CubitVector& coord2,
		    CubitVector& mid_point) const;
  
  void mid_point(CubitVector& mp) const;
  void mid_point(const double param1, const double param2,
		 CubitVector& mp) const;
  void mid_point(const CubitVector& coord1, const CubitVector& coord2,
		 CubitVector& mp) const;
  
  double TVT() const;
  double TVT(const double param1, const double param2) const;  
  double TVT(const CubitVector& coord1, const CubitVector& coord2) const;
  
  void interior_extrema(DLIList<CubitVector*>& point_list, 
			CubitSense return_sense) const;
  
  double X(const int pt) const;
  double Y(const int pt) const;
  double Z(const int pt) const;
  void get_coord(const int pt, CubitVector& coord) const;
  const CubitVector& get_coord(const int pt) const;

  //Sets Points[iPt] to new location
  void set_X(const int pt, const double X);
  void set_Y(const int pt, const double Y);
  void set_Z(const int pt, const double Z);
  void set_coord(const int pt, const CubitVector& new_coord); 

  bool qValid() const;  

};



class GRArc {
  
 private:

/*   //Disabling copy construction and operator=. */
/*   GRCurveGeom(const GRCurveGeom& C)  */
/*     { assert(0); n_pts = 0; points = NULL; }   */
/*   GRCurveGeom& operator=(const GRCurveGeom& C)  */
/*     { assert(0); n_pts = 0; points = NULL; return *this; } */
 
 public:
  
/*   GRCurveGeom() : n_pts(0), points(NULL), n_dim(0) {} */
/*   virtual ~GRCurveGeom() {} */

/*   //Returns the curve number of topological dimensions. */
/*   int topo_dim() const { return 1; } */

/*   //Returns the curve number of geometric dimensions. */
/*   virtual int geom_dim() const = 0; */

/*   //Returns the type of geometric entity. */
/*   virtual curve_type get_curve_type() const = 0; */

/*   //Returns the min and max param for this curve; */
/*   virtual double min_param() const = 0; */
/*   virtual double max_param() const = 0; */

/*   //Calculates the bounding box for this curve. */
/*   virtual CubitBox bounding_box() const = 0; */

/*   //Returns parameter at a given coordinate. */
/*   virtual double param_at_coord(const CubitVector& coord) const = 0; */

/*   //Gets the parameter located at arc_length from param_start  */
/*   virtual double param_at_arc_length(const double param_start, const double arc_length) const = 0; */

/*   //Returns coordinate at a given parameter. */
/*   virtual void coord_at_param(const double param, CubitVector& coord) const = 0; */

/*   //Given a coordinate near the curve, computes neareast coord on the curve. */
/*   virtual void closest_coord_on_curve(const CubitVector& coord, */
/* 				      CubitVector& closest) const = 0; */

/*   //Given a coordinate near the curve, computes neareast parameter on the curve. */
/*   virtual double closest_param_on_curve(const CubitVector& coord) const = 0; */

/*   //Returns true if coord is directly on the curve. */
/*   virtual bool coord_on_curve(const CubitVector& coord, */
/* 			      double* const param_X = NULL, */
/* 			      double* const param_Y = NULL, */
/* 			      double* const param_Z = NULL) const = 0;  */

/*   //Returns the first and second derivative with respect to the parameter */
/*   virtual void first_deriv(const double param, CubitVector& FD) const = 0; */
/*   virtual void first_deriv(const CubitVector& coord, CubitVector& FD) const = 0; */

/*   virtual void second_deriv(const double param, CubitVector& SD) const = 0; */
/*   virtual void second_deriv(const CubitVector& coord, CubitVector& SD) const = 0; */

/*   //Curvature at param and at coord (located on the curve). */
/*   virtual double curvature(const double param) const = 0; */
/*   virtual double curvature(const CubitVector& coord) const = 0; */

/*   //Curvature vector */
/*   virtual void curvature(const double param, CubitVector& curvature) const = 0; */
/*   virtual void curvature(const CubitVector& coord, CubitVector& Curvature) const = 0; */

/*   //Unit normal points to the left side of the curve. */
/*   //ONLY WORKS IN 2D!!! */
/*   virtual void unit_normal(const double param, CubitVector& normal) const = 0; */
/*   virtual void unit_normal(const CubitVector& coord, CubitVector& normal) const = 0; */
  
/*   //Tangent vector has the same "sense" as the curve. */
/*   virtual void unit_tangent(const double param, CubitVector& tangent) const = 0; */
/*   virtual void unit_tangent(const CubitVector& coord, CubitVector& tangent) const = 0; */

/*   //Total length, length between dParam1 and dParam2, between adCoord1 and adCoord2. */
/*   virtual double arc_length() const = 0; */
/*   virtual double arc_length(const double param1, const double param2) const = 0; */
/*   virtual double arc_length(const CubitVector& coord1, const CubitVector& coord2) const = 0; */

/*   //Rerturns the center point based on arc length. Either based on total */
/*   //arc length or between params or points. */
/*   virtual void center_point(CubitVector& mid_point) const = 0; */
/*   virtual void center_point(const double param1, const double param2, */
/* 			    CubitVector& mid_point) const = 0; */
/*   virtual void center_point(const CubitVector& coord1, const CubitVector& coord2, */
/* 			    CubitVector& mid_point) const = 0; */

/*   //Total variation of the tangent angle */
/*   virtual double TVT() const = 0; */
/*   virtual double TVT(const double param1, const double param2) const = 0;   */
/*   virtual double TVT(const CubitVector& coord1, const CubitVector& coord2) const = 0; */

/*   //Returns the local interior extrema for the curve. */
/*   virtual void interior_extrema(DLIList<CubitVector*>& point_list,  */
/* 				CubitSense return_sense) const = 0; */

/*   //Gets points[pt] in some form or another. */
/*   virtual double X(const int pt) const = 0; */
/*   virtual double Y(const int pt) const = 0; */
/*   virtual double Z(const int pt) const = 0; */
/*   virtual void get_coord(const int pt, CubitVector& coord) const = 0; */
/*   virtual const CubitVector& get_coord(const int pt) const = 0; */

/*   //Sets Points[iPt] to new location */
/*   virtual void set_X(const int pt, const double X) = 0; */
/*   virtual void set_Y(const int pt, const double Y) = 0; */
/*   virtual void set_Z(const int pt, const double Z) = 0; */
/*   virtual void SetCoord(const int pt, const CubitVector& new_coord) = 0;  */

  void get_center(CubitVector& center) { assert(0); }
  double get_radius() { assert(0); return -1.; }

/*   virtual bool qValid() const = 0;   */

};

///INLINE DEFINITIONS FOR BASE CLASS///

inline void GRCurveGeom::mid_point(CubitVector& mp) const {
  assert(qValid()); coord_at_param(0.5 * (max_param() + min_param()), mp);
}

inline void GRCurveGeom::mid_point(const double param1, 
				   const double param2, 
				   CubitVector& mp) const {
  assert(qValid());
  assert(iFuzzyComp(param1, min_param()) >= 0 &&
	 iFuzzyComp(param1, max_param()) <= 0);
  assert(iFuzzyComp(param2, min_param()) >= 0 &&
	 iFuzzyComp(param2, max_param()) <= 0);
  coord_at_param(0.5 * (param1 + param2), mp);
}

inline void GRCurveGeom::mid_point(const CubitVector& coord1, 
				   const CubitVector& coord2, 
				   CubitVector& mp) const {
  assert(qValid());
  assert(coord_on_curve(coord1));
  assert(coord_on_curve(coord2));
  coord_at_param(0.5 * (param_at_coord(coord1) + param_at_coord(coord2)), mp);
}

///INLINE DEFINITIONS FOR GRLine CLASS///

inline void GRLine::second_deriv(const double param,
				 CubitVector& SD) const { SD.set(0., 0., 0.); }
inline void GRLine::second_deriv(const CubitVector& coord,
				 CubitVector& SD) const { SD.set(0., 0., 0.); }

inline double GRLine::curvature(const double param) const { return 0.; }
inline double GRLine::curvature(const CubitVector& coord) const { return 0.; }

inline void GRLine::curvature(const double param,
			       CubitVector& curvature) const { curvature.set(0., 0., 0.); }

inline void GRLine::curvature(const CubitVector& coord,
			       CubitVector& curvature) const { curvature.set(0., 0., 0.); }

inline double GRLine::TVT() const { return 0.; }
inline double GRLine::TVT(const double param1,
			  const double param2) const { return 0.; }
inline double GRLine::TVT(const CubitVector& coord1,
			  const CubitVector& coord2) const { return 0.; }

//No interior extremum for a line.
inline void GRLine::interior_extrema(DLIList<CubitVector*>& point_list, 
				     CubitSense return_sense) const { point_list.clean_out(); }

inline double GRLine::X(const int pt) const { 
  assert(pt >= 0 && pt < n_pts);
  return points[pt].x();
}
inline double GRLine::Y(const int pt) const {
  assert(pt >= 0 && pt < n_pts);
  return points[pt].y();
}
inline double GRLine::Z(const int pt) const {
  assert(pt >= 0 && pt < n_pts);
  return points[pt].z();
}
inline void GRLine::get_coord(const int pt,
			      CubitVector& coord) const {
  assert(pt >= 0 && pt < n_pts);
  coord.set(points[pt].x(), points[pt].y(), points[pt].z());  
}
inline const CubitVector& GRLine::get_coord(const int pt) const {
  assert(pt >= 0 && pt < n_pts);
  return points[pt];
}

inline void GRLine::set_X(const int pt,
			  const double X) {
  assert(n_pts == 2);
  assert(points != NULL);
  assert(pt >= 0 && pt < n_pts);
  
  points[pt].x(X);

  if(iFuzzyComp(points[0].x(), points[1].x()) == 0 && 
     iFuzzyComp(points[0].y(), points[1].y()) == 0 &&
     iFuzzyComp(points[0].z(), points[1].z()) == 0) {
    vFatalError("Both points of a line are at the same coordinate.",
		"GRLine::set_x(const int, const double)");    
  }
}
inline void GRLine::set_Y(const int pt,
			  const double Y) {
  assert(n_pts == 2);
  assert(points != NULL);
  assert(pt >= 0 && pt < n_pts);
  
  points[pt].y(Y);

  if(iFuzzyComp(points[0].x(), points[1].x()) == 0 && 
     iFuzzyComp(points[0].y(), points[1].y()) == 0 &&
     iFuzzyComp(points[0].z(), points[1].z()) == 0) {
    vFatalError("Both points of a line are at the same coordinate.",
		"GRLine::set_y(const int, const double)");    
  }
}
inline void GRLine::set_Z(const int pt,
			  const double Z) {
  assert(n_pts == 2);
  assert(points != NULL);
  assert(pt >= 0 && pt < n_pts);

  if(n_dim == 2) { assert(0); }
  points[pt].z(Z);

  if(iFuzzyComp(points[0].x(), points[1].x()) == 0 && 
     iFuzzyComp(points[0].y(), points[1].y()) == 0 &&
     iFuzzyComp(points[0].z(), points[1].z()) == 0) {
    vFatalError("Both points of a line are at the same coordinate.",
		"GRLine::vSetY(const int, const double)");    
  }
}
inline void GRLine::set_coord(const int pt,
			      const CubitVector& new_coord) {
  assert(n_pts == 2);
  assert(points != NULL);
  assert(pt >= 0 && pt < n_pts);
  if(n_dim == 2) { assert(iFuzzyComp(new_coord.z(), 0.) == 0); }

  points[pt].set(new_coord.x(), new_coord.y(), new_coord.z());

  if(iFuzzyComp(points[0].x(), points[1].x()) == 0 && 
     iFuzzyComp(points[0].y(), points[1].y()) == 0 &&
     iFuzzyComp(points[0].z(), points[1].z()) == 0) {
    vFatalError("Both points of a line are at the same coordinate.",
		"GRLine::set_coord(const int, const CubitVector&).");    
  }
}



#endif 
