%% Sample LaTeX file for ePiX                July, 2002
%%
%% Andrew D. Hwang            ahwang@mathcs.holycross.edu
%%
\documentclass[12pt]{article}
\usepackage[leqno]{amsmath}
\usepackage{latexsym,pstcol,epic,eepic}

\def\ePiX{{\tt ePiX}}
\def\C{{\tt C}}
\def\R{{\bf R}}

\DeclareMathOperator{\re}{Re}
\DeclareMathOperator{\im}{Im}

\setlength{\textheight}{9.5in}
\setlength{\textwidth}{6.5in}
\setlength{\oddsidemargin}{0in}
\setlength{\topmargin}{0in}

\begin{document}

\section{Overview}

\ePiX\ is a powerful, flexible, lightweight runtime utility for
creating mathematically accurate 2-dimensional \LaTeX\ plots and
figures from simple, mnemonic commands. A detailed user's manual is
available in several formats from
\begin{verbatim}
http://mathcs.holycross.edu/~ahwang/current/ePiX.html
\end{verbatim}
This sample document demonstrates what \ePiX\ can do, with
side-by-side comparisons of input files\footnote{The sample files
distributed with \ePiX\ have {\tt offset} lines to place the figures
next to their source code.} and corresponding output.

\subsection{Line Figures}

\ePiX\ was originally designed for function plotting, but contains
provisions for creation and manipulation of ``standard'' plot objects,
including arbitrary triangles, coordinate rectangles, splines, and
half-ellipses.

Figure~\ref{fig:golden_rect} is very basic, but nonetheless illustrates
useful features; particularly, sizes and locations are given
symbolically, in terms of the golden ratio, and the arc is specified
by center and angle, so the ends are placed precisely.

\begin{figure}[hbt]
\begin{center}
\input{golden_rect.eepic}
\vspace*{-2.00in}
\begin{small}
\begin{verbatim}
#include "epix.h"

main() 
{
  double tau = (1+sqrt(5))/2;

  bounding_box(P(0,0), P(tau,1));
  picture(P(100*(1+tau), 100*tau));
  unitlength("0.0125in");

  begin();

  dashed;
  line(P(0.5,0), P(0.5,1));
  solid;
  line(P(1,0), P(1,1));
  line(P(0.5,0), P(1,1));

  arc(P(0.5,0), tau-0.5, 0, atan(2));
  boldrect(P(0,0), P(tau, 1));

  end();
}
\end{verbatim}
\end{small}
\caption{Constructing a golden rectangle.}
\label{fig:golden_rect}
\end{center}
\end{figure}

\clearpage

With variables and {\tt for} loops, objects can be placed regularly,
accurately, and in an easily-modified way. For example, to change the
number of rectangles in the figure below, only the definition of~{\tt
N} needs to be changed. Variables in \ePiX\ can be used like macros
and user-defined commands in \LaTeX, to express a figure in terms of
its logical structure.

\begin{figure}[hbt]
\begin{center}
\input{geomsum.eepic}
\vspace*{-192pt}
%%\begin{small}
\begin{verbatim}
#include <epix.h>
#define N 8

main() 
{
  double t;

  unitlength("1.5pt");
  bounding_box(P(0,0), P(1,1));
  picture(P(128,128));

  begin();   

  for(int i=0; i<8; ++i)
    {
      t = pow(0.5, i+1); // 2^-(i+1)
      line(P(1-t, 1-2*t), P(1-t, 1));
      line(P(1-t, 1-t), P(1, 1-t)); 
    }
  bold;
  rect(P(x_min, y_min), P(x_max, y_max));

  end();
}
\end{verbatim}
%\end{small}
\caption{A geometric series with ratio $1/2$.}
\label{fig:geomsum}
\end{center}
\end{figure}

\ePiX\ now uses the \texttt{C++} compiler, so variables need not be
declared at the start of a file; the loop index is declared in the
loop itself. It is a good idea to declare ``tweakable'' variables in
the preamble, but indices can generally be hidden safely.

\clearpage

The ``denominator'' function~$f$ is defined by
$$
f(x)=\begin{cases}
\frac{1}{q} & \text{if $x=\frac{p}{q}$ in lowest terms} \\
0 & \text{if $x$~is irrational}
\end{cases}
$$

