/*****
 * runpath3.in
 *
 * Runtime functions for path3 operations.
 *
 *****/

pair     => primPair()
triple   => primTriple()
path3     => primPath3()
boolarray* => booleanArray()
realarray* => realArray()
realarray2* => realArray2()
triplearray* => tripleArray()
triplearray2* => tripleArray2()

#include "path3.h"
#include "array.h"
#include "drawsurface.h"
#include "predicates.h"

using namespace camp;
using namespace vm;

typedef array boolarray;
typedef array realarray;
typedef array realarray2;
typedef array triplearray;
typedef array triplearray2;

using types::booleanArray;
using types::realArray;
using types::realArray2;
using types::tripleArray;
using types::tripleArray2;

// Autogenerated routines:


path3 path3(triplearray *pre, triplearray *point, triplearray *post,
            boolarray *straight, bool cyclic)
{
  size_t n=checkArrays(pre,point);
  checkEqual(n,checkArray(post));
  checkEqual(n,checkArray(straight));
  mem::vector<solvedKnot3> nodes(n);
  for(size_t i=0; i < n; ++i) {
    nodes[i].pre=read<triple>(pre,i);
    nodes[i].point=read<triple>(point,i);
    nodes[i].post=read<triple>(post,i);
    nodes[i].straight=read<bool>(straight,i);
  }

  return path3(nodes,(Int) n,cyclic);
}

path3 :nullPath3()
{
  return nullpath3;
}

bool ==(path3 a, path3 b)
{
  return a == b;
}

bool !=(path3 a, path3 b)
{
  return !(a == b);
}

triple point(path3 p, Int t)
{
  return p.point((Int) t);
}

triple point(path3 p, real t)
{
  return p.point(t);
}

triple precontrol(path3 p, Int t)
{
  return p.precontrol((Int) t);
}

triple precontrol(path3 p, real t)
{
  return p.precontrol(t);
}

triple postcontrol(path3 p, Int t)
{
  return p.postcontrol((Int) t);
}

triple postcontrol(path3 p, real t)
{
  return p.postcontrol(t);
}

triple dir(path3 p, Int t, Int sign=0, bool normalize=true)
{
  return p.dir(t,sign,normalize);
}

triple dir(path3 p, real t, bool normalize=true)
{
  return p.dir(t,normalize);
}

triple accel(path3 p, Int t, Int sign=0)
{
  return p.accel(t,sign);
}

triple accel(path3 p, real t)
{
  return p.accel(t);
}

real radius(path3 p, real t)
{
  triple v=p.dir(t,false);
  triple a=p.accel(t);
  real d=dot(a,v);
  real v2=v.abs2();
  real a2=a.abs2();
  real denom=v2*a2-d*d;
  real r=v2*sqrt(v2);
  return denom > 0 ? r/sqrt(denom) : 0.0;
}

real radius(triple z0, triple c0, triple c1, triple z1, real t)
{
  triple v=(3.0*(z1-z0)+9.0*(c0-c1))*t*t+(6.0*(z0+c1)-12.0*c0)*t+3.0*(c0-z0);
  triple a=6.0*(z1-z0+3.0*(c0-c1))*t+6.0*(z0+c1)-12.0*c0;
  real d=dot(a,v);
  real v2=v.abs2();
  real a2=a.abs2();
  real denom=v2*a2-d*d;
  real r=v2*sqrt(v2);
  return denom > 0 ? r/sqrt(denom) : 0.0;
}

path3 reverse(path3 p)
{
  return p.reverse();
}

path3 subpath(path3 p, Int a, Int b)
{
  return p.subpath((Int) a, (Int) b);
}

path3 subpath(path3 p, real a, real b)
{
  return p.subpath(a,b);
}

Int length(path3 p)
{
  return p.length();
}

bool cyclic(path3 p)
{
  return p.cyclic();
}

bool straight(path3 p, Int t)
{
  return p.straight(t);
}

path3 unstraighten(path3 p)
{
  return p.unstraighten();
}

// return the maximum distance squared of points c0 and c1 from
// the respective internal control points of z0--z1.
real straightness(triple z0, triple c0, triple c1, triple z1)
{
  return Straightness(z0,c0,c1,z1);
}

// return the straightness of segment i of path3 g.
real straightness(path3 p, Int t)
{
  if(p.straight(t)) return 0;
  return Straightness(p.point(t),p.postcontrol(t),p.precontrol(t+1),
                      p.point(t+1));
}

bool piecewisestraight(path3 p)
{
  return p.piecewisestraight();
}

real arclength(path3 p)
{
  return p.arclength();
}

real arclength(triple z0, triple c0, triple c1, triple z1)
{
  return arcLength(z0,c0,c1,z1);
}

real arctime(path3 p, real dval)
{
  return p.arctime(dval);
}

realarray* intersect(path3 p, path3 q, real fuzz=-1)
{
  bool exact=fuzz <= 0.0;
  if(fuzz < 0)
    fuzz=BigFuzz*::max(::max(length(p.max()),length(p.min())),
                       ::max(length(q.max()),length(q.min())));

  std::vector<real> S,T;
  real s,t;
  if(intersections(s,t,S,T,p,q,fuzz,true,exact)) {
    array *V=new array(2);
    (*V)[0]=s;
    (*V)[1]=t;
    return V;
  } else
    return new array(0);
}

