/* -*-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(72);
int N2(24);

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.0); // fact shrink factor
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);
}

// facet-like class with spot reflection
class element
{
private:
  P pt1;
  P pt2;
  P pt3;
  P pt4;
  
  double distance;

public:
  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.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)));

    fill(RGB(dens_r, dens_g, dens_b));
    ePiX::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);
  picture(P(-MAX,-MAX),P(MAX,MAX), "5x5in");

  begin();
  pen(Neutral(), 0); // no grid lines

  // draw bounding square for uniform frame size
  backing(Black());
  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);
  */

  camera.range(20);

  // build surface
  std::vector<element> mesh;

  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());

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

  pst_format();
  end();
}