/* -*-ePiX-*- */
/* levelset.xp -- level sets */
#include "epix.h"
using namespace ePiX;

const int MAX(3);
const double k(0.25);

P Fp(double u, double v) 
{
  double r(sqrt(fabs(v)));
  if (v>=0)
    return P(r*cosh(u), r*sinh(u), k*v);

  else
    return P(-r*sinh(u), -r*cosh(u), k*v);

}

P Fm(double u, double v) 
{
  double r(sqrt(fabs(v)));
  if (v>=0)
    return P(-r*cosh(u), r*sinh(u), k*v);

  else
    return P(-r*sinh(u), r*cosh(u), k*v);

}

double f(double x, double y) 
{ 
  return k*(x-y)*(x+y); 
}

int main()
{
  picture(P(-MAX,-MAX), P(MAX, MAX), "2x2in");

  begin();

  camera.at(P(1,-3,5));
  camera.range(40);

  clip_box(P(MAX,MAX,2*MAX));
  clip();

  red();
  plot(f, P(-MAX,-MAX), P(MAX, MAX), mesh(4,4), mesh(40,40));

  domain R(P(-4,-9), P(4,9), mesh(1,12), mesh(90,1));

  bold(); rgb(0.5,0.5,1); 
  line(P(MAX,-MAX,0), P(-MAX, MAX,0));
  line(P(MAX, MAX,0), P(-MAX,-MAX,0));
  plot(Fp, R.slices2());
  plot(Fm, R.slices2());

  end(); 
}