/* -*-ePiX-*- */
// ripple_sm.xp -- October 04, 2004
#include "epix.h"

#define RIPPLE
//#define SADDLE

const int N=20; // (2*N)^2 mesh elements, cannot be made much larger
const int MIN=4; // fabs(x_min/max)

const P VIEWPT=P(18,6,7); // right
// const P VIEWPT=P(18,8,7); // left

const double EPS=EPIX_EPSILON;

// plot this
#ifdef RIPPLE // orange juicer
double f(double u, double v) 
{ 
  return exp(-0.125*(u*u+v*v))*Cos(3*sqrt(v*v+u*u)); 
}

const P LIGHT(0,1,1); // location of "light source"
#endif // RIPPLE

#ifdef SADDLE // dog saddle
double f(double u, double v) 
{ 
  return 0.04*u*v*(-u*u+v*v); 
}

const P LIGHT(0,1,1); // location of "light source"
#endif // SADDLE

// scale [-N,N] to [-MIN, MIN]
double g(int i) { return MIN*i*1.0/N; }

/*
 * gray density of mesh element; use cosine^2 of angle between the surface
 * normal and the light source (0=white, 1=black, so subtract from 1)
 */
double density(int i, int j)
{
  // compute partial derivatives
  double f_x = (f(g(i)+0.5*EPS, g(j)) - f(g(i)-0.5*EPS, g(j)))/EPS;
  double f_y = (f(g(i), g(j)+0.5*EPS) - f(g(i), g(j)-0.5*EPS))/EPS;

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

  // scale [0,1] to [0.125, 0.875]
  return 0.0125+0.625*(1-pow(normal|LIGHT, 2)/(LIGHT|LIGHT));
}

main() {
  picture(6,6);
  bounding_box(P(-4,-4), P(4,4));
  unitlength("0.75in");
  use_pstricks();

  begin();
  viewpoint(VIEWPT);
  use_pstricks(false);

  //  camera.range(20);
  fill();

  // back half of mesh (x < 0)
  for (int i=-N; i < 0; ++i)
    for (int j=-N; j < N; ++j)
      {
	gray(density(i,j));
	quad(P(g(i),   g(j),   f(g(i),   g(j))), 
	     P(g(i+1), g(j),   f(g(i+1), g(j))),
	     P(g(i+1), g(j+1), f(g(i+1), g(j+1))), 
	     P(g(i),   g(j+1), f(g(i),   g(j+1))));
      }

  // left rear quarter of mesh (x > 0 and y < 0)
  for (int i=0; i < N; ++i)
    for (int j=-N; j < 0; ++j)
      {
	gray(density(i,j));
	quad(P(g(i),   g(j),   f(g(i),   g(j))), 
	     P(g(i+1), g(j),   f(g(i+1), g(j))),
	     P(g(i+1), g(j+1), f(g(i+1), g(j+1))), 
	     P(g(i),   g(j+1), f(g(i),   g(j+1))));
      }

  // axes
  bold();
  use_pstricks();
  psset("linecolor=red, linewidth=1pt");
  dart(P(0,0,0), P(1+MIN,0,0));
  dart(P(0,0,0), P(0,1+MIN,0));
#ifdef SADDLE
  dart(P(0,0,0), P(0,0,1+MIN));
#endif // SADDLE

#ifdef RIPPLE
  dart(P(0,0,0), P(0,0,f(0,0)+0.5));
#endif // RIPPLE

  use_pstricks(false);

  plain();

  // from right mesh (x, y > 0)
  for (int i=0; i < N; ++i)
    for (int j=0; j < N; ++j)
      {
	gray(density(i,j));
	quad(P(g(i),   g(j),   f(g(i),   g(j))), 
	     P(g(i+1), g(j),   f(g(i+1), g(j))),
	     P(g(i+1), g(j+1), f(g(i+1), g(j+1))), 
	     P(g(i),   g(j+1), f(g(i),   g(j+1))));
      }

  // axis labels
#ifdef SADDLE
  white();
#endif // SADDLE

  red();
  label(P(1+MIN,0,0), P(-2,2), "$x$", bl);
  label(P(0,1+MIN,0), P(6,-2), "$y$", r);

#ifdef SADDLE
  black();
  label(P(0, 1.5*MIN, -MIN), P(0,0), "$z=-0.04xy(x^2-y^2)$", l);
  red();
  label(P(0,0,1+MIN), P(0,6), "$z$", t);
#endif // SADDLE

#ifdef RIPPLE
  label(P(0,0,f(0,0)+0.5), P(0,6), "$z$", t);
#endif // RIPPLE

  end();
}