#include	"main.h"


int edge_code(n0, n1)
        Vector n0, n1;
{
  if (g3_dot(n0, n1) < dot_tol)
    return 0;
  else
    return 1;
}


#define MIDPOINT(e, p0, p1, pm, nm) \
   pm = g3_scalev(0.5, g3_sumv(p0, p1)); \
   nm = f_ngrad(pm); \
   if (e == 0) \
      pm = project(STEP_0, pm, nm, f_value(pm), MAX_REC);


adapt_tri(p0, p1, p2, n0, n1, n2)
     Vector p0, p1, p2;
     Vector n0, n1, n2;
{
  int e1, e2, e3;
  Vector pm1, pm2, pm3;
  Vector nm1, nm2, nm3;

  e1 = edge_code(n0, n1);
  e2 = edge_code(n1, n2);
  e3 = edge_code(n2, n0);
	
  if (e1 && e2 && e3)  {
    output_tri(p0, p1, p2, n0, n1, n2);
  } else {
    MIDPOINT(e1, p0, p1, pm1, nm1);
    MIDPOINT(e2, p1, p2, pm2, nm2);
    MIDPOINT(e3, p2, p0, pm3, nm3);
	
    adapt_tri(p0, pm1, pm3, n0, nm1, nm3);
    adapt_tri(p1, pm2, pm1, n1, nm2, nm1);
    adapt_tri(p2, pm3, pm2, n2, nm3, nm2);
    adapt_tri(pm1, pm2, pm3, nm1, nm2, nm3);
  }
}


Vector project(s, p, n, v0, k)
     double s;
     Vector p, n;
     double v0;
     int k;
{
  double v1;
  
  p = g3_sumv(p, g3_scalev(s * sign(v0) * -1.0, n));
  v1 = f_value(p);
	
  if ((fabs(v1) < eps_dist) || k-- == 0)
    return p;
	
  if (v0*v1 < 0.0)
    s = s * 0.5;

  return project(s, p, n, v1, k);
}