\begin{figure}[hbt]
\begin{center}
\input{denom.eepic}
\vfil

\begin{verbatim}
#include "epix.h"

#define N 20

main() 
{
  unitlength("0.0125in");
  bounding_box(P(-2,0), P(2,1));
  picture(P(400,100));

  begin();   

  h_axis(P(x_min, 0), P(x_max, 0), 2*x_size);
  v_axis(P(0, 0), P(0, y_max), 2);

  h_axis_labels(P(x_min, 0), P(x_max, 0), x_size, P(-12, -12));

  for (int i=1; i< N; ++i) 
    for (int j=i*x_min; j <= i*x_max; ++j) 
      if (gcd(i, j) == 1) 
        marker(P(((double)j)/i, 1.0/i), DOT);

  end();
}
\end{verbatim}
\caption{The graph of the denominator function.}
\label{fig:denom}
\end{center}
\end{figure}

\clearpage

In \ePiX, points in the plane are specified by an ordered pair data
structure; \ePiX\ provides functions for adding/subtracting, scalar
multiplying, and complex multiplying pairs, and for creating pairs in
both Cartesian and polar coordinates. 

\begin{figure}[hbt]
\begin{center}
\input{demoivre.eepic}
\vspace*{-75pt}
\begin{verbatim}
#include "epix.h"
int n=24;

