/* -*-ePiX-*- */
#include "epix.h"
using namespace ePiX;

const double MAX=3;

const double R_0=4; // radius of strip
const double r_0=1; // width of strip

double param(double u) 
{ 
  double temp = Acos(Cos(u));
  
  return pow(temp*(1-temp), 20); 
}

I twist(param); // definite integral class

const double c_0=1/twist.eval(1);

P moebius(double u, double v)
{
  double th=0.5+c_0*twist.eval(u);

  return polar(R_0+v*Cos(0.5*th), u) + v*P(0,0,Sin(0.5*th));
}

// for erasure
domain R0(P(0.25,-r_0), P(0.75,r_0), mesh(36,48), mesh(72,1));

// back and front halves of strip
domain R1(P(-0.25,-r_0), P(0.25,r_0), mesh(36,8), mesh(72,1));
domain R2(P( 0.25,-r_0), P(0.75,r_0), mesh(36,8), mesh(72,1));

int main() {
  bounding_box(P(-MAX,-MAX), P(MAX,MAX));
  unitlength("1in");
  picture(4,4);

  begin();
  revolutions();

  viewpoint(P(-100,0,20));

  plot(moebius, R1);

  bold();
  red();
  ellipse(P(0,0), R_0*E_1, R_0*E_2, -0.25, 0.25, 72);

  // hidden line removal
  pen(3);
  white();
  plot(moebius, R0);

  pen(2);
  red();
  ellipse(P(0,0), R_0*E_1, R_0*E_2, 0.25, 0.75, 72);

  bold();
  black();
  plot(moebius, R2);

  end();
}