/* -*-flix-*- */
/*
 * Animation depicting the local isometric deformation of a
 * catenoid to a helicoid through minimal immersions
 */
#include "epix.h"
using namespace ePiX;

// number of meridians and latitudes
int N1=48;
int N2=18;

P VIEWPT(4, 3, 4);

// locations of "lights"; see mesh::draw() below for implementation
P LIGHT_R(10,0,10);
P LIGHT_G(0,10,10);
P LIGHT_B(0,-10,10);

// internal constants
double r_0=1.0;
double EPS=0.00125;
double du=2.0/N1, dv=4.0/N2;

// parametrized surfaces
P helicoid(double u, double v)
{
  return P(sinh(v)*Cos(M_PI*u), sinh(v)*Sin(M_PI*u), 2*M_PI*u);
}

P catenoid(double u, double v)
{
  return P(cosh(v)*Sin(M_PI*u), -cosh(v)*Cos(M_PI*u), -v);
}

P morph(double u, double v)
{
  return Cos(2*M_PI*tix)*helicoid(u,v) + Sin(2*M_PI*tix)*catenoid(u,v);
}

// generalities
class element
{
private:
  P pt1;
  P pt2;
  P pt3;
  P pt4;
  
  double distance;

public:
  element() 
  {
    pt1=P();
    pt2=P();
    pt3=P();
    pt4=P();
    distance=camera.get_range();
  }

  element(P f(double u, double v), double u0, double v0)
  {
    pt1=f(u0+EPS,v0+EPS);
    pt2=f(u0+du-EPS,v0+EPS);
    pt3=f(u0+du-EPS,v0+dv-EPS);
    pt4=f(u0+EPS,v0+dv-EPS);
    
    P center = 0.25*(pt1 + (pt2 + (pt3 + pt4)));
    P temp = camera.get_viewpt();

    distance = norm(center-temp);
  }

  double how_far() const { return distance; }

  void draw() 
  { 
    P normal = (pt2 - pt1)*(pt4 - pt1);
    normal *= 1/norm(normal);

    double dens_r  = 0.75*(pow(normal|LIGHT_R, 2)/(LIGHT_R|LIGHT_R));
    double dens_g  = 0.75*(pow(normal|LIGHT_G, 2)/(LIGHT_G|LIGHT_G));
    double dens_b  = 0.75*(pow(normal|LIGHT_B, 2)/(LIGHT_B|LIGHT_B));

    // white-filled quad
    fill();
    gray(0);
    quad(pt1, pt2, pt3, pt4); 
    fill(false);

    // with boundary colored according to normal shading
    rgb(dens_r, dens_g, dens_b);
    quad(pt1, pt2, pt3, pt4); 
  }

};

class By_distance {
public:
  bool operator() (const element& arg1, const element& arg2)
  {
    return arg1.how_far() > arg2.how_far(); 
  }
};


int main(int argc, char* argv[])
{
  if (argc == 3)
    {
      char* arg;
      double temp1, temp2;
      temp1=strtod(argv[1], &arg);
      temp2=strtod(argv[2], &arg);

      tix=temp1/temp2;
    }

  double MAX=8;
  bounding_box(P(-MAX,-MAX),P(MAX,MAX));
  unitlength("1in");
  picture(P(5,5));

  begin();

  // draw bounding square for uniform frame size
  grid();
  viewpoint(VIEWPT);

  /* rotating lights
  LIGHT_R=cis(2*M_PI*tix);
  LIGHT_G=cis(2*M_PI*(tix+0.25));
  LIGHT_B=P(-2,-2,6);
  */

  crop();
  camera.range(20);

  std::vector<element> mesh(0);

  for (int i=0; i<N1; ++i)
    for (int j=0; j<N2; ++j)
      mesh.push_back(element(morph, -1+du*i, -2+dv*j));

  sort(mesh.begin(), mesh.end(), By_distance());

  pen(1.5);
  for (unsigned int i=0; i<mesh.size(); ++i)
    mesh.at(i).draw();

  end();
}