main()
{
  pair power = e_1; // (1,0)
  pair alpha = e_1 + (M_PI/n)*e_2;

  bounding_box(P(-1.5,0), P(1,1.25));
  picture(P(200,100));
  unitlength("1pt");

  begin();
  
  for(int i=0; i<n; ++i)
    {
      line(P(0,0), alpha*power); // draw two sides
      line(power, alpha*power);  // of triangle
      power *= alpha;            // and increment power
    }

  label(alpha, P(2,0), "$\\alpha=1+\\frac{i\\pi}{n}$", r);
  label(power, P(-1,0), "$\\alpha^n$", l);

  bold;
  red;
  triangle(P(0, 0), e_1, alpha);

  end();
}
\end{verbatim}
\caption{A geometric interpretation of Euler's formula $e^{i\pi}+1=0$.}
\label{fig:demoivre}
\end{center}
\end{figure}

\clearpage

Figure~\ref{fig:contour} uses variables to control the appearance of a
figure.  Inner and outer radii, arrow and label locations, and points
of intersection are all computed in the figure itself, rather than
being hard-wired. Reverse \LaTeX-style label alignment is used, which
is easier than plain \texttt{offset}s in this example. To rescale the
figure, only the radius and/or width parameters need be changed.

\begin{figure}[hbt]
\begin{center}
\input{contour.eepic}
\vspace*{-6cm}
\begin{footnotesize}
\begin{verbatim}
#include "epix.h"

double pi=M_PI, theta=1;
double Rad=4.5; // outer radius
double rad=0.5; // inner radius
double D=0.125;  // half-width of slit

main() 
{
  bounding_box(P(-5, -5), P(5,5));
  picture(P(160, 160));
  unitlength("0.35mm");

  begin();

  spot(P(0,0)); 

  // Circular arcs
  arc(P(0,0), Rad, asin(D/Rad), 2*M_PI-asin(D/Rad));
  arc(P(0,0), rad, asin(D/rad), 2*M_PI-asin(D/rad));

  // Slit
  line(P(rad*cos(asin(D/rad)),  D), P(Rad*cos(asin(D/Rad)),  D));
  line(P(rad*cos(asin(D/rad)), -D), P(Rad*cos(asin(D/Rad)), -D));

  // Arrows
  arrow(P(  Rad/4,  rad), P(3*Rad/4,  rad));
  arrow(P(3*Rad/4, -rad), P(  Rad/4, -rad));

  arc_arrow(P(0,0), 0.9*Rad, pi-theta, pi+theta);

  label(P(Rad,D), P(2,4), "$R\\to\\infty$", tr);
  label(P(0,rad), P(0,4), "$\\delta\\to0$", tl);
  label(polar(Rad, M_PI_4), P(0,0), "$\\gamma$", tr);

  end();
}
\end{verbatim}
\end{footnotesize}
\caption{A keyhole contour for $\displaystyle\int_0^\infty f(z)\,dz$.}
\label{fig:contour}
\end{center}
\end{figure}

\clearpage

\subsection{Function Plotting}

\ePiX\ provides Cartesian and polar plotting, data plotting from a
file, and parametric plotting with a few choices of plot attributes,
such as clipping the plot to the bounding box or plotting one or more
curves in a family. Figure~\ref{fig:clipplot} illustrates {\tt
clipplot} and reverse \LaTeX\ alignment of axis labels. (\texttt{tl}
and \texttt{lt} align the labels above and left of their Cartesian
location.) 

\begin{figure}[hbt]
\begin{center}
\input{clipplot.eepic}
\vspace*{-200pt}
\begin{footnotesize}
\begin{verbatim}
#include <epix.h>

double f(double t)
{
  return 2*t*(1-t)*(1-t);
}

double g(double t)
{
  return 1/(1-t*t);
}

main() {

  bounding_box(P(-2,-4), P(2,4));
  picture(P(200,200));
  unitlength("1pt");

  begin();

  /* Vertical asymptotes */
  dashed;
  line(P(-1, y_min), P(-1, y_max));
  line(P( 1, y_min), P( 1, y_max));

  /* Axes */
  solid;
  h_axis(P(x_min, 0), P(x_max, 0), 8);
  v_axis(P(0, y_min), P(0, y_max), 8);

  h_axis_labels(P(x_min, 0), P(x_max, 0), 4, P(-1, 2), tl);
  v_axis_labels(P(0, y_min), P(0, y_max), 4, P(-1, 2), lt);

  /* Graphs */
  clipplot(f, x_min, x_max, 40);
  boldclipplot(g, x_min, x_max, 80);

  end();
}
\end{verbatim}
\end{footnotesize}
\caption{{\tt clipplot} can be used to graph functions with poles.}
\label{fig:clipplot}
\end{center}
\end{figure}

\clearpage

\ePiX\ contains primitives for Cartesian and polar grids, as well as
for polar plotting. Note that angles in polar plots are measured in
revolutions rather than in radians. It is prudent to make the number
of points plotted (120~here) divisible by the order of symmetry~(3).

\begin{figure}[hbt]
\begin{center}
\input{polar.eepic}
\vspace*{-150pt}
\begin{verbatim}
#include <epix.h>

double f(double t)
{
  return 2*cos(3*t);
}

main() 
{
  unitlength("1pt");
  bounding_box(P(-2,-2), P(2,2));
  picture(P(180, 180));

  begin();

  h_axis(P(x_min, y_min), P(x_max, y_min), 8);
  v_axis(P(x_min, y_min), P(x_min, y_max), 8);
  
  polar_grid(2, 4, 24); // radius, rings, sectors

  h_axis_labels(P(x_min,y_min), P(x_max,y_min), x_size, P(-12,-14));
  v_axis_labels(P(x_min,y_min), P(x_min,y_max), y_size, P(-4,0), l);

  bold;
  polarplot(f, 0, 0.5, 120);

  end();
}
\end{verbatim}
\caption{The polar graph $r=2\cos 3\theta$ for $0\leq\theta\leq\pi$.}
\label{fig:polar}
\end{center}
\end{figure}

\clearpage

A loop index may be used to control an entire figure, generating a
sequence of snapshots of a time-varying picture, such as a rolling
wheel, or the flow of an ODE.

\begin{footnotesize}
\begin{verbatim}
#include "epix.h"

pair F(double t, double r)
{
  return P(t,0) - r*P(sin(t), cos(t));
}

main() 
{
  double dt = 5*M_PI/11;
  double t = 0;

  picture(P(420, 120));
  bounding_box(P(0,-1), P(7,1));
  unitlength("0.005in");

  for(int i=0; i < 9; ++i)
    {
      t += dt;
      begin(); // Entire picture inside loop body

      line(P(x_min - 2, y_min), P(x_max + 6, y_min)); // the ground
      
      ellipse(P(t,0), P(1,1)); // the wheel

      /* the paths */
      bold;
      for (int j=0; j < 6; ++j)
	{
	  rgb(1-0.125*j, 0.125*j, 0.5+0.25*j);
	  sliceplot2(F, 0, t, 0.2*j, (int) ceil(1+10*t));
	}

      /* the spoke */
      green;
      line(P(t,0), F(t,1));

      end();
    
      /* figure separator and vertical space */
      printf("\n\\vspace*{3ex}\n%%%%");
    }
}
\end{verbatim}
\end{footnotesize}

\vspace*{-100pt}
\begin{figure}[hb]
\begin{center}
\input{wheel.eepic}
\caption{Snapshots of cycloids.}
\label{fig:wheel}
\end{center}
\end{figure}

\clearpage

\noindent Because an \ePiX\ file is a \C\ program, it is easy to plot
finite sums of series, as in Figure~\ref{fig:weierstrass}. The
function being scaled and summed is ${\rm cb}$, the ``Charlie Brown''
function (blue). Nowhere differentiability is essentially obvious from
the picture, because the graph is self-similar (red), and is not a
line.

\begin{figure}[hbt]
\begin{center}
\input{weierstrass.eepic}
\vspace*{-1.5in}
\begin{footnotesize}
\begin{verbatim}
#include "epix.h"
#define N 8 // Number of summands

/* "cb" is the Charlie Brown function, see functions.c */

double weierstrass(double t)
{
  double y=0;
  for(int i=0; i < N; ++i)
    y += pow(2,-i)*cb(pow(2,i)*t);

  return y;
}

main() {

  bounding_box(P(-2, 0), P(2, 1.5));
  picture(P(320, 120));
  unitlength("0.01in");

  begin();

  h_axis(P(x_min,0), P(x_max,0), 2*x_size);
  v_axis(P(0,y_min), P(0,y_max), 2*y_size);
  h_axis_labels(P(x_min,0), P(x_max,0), x_size, P(-12, -14));
  
  blue;
  plot(cb, x_min - 0.25, x_max+0.25, 4*x_size + 2);

  bold;
  black;
  plot(weierstrass, x_min, x_max, pow(2, N));

  red;
  plot(weierstrass, 0.5, 1.5, pow(2, N-2));

  end();
}
\end{verbatim}
\end{footnotesize}
\caption{A Weierstrass nowhere-differentiable function.}
\label{fig:weierstrass}
\end{center}
\end{figure}

\clearpage

\ePiX\ can approximate the extreme values of a function on an
interval, which is useful for drawing inscribed or circumscribed
rectangles in a graph. The sine function is symbolized by~\texttt{f}
throughout the body (except the label, of course), so re-drawing the
figure with a different integrand requires changing only the
\texttt{define} line and re-sizing the \verb+bounding_box+.  The
number of rectangles is similarly ``parametrized''.

\begin{figure}[hbt]
\begin{center}
\input{uppersum.eepic}
\vspace*{-1in}
\begin{footnotesize}
\begin{verbatim}
#include "epix.h"
#define N 12.0 // Number of rectangles
#define f sin  // gather references to integrand

main() 
{
  picture(P(300, 100));
  bounding_box(P(0,0), P(3,1));
  unitlength("0.01in");

  begin();

  double dx = x_size/N;

  for(int i=0; i < N; ++i)
    {
      double ai=x_min + i*dx;
      double bi=x_min + (i+1)*dx;

      rect(P(ai, 0), P(bi, sup(f, ai, bi)));
    }

  h_axis(P(x_min, y_min), P(x_max, y_min), x_size);
  v_axis(P(x_min, y_min), P(x_min, y_max), 2*y_size);

  h_axis_labels(P(x_min, 0), P(x_max, 0), x_size, P(0,-4), b);
  label(P(x_max, f(x_max)), P(4,0), "$y=\\sin x$", r);

  bold;
  blue;
  plot(f, x_min, x_max, 40);

  end();
}
\end{verbatim}
\end{footnotesize}
\caption{Upper sums for $\displaystyle\int_0^3\sin x\,dx$.}
\label{fig:uppersum}
\end{center}
\end{figure}

\clearpage

Simple level sets may be plotted as parametrized paths using {\tt
sliceplot}, the restriction of a mapping $T:\R^2\to\R^2$ to horizontal
or vertical lines; \ePiX\ does not do implicit plotting, and probably
never will because of inherent instability near critical points and
the resulting quality of output. 

\begin{figure}[hbt]
\begin{center}
\input{levelset.eepic}
\vspace*{-2.75in}
\begin{footnotesize}
\begin{verbatim}
#include "epix.h"
#define MAX 3

pair F1(double t, double r) {
  return r*P(sinh(t), cosh(t));
}
pair F2(double t, double r) {
  return r*P(cosh(t), sinh(t));
}

main()
{
  int i;

  bounding_box(P(-(MAX),-(MAX)), P((MAX), (MAX)));
  unitlength("0.625in");
  picture(P(4,4));

  begin();

  bold;
  rgb(0.5,0.5,1);
  line(P(x_min, y_min), P(x_max, y_max));
  line(P(x_min, y_max), P(x_max, y_min));

  for (i=1; i <= (MAX)*(MAX); ++i) 
    {
      double t0 = acosh((MAX)/sqrt(i));
      double D = i*i/pow(MAX, 4); //((MAX)*(MAX));

      rgb(0.5+0.75*D, 0.5-0.25*D, 1-D);
      sliceplot2(F1, -t0, t0,  sqrt(i), 40);
      sliceplot2(F1, -t0, t0, -sqrt(i), 40);

      rgb(0.5-0.5*D, 0.5+0.25*D, 1-0.25*D);
      sliceplot2(F2, -t0, t0,  sqrt(i), 40);
      sliceplot2(F2, -t0, t0, -sqrt(i), 40);
  }

  plain;
  red;
  rect(P(x_min, y_min), P(x_max, y_max));
  end(); 
}
\end{verbatim}
\end{footnotesize}
\end{center}
\end{figure}

\clearpage

In addition to handling line objects in space, \ePiX\ plots space
curves; even text labels may be placed in their correct spatial
location. Pairing two figures with slightly different {\tt viewpoint}s
gives a cross-your-eyes stereogram. The source code for the left frame
is shown; the only essential difference between the frames is their
respective {\tt viewpoint}s. Because PostScript builds figures in
layers, and because \ePiX\ does not do hidden line removal, the source
file must be constructed with some thought.

\begin{figure}[hbt]
\begin{center}
\input{stereo-right.eepic}\qquad
\input{stereo-left.eepic}
\caption{The twisted cubic and its coordinate plane projections.}
\label{fig:stereo}
\end{center}
\end{figure}

\begin{verbatim}
#include <epix.h>

double   u1(double t) {return t;}
double   u2(double t) {return t*t;}
double   u3(double t) {return t*t*t;}
double zero(double t) {return 0;}

int N=60;

main() 
{
  int i;
  double t;

  unitlength("1pt");
  bounding_box(P(-1,-1), P(1,1));
  picture(P(150,150));

  viewpoint(1, 1.8, 0.5);

  begin();

  // coordinate axes
  boldarrow(V(-1,0,0), V(1,0,0));
  boldarrow(V(0,0,-1), V(0,0,1));

  bold;
  blue;
  plot(u1, zero, u3, -1, 1, N);
  for (i=0; i<=20; ++i)
    {
      t = -1 + i/10.0;
      line(V(u1(t), 0, u3(t)), V(u1(t), u2(t), u3(t)));
    }

  label(V(1,0,1), P(2,2), "$x$-$z$ plane");

  red;
  plot(zero, u2, u3, -1, 1, N);
  for (i=0; i<=20; ++i)
    {
      t = -1 + i/10.0;
      line(V(0, u2(t), u3(t)), V(u1(t), u2(t), u3(t)));
    }

  label(V(0,1,1), P(2,2), "$y$-$z$ plane");

  black;
  boldarrow(V(0,0,0), V(0,1,0));

  magenta;
  plot(u1, u2, zero, -1, 1, N);
  for (i=0; i<=20; ++i)
    {
      t = -1 + i/10.0;
      line(V(u1(t), u2(t), 0), V(u1(t), u2(t), u3(t)));
    }

  label(V(1,1,0), P(-45,-12), "$x$-$y$ plane");
  black;

  /* The cubic */
  plot(u1, u2, u3, -1, 1, N);

  end();
}
\end{verbatim}

\clearpage

Figure~\ref{fig:sphere} demonstrates the spherical geometry abilities
of \ePiX: Space curves can be projected to the unit sphere radially,
and plane curves can be projected stereographically from the north or
south pole, with or without hidden line removal.

\begin{figure}[hbt]
\begin{center}
\input{sphere.eepic}
\vspace*{-220pt}
\begin{footnotesize}
\begin{verbatim}
#include "epix.h"
static double zero (double t) {return 0;}
triple phi1(double t) {return P( sin(t), cos(t), t);}
triple phi2(double t) {return P(-cos(t), sin(t), t);}

main() 
{
  bounding_box(P(-1, -1), P(1,1));
  picture(P(160, 160));
  unitlength("1.25pt");
  viewpoint(1.5, 2, 1);

  begin();
  ellipse(P(0,0), P(1,1));

  dotted;
  backplot_R(cos, sin, zero, 0, 2*M_PI, 60);
  backplot_R(zero, cos, sin, 0, 2*M_PI, 60);
  backplot_R(sin, zero, cos, 0, 2*M_PI, 60);

  solid;
  frontplot_R(cos, sin, zero, 0, 2*M_PI, 60);
  frontplot_R(zero, cos, sin, 0, 2*M_PI, 60);
  frontplot_R(sin, zero, cos, 0, 2*M_PI, 60);

  red;
  pen(0.15);
  dashed;
  backplot_R(phi1, -4*M_PI, 6*M_PI, 240);
  backplot_R(phi2, -4*M_PI, 6*M_PI, 240);

  bold;
  solid;
  frontplot_R(phi1, -4*M_PI, 6*M_PI, 240);
  frontplot_R(phi2, -4*M_PI, 6*M_PI, 240);

  end();
}
\end{verbatim}
\end{footnotesize}
\caption{Helices project to loxodromes on the unit sphere.}
\label{fig:sphere}
\end{center}
\end{figure}

\clearpage

\subsection{Calculus Plots}

\ePiX\ can graph derivatives, integrals, slope and vector fields, and
can solve ODEs.

\begin{figure}[hbt]
\begin{center}
\input{calculus.eepic}
\vspace*{-120pt}
\begin{footnotesize}
\begin{verbatim}
#include "epix.h"
double MAX=2*M_PI;
double f(double t)
{
  return t*sin(t);
}

main() {

  unitlength("1pt");
  picture(P(240, 120));
  bounding_box(P(-MAX,-MAX), P(MAX,MAX));

  begin();

  // Coordinate axes and labels
  h_axis(P(x_min,0), P(x_max,0), 8);
  v_axis(P(0,y_min), P(0,y_max), 4);

  label(P(0,y_max), P(-2,0), "$\\phantom{-}2\\pi$", l);
  label(P(0,y_min), P(-2,0), "$-2\\pi$", l);

  label(P(x_min,0), P(-4,2), "$-2\\pi$", t);
  label(P(x_max,0), P( 0,2), "$2\\pi$", t);

  bold;
  plot(f, x_min, x_max, 90);
  green;
  plot_deriv(f, x_min, x_max, 90);

  /* Definite integral from x=0 */
  blue;
  plot_int(f, 0, x_min, x_max, 90);

  end();
}
\end{verbatim}
\end{footnotesize}
\caption{$y=x\sin x$ (black), its derivative (green) and integral from~$0$
(blue).}
\label{fig:calculus}
\end{center}
\end{figure}

\clearpage

Solutions of planar ODEs are computed with Euler's method. Curves
terminate if they leave the \verb+bounding_box+. 

\begin{figure}[ht]
\begin{center}
\input{slopefield.eepic}
\vspace*{-180pt}
%%\begin{footnotesize}
\begin{verbatim}
#include "epix.h"
pair F(double s, double t) 
{
  pair temp = P(s,t);
  return (P(0.1,0.025)&temp) 
    + (1/(0.01+s*s+t*t))*J(temp);
}

main() 
{
  unitlength("1pt");
  bounding_box(P(-4, -3), P(2,2));
  picture(P(234, 195));

  begin();
  dart_field(F, P(x_min, y_min), P(x_max, y_max), 4*x_size, 4*y_size);

  bold;
  for (int i=0; i<7; ++i)
    {
      rgb(0.25 - 0.05*i, 1 - 0.1*i, 0.2*i);
      ode_plot(F, P(-0.9-0.025*i,0), 20, 200);
    }
  end();
}
\end{verbatim}
%%\end{footnotesize}
\caption{A slope field and six solutions of an ODE.}
\label{fig:slopefield}
\end{center}
\end{figure}

\clearpage

\noindent \verb+tan_field+ graphs the field of tangent vectors along a
parametrized path; the path itself is drawn separately. The choice of
contrasting colors emphasizes places where the curvature is
large. This figure again demonstrates how PostScript builds up a
figure in layers; in particular, the axis labels are ``masked''. Note
that a \texttt{grid} is by default the size of the \verb+bounding_box+.

\begin{figure}[hbt]
\begin{center}
\input{lissajous.eepic}
\vspace*{-3in}
\begin{small}
\begin{verbatim}
#include "epix.h"

pair f(double t) 
{
  return P(sin(3*M_PI*t+0.5), sin(4*M_PI*t+0.5));
}

main() {

  bounding_box(P(-1,-1), P(1,1));
  picture(P(300, 300));
  unitlength("0.01in");

  begin();

  grid(4*x_size, 4*y_size);

  bold;
  line(P(x_min, 0), P(x_max, 0));
  line(P(0, y_min), P(0, y_max));

  red;
  plot(f, 0, 2, 200);

  plain;
  blue;
  tan_field(f, 0, 2, 81);

  black;
  v_axis_masklabels(P(0, y_min), P(0, y_max), 2*y_size, P(0,0), c);

  end();
}
\end{verbatim}
\end{small}
\caption{The tangent field of a Lissajous curve.}
\label{fig:lissajous}
\end{center}
\end{figure}

\clearpage

\ePiX\ plots fractal curves generated by ``seeds''. A seed is a
(usually short) piecewise-linear curve whose equal-length segments
point in one of~$n$ equally-spaced directions
(``spokes''). Recursively up to a pre-specified depth~$D$, each
segment is replaced by a scaled, rotated, translated copy of the seed
(shown in red). The sequence \{\texttt{6,4,0,1,-1,0}\} generates the
Koch snowflake curve: There are 6~spokes, the seed has length~4, and
the remaining entries are the seed proper.

\begin{figure}[hbt]
\begin{center}
\input{koch.eepic}
\caption{Approximations of the Koch snowflake curve.}
\label{fig:koch}
\end{center}
\end{figure}

%%\begin{scriptsize}
\begin{verbatim}
#include <epix.h>

const int koch_seed[] = {6, 4, 0, 1, -1, 0};

main() 
{
  bounding_box(P(-1,-1), P(1, 1));
  picture(P(75, 75));
  unitlength("1pt");

  for (int i=1; i <= 5; ++i)
    {
      begin();   
      
      fractal(cis(5*M_PI/6), cis(  M_PI/6), i, koch_seed);
      fractal(cis(  M_PI/6), cis( -M_PI_2), i, koch_seed);
      fractal(cis( -M_PI_2), cis(5*M_PI/6), i, koch_seed);

      if (i==1)
        {
          bold;
          red;
          fractal(cis(5*M_PI/6), cis(M_PI/6), i, koch_seed);
          black;
          end_bold;
        }
        end();
    }
}
\end{verbatim}
%%\end{scriptsize}

\clearpage
One way to visualize the planar flow of a vector field is to put a
``test object'' into the flow and follow it as time
passes. Figure~\ref{fig:flow-plot} follows three circles under a flow
that simulates water spiralling down a drain. The curves progress from
blue to red as time passes.

\begin{figure}[hbt]
\begin{center}
\input{flow-plot.eepic}
\caption{Test curves in a draining sink.}
\label{fig:flow-plot}
\end{center}
\end{figure}

\clearpage

\begin{verbatim}
#include "epix.h"
pair F(double s, double t)
{
  pair point = P(s,t);
  return (1/(0.01 + raw_len(point)))*J(point);
}
 
pair curve1(double t)
{
  return P(0.5,0)+polar(0.05, 2*M_PI*t);
}

pair curve2(double t)
{
  return P(0.75,0)+polar(0.05, 2*M_PI*t);
}

pair curve3(double t)
{
  return P(1,0)+polar(0.05, 2*M_PI*t);
}
 
triple color(double t)
{
  return P(t, t*(1-t), 1-t);
}

main() {
  unitlength("0.01in");
  bounding_box(P(-1,-1), P(1,1));
  picture(P(550,550));
 
  begin();
 
  bold;
  flow_plot(F, curve1, color, 3, 16, 80);
  flow_plot(F, curve2, color, 3, 16, 240);                                  
  flow_plot(F, curve3, color, 3, 16, 240);                                  

  end();
}
\end{verbatim}

\end{document}
