///////////////// Code for functions of famous surfaces /////////////////


// Sphere (r=1)

double sphere_function (double x, double y, double z) {
  double x2=x*x, y2=y*y, z2=z*z;

  return x2+y2+z2-1;
}


// Ellipsoid (r=1)

double ellipsoid_function (double x, double y, double z) {
  double x2=x*x, y2=y*y, z2=z*z;

  return x2+2*y2+4*z2-1;
}


// Torus (r=2)

double torus_function (double x, double y, double z) {
  double x2=x*x, y2=y*y, z2=z*z;
  double x4=x2*x2, y4=y2*y2, z4=z2*z2;

  return x4  + y4  + z4  + 2 *x2*  y2  + 2*
    x2*z2  + 2*y2*  z2  - 5 *x2  + 4* y2  - 5*z2+4;
}


// "Chair" (r=6)

double chair_function (double x, double y, double z) {
  double x2=x*x, y2=y*y, z2=z*z;
  double x4=x2*x2, y4=y2*y2, z4=z2*z2;

  return x4-1.2*x2*y2+3.6*x2*z2-7.50*x2+y4+3.6*y2*z2-7.50*y2+.2*z4-7.50*z2+64.0625-16.0*z*y2+16.0*x2*z;
}


// "Tanglecube" (r=4)

double tanglecube_function (double x, double y, double z) {
  double x2=x*x, y2=y*y, z2=z*z;
  double x4=x2*x2, y4=y2*y2, z4=z2*z2;

  return x4 - 5*x2 + y4 - 5*y2 + z4 - 5*z2 + 11.8;
}

double cube_function (double x, double y, double z){
  if( x < 1 && -x < 1 &&
      y < 1 && -y < 1 &&
      z < 1 && -z < 1 )
    return -1.;
  return 1.;
}

// Barth's octic surface (degree 8)

double octic_function (double x, double y, double z) {  // r=2
  double x2=x*x, y2=y*y, z2=z*z;
  double x4=x2*x2, y4=y2*y2, z4=z2*z2;
  double x6=x4*x2, y6=y4*y2, z6=z4*z2;
  double x8=x4*x4, y8=y4*y4, z8=z4*z4;

  return 43.30495169 *x2*  y2  + 43.30495169 *x2*  z2  + 43.30495169 *y2 *    z2 + 44.36067980 *x6*  y2     + 44.36067980* x6 * z2  + 66.54101970* x4*  y4  + 66.54101970* x4 * z4    + 44.36067980 *x2 * y6  - 11.70820393* x2  - 11.70820393* y2  - 11.70820393* z2     + 37.65247585 *x4 + 37.65247585 *y4  + 37.65247585* z4  + 11.09016995* x8     + 11.09016995  *y8  + 11.09016995* z8  + 133.0820394 *x2*  y4*  z2     + 133.0820394   *x2 * y2 * z4  + 44.36067980* x2 * z6  + 44.36067980 *y6 * z2     + 66.54101970 *y4 * z4  + 44.36067980 *y2 * z6  + 133.0820394*    x4*  y2 * z2     - 91.95742756 *x4 * y2  - 91.95742756 *x4  *z2  - 91.95742756* x2 * y4     - 91.95742756 *x2  *z4  - 91.95742756* y4 *   z2  - 183.9148551* x2  *y2  *z2     - 30.65247585 *x6  - 30.65247585* y6  - 91.95742756 *y2 * z4  - 30.65247585* z6   + 1.618033988;
}


// "Heart"

double heart_function (double x, double y, double z) {  // radius = 2

  return (2*x*x+y*y+z*z-1)*(2*x*x+y*y+z*z-1)*(2*x*x+y*y+z*z-1) - (0.1*x*x+y*y)*z*z*z;
}


// Klein's bottle

double klein_function (double x, double y, double z) {  // radius = 4

  return (x*x+y*y+z*z+2*y-1)*((x*x+y*y+z*z-2*y-1)*(x*x+y*y+z*z-2*y-1)-8*z*z)+16*x*z*(x*x+y*y+z*z-2*y-1);
}


// Ring

double ring_function (double x, double y, double z) {  // radius = ?
  double e=0.1;

  double f1 = x*x+y*y+z*z-1;
  double f2 = x;

  f1 = f1*f1-e*e;
  f2 = f2*f2-e*e;

  if (f1 < 0 && f2 < 0)
    return -1.;
  else if (f1 > 0 || f2 > 0)
    return 1.;
  else
    return 0.;
}


// False knot

double false_knot_function (double x, double y, double z) {  // radius = 1
  double d=1.2, e=0.1;

  double f1 = x*(x*x-z*z)-2*x*z*z-y*y+d*d-x*x-z*z-y*y;

  double m2 = z*(x*x-z*z)+2*x*x*z;
  double f2 = 4*y*y*(d*d-x*x-z*z-y*y) - m2*m2;

  f1 = f1*f1-e*e;
  f2 = f2*f2-e*e;

  if (f1 < 0 && f2 < 0)
    return -1.;
  else if (f1 > 0 || f2 > 0)
    return 1.;
  else
    return 0.;
}


