/*******************************************************************************
+
+  LEDA 3.5
+
+  _circle.c
+
+  This file is part of the LEDA research version (LEDA-R) that can be 
+  used free of charge in academic research and teaching. Any commercial
+  use of this software requires a license which is distributed by the
+  LEDA Software GmbH, Postfach 151101, 66041 Saarbruecken, FRG
+  (fax +49 681 31104).
+
+  Copyright (c) 1991-1997  by  Max-Planck-Institut fuer Informatik
+  Im Stadtwald, 66123 Saarbruecken, Germany     
+  All rights reserved.
+ 
*******************************************************************************/

#include <LEDA/circle.h>
#include <math.h>


//------------------------------------------------------------------------------
// circles
//
// S. N"aher (1996)
//------------------------------------------------------------------------------


circle_rep::circle_rep(const point& p1, const point& p2, const point& p3) : 
a(p1), b(p2), c(p3), cp(0), rp(0), first_side_of(true)
{ 
  orient = orientation(p1,p2,p3);

  if (orient == 0 && (p1 == p2 || p1 == p3 || p2 == p3) )
  { if (p1 == p2 && p1 == p3)
      { cp = new point(p1);
        rp = new double(0);
       }
    else
       error_handler(1,"circle::circle: non-admissable triple");
  }
} 
                              

circle_rep::~circle_rep() 
{ if (cp) delete cp;
  if (rp) delete rp;
 }



circle::circle(const point& a, const point& b, const point& c)
{ PTR = new circle_rep(a,b,c); }

circle::circle(const point& a, const point& b0)
{ // center a and b0 on circle
  point b1 = b0.rotate90(a);
  point b2 = b1.rotate90(a);
  PTR = new circle_rep(b0,b1,b2);
  ptr()->cp = new point(a);
  ptr()->rp = new double(a.distance(b0));
 }

circle::circle(const point& a) 
{ PTR = new circle_rep(a,a,a); }


circle::circle() 
{ point p(0,0);
  PTR = new circle_rep(p,p,p);
 }

circle::circle(const point& m, double r)       
{ point a(m.xcoord(),  m.ycoord()+r); 
  point b(m.xcoord(),  m.ycoord()-r); 
  point c(m.xcoord()+r,m.ycoord()); 
  PTR = new circle_rep(a,b,c); 
  ptr()->cp = new point(m);
  ptr()->rp = new double(r);
}

circle::circle(double x, double y, double r)  
{ PTR = new circle_rep(point(x,y+r),point(x+r,y),point(x,y-r)); 
  ptr()->cp = new point(x,y);
  ptr()->rp = new double(r);
}



point circle::center() const
{ if (ptr()->cp == 0)
  { if (collinear(ptr()->a,ptr()->b,ptr()->c)) 
       error_handler(1,"circle::center(): points are collinear.");
    line l1 = p_bisector(ptr()->a,ptr()->b);
    line l2 = p_bisector(ptr()->b,ptr()->c);
    point m;
    l1.intersection(l2,m);
    ptr()->cp = new point(m);
   }
  return *(ptr()->cp);
 }


double circle::radius() const
{ if (ptr()->rp == 0)
    ptr()->rp = new double(center().distance(ptr()->a));
 return *(ptr()->rp);
}

  

bool circle::operator==(const circle& c)  const
{ if (!contains(c.ptr()->a)) return false;
  if (!contains(c.ptr()->b)) return false;
  if (!contains(c.ptr()->c)) return false;
  return true;
}


double circle::distance(const point& p) const
{ double d = p.distance(center());
  return (d - radius());
 }


double circle::distance(const line& l) const
{ double d = l.distance(center());
  return (d - radius());
 }

double circle::distance(const circle& c) const
{ double d = center().distance(c.center());
  return (d - radius() - c.radius());
 }

int circle::side_of(const point& p) const
{ 
  double ax = ptr()->a.xcoord();
  double ay = ptr()->a.ycoord();

  // compute D1 D2 D3 (if not done before)

  if (ptr()->first_side_of) 
  { double bx = ptr()->b.xcoord() - ax;
    double by = ptr()->b.ycoord() - ay;
    double bw = bx*bx + by*by;

    double cx = ptr()->c.xcoord() - ax;
    double cy = ptr()->c.ycoord() - ay;
    double cw = cx*cx + cy*cy;

    ptr()->D1 = cy*bw - by*cw;
    ptr()->D2 = bx*cw - cx*bw;
    ptr()->D3 = by*cx - bx*cy;

    ptr()->first_side_of = false;
   }

   double px = p.xcoord() - ax;
   double py = p.ycoord() - ay;
   double pw = px*px + py*py;
 
   double D  = px*ptr()->D1  + py*ptr()->D2 + pw*ptr()->D3;

   if (D != 0)
      return (D > 0) ? 1 : -1;
   else
      return 0;
}



bool circle::outside(const point& p) const 
{ return (ptr()->orient * side_of(p)) < 0; }


bool circle::inside(const point& p) const 
{ return (ptr()->orient * side_of(p)) > 0; }


bool circle::contains(const point& p) const 
{ return side_of(p) == 0; }


circle circle::translate_by_angle(double alpha, double d) const
{ point a1 = ptr()->a.translate_by_angle(alpha,d);
  point b1 = ptr()->b.translate_by_angle(alpha,d);
  point c1 = ptr()->c.translate_by_angle(alpha,d);
  return circle(a1,b1,c1);
 }


circle circle::translate(double dx, double dy) const
{ point a1 = ptr()->a.translate(dx,dy);
  point b1 = ptr()->b.translate(dx,dy);
  point c1 = ptr()->c.translate(dx,dy);
  return circle(a1,b1,c1);
}


circle circle::translate(const vector& v) const
{ point a1 = ptr()->a.translate(v);
  point b1 = ptr()->b.translate(v);
  point c1 = ptr()->c.translate(v);
  return circle(a1,b1,c1);
 }

circle  circle::rotate(const point& o, double alpha) const
{ point a1 = ptr()->a.rotate(o,alpha);
  point b1 = ptr()->b.rotate(o,alpha);
  point c1 = ptr()->c.rotate(o,alpha);
  return circle(a1,b1,c1);
 }

circle  circle::rotate(double alpha)  const
{ point o(0,0);
  point a1 = ptr()->a.rotate(o,alpha);
  point b1 = ptr()->b.rotate(o,alpha);
  point c1 = ptr()->c.rotate(o,alpha);
  return circle(a1,b1,c1);
}


circle  circle::rotate90(const point& o)  const
{ point a1 = ptr()->a.rotate90(o);
  point b1 = ptr()->b.rotate90(o);
  point c1 = ptr()->c.rotate90(o);
  return circle(a1,b1,c1);
}


circle  circle::reflect(const point& p, const point& q)  const
{ point a1 = ptr()->a.reflect(p,q);
  point b1 = ptr()->b.reflect(p,q);
  point c1 = ptr()->c.reflect(p,q);
  return circle(a1,b1,c1);
}

circle  circle::reflect(const point& p)  const
{ point a1 = ptr()->a.reflect(p);
  point b1 = ptr()->b.reflect(p);
  point c1 = ptr()->c.reflect(p);
  return circle(a1,b1,c1);
}



list<point> circle::intersection(const line& l) const
{ list<point> result;

  segment s = l.perpendicular(center());
  double  d = s.length();
  double  r = radius();

  if (d==r) result.append(s.end());
  
  if (d < r)
  { segment sl = l.ptr()->seg;
    double t = sqrt(r*r - d*d)/sl.length();
    double dx = t*sl.dx();
    double dy = t*sl.dy();
    result.append(s.end().translate(vector(dx,dy)));
    result.append(s.end().translate(vector(-dx,-dy)));
  }

  return result;
}


list<point> circle::intersection(const segment& s) const
{ list<point> result,L;

  L = intersection(line(s));

  point p;
  double  d  = s.length();

  forall(p,L)
  { double d1 = s.ptr()->start.distance(p);
    double d2 = s.ptr()->end.distance(p);
    if (d1 <= d && d2 <= d) result.append(p);
   }

  return result;
}

list<point> circle::intersection(const circle& c) const
{ list<point> result;
  segment s(center(), c.center());
  double d  = s.length();
  double r1 = radius();
  double r2 = c.radius();

  if (d > (r1+r2) || (d+r2) < r1 || (d+r1) < r2) return result;


  double x = (d*d + r1*r1 - r2*r2)/(2*d);
  double alpha = acos(x/r1);
  double beta  = s.angle() + alpha;
  double gamma = s.angle() - alpha;

  result.append(center().translate_by_angle(beta,r1));
  if (alpha!=0) result.append(center().translate_by_angle(gamma,r1));

  return result;
 }


segment circle::left_tangent(const point& p) const
{ if (inside(p)) error_handler(1,"left_tangent:: point inside circle");
  segment s(p,center());
  double d = s.length();
  double alpha = asin(radius()/d) + s.angle();
  point touch = p.translate_by_angle(alpha,sqrt(d*d - radius()*radius()));
  return segment(p,touch);
}

segment circle::right_tangent(const point& p) const
{ if (inside(p)) error_handler(1,"right_tangent:: point inside circle");
  segment s(p,center());
  double d = s.length();
  double alpha = s.angle() - asin(radius()/d);
  point touch = p.translate_by_angle(alpha,sqrt(d*d - radius()*radius()));
  return segment(p,touch);
}


point circle::point_on_circle(double phi,double) const
{ point  c = center();
  double r = radius();
  double dx = cos(phi) * r;
  double dy = sin(phi) * r;
  return point(c.xcoord()+dx,c.ycoord()+dy);
}



ostream& operator<<(ostream& out, const circle& c) 
{ out << c.point1() << " " << c.point2() << " " << c.point3();
  return out;
 } 

istream& operator>>(istream& in, circle& c) 
{ point p,q,r;
  in >> p >> q >> r;
  c = circle(p,q,r);
  return in;
}



