#include <math.h>
#include <cassert>

#ifndef _geom_include
#define _geom_include

extern long ccwcount;

inline double ccw(double x1,double y1,double x2,double y2,double x3,double y3)
{
  ++ccwcount;
  return(x1-x3)*(y2-y3)-(x2-x3)*(y1-y3);
}

inline double ccw(const double*p1,const double*p2,const double*p3)
{
  return ccw(p1[0],p1[1],p2[0],p2[1],p3[0],p3[1]);
}

inline double sqrtdist(double p1[2], double p2[2])
{
   return ((p2[0]-p1[0])*(p2[0]-p1[0]))+((p2[1]-p1[1])*(p2[1]-p1[1]));
}

inline double dist(double p1[2], double p2[2])
{
  return sqrt(sqrtdist(p1, p2));
}

inline double sqrsum(double x,double y)
{
  return x*x+y*y;
}

inline int intersect(double e11[2], double e12[2], double e21[2], double e22[2])
{
// 	printf(" ccw1 = %lf, ccw2 = %lf\n", ccw (e11, e22, e21), ccw (e21, e22, e12));
  	if (ccw (e11, e22, e21) > 0 && ccw (e21, e22, e12) > 0)
  		return 1;
  	else
  		return 0;
}

inline double linear_normalized_length(double p1[2], double h1, double p2[2], double h2)
{
   assert(0.0<h1);
   assert(0.0<h2);

   double d = dist(p1, p2);
   double dh = h2 - h1;
   double lenght;

   if(h1/h2 < 0.999 || h2/h1 < 0.999)
   {
	  double dln = log(h2) - log(h1);
	  lenght = (d*dln)/dh;
   }
   else lenght = d/h1;

   return lenght;
}

// inline double geom_normalized_length(double p1[2], double h1, double p2[2], double h2)
// {
//    assert(0.0<h1);
//    assert(0.0<h2);
// 
//    double d = dist(p1, p2);
//    double dh = h2 - h1;
//    double lenght;
// 
//    /** geometric */
//    if(fabs(dh) <= 0.001) lenght = d/a->h;
//    else
//    {
// 	  double dln = log(b->h) - log(a->h);
// 	  lenght = (d*dh)/(a->h*b->h*dln);
//    }
// 
//    return lenght;
// }

inline double get_abs_curv(double l, double lt, double h1, double h2)
{
	double abs;
	double dh = h2-h1;
	if (fabs(dh) < 0.001)
		abs = (l/lt)*h1;
	else
		abs = ((exp(l*dh/lt)*h1)-h1)/dh;
	return abs;
}

#endif