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

double f(double x, double y)
{ 
  return 2+0.075*(Sin(x*x-y) + 0.25*Cos(x*x+y*y)); 
}

// Half-size of bounding box
const int N=4;
const int blocks = 4;
const int coarse = 20;
const int fine = 120;

main() {
  bounding_box(P(-N,-N), P(N,N));
  unitlength("1in");
  picture(4,4);

  begin();
  viewpoint(-10,-30,20);

  camera.range(10);

  domain R(P(x_min,y_min), P(x_max,y_max), 
	   mesh(coarse,coarse), mesh(fine,fine));

  // bottom grid
  green();
  pen(0.15); // fine mesh
  grid(P(x_min,y_min,-1), P(x_max, y_max, -1), coarse, coarse);

  plain();
  blue(); // coarse mesh
  grid(P(x_min,y_min,-1), P(x_max, y_max, -1), 
       mesh(blocks, blocks), mesh(coarse, coarse));

  // vertical segments
  red();
  for (int i=0; i <= blocks; ++i)
    {
      double xi = x_min + i*x_size/blocks;
      for (int j=0; j <= blocks; ++j)
	{
	  double yj = y_min + j*y_size/blocks;

	  line(P(xi, yj, -1), P(xi, yj, f(xi, yj)));
	}
    }

  // erase layer...
  white();
  pen(3);
  plot(f, R);

  black();
  h_axis_labels(P(x_min,y_min,-1), P(x_max,y_min,-1), blocks, P(-6,-12));
  v_axis_labels(P(x_min,y_min,-1), P(x_min,y_max,-1), blocks, P(-18,-4));

  // function graph
  plain();
  plot(f, R);

  magenta();
  pen(0.5);
  plot(f, R);

  end();
}