/* -*-ePiX-*- */
#include <iostream>
#include <vector>
#include "epix.h"
using namespace ePiX;

const int N1=24;
const int N2=48;

const double du=1.0/N1, dv=1.0/N2;
const P LIGHT(1,2,2);
// const P VIEWPT(4, 3, 4);

const double r_0=0.95;

const double EPS=0.00125;

const double k = 2*M_PI/(360*sqrt(3)); // assume "degrees" mode

class mesh_quad
{
private:
  P pt1;
  P pt2;
  P pt3;
  P pt4;
  
  double distance;

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

  mesh_quad(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  = 0.125+0.75*(1-pow(normal|LIGHT, 2)/(LIGHT|LIGHT));
    double dens  = 0.75*(1-pow(normal|LIGHT, 2)/(LIGHT|LIGHT));

    gray(dens); 
    quad(pt1, pt2, pt3, pt4); 
  }

};

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

P f1(double u, double v)
{
  return polar(2+r_0*Cos(u), v) + P(0,0,r_0*Sin(u)); // torus
}

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;
    }

  bounding_box(P(-3,-3),P(3,3));
  unitlength("1in");
  picture(4,4);

  begin();
  revolutions();

  rgb(0.95,0.95,0.95);
  grid();

  //  viewpoint(VIEWPT);

  camera.range(10);
  camera.rotate_eye(0.25);
  camera.rotate_sea(0.0625);

  //  crop();
  clip_box(P(-3,-3,-1), P(3*Sin(tix),3,1));
  clip();
  fill();
  //  camera.range(20);

  std::vector<mesh_quad> torus(0);

  for (int i=0; i<N1; ++i)
    for (int j=0; j<N2; ++j)
      torus.push_back(mesh_quad(f1, du*i, dv*j));

  sort(torus.begin(), torus.end(), by_distance());

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

  fill(false);
  clip(false);
  pen(0.15);
  red();
  plot(f1, domain(P(0,0), P(1,1), mesh(24,48), mesh(24,48)));

  end();
}