/* -*-ePiX-*- */
/* ePiX surface demo image, May 12, 2004 */
#include "epix.h"
using namespace ePiX;

#define COLOR
#define SHOW_YSLICE

// Style parameters
const triple LIGHT(1,0,4); //(-1,-1,-2)

#ifdef COLOR
inline void color_coord(void)  { rgb(0.7, 0.9, 0.5); }
inline void color_axis(void)   { rgb(0.8, 0.2, 0.9); } 
inline void color_xslice(void) { red(); }
inline void color_yslice(void) { blue(); }
#else
inline void color_coord(void)  { rgb(0.7, 0.7, 0.7); }
inline void color_axis(void)   { black(); } 
inline void color_xslice(void) { black(); }
inline void color_yslice(void) { rgb(0.4,0.4,0.4); }
#endif

const int MESH=8; // number of coordinate grid squares

// location of tangency point
const double x_0=0.625;
const double y_0=0.5;

const double z_0=0.25; // height of top of slicing planes

// function to be graphed
double f(double x, double y) { return 0.5*y*(y*y-3*x*x); }

// Auxiliary internal parameters
const int N=2*MESH; // number of surface mesh subdivisions
const int MAX=1; // maximum coordinate

// constant y slice of the graph of f
inline triple Df_x(double t) { return P(t+x_0, y_0, f(t+x_0, y_0)); }

// constant x slice of the graph of f
inline triple Df_y(double t) { return P(x_0, t+y_0, f(x_0, t+y_0)); }

// shading density function and helpers
const double EPS=0.001; // for approximating derivatives
double g(int i) { return MAX*i*1.0/N; } // int -> grid coordinate

// 0.5*(1-cos^2(t)), t=angle from camera to surface normal
double density(int i, int j)
{
  // compute unit surface normal
  double f_x = 0.5*(f(g(i)+EPS, g(j)) - f(g(i)-EPS, g(j)))/EPS;
  double f_y = 0.5*(f(g(i), g(j)+EPS) - f(g(i), g(j)-EPS))/EPS;

  triple normal = P(-f_x, -f_y, 1);
  normal *= 1/sqrt(normal|normal);

  // convert cos^2 to a gray density
  //  return  0.5*(1-pow(normal|VIEWPT, 2)/(VIEWPT|VIEWPT));
  return  0.625*sgn(normal|LIGHT)*(1-pow(normal|LIGHT, 2)/(LIGHT|LIGHT));
}

// mesh location -> draw surface element
inline void draw_mesh_quad(int i, int j)
{
  gray(density(i,j));
  P p1 = P(g(i),   g(j),   f(g(i),   g(j)));
  P p2 = P(g(i+1), g(j),   f(g(i+1), g(j)));
  P p3 = P(g(i+1), g(j+1), f(g(i+1), g(j+1)));
  P p4 = P(g(i),   g(j+1), f(g(i),   g(j+1)));

  path temp = polygon(4, &p1, &p2, &p3, &p4);
  temp.set_fill();
  temp.draw();
}


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

  begin();
  // camera location (in spherical coordinates); must be in first orthant
  const triple VIEWPT=sph(4, M_PI/6, M_PI/6);

  viewpoint(VIEWPT);
  clip_box(P(-MAX,-MAX,-0.75*MAX), P(MAX,MAX,0.75*MAX));

  // coordinate grids
  color_coord();
  grid(P(-MAX,-MAX,-MAX), P(-MAX, MAX, MAX), mesh(MESH,MESH), mesh(1,1));
  grid(P(-MAX,-MAX,-MAX), P( MAX,-MAX, MAX), mesh(MESH,MESH), mesh(1,1));
  grid(P(-MAX,-MAX,-MAX), P( MAX, MAX,-MAX), mesh(MESH,MESH), mesh(1,1));

  color_axis();

  label(P(0,0.25+MAX,0), P(4,0), "$y$", r);
  label(P(0,0,0.25+MAX), P(0,4), "$z$", t);

  clip();
  // back 1/2 of surface (x<0), drawn back-to-front
  for (int i=-N; i < 0; ++i)
    for (int j=-N; j < N; ++j)
      draw_mesh_quad(i,j);
  clip(false);

  // coordinate axes
  bold();
  color_axis();
  line(P(0,-MAX,0), P(0,0.25+MAX,0));
  line(P(0,0,0), P(0,0,0.25+MAX));

  const P pt1 = P(x_0, 0, z_0);
  const P pt2 = P(x_0, 0, -MAX);
  const P pt3 = P(x_0, MAX, -MAX);

  const P pt4 = P(0, y_0, z_0);
  const P pt5 = P(0, y_0, -MAX);
  const P pt6 = P(MAX, y_0, -MAX);

  const P pt7 = P(x_0, 0, 0);
  const P pt8 = P(x_0, MAX, z_0);

  const P pt10 = P(0, y_0, f(0,y_0));
  const P pt11 = P(MAX, y_0, z_0);

  // lower outlines of slicing planes (partially hidden by surface)
  color_xslice();
  polyline(3, &pt1, &pt2, &pt3).draw();
#ifdef SHOW_YSLICE
  color_yslice();
  polyline(3, &pt4, &pt5, &pt6).draw();
#endif

  plain();
  clip();
  // front 1/2 of surface (x>0)
  for (int i=0; i < N; ++i)
    for (int j=-N; j < N; ++j)
      draw_mesh_quad(i,j);
  clip(false);

  // last coordinate axis
  bold();
  color_axis();
  line(P(-MAX,0,0), P(0.25+MAX,0,0));
  label(P(0.25+MAX,0,0), P(-2,-2), "$x$", bl);

  // slices of the graph and labels
  pen(1.5);
#ifdef SHOW_YSLICE
  color_yslice();
  plot(Df_x, -x_0, MAX-x_0, 20);
  label(Df_x(MAX-x_0), P(-4,-8), 
	"$\\displaystyle\\frac{\\partial f}{\\partial x}$: $y$ constant", l);
#endif

  color_xslice();
  plot(Df_y, -y_0, MAX-y_0, 20);
  label(Df_y(MAX-y_0), P(4,-8), 
	"$\\displaystyle\\frac{\\partial f}{\\partial y}$: $x$ constant", r);

  // overall label
  black();
  label(P(-MAX, 0, MAX), P(4,8), "$z=\\frac{1}{2}y^3-3x^2y$", r);

  // fake transparency: redraw slicing planes thinly
  pen(0.15);
  color_xslice();
  polyline(3, &pt1, &pt2, &pt3).draw();
#ifdef SHOW_YSLICE
  color_yslice();
  polyline(3, &pt4, &pt5, &pt6).draw();
#endif

  // upper outlines of slicing planes (not hidden by surface) 
  bold();
  color_xslice();
  polyline(4, &pt7, &pt1, &pt8, &pt3).draw();

#ifdef SHOW_YSLICE
  color_yslice();
  polyline(4, &pt10, &pt4, &pt11, &pt6).draw();
#endif

  end();
}