// Knot 1

void puiss(double& x, double& y, int n) {

  double xx = 1, yy = 0;

  while(n>0) {
    if (n&1) {
      double xxx = xx, yyy = yy;
      xx = xxx*x - yyy*y;
      yy = xxx*y + yyy*x;
    }

    double xxx = x, yyy = y;
    x=xxx*xxx-yyy*yyy;
    y=2*xxx*yyy;

    n/=2;
  }

  x = xx;
  y = yy;
}



double knot1_function (double a, double b, double c) {  // radius = 4
  double e=0.1;

  double x, y, z, t, den;
  den=1+a*a+b*b+c*c;
  x=2*a/den;
  y=2*b/den;
  z=2*c/den;
  t=(1-a*a-b*b-c*c)/den;

  double x3=x, y3=y, z2=z, t2=t;
  puiss(x3,y3,3);
  puiss(z2,t2,2);

  double f1 = z2-x3;
  double f2 = t2-y3;

  f1 = f1*f1;
  f2 = f2*f2;
  e=e*e/(den-1);

  if (f1+f2-e < 0)
    return -1.;
  else if (f1+f2-e > 0)
    return 1.;
  else
    return 0.;
}

/*
double knot1_function (double x, double y, double z) {  // radius = 4
  double e=0.1;

  double c1, c2, c3, c4, den;
  den=1+x*x+y*y+z*z;
  c1=2*x/den;
  c2=2*y/den;
  c3=2*z/den;
  c4=(1-x*x-y*y-z*z)/den;

  double f1 = c1*(c1*c1-c2*c2)-2*c1*c2*c2-c3*c3+c4*c4;
  double f2 = c2*(c1*c1-c2*c2)+2*c1*c1*c2-2*c3*c4;

  f1 = f1*f1;
  f2 = f2*f2;
  e=e*e/(den-1);

  if (f1+f2-e < 0)
    return -1.;
  else if (f1+f2-e > 0)
    return 1.;
  else
    return 0.;
}
*/

double knot1_surf1_function (double x, double y, double z) {  // radius = 4
  double c1, c2, c3, c4, den;
  den=1+x*x+y*y+z*z;
  c1=2*x/den;
  c2=2*y/den;
  c3=2*z/den;
  c4=(1-x*x-y*y-z*z)/den;

  return c1*(c1*c1-c2*c2)-2*c1*c2*c2-c3*c3+c4*c4;
}


double knot1_surf2_function (double x, double y, double z) {  // radius = 4
  double c1, c2, c3, c4, den;
  den=1+x*x+y*y+z*z;
  c1=2*x/den;
  c2=2*y/den;
  c3=2*z/den;
  c4=(1-x*x-y*y-z*z)/den;

  return c2*(c1*c1-c2*c2)+2*c1*c1*c2-2*c3*c4;
}


// Knot 2

double knot2_function (double a, double b, double c) {  // radius = 4
  double e=0.025;

  double x, y, z, t, den;
  den=1+a*a+b*b+c*c;
  x=2*a/den;
  y=2*b/den;
  z=2*c/den;
  t=(1-a*a-b*b-c*c)/den;

  double x7=x, y7=y, x13=x, y13=y;
  puiss(x7,y7,7);
  puiss(x13,y13,13);

  double z3=z, t3=t;
  puiss(z3,t3,3);

  double f1t = (z3-x7)*(z3-x7+100*x13) - (t3-y7)*(t3-y7+100*y13);
  double f2t = (z3-x7)*(t3-y7+100*y13) + (t3-y7)*(z3-x7+100*x13);
  double f1 = f1t*(z3-x7-100*x13) - f2t*(t3-y7-100*y13);
  double f2 = f1t*(t3-y7-100*y13) + f2t*(z3-x7-100*x13);

  f1 = f1*f1;
  f2 = f2*f2;
  e=e*e/(den-1);

  if (f1+f2-e < 0)
    return -1.;
  else if (f1+f2-e > 0)
    return 1.;
  else
    return 0.;
}


// Knot 3


double knot3_function (double a, double b, double c) {  // radius = 4
  double e=0.025;

  double x, y, z, t, den;
  den=1+a*a+b*b+c*c;
  x=2*a/den;
  y=2*b/den;
  z=2*c/den;
  t=(1-a*a-b*b-c*c)/den;

  double x19=x, y19=y, z17=z, t17=t;
  puiss(x19,y19,19);
  puiss(z17,t17,17);

  double f1 = z17-x19;
  double f2 = t17-y19;

  f1 = f1*f1;
  f2 = f2*f2;
  e=e*e/(den-1);

  if (f1+f2-e < 0)
    return -1.;
  else if (f1+f2-e > 0)
    return 1.;
  else
    return 0.;
}