/*******************************************************************************
+
+  LEDA 3.5
+
+  _d3_rat_plane.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/d3_rat_plane.h>
#include <math.h>
#include <ctype.h>

//------------------------------------------------------------------------------
// d3_rat_plane 
//
// by S. Naeher (1997)
//------------------------------------------------------------------------------


leda_mutex d3_rat_plane_rep::mutex_id_counter;

unsigned long d3_rat_plane_rep::id_counter = 0;


d3_rat_plane_rep::d3_rat_plane_rep(const d3_rat_point& p, const integer& x,
                                                          const integer& y,
                                                          const integer& z)
: a(p),nx(x),ny(y),nz(z)

{ if (nx == 0 && ny == 0 && nz == 0)
     error_handler(1,"d3_rat_plane::d3_rat_plane(): cannot construct plane.");
  mutex_id_counter.lock();
  id  = id_counter++; 
  mutex_id_counter.unlock();
}


d3_rat_plane::d3_rat_plane(const d3_rat_point& a, const d3_rat_point& b, const d3_rat_point& c)
{ integer X1 = b.X() - a.X();
  integer Y1 = b.Y() - a.Y();
  integer Z1 = b.Z() - a.Z();
  integer X2 = c.X() - a.X();
  integer Y2 = c.Y() - a.Y();
  integer Z2 = c.Z() - a.Z();
  integer nx = Z1*Y2 - Y1*Z2;
  integer ny = X1*Z2 - Z1*X2;
  integer nz = Y1*X2 - X1*Y2;
  PTR = new d3_rat_plane_rep(a,nx,ny,nz); 
}


d3_rat_plane::d3_rat_plane(const d3_rat_point& p, const rat_vector& v) 
{ PTR = new d3_rat_plane_rep(p,v.hcoord(0),v.hcoord(1),v.hcoord(2)); }


d3_rat_plane d3_rat_plane::translate(integer dx, integer dy, integer dz, integer dw) const
{ d3_rat_point p = point1().translate(dx,dy,dz,dw);
  return d3_rat_plane(p,normal());
}

d3_rat_plane d3_rat_plane::translate(const rational& dx, 
                                     const rational& dy, 
                                     const rational& dz) const
{ integer W = dx.denominator()*dy.denominator()*dz.denominator(); 
  return translate(dx.numerator(),dy.numerator(),dz.numerator(), W); 
 }


d3_rat_plane d3_rat_plane::translate(const rat_vector& v) const 
{ d3_rat_point p = point1().translate(v);
  return d3_rat_plane(p,normal());
}


// Distances

d3_rat_plane d3_rat_plane::reflect(const d3_rat_point& q) const
{ // reflect plane across point q
  d3_rat_point a = point1();
  d3_rat_point b = a.translate(normal());
  d3_rat_point a1 = a.reflect(q);
  d3_rat_point b1 = b.reflect(q);
  return d3_rat_plane(a1,b1-a1);
 }


rat_vector d3_rat_plane::normal_project(const d3_rat_point& p) const
{ rat_vector v = p - point1();
  integer nx = ptr()->nx;
  integer ny = ptr()->ny;
  integer nz = ptr()->nz;

  integer A = -(nx*v.hcoord(0)+ny*v.hcoord(1)+nz*v.hcoord(2));
  integer W = nx*nx+ny*ny+nz*nz*v.hcoord(3);
  return rat_vector(A*nx,A*ny,A*nz,W);
}

d3_rat_point d3_rat_plane::reflect_point(const d3_rat_point& p) const
{ return p.translate(2*normal_project(p)); } 


d3_rat_plane  d3_rat_plane::reflect(const d3_rat_plane& Q) const
{ // reflect plane across plane q
  d3_rat_point a = point1();
  d3_rat_point b = a.translate(normal());
  d3_rat_point a1 = Q.reflect_point(a);
  d3_rat_point b1 = Q.reflect_point(b);
  return d3_rat_plane(a1,b1-a1);
 }


rational d3_rat_plane::sqr_dist(const d3_rat_point& p)  const
{ rat_vector vec = normal_project(p);
  return vec.sqr_length();
 }


bool d3_rat_plane::parallel(const d3_rat_plane& Q) const
{ integer nx1 = ptr()->nx;
  integer ny1 = ptr()->ny;
  integer nz1 = ptr()->nz;
  integer nx2 = Q.ptr()->nx;
  integer ny2 = Q.ptr()->ny;
  integer nz2 = Q.ptr()->nz;
  return nz1*ny2 == ny1*nz2 && nx1*nz2 == nz1*nx2 && ny1*nx2 == nx1*ny2;
}


int d3_rat_plane::side_of(const d3_rat_point& q) const
{ 
  d3_rat_point a = point1();

  integer qx = q.X()*a.W() - a.X()*q.W();
  integer qy = q.Y()*a.W() - a.Y()*q.W();
  integer qz = q.Z()*a.W() - a.Z()*q.W();

  integer nx = ptr()->nx;
  integer ny = ptr()->ny;
  integer nz = ptr()->nz;

  integer D = nx*qx + ny*qy + nz*qz;

  if (D > 0) return +1;
  if (D < 0) return -1;
  return 0;
}


d3_plane d3_rat_plane::to_d3_plane() const
{ return d3_plane(point1().to_d3_point(),normal().to_vector()); }


int d3_rat_plane::operator==(const d3_rat_plane& P) const 
{ return parallel(P) && contains(P.point1()); }

   
ostream& operator<<(ostream& out, const d3_rat_plane& P)
{ cout << P.point1() << " " << P.normal(); return out; } 

istream& operator>>(istream& in, d3_rat_plane& P) 
{ d3_rat_point p;
  rat_vector n(3);
  cin >> p >> n;
  P = d3_rat_plane(p,n);
  return in; 
} 
