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

const int M(8);
double MAX(1.5), Bd(2);

P F(double x, double y) // fake pole
{
  return P(x, y, x/(0.01+x*x+y*y));
}

P pt1(-MAX,-MAX), pt2(MAX, MAX);
domain R(pt1, pt2, mesh(2*M,2*M), mesh(6*M,6*M));

int main()
{
  picture(P(-Bd,-Bd), P(Bd,Bd), "2.5x2.5in");

  begin();

  label(P(xmax(), 2.5), P(-2,-2),
        "$z=\\mathrm{Re}\\,\\displaystyle\\frac{1}{x+iy}$", bl);

  camera.at(P(6,8,5));

  clip_box(P(-Bd,-Bd,-Bd), P(Bd,Bd,0));

  plain(Magenta());
  fill(White());
  surface(F, R); // bottom half

  grid(pt1, pt2, M, M); // axes

  pen(White(0.75), "0.15pt");
  plot(F,R);    // simulate transparency

  clip_box(P(-Bd,-Bd,-Bd), P(Bd,Bd,2*Bd)); 
  plain(Magenta());

  // axes and labels
  arrow(P(-Bd,0,0), P(Bd,0,0));
  label(P(Bd,0), P(-4,-2), "$x$", l); 

  arrow(P(0,-Bd,0), P(0,Bd,0));
  label(P(0,Bd), P( 2,-2), "$y$", br);

  arrow(P(0, 0, 0), P(0,0,2.5));
  label(P(0,0,2.5), P(0,4), "$z$", t);

  clip_box(P(-Bd,-Bd,0),P(Bd,Bd,Bd));

  surface(F,R); // top half

  bbold(Red());
  plot(F, R.resize1(0.25,MAX).slice2(0));

  pst_format();
  end();
}