realarray2* intersections(path3 p, path3 q, real fuzz=-1)
{
  bool exact=fuzz <= 0.0;
  if(fuzz < 0)
    fuzz=BigFuzz*::max(::max(length(p.max()),length(p.min())),
                       ::max(length(q.max()),length(q.min())));
  bool single=!exact;

  real s,t;
  std::vector<real> S,T;
  bool found=intersections(s,t,S,T,p,q,fuzz,single,exact);
  if(!found) return new array(0);
  array *V;
  if(single) {
    V=new array(1);
    array *Vi=new array(2);
    (*V)[0]=Vi;
    (*Vi)[0]=s;
    (*Vi)[1]=t;
  } else {
    size_t n=S.size();
    V=new array(n);
    for(size_t i=0; i < n; ++i) {
      array *Vi=new array(2);
      (*V)[i]=Vi;
      (*Vi)[0]=S[i];
      (*Vi)[1]=T[i];
    }
  }
  stable_sort(V->begin(),V->end(),run::compare2<real>());
  return V;
}

realarray* intersect(path3 p, triplearray2 *P, real fuzz=-1)
{
  triple *A;
  copyArray2C(A,P,true,4);
  if(fuzz <= 0) fuzz=BigFuzz*::max(::max(length(p.max()),length(p.min())),
                                   norm(A,16));
  std::vector<real> T,U,V;
  bool found=intersections(T,U,V,p,A,fuzz,true);
  delete[] A;
  if(found) {
    array *W=new array(3);
    (*W)[0]=T[0];
    (*W)[1]=U[0];
    (*W)[2]=V[0];
    return W;
  } else
    return new array(0);
}

realarray2* intersections(path3 p, triplearray2 *P, real fuzz=-1)
{
  triple *A;
  copyArray2C(A,P,true,4);
  if(fuzz <= 0) fuzz=BigFuzz*::max(::max(length(p.max()),length(p.min())),
                                   norm(A,16));
  std::vector<real> T,U,V;
  intersections(T,U,V,p,A,fuzz,false);
  delete[] A;
  size_t n=T.size();
  array *W=new array(n);
  for(size_t i=0; i < n; ++i) {
    array *Wi=new array(3);
    (*W)[i]=Wi;
    (*Wi)[0]=T[i];
    (*Wi)[1]=U[i];
    (*Wi)[2]=V[i];
  }
  return W; // Sorting will done in asy.
}

Int size(path3 p)
{
  return p.size();
}

path3 &(path3 p, path3 q)
{
  return camp::concat(p,q);
}

triple min(path3 p)
{
  return p.min();
}

triple max(path3 p)
{
  return p.max();
}

realarray *mintimes(path3 p)
{
  array *V=new array(3);
  triple v=p.mintimes();
  (*V)[0]=v.getx();
  (*V)[1]=v.gety();
  (*V)[2]=v.getz();
  return V;
}

realarray *maxtimes(path3 p)
{
  array *V=new array(3);
  triple v=p.maxtimes();
  (*V)[0]=v.getx();
  (*V)[1]=v.gety();
  (*V)[2]=v.getz();
  return V;
}

path3 Operator *(realarray2 *t, path3 g)
{
  return transformed(*t,g);
}

pair minratio(path3 g)
{
  return g.ratio(::min);
}

pair maxratio(path3 g)
{
  return g.ratio(::max);
}

// Return a negative (positive) value if a--b--c--cycle is oriented
// counterclockwise (clockwise) when viewed from d or zero if all four
// points are coplanar.
// The value returned is the determinant
// |a.x a.y a.z 1|
// |b.x b.y b.z 1|
// |c.x c.y c.z 1|
// |d.x d.y d.z 1|
real orient(triple a, triple b, triple c, triple d)
{
  real A[]={a.getx(),a.gety(),a.getz()};
  real B[]={b.getx(),b.gety(),b.getz()};
  real C[]={c.getx(),c.gety(),c.getz()};
  real D[]={d.getx(),d.gety(),d.getz()};
  return orient3d(A,B,C,D);
}

// Return a positive (negative) value if e lies inside (outside)
// the sphere passing through the points a,b,c,d oriented so that
// a--b--c--cycle appears in clockwise order when viewed from d
// or zero if all five points are cospherical.
// The value returned is the determinant
// |a.x a.y a.z a.x^2+a.y^2+a.z^2 1|
// |b.x b.y b.z b.x^2+b.y^2+b.z^2 1|
// |c.x c.y c.z c.x^2+c.y^2+c.z^2 1|
// |d.x d.y d.z d.x^2+d.y^2+d.z^2 1|
// |e.x e.y e.z e.x^2+e.y^2+e.z^2 1|
real insphere(triple a, triple b, triple c, triple d, triple e)
{
  real A[]={a.getx(),a.gety(),a.getz()};
  real B[]={b.getx(),b.gety(),b.getz()};
  real C[]={c.getx(),c.gety(),c.getz()};
  real D[]={d.getx(),d.gety(),d.getz()};
  real E[]={e.getx(),e.gety(),e.getz()};
  return insphere(A,B,C,D,E);
}