// GeoElt2d.cpp
//
// J.Liu, R.Cali
// ColoState, Jan.-Jun. 2007


#include <cmath>
#include "GeoElt2d.h"
#include "PtVec2d.h"


// Line2d: copy assignment 

Line2d &Line2d::operator=(const Line2d &L) 
{
   if (this != &L) {
      P1 = L.P1;
      P2 = L.P2;
   }
   return *this;
}


// Line2d: a normal vector

PtVec2d Line2d::normal() const 
{
   double x = P2.xCrd() - P1.xCrd();
   double y = P2.yCrd() - P1.yCrd();
   double tmp = sqrt(x*x+y*y);
   return PtVec2d(y/tmp, -x/tmp);
}


// Rect2d: copy assignment 

Rect2d &Rect2d::operator=(const Rect2d &R) 
{
   if (this != &R) {
      P1 = R.P1;
      P2 = R.P2;
   }
   return *this;
}


// Rect2d: get the coordinates

void Rect2d::getCrds(double &x1, double &y1, double &x2, double &y2) const 
{
   x1 = P1.xCrd();
   y1 = P1.yCrd();
   x2 = P2.xCrd();
   y2 = P2.yCrd();
}


// Rect2d: maximal and minimal x-, y- coordinates

void Rect2d::getMaxMinCrds(double &xmax, double &xmin, 
   double &ymax, double &ymin) const 
{
   xmax = P2.xCrd();
   xmin = P1.xCrd();
   ymax = P2.yCrd();
   ymin = P1.yCrd();
   return;
}


// Rect2d: maximal and minimal Edge lengths

void Rect2d::getMaxMinEdges(double &maxEdge, double &minEdge) const 
{
   double sx = fabs(P2.xCrd()-P1.xCrd());
   double sy = fabs(P2.yCrd()-P1.yCrd());

   if (sx>=sy) {
      maxEdge = sx;
      minEdge = sy;
   } else {
      maxEdge = sy;
      minEdge = sx;
   }

   return;
}


// Rect2d: mapping from the reference rectangle [0,1]^2 
// with natural coordinates a, b

PtVec2d Rect2d::mapping(double a, double b) const 
{
   double x1 = P1.xCrd();
   double y1 = P1.yCrd();
   double x2 = P2.xCrd();
   double y2 = P2.yCrd();
   double x = (1-a)*x1 + a*x2;
   double y = (1-b)*y1 + b*y2;
   return PtVec2d(x, y);
}


// Rect2d: mapping from the reference rectangle [-1,1]^2 
// with natural coordinates a, b

PtVec2d Rect2d::mapping1(double a, double b) const 
{
   double x1 = P1.xCrd();
   double y1 = P1.yCrd();
   double x2 = P2.xCrd();
   double y2 = P2.yCrd();
   double x = (1-a)/2*x1 + (1+a)/2*x2;
   double y = (1-b)/2*y1 + (1+b)/2*y2;
   return PtVec2d(x, y);
}


// Rect2d: is a given point in the rectangle 

bool Rect2d::isInRect2d(PtVec2d P) const 
{
   bool b = false;
   if (P1<=P && P<=P2)  b = true;
/*
   double x = P.xCrd();
   double y = P.yCrd();
   double x1 = P1.xCrd();
   double y1 = P1.yCrd();
   double x2 = P2.xCrd();
   double y2 = P2.yCrd();
   bool b = false;
   if (x1<=x && x<=x2 && y1<=y && y<=y2)  b = true;
*/
   return b;
}


// Trig: copy assignment 

Trig &Trig::operator=(const Trig &T) 
{
   if (this != &T) {
      P1 = T.P1;
      P2 = T.P2;
      P3 = T.P3;
   }
   return *this;
}


// Trig: get the coordinates

void Trig::getCrds(double &x1, double &y1, double &x2, double &y2, 
   double &x3, double &y3) const 
{
   x1 = P1.xCrd();  y1 = P1.yCrd();
   x2 = P2.xCrd();  y2 = P2.yCrd();
   x3 = P3.xCrd();  y3 = P3.yCrd();
   return;
}


// Trig: maximal and minimal x-, y- coordinates

void Trig::getMaxMinCrds(double &xmax, double &xmin, 
   double &ymax, double &ymin) const 
{
   double x1 = P1.xCrd();
   double y1 = P1.yCrd();
   double x2 = P2.xCrd();
   double y2 = P2.yCrd();
   double x3 = P3.xCrd();
   double y3 = P3.yCrd();

   xmax = x1;
   xmin = x1;
   if (x2>xmax) xmax = x2;
   if (x3>xmax) xmax = x3;
   if (x2<xmin) xmin = x2;
   if (x3<xmin) xmin = x3;

   ymax = y1;
   ymin = y1;
   if (y2>ymax) ymax = y2;
   if (y3>ymax) ymax = y3;
   if (y2<ymin) ymin = y2;
   if (y3<ymin) ymin = y3;

   return;
}


// Trig: maximal and minimal edges

void Trig::getMaxMinEdges(double &maxEdge, double &minEdge) const 
{
   double s1 = dist(P2, P3);
   double s2 = dist(P3, P1);
   double s3 = dist(P1, P2);
// cout << "s1, s2, s3= " << s1 << "," << s2 << "," << s3 << "\n";
/*
   double x1 = P1.xCrd();
   double y1 = P1.yCrd();
   double x2 = P2.xCrd();
   double y2 = P2.yCrd();
   double x3 = P3.xCrd();
   double y3 = P3.yCrd();
   double s1 = sqrt((x2-x3)*(x2-x3)+(y2-y3)*(y2-y3));
   double s2 = sqrt((x3-x1)*(x3-x1)+(y3-y1)*(y3-y1));
   double s3 = sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
*/
   maxEdge = s1;
   minEdge = s1;
   if (s2>maxEdge)  maxEdge = s2;
   if (s3>maxEdge)  maxEdge = s3;
   if (s2<minEdge)  minEdge = s2;
   if (s3<minEdge)  minEdge = s3;

   return;
}


// Trig: area

double Trig::area() const 
{
   double x1 = P1.xCrd();
   double y1 = P1.yCrd();
   double x2 = P2.xCrd();
   double y2 = P2.yCrd();
   double x3 = P3.xCrd();
   double y3 = P3.yCrd();
   return 0.5*fabs((x2-x1)*(y3-y1)-(x3-x1)*(y2-y1));
}


// Trig: inverse mapping to the reference triangle

PtVec2d Trig::invMapping(PtVec2d P) const
{
   double a, b, det1;

   double x = P.xCrd();
   double y = P.yCrd();
   double x1 = P1.xCrd();
   double y1 = P1.yCrd();
   double x2 = P2.xCrd();
   double y2 = P2.yCrd();
   double x3 = P3.xCrd();
   double y3 = P3.yCrd();

   det1 = 1.0/((x2-x1)*(y3-y1)-(x3-x1)*(y2-y1));
   a = det1*((x-x1)*(y3-y1)+(x1-x3)*(y-y1));   
   b = det1*((x-x1)*(y1-y2)+(x2-x1)*(y-y1));

   return PtVec2d(a, b);
} 


// Trig: is a given point inside the triangle

bool Trig::isInTrig(PtVec2d P) const 
{
   PtVec2d Q = invMapping(P);
   double a = Q.xCrd();
   double b = Q.yCrd();

   if (a>=0.0 && b>=0.0 && a+b<=1.0) {
      return true;
   }
   else {
      return false;
   }
}


// Trig: hyp intersection 

PtVec2d Trig::hypInterSxn(const PtVec2d &Pi, const PtVec2d &Po) const 
{
   double a, b, s;

   PtVec2d Qi = invMapping(Pi);
   PtVec2d Qo = invMapping(Po);
   double xi = Qi.xCrd();
   double yi = Qi.yCrd();
   double xo = Qo.xCrd();
   double yo = Qo.yCrd();

   if (xi!=xo) {
      s = (yi-yo)/(xi-xo);
      a = (s*xo+1-yo)/(1+s);
      b = 1-a;
   }

   if (xi==xo) {
      a = xi;
      b = 1-a;
   }

   return mapping(a, b);
}


PtVec2d Trig::normalEdge1() const
{
   double x1 = P1.xCrd();
   double y1 = P1.yCrd();
   double x2 = P2.xCrd();
   double y2 = P2.yCrd();
   double x3 = P3.xCrd();
   double y3 = P3.yCrd();
   double a, tmp, x, y;

   a = (-x2*(x3-x2)+x1*(x3-x2)-y2*(y3-y2)+y1*(y3-y2))/(-x2*(x3-x2)+x3*(x3-x2)-y2*(y3-y2)+y3*(y3-y2));
   x = (1-a)*x2+a*x3-x1;
   y = (1-a)*y2+a*y3-y1;
   tmp = sqrt(x*x+y*y);
   
   return PtVec2d(x/tmp, y/tmp);
}


PtVec2d Trig::normalEdge2() const
{
   double x1 = P1.xCrd();
   double y1 = P1.yCrd();
   double x2 = P2.xCrd();
   double y2 = P2.yCrd();
   double x3 = P3.xCrd();
   double y3 = P3.yCrd();
   double a, tmp, x, y;

   a = (-x1*(x3-x1)+x2*(x3-x1)-y1*(y3-y1)+y2*(y3-y1))/(-x1*(x3-x1)+x3*(x3-x1)-y1*(y3-y1)+y3*(y3-y1));
   x = (1-a)*x1+a*x3-x2;
   y = (1-a)*y1+a*y3-y2;
   tmp = sqrt(x*x+y*y);

   return PtVec2d(x/tmp, y/tmp);
}


PtVec2d Trig::normalEdge3() const
{
   double x1 = P1.xCrd();
   double y1 = P1.yCrd();
   double x2 = P2.xCrd();
   double y2 = P2.yCrd();
   double x3 = P3.xCrd();
   double y3 = P3.yCrd();
   double a, tmp, x, y;
 
   a = (-x1*(x2-x1)+x3*(x2-x1)-y1*(y2-y1)+y3*(y2-y1))/(-x1*(x2-x1)+x2*(x2-x1)-y1*(y2-y1)+y2*(y2-y1));
   x = (1-a)*x1+a*x2-x3;
   y = (1-a)*y1+a*y2-y3;
   tmp = sqrt(x*x+y*y);

   return PtVec2d(x/tmp, y/tmp);
}