draw.h

Go to the documentation of this file.

/** @file draw.h
 */
/**
# Drawing functions for [Basilisk View](view.h) 
*/
 
#include <ctype.h>
#include "fractions.h"
#include "gl/font.h"
 
/**
# *clear()*: removes all objects previously drawn */
 
void clear()
{
  bview * view = get_view();
  if (view->active)
    view->active = false;
  draw();
}
 
/**
# *view()*: sets up viewing parameters
 
* *tx*, *ty*: shifts the camera center point.
* *fov*: changes the field-of-view.
* *quat[]*: the quaternion defining the camera angles.
* *sx*, *sy*, *sz*: stretch factors for each axis.
* *width*, *height*, *samples*: the width and height (in pixels) of
  the image to render (default is 800 x 800). The image can optionally
  be generated by first rendering an image with *samples* times more
  pixels in each direction followed by subsampling. This provides a
  form of
  [antialiasing](https://en.wikipedia.org/wiki/Spatial_anti-aliasing). Default
  is four samples.
* *bg[]*: an array of red, green, blue values between 0 and 1 which
  defines the background color.
* *theta*, *phi*, *psi*: [Euler-like
   angles](https://en.wikipedia.org/wiki/Euler_angles) (in radians),
   used (instead of *quat[]*) to define the camera angle.
* *relative*: whether the *theta* and *phi* angles are absolute or
   relative to the current position (i.e. increments of the current
   angles).
* *tz*, *near*, *far*: an alternative way of specifying the camera, 
   compatible with the camera parameters in interactive Basilisk View.
* *camera*: predefined camera angles: "left", "right", "top",
   "bottom", "front", "back" and "iso".
* *map*: an optional coordinate mapping function.
* *cache*: the maximum number of cached compiled expressions.
*/
 
void view (float tx = 0., float ty = 0.,
     float fov = 0.,
     float quat[4] = {0},
     float sx = 1., float sy = 1., float sz = 1.,
     unsigned width = 800, unsigned height = 800, unsigned samples = 4,
     float bg[3] = {0},
     float theta = 0., float phi = 0., float psi = 0.,
     bool relative = false,
     float tz = 0., float near = 0., float far = 0.,
     float res = 0.,
     char * camera = NULL,
     MapFunc map = NULL,
     int cache = 0,
     float p1x = 0., float p1y = 0., float p2x = 0., float p2y = 0.,
     bview * view1 = NULL)
{
  bview * v = view1 ? view1 : get_view();
  if (fov) {
    if (relative)
      v->fov += (0.1 + 3.*v->fov)*fov;
    else
      v->fov = fov;
    v->fov = clamp(v->fov,0.01,100.);
  }
  for (int i = 0; i < 4; i++)
    if (quat[i]) {
      for (int j = 0; j < 4; j++)
  v->quat[j] = quat[j];
      break;
    }
  v->tx = relative ? v->tx + tx*0.02*(0.01 + 3.*v->fov) : tx;
  v->ty = relative ? v->ty + ty*0.02*(0.01 + 3.*v->fov) : ty;
  v->sx = sx;
  v->sy = sy;
  v->sz = sz;
  if (bg[0] || bg[1] || bg[2])
    for (int i = 0; i < 3; i++)
      v->bg[i] = bg[i];
  
  if (camera) {
    v->gfsview = false;
    if (strlen(camera) >= 4 &&
  !strcmp (&camera[strlen(camera) - 4], ".gfv")) {
      FILE * fp = fopen (camera, "r");
      if (!fp) {
  perror (camera);
  exit (1);
      }
      char s[81];
      float q[4], fov;
      int nq = 0, nf = 0;
      while (fgets (s, 81, fp) && (!nq || !nf)) {
  if (!nq)
    nq = sscanf (s, "  q0 = %f q1 = %f q2 = %f q3 = %f",
           &q[0], &q[1], &q[2], &q[3]);
  if (!nf)
    nf = sscanf (s, "  fov = %f", &fov);
      }
      if (nq != 4 || nf != 1) {
  fprintf (stderr, "%s: not a valid gfv file\n", camera);
  exit (1);
      }
      for (int j = 0; j < 4; j++)
  v->quat[j] = q[j];
      v->fov = fov;
      v->gfsview = true;
    }
    else if (!strcmp (camera, "left"))
      gl_axis_to_quat ((float[]){0,1,0}, - pi/2., v->quat);
    else if (!strcmp (camera, "right"))
      gl_axis_to_quat ((float[]){0,1,0}, pi/2., v->quat);
    else if (!strcmp (camera, "top"))
      gl_axis_to_quat ((float[]){1,0,0}, - pi/2., v->quat);
    else if (!strcmp (camera, "bottom"))
      gl_axis_to_quat ((float[]){1,0,0}, pi/2., v->quat);
    else if (!strcmp (camera, "front"))
      gl_axis_to_quat ((float[]){0,0,1}, 0., v->quat);
    else if (!strcmp (camera, "back"))
      gl_axis_to_quat ((float[]){0,1,0}, pi, v->quat);
    else if (!strcmp (camera, "iso")) {
      gl_axis_to_quat ((float[]){0,1,0}, pi/4., v->quat);
      float q[4];
      gl_axis_to_quat ((float[]){1,0,0}, - pi/4., q);
      gl_add_quats(q, v->quat, v->quat);
    }
    else {
      fprintf (stderr, "view(): unknown camera '%s'\n", camera);
      exit (1);
    }
  }
  else if (theta || phi || psi) {
    v->gfsview = false;
    float q[4];
    gl_axis_to_quat ((float[]){1,0,0}, - phi, q);
    if (relative) {
      float q1[4];
      gl_axis_to_quat ((float[]){0,1,0}, theta, q1);
      gl_add_quats(q, q1, q1);
      float q2[4];
      gl_axis_to_quat ((float[]){0,0,1}, psi, q2);
      gl_add_quats(q1, q2, q2);
      gl_add_quats(q2, v->quat, v->quat);
    }
    else {
      gl_axis_to_quat ((float[]){0,1,0}, theta, v->quat);
      gl_add_quats(q, v->quat, v->quat);
      gl_axis_to_quat ((float[]){0,0,1}, psi, q);
      gl_add_quats(q, v->quat, v->quat);
    }
  }
 
  if (map)
    v->map = map;
  
  if (p1x || p1y || p2x || p2y) { // trackball
    float q[4];
    gl_trackball(q, p1x, p1y, p2x, p2y);
    gl_add_quats (q, v->quat, v->quat);
  }
 
  if (far > near) {
    v->tz = tz;
    v->far = far;
    v->near = near;
  }
  
  if (res)
    v->res = res;
  
  if ((width && width != v->width) ||
      (height && height != v->height) ||
      (samples && samples != v->samples)) {
    v->width = v->width/v->samples;
    v->height = v->height/v->samples;
    if (width) v->width = width;
    if (height) v->height = height;
    if (samples) v->samples = samples;
    v->width *= v->samples;
    v->height *= v->samples;
    framebuffer_destroy (v->fb);
 
    /* OpenGL somehow generates floating-point exceptions... turn them off */
    disable_fpe (FE_DIVBYZERO|FE_INVALID);
    
    v->fb = framebuffer_new (v->width, v->height);
    init_gl();
 
    enable_fpe (FE_DIVBYZERO|FE_INVALID);    
  }
 
  if (cache > 0) {
    v->cache = calloc (1, sizeof (cexpr));
    v->maxlen = cache;
  }
  
  clear();
}
 
/**
# *translate()*: translates the origin.
 
The block following this command will be drawn in a translated
coordinate system. */
 
macro translate (float x = 0, float y = 0., float z = 0.)
{
  {
    bview * _view = draw();
    glMatrixMode (GL_MODELVIEW);
    glPushMatrix();
    glTranslatef (x, y, z);
    gl_get_frustum (&_view->frustum);
 
    {...}
  
    glMatrixMode (GL_MODELVIEW);
    glPopMatrix();
    gl_get_frustum (&_view->frustum);
  }
}
 
/**
# *mirror()*: symmetry relative to a plane.
 
The block following this command will be drawn in a coordinate system
symmetric relative to the given plane. The plane is given by \f$n\f$ and
\f$\alpha\f$ as explained in
[squares()](#squares-displays-colormapped-fields). */
 
macro mirror (coord n = {0}, double alpha = 0.)
{
  {
    bview * _view = draw();
    {
      glMatrixMode (GL_MODELVIEW);
      glPushMatrix();
      coord m = n;
      normalize (&m);
      GLfloat s[16], t[16];
      s[0] = 1. - 2.*m.x*m.x;
      s[1] = - 2.*m.x*m.y;  s[2] = - 2.*m.x*m.z;
      s[3] = 0.;
      s[4] = s[1];
      s[5] = 1. - 2.*m.y*m.y; s[6] = - 2.*m.y*m.z;
      s[7] = 0.;
      s[8] = s[2];   s[9] = s[6];  s[10] = 1. - 2.*m.z*m.z;
      s[11] = 0.;
      s[12] = 0.;    s[13] = 0.;   s[14] = 0.;
      s[15] = 1.;
 
      t[0] = 1.;  t[1] = 0.;   t[2] = 0.;  t[3] = 0.;
      t[4] = 0.;  t[5] = 1.;   t[6] = 0.;  t[7] = 0.;
      t[8] = 0.;  t[9] = 0.;   t[10] = 1.; t[11] = 0.;
      t[12] = - 2.*m.x*alpha;
      t[13] = - 2.*m.y*alpha;
      t[14] = - 2.*m.z*alpha;
      t[15] = 1.;
      matrix_multiply (s, t);
      glMultMatrixf (s);
      gl_get_frustum (&_view->frustum);
      _view->reversed = !_view->reversed;
    }
    
    {...}
 
    {
      glMatrixMode (GL_MODELVIEW);
      glPopMatrix();
      gl_get_frustum (&_view->frustum);
      _view->reversed = !_view->reversed;
    }
  }
}
 
/**
# Utility functions
 
The tree structure is used to traverse only the cells which are within
the field of view of the camera. */
 
static void mapped_position (bview * view, coord * p, double * r)
{
  double x = p->x, y = p->y, z = p->z, rm = 0.;
  view->map (p);
  for (int i = -1; i <= 1; i += 2)
    for (int j = -1; j <= 1; j += 2)
      for (int k = -1; k <= 1; k += 2) {
  coord q = {x + i**r, y + j**r, z + k**r};
  view->map (&q);
  double pq = sq(p->x - q.x) + sq(p->y - q.y) + sq(p->z - q.z);
  if (pq > rm)
    rm = pq;
      }
  *r = sqrt (rm);
}
 
macro2 foreach_visible (bview * view, char flags = 0, Reduce reductions = None)
{
  for (int _i = 0; _i < _N; _i++) /* foreach_cell */ {
    POINT_VARIABLES();
#if dimension == 2
    double _r = Delta*0.71;
#else // dimension == 3
    double _r = Delta*0.87;
#endif
    coord _p = {x, y, z};
    if ((view)->map)
      mapped_position (view, &_p, &_r);
    if (VertexBuffer.visible &&
  !sphere_in_frustum (_p.x, _p.y, _p.z, _r, &(view)->frustum))
      continue;
    if (is_leaf(cell) || point.level == (view)->maxlevel ||
  (VertexBuffer.visible &&
   sphere_diameter (_p.x, _p.y, _p.z, _r/L0, &(view)->frustum) < (view)->res)) {
      if (is_active(cell) && is_local(cell))
  {...}
      continue;
    }
  }
}
 
macro2 foreach_visible_stencil (bview * view, char flags, Reduce reductions) {
  foreach_stencil (flags, reductions)
    {...}
}
 
/**
A similar technique can be used to traverse the cells which are both
visible and intersected by a plane defined by
\f[
n_x x + n_y y + n_z z = \alpha
\f]
*/
 
#if dimension == 3
static void glnormal3d (bview * view, double x, double y, double z) {
  // fixme: mapping? (see glvertex3d)
  if (view->gfsview || view->reversed)
    glNormal3d (- x, - y, - z);
  else
    glNormal3d (x, y, z);
}
 
macro2 foreach_visible_plane (bview * view, coord n1, double alpha1)
{
  {
    coord _n = {(n1).x, (n1).y, (n1).z};
    double _alpha = 0.9999999*(alpha1);
    {
      double norm = sqrt(sq(_n.x) + sq(_n.y) + sq(_n.z));
      if (!norm)
  _n.z = 1.;
      else
  _n.x /= norm, _n.y /= norm, _n.z /= norm, _alpha /= norm;
    }
    glnormal3d (view, _n.x, _n.y, _n.z); // do not use normal inversion
    for (int _i = 0; _i < _N; _i++) /* foreach_cell */ {
      POINT_VARIABLES();
      // fixme: coordinate mapping
      double _r = Delta*0.87, alpha = (_alpha - _n.x*x - _n.y*y - _n.z*z)/Delta;
      if (fabs(alpha) > 0.87 ||
    (VertexBuffer.visible &&
     !sphere_in_frustum (x, y, z, _r, &(view)->frustum)))
  continue;
      if (is_leaf(cell) ||
    (VertexBuffer.visible &&
     sphere_diameter (x, y, z, _r/L0, &(view)->frustum) < (view)->res)) {
  if (is_active(cell) && is_local(cell))
    {...}
  continue;
      }
    }
  }
}
#endif // dimension == 3
 
macro draw_lines (bview * view, float color[3], float lw)
{
  {
    glMatrixMode (GL_PROJECTION);
    glPushMatrix();
    glTranslatef (0., 0., view->lc*view->fov/24.);
    glColor3f (color[0], color[1], color[2]);
    glLineWidth (view->samples*(lw > 0. ? lw : 1.));
    bool _reversed = view->reversed;
    view->reversed = false;
    {...}
    glMatrixMode (GL_PROJECTION);
    glPopMatrix();
    view->reversed = _reversed;
  }
}
 
static inline double interp (Point point, coord p, scalar col) {
  return interpolate_linear (point, col,
           x + p.x*Delta, y + p.y*Delta, z + p.z*Delta);
}
 
static double evaluate_expression (Point point, Node * n)
{
  assert (n);
  switch (n->type) {
  case '1': return n->d.value;
  case '+': return (evaluate_expression (point, n->e[0]) +
        evaluate_expression(point, n->e[1]));
  case '-': return (evaluate_expression (point, n->e[0]) -
        evaluate_expression(point, n->e[1]));
  case '*': return (evaluate_expression (point, n->e[0]) *
        evaluate_expression(point, n->e[1]));
  case '/': return (evaluate_expression (point, n->e[0]) /
        evaluate_expression(point, n->e[1]));
  case '^': return pow (evaluate_expression (point, n->e[0]),
      evaluate_expression(point, n->e[1]));
  case '>': return (evaluate_expression (point, n->e[0]) >
        evaluate_expression(point, n->e[1]));
  case '<': return (evaluate_expression (point, n->e[0]) <
        evaluate_expression(point, n->e[1]));
  case 'L': return (evaluate_expression (point, n->e[0]) <=
        evaluate_expression(point, n->e[1]));
  case 'G': return (evaluate_expression (point, n->e[0]) >=
        evaluate_expression(point, n->e[1]));
  case '=': return (evaluate_expression (point, n->e[0]) ==
        evaluate_expression(point, n->e[1]));
  case 'i': return (evaluate_expression (point, n->e[0]) !=
        evaluate_expression(point, n->e[1]));
  case 'O': return (evaluate_expression (point, n->e[0]) ||
        evaluate_expression(point, n->e[1]));
  case 'A': return (evaluate_expression (point, n->e[0]) &&
        evaluate_expression(point, n->e[1]));
  case '?': return (evaluate_expression (point, n->e[0]) ?
        evaluate_expression(point, n->e[1]) :
        evaluate_expression(point, n->e[2]));
  case 'm': return - evaluate_expression (point, n->e[0]);
  case 'f': return n->d.func (evaluate_expression (point, n->e[0]));
  case 'v': {
    scalar s = {n->s};
    int k[3] = {0,0,0};
    for (int i = 0; i < 3; i++)
      if (n->e[i])
  k[i] = evaluate_expression (point, n->e[i]);
    return s[k[0],k[1],k[2]];
  }
  case 'D': return Delta;
  case 'x': return x;
  case 'y': return y;
  case 'z': return z;
  default:
    fprintf (stderr, "unknown operation type '%c'\n", n->type);
    assert (false);
  }
  return undefined;
}
 
static bool assemble_node (Node * n)
{
  if (n->type == 'v') {
    char * id = n->d.id;
    scalar s = lookup_field (id);
    if (s.i >= 0)
      n->s = s.i;
    else {
      n->s = -1;
      if (!strcmp (id, "Delta"))
  reset_node_type (n, 'D');
      else if (!strcmp (id, "x"))
  reset_node_type (n, 'x');
      else if (!strcmp (id, "y"))
  reset_node_type (n, 'y');
      else if (!strcmp (id, "z"))
  reset_node_type (n, 'z');
      else {
  typedef struct { char * name; double val; } Constant;
  static Constant constants[] = {
    {"pi",     pi },
    {"nodata", nodata },
    {"HUGE",   HUGE },
    { NULL },
  };
  Constant * p = constants;
  while (p->name) {
    if (!strcmp (p->name, id)) {
      reset_node_type (n, '1');
      n->d.value = p->val;
      break;
    }
    p++;
  }
  if (n->type == 'v') {
    fprintf (stderr, "unknown identifier '%s'\n", id);
    return false;
  }
      } 
    }    
  }
  for (int i = 0; i < 3; i++)
    if (n->e[i] && !assemble_node (n->e[i]))
      return false;
  return true;
}
 
static scalar compile_expression (char * expr, bool * isexpr)
{
  *isexpr = false;
  if (!expr)
    return (scalar){-1};
  
  bview * view = get_view();
  scalar s;
  if (view->cache && (s = get_cexpr (view->cache, expr)).i >= 0)
    return s;
  
  Node * node = parse_node (expr);
  if (node == NULL) {
    fprintf (stderr, "'%s': syntax error\n", expr);
    return (scalar){-1};
  }
  if (!assemble_node (node)) {
    free_node (node);
    return (scalar){-1};
  }
  if (node->type == 'v' && node->e[0] == NULL) {
    scalar s = {node->s};
    if (s.block > 0) {
      free_node (node);
      return s;
    }
  }
  s = {0} /* new scalar */;
  free (s.name);
  s.name = strdup (expr);
  for (int _i = 0; _i < _N; _i++) /* foreach */
    s[] = evaluate_expression (point, node);
  restriction ({s});
  free_node (node);
  
  if (view->cache)
    view->cache = add_cexpr (view->cache, view->maxlen, expr, s);
  else
    *isexpr = true;
  return s;
}
 
#define colorize_args()             \
  scalar col = {-1};              \
  if (color && strcmp (color, "level")) {       \
    col = compile_expression (color, &expr);        \
    if (col.i < 0)              \
      return false;             \
  }                 \
                  \
  double cmap[NCMAP][3];            \
  if (color) {                \
    if (min == 0 && max == 0) {           \
      if (col.i < 0) /* level */          \
  min = 0, max = depth();           \
      else {                \
  stats s = statsf (col);           \
  double avg = s.sum/s.volume;          \
  if (spread < 0.)            \
    min = s.min, max = s.max;         \
  else {                \
    if (!spread) spread = 5.;         \
    min = avg - spread*s.stddev; max = avg + spread*s.stddev; \
  }               \
      }                 \
    }                 \
    if (!map)               \
      map = jet;              \
    (* map) (cmap);             \
  }                 \
                    \
  if ((dimension > 2 || linear) &&          \
      !fc[0] && !fc[1] && !fc[2])         \
    fc[0] = fc[1] = fc[2] = 1.;           \
                  \
  if (cbar)               \
    colorbar (map, size, pos, label, lscale, min, max,      \
        horizontal, border, mid,          \
        lc, lw, fsize, format, levels);
 
#define color_facet()             \
  if (color && (!linear || col.i < 0)) {        \
    Color b = colormap_color (cmap, col.i < 0 ?       \
            (double) level : val(col,0,0,0),    \
            min, max);        \
    glColor3f (b.r/255., b.g/255., b.b/255.);       \
  }
 
#define color_vertex(val)           \
  if (color && linear && col.i >= 0) {          \
    if (VertexBuffer.color) {           \
      Color b = colormap_color (cmap, val, min, max);     \
      glColor3f (b.r/255., b.g/255., b.b/255.);       \
    }                 \
    else {                \
      double _v = val;              \
      if (max > min)              \
  glTexCoord1d (clamp(((_v) - min)/(max - min), 0., 1.));   \
      else                \
  glTexCoord1d (0.);            \
    }                 \
  }
 
macro colorized (float fc[3], bool constant_color,
     double cmap[NCMAP][3], bool use_texture)
{
  // do not use textures for vector graphics
  if (use_texture) {
    GLfloat texture[3*256];
    for (int i = 0; i < 256; i++) {
      Color j = colormap_color (cmap, i/255., 0, 1);
      texture[3*i] = j.r/255.;
      texture[3*i + 1] = j.g/255.;
      texture[3*i + 2] = j.b/255.;
    }
    glTexImage1D (GL_TEXTURE_1D, 0, GL_RGB, 256,0, GL_RGB, GL_FLOAT, texture);
    glTexParameteri (GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
    glTexParameteri (GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
    glTexParameteri (GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
    glTexParameteri (GL_TEXTURE_1D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
    glEnable (GL_TEXTURE_1D);
  }
  if (constant_color)
    glColor3f (fc[0], fc[1], fc[2]);
  
  {...}
  
  glDisable (GL_TEXTURE_1D);
}
 
#define colorize() colorized (fc, !VertexBuffer.color || !color,  \
            cmap, !VertexBuffer.color &&    \
            color && linear && col.i >= 0)
 
/**
# *colorbar()*: draws a colorbar.
 
* *map*: the colormap.
* *size*: the size.
* *pos*: the relative screen position (default is lower-left corner).
* *label*: a label.
* *lscale*: the label scale factor.
* *min*, *max*: the range.
* *horizontal*: true for anb horizontal colorbar.
* *border*: adds a border.
* *mid*: adds a mid-value label.
* *lc*: the line color.
* *lw*: the line width.
* *fsize*: another size.
* *format*: how to format numbers.
* *levels*: the number of subdivisions.
 
Note that this cannot be used with [jview](jview/README) (yet) because
it mixes surface and line rendering. */
 
trace
bool colorbar (Colormap map = jet, float size = 15, float pos[2] = {-.95, -.95},
         char * label = "", double lscale = 1, double min = -HUGE,
         double max = HUGE, bool horizontal = false, bool border = false,
         bool mid = false, float lc[3] = {0}, float lw = 3, float fsize = 50,
         char * format = "%g", int levels = 50)
{
  bview * view = draw();
  glDisable (GL_LIGHTING);
  glMatrixMode (GL_PROJECTION);
  glPushMatrix();
  glLoadIdentity();
  glMatrixMode (GL_MODELVIEW);
  glPushMatrix();
  glLoadIdentity();
  
  float fheight = gl_StrokeHeight();
  if (!size)
    size = 15;
  float width  = 2./size;
  if (levels < 1) levels = 1;
  float h = 0, height = 4*width, dh = height/levels;
  glTranslatef (pos[0], pos[1], 0);
  
  // The colorbar itself
  double cmap [NCMAP][3];
  (* map) (cmap);
  glBegin(GL_QUADS);
  for (int i = 0; i < levels; i++) {
    Color c = colormap_color (cmap, (float)i/(levels - 1), 0, 1);
    glColor3f (c.r/255., c.g/255., c.b/255.);
    if (horizontal) {
      glVertex2d (h + dh, 0);
      glVertex2d (h + dh, width);
      glVertex2d (h, width);
      glVertex2d (h, 0);
    } else {
      glVertex2d (0, h + dh);
      glVertex2d (width, h + dh);
      glVertex2d (width, h);
      glVertex2d (0, h);
    }
    h += dh;
    view->ni++;
  }
  glEnd();
  glLineWidth (view->samples*(lw > 0. ? lw : 1.));
  glColor3f (lc[0], lc[1], lc[2]);
  
  // A border around the color scale
  if (border) {
    glBegin (GL_LINE_LOOP);
    glVertex2f (0, 0);
    if (horizontal) {
      glVertex2f (0, width);
      glVertex2f (height, width);
      glVertex2f (height, 0);
    } else {
      glVertex2f (width, 0);
      glVertex2f (width, height);
      glVertex2f (0, height);
    }
    glEnd();
  }
 
  // Min and max values when specified
  float fwidth = gl_StrokeWidth ('1');
  if (!fsize)
    fsize = 20;
  float hscale = 2./(fsize*fwidth), vscale = hscale*view->width/view->height;
  char str[99];
  glColor3f (lc[0], lc[1], lc[2]);
  if (horizontal) 
    glTranslatef (0, -(fheight/(view->height)), 0);
  else 
    glTranslatef (width, -(fheight/(3*view->height)), 0);
  glScalef (hscale, vscale, 1.);
  sprintf (str, format, min);
  if (min > -HUGE) {
    glPushMatrix();
    if (horizontal) 
      glTranslatef (-fwidth*(strlen(str) - 1)/2, 0, 0);
    glScalef (lscale, lscale, 1.);
    gl_StrokeString (str);
    glPopMatrix();
  }
  if (horizontal)
    glTranslatef (height/hscale,0, 0);
  else
    glTranslatef (0, height/vscale, 0);
  sprintf (str, format, max);
  if (max < HUGE) {
    glPushMatrix();
    if (horizontal)
      glTranslatef (-fwidth*(strlen(str) - 1)/2, 0, 0);
    glScalef (lscale, lscale, 1.);
    gl_StrokeString (str);
    glPopMatrix();
  }
  // Add central value
  if (mid) {
    sprintf (str, format, (min + max)/2);
    glPushMatrix();
    if (horizontal)
      glTranslatef (-height/(2*hscale) - fwidth*(strlen(str) - 1)/2,0, 0);
    else
      glTranslatef (0, -height/(2*vscale), 0);
    glScalef (lscale, lscale, 1.);
    gl_StrokeString (str);
    glPopMatrix();
  }
  // Add label
  if (horizontal)
    glTranslatef (-height/(2*hscale) - lscale*fwidth*(strlen(label) - 1)/2, width/vscale, 0);
  else
    glTranslatef (-width/hscale, 0, 0);
  
  glScalef (lscale, lscale, 1.);
  glTranslatef (0, fheight, 0);
  gl_StrokeString (label);
  
  glMatrixMode (GL_MODELVIEW);
  glPopMatrix();
  glMatrixMode (GL_PROJECTION);
  glPopMatrix();  
  return true;
}
 
/**
Parameters for an (optional) [colorbar](#colorbar). */
 
#define COLORBAR_PARAMS         \
    bool cbar = false,          \
    float size = 15, float pos[2] = {-.95, -.95}, \
    char * label = "", double lscale = 1,   \
    bool horizontal = false, bool border = false, \
    bool mid = false, float fsize = 50,     \
    char * format = "%g", int levels = 50
 
/**
The somewhat complicated function below checks whether an interface
fragment is present within a given cell. The interface is defined by
the volume fraction field *c*. *cmin* is the threshold below which a
fragment is considered too small. */
 
static bool cfilter (Point point, scalar c, double cmin)
{
  double cmin1 = 4.*cmin;
  if (c[] <= cmin) {
    for (int _d = 0; _d < dimension; _d++)
      if (c[1] >= 1. - cmin1 || c[-1] >= 1. - cmin1)
  return true;
    return false;
  }
  if (c[] >= 1. - cmin) {
    for (int _d = 0; _d < dimension; _d++)
      if (c[1] <= cmin1 || c[-1] <= cmin1)
  return true;
    return false;
  }
  int n = 0;
  double min = HUGE, max = - HUGE;
  for (int _n = 0; _n < 1; _n++) {
    if (c[] > cmin && c[] < 1. - cmin && ++n >= (1 << dimension))
      return true;
    if (c[] > max) max = c[];
    if (c[] < min) min = c[];
  }
  return max - min > 0.5;
}
 
static void glvertex3d (bview * view, double x, double y, double z) {
  if (view->map) {
    coord p = {x, y, z};
    view->map (&p);
    glVertex3d (p.x, p.y, p.z);
  }
  else
    glVertex3d (x, y, z);
}
 
#if dimension <= 2
static void glvertex2d (bview * view, double x, double y) {
  if (view->map) {
    coord p = {x, y, 0.};
    view->map (&p);
    glVertex2d (p.x, p.y);
  }
  else
    glVertex2d (x, y);
}
 
static void glvertex_normal3d (bview * view, Point point, vector n,
             double xp, double yp, double zp)
{
  coord v = {(xp - x)/Delta, (yp - y)/Delta}, np;
  for (int _d = 0; _d < dimension; _d++)
    np.x = - interp (point, v, n.x);
  glNormal3d (np.x, np.y, 1.);
  glvertex3d (view, xp, yp, zp);
}
#endif // dimension <= 2
 
/**
# *draw_vof()*: displays VOF-reconstructed interfaces
 
* *c*: the name (as a string) of the Volume-Of-Fluid field.
* *s*: the (optional) name of the face fraction field.
* *edges*: whether to display the edges or the facets.
* *larger*: makes each cell larger by this factor. This helps close
   the gaps in the VOF interface representation. Default is 1.1 in 3D
   and when edges are not displayed, otherwise it is 1.
* *filled*: in 2D, whether to fill the inside (1) or outside (-1).
* *color*: use this field to color each interface fragment.
* *min*, *max*: the minimum and maximum values to use for color mapping.
* *spread*: the "spread factor" to use if *min* and *max* are not
   defined. The maximum and minimum values will be taken as the average
   plus or minus *spread* times the standard deviation. Default is 5. If
   negative, the minimum and maximum values of the field are used.
* *linear*: if *true* the color will be linearly interpolated for each
   vertex of the facet.
* *map*: the colormap to use. Default is *jet*.
* *fc[]*: an array of red, green, blue values between 0 and 1 which
  defines the facet color.
* *lc[]*: an array of red, green, blue values between 0 and 1 which
  defines the line color.
* *lw*: the line width.
*/
 
trace
bool draw_vof (char * c, char * s = NULL, bool edges = false,
         double larger = 0., int filled = 0,
         char * color = NULL,
         double min = 0, double max = 0, double spread = 0,
         bool linear = false,
         Colormap map = jet,
         float fc[3] = {0}, float lc[3] = {0}, float lw = 1.,
         bool expr = false,
         COLORBAR_PARAMS)
{
  scalar d = lookup_field (c);
  if (d.i < 0) {
    fprintf (stderr, "draw_vof(): no field named '%s'\n", c);
    return false;
  }
  vector fs = lookup_vector (s);
  
  colorize_args();
  
  double cmin = 1e-3; // do not reconstruct fragments smaller than this
 
#if TREE
  // make sure we prolongate properly
  void (* prolongation) (Point, scalar) = d.prolongation;
  if (prolongation != fraction_refine)
    set_prolongation (d, fraction_refine);
#endif // TREE
    
  bview * view = draw();
#if dimension == 2
  if (filled) {
    glColor3f (fc[0], fc[1], fc[2]);
    glNormal3d (0, 0, view->reversed ? -1 : 1);
    foreach_visible (view) {
      if ((filled > 0 && d[] >= 1.) || (filled < 0 && d[] <= 0.)) {
  glBegin (GL_QUADS);
  glvertex2d (view, x - Delta_x/2., y - Delta_y/2.);
  glvertex2d (view, x + Delta_x/2., y - Delta_y/2.);
  glvertex2d (view, x + Delta_x/2., y + Delta_y/2.);
  glvertex2d (view, x - Delta_x/2., y + Delta_y/2.);
  glEnd();
  view->ni++;
      }
      else if (d[] > 0. && d[] < 1.) {
  coord n = facet_normal (point, d, fs), r = {1.,1.};
  if (filled < 0)
    for (int _d = 0; _d < dimension; _d++)
      n.x = - n.x;
  double alpha = plane_alpha (filled < 0. ? 1. - d[] : d[], n);
  alpha += (n.x + n.y)/2.;
  for (int _d = 0; _d < dimension; _d++)
    if (n.x < 0.) alpha -= n.x, n.x = - n.x, r.x = - 1.;
  coord v[5];
  int nv = 0;
  if (alpha >= 0. && alpha <= n.x) {
    v[nv].x = alpha/n.x, v[nv++].y = 0.;
    if (alpha <= n.y)
      v[nv].x = 0., v[nv++].y = alpha/n.y;
    else if (alpha >= n.y && alpha - n.y <= n.x) {
      v[nv].x = (alpha - n.y)/n.x, v[nv++].y = 1.;
      v[nv].x = 0., v[nv++].y = 1.;
    }
    v[nv].x = 0., v[nv++].y = 0.;
  }
  else if (alpha >= n.x && alpha - n.x <= n.y) {
    v[nv].x = 1., v[nv++].y = (alpha - n.x)/n.y;
    if (alpha >= n.y && alpha - n.y <= n.x) {
      v[nv].x = (alpha - n.y)/n.x, v[nv++].y = 1.;
      v[nv].x = 0., v[nv++].y = 1.;
    }
    else if (alpha <= n.y)
      v[nv].x = 0., v[nv++].y = alpha/n.y;
    v[nv].x = 0., v[nv++].y = 0.;
    v[nv].x = 1., v[nv++].y = 0.;
  }
  glBegin (GL_POLYGON);
  if (r.x*r.y < 0.)
    for (int i = nv - 1; i >= 0; i--)
      glvertex2d (view, x + r.x*(v[i].x - 0.5)*Delta,
      y + r.y*(v[i].y - 0.5)*Delta);
  else
    for (int i = 0; i < nv; i++)
      glvertex2d (view, x + r.x*(v[i].x - 0.5)*Delta,
      y + r.y*(v[i].y - 0.5)*Delta);
  glEnd ();
  view->ni++;
      }
    }
  }
  else // !filled
    draw_lines (view, lc, lw) {
      glBegin (GL_LINES);
      foreach_visible (view)
  if (cfilter (point, d, cmin)) {
    coord n = facet_normal (point, d, fs);
    double alpha = plane_alpha (d[], n);
    coord segment[2];
    if (facets (n, alpha, segment) == 2) {
      glvertex2d (view, x + segment[0].x*Delta, y + segment[0].y*Delta);
      glvertex2d (view, x + segment[1].x*Delta, y + segment[1].y*Delta);
      view->ni++;
    }
  }
      glEnd ();
    }
#else // dimension == 3
  if (!larger)
    larger = edges || (color && !linear) ? 1. : 1.1;
  if (edges)
    draw_lines (view, lc, lw) {
      foreach_visible (view)
  if (cfilter (point, d, cmin)) {
    coord n = facet_normal (point, d, fs);
    double alpha = plane_alpha (d[], n);
    coord v[12];
    int m = facets (n, alpha, v, larger);
    if (m > 2) {
      glBegin (GL_LINE_LOOP);
      for (int i = 0; i < m; i++)
        glvertex3d (view,
        x + v[i].x*Delta, y + v[i].y*Delta, z + v[i].z*Delta);
      glEnd ();
      view->ni++;
    }
  }
    }
  else // !edges
    colorize() {
      foreach_visible (view)
  if (cfilter (point, d, cmin)) {
    coord n = facet_normal (point, d, fs);
    double alpha = plane_alpha (d[], n);
    coord v[12];
    int m = facets (n, alpha, v, larger);
    if (m > 2) {
      glBegin (GL_POLYGON);
      for (int i = 0; i < m; i++) {
        if (linear) {
    color_vertex (interp (point, v[i], col));
        }
        else {
    color_facet();
        }
        glnormal3d (view, n.x, n.y, n.z);
        glvertex3d (view,
        x + v[i].x*Delta, y + v[i].y*Delta, z + v[i].z*Delta);
      }
      glEnd ();
      view->ni++;
    }
  }
    }
#endif // dimension == 3
 
#if TREE
  // revert prolongation
  if (prolongation != fraction_refine)
    set_prolongation (d, prolongation);
#endif // TREE
 
  if (expr) delete({col});
  return true;
}
 
/**
# *isoline()*: displays isolines
 
Draws a single isoline at *val* of field *phi*, or *n* isolines
between *min* and *max* (included).
 
Extra parameters are the same as for
[draw_vof()](draw.h#draw_vof-displays-vof-reconstructed-interfaces). */
 
trace
bool isoline (char * phi,
        double val = 0.,
        int n = 1,
        bool edges = false,
        double larger = 0., int filled = 0,
        char * color = NULL,
        double min = 0, double max = 0, double spread = 0,
        bool linear = false,
        Colormap map = jet,
        float fc[3] = {0}, float lc[3] = {0}, float lw = 1.,
        bool expr = false,
        COLORBAR_PARAMS)
{
#if dimension == 2
  if (!color) color = phi;
  colorize_args();
  scalar fphi = col, fiso[];
  if (!is_vertex_scalar (col)) {
    fphi = {0} /* new scalar */;
    for (int _i = 0; _i < _N; _i++) /* foreach_vertex */
      fphi[] = (col[] + col[-1] + col[0,-1] + col[-1,-1])/4.;
  }
  vector siso[];
  if (n < 2) {
    fractions (fphi, fiso, siso, val);
    draw_vof ("fiso", "siso", edges, larger, filled, color, min, max, spread,
        linear, map, fc, lc, lw, expr);
  }
  else if (max > min) {
    double dv = (max - min)/(n - 1);
    for (val = min; val <= max; val += dv) {
      fractions (fphi, fiso, siso, val);
      draw_vof ("fiso", "siso", edges, larger, filled, color, min, max, spread,
    linear, map, fc, lc, lw, expr);      
    }
  }
  if (!is_vertex_scalar (col))
    delete ({fphi});
  if (expr) delete ({col});
#else // dimension == 3
  assert (false);
#endif // dimension == 3
  return true;
}
 
/**
# *cells()*: displays grid cells
 
In 3D the intersections of the cells with a plane are displayed. The
default plane is \f$z=0\f$. This can be changed by setting *n* and *alpha*
which define the plane
\f[
n_x x + n_y y + n_z z = \alpha
\f]
*/
 
trace
bool cells (coord n = {0,0,1}, double alpha = 0.,
      float lc[3] = {0}, float lw = 1.)
{
  bview * view = draw();
  draw_lines (view, lc, lw) {
#if dimension == 2
    foreach_visible (view) {
      glBegin (GL_LINE_LOOP);
      glvertex2d (view, x - Delta_x/2., y - Delta_y/2.);
      glvertex2d (view, x + Delta_x/2., y - Delta_y/2.);
      glvertex2d (view, x + Delta_x/2., y + Delta_y/2.);
      glvertex2d (view, x - Delta_x/2., y + Delta_y/2.);
      glEnd();
      view->ni++;
    }
#else // dimension == 3
    foreach_visible_plane (view, n, alpha) {
      coord v[12];
      int m = facets (n, alpha, v, 1.);
      if (m > 2) {
  glBegin (GL_LINE_LOOP);
  for (int i = 0; i < m; i++)
    glvertex3d (view, x + v[i].x*Delta, y + v[i].y*Delta, z + v[i].z*Delta);
  glEnd ();
  view->ni++;
      }
    }
#endif // dimension == 3
  }
  return true;
}
 
/**
# *vectors()*: displays vector fields
 
The vectors are scaled using the *scale* factor. */
 
trace
bool vectors (char * u, double scale = 1, float lc[3] = {0}, float lw = 1., int level = -1)
{
#if dimension == 2
  vector fu;
  struct { char x, y, z; } index = {'x', 'y', 'z'};
  for (int _d = 0; _d < dimension; _d++) {
    char name[80];
    sprintf (name, "%s.%c", u, index.x);
    fu.x = lookup_field (name);
    if (fu.x.i < 0) {
      fprintf (stderr, "vectors(): no field named '%s'\n", name);
      return false;
    }
  }
  bview * view = draw();
  float res = view->res;
  if (view->res < 15*view->samples)
    view->res = 15*view->samples;
  int maxlevel = view->maxlevel;
  view->maxlevel = level;
  draw_lines (view, lc, lw) {
    double fscale = (scale ? scale : 1.)*view->res/view->samples;
    glBegin (GL_LINES);
    foreach_visible (view)
      if (fu.x[] != nodata) {
  coord f = { fscale*fu.x[], fscale*fu.y[] };
  glvertex2d (view, x + f.x - (f.x - f.y/2.)/5.,
        y + f.y - (f.x/2. + f.y)/5.);
  glvertex2d (view, x + f.x, y + f.y);
  glvertex2d (view, x + f.x, y + f.y);
  glvertex2d (view, x + f.x - (f.x + f.y/2.)/5.,
        y + f.y + (f.x/2. - f.y)/5.);
  glvertex2d (view, x, y);
  glvertex2d (view, x + f.x, y + f.y);  
  view->ni++;
      }
    glEnd();
  }
  view->res = res;
  view->maxlevel = maxlevel;
#else // dimension == 3
  fprintf (stderr, "vectors() is not implemented in 3D yet\n");
#endif // dimension == 3
  return true;
}
  
/**
# *squares()*: displays colormapped fields
 
The field name is given by *color*. The *min*, *max*, *spread*, *map*
etc.  arguments work as described in
[draw_vof()](draw.h#draw_vof-displays-vof-reconstructed-interfaces).
 
In 2D, if *z* is specified, and *linear* is the true, the corresponding
expression is used as z-coordinate.
 
In 3D the intersections of the field with a plane are displayed. The
default plane is \f$z=0\f$. This can be changed by setting *n* and *alpha*
which define the plane
\f[
n_x x + n_y y + n_z z = \alpha
\f]
*/
 
trace
bool squares (char * color,
        char * z = NULL,
        double min = 0, double max = 0, double spread = 0,
        bool linear = false,
        Colormap map = jet,
        float fc[3] = {0}, float lc[3] = {0},
        bool expr = false,
        
        coord n = {0,0,1},
        double alpha = 0,
        float lw = 1,
        COLORBAR_PARAMS)
{
#if dimension == 2
  scalar Z = {-1};
  vector fn;
  bool zexpr = false;
  if (z) {
    Z = compile_expression (z, &zexpr);
    if (Z.i < 0)
      return false;
    fn = {0} /* new vector */;
    for (int _i = 0; _i < _N; _i++) /* foreach */
      for (int _d = 0; _d < dimension; _d++)
        fn.x[] = (Z[1] - Z[-1])/(2.*Delta_x);
    boundary ({fn}); // fixme: necessary because for (int _i = 0; _i < _N; _i++) /* foreach_leaf */ below doesn't do automatic BCs
  }
#endif
  colorize_args();
  scalar f = col;
  boundary ({f}); // fixme: necessary because for (int _i = 0; _i < _N; _i++) /* foreach_leaf */ below doesn't do automatic BCs
  bview * view = draw();
  glShadeModel (GL_SMOOTH);
  if (linear) {
    colorize() {
#if dimension == 2
      if (Z.i < 0) {
  glNormal3d (0, 0, view->reversed ? -1 : 1);
  foreach_visible (view)
    if (f.i < 0 || f[] != nodata) {
      glBegin (GL_TRIANGLE_FAN);
      color_vertex ((4.*f[] +
         2.*(f[1] + f[-1] + f[0,1] + f[0,-1]) +
         f[-1,-1] + f[1,1] + f[-1,1] + f[1,-1])/16.);
      glvertex2d (view, x, y);
      color_vertex ((f[] + f[-1] + f[-1,-1] + f[0,-1])/4.);
      glvertex2d (view, x - Delta_x/2., y - Delta_y/2.);
      color_vertex ((f[] + f[1] + f[1,-1] + f[0,-1])/4.);
      glvertex2d (view, x + Delta_x/2., y - Delta_y/2.);
      color_vertex ((f[] + f[1] + f[1,1] + f[0,1])/4.);
      glvertex2d (view, x + Delta_x/2., y + Delta_y/2.);
      color_vertex ((f[] + f[-1] + f[-1,1] + f[0,1])/4.);
      glvertex2d (view, x - Delta_x/2., y + Delta_y/2.);
      color_vertex ((f[] + f[-1] + f[-1,-1] + f[0,-1])/4.);
      glvertex2d (view, x - Delta_x/2., y - Delta_y/2.);
      glEnd();
      view->ni++;
    }
      }
      else // Z.i > 0
  for (int _i = 0; _i < _N; _i++) /* foreach_leaf */ // fixme: foreach_visible() would be better
    if (f.i < 0 || f[] != nodata) {
      glBegin (GL_TRIANGLE_FAN);
      color_vertex ((4.*f[] +
         2.*(f[1] + f[-1] + f[0,1] + f[0,-1]) +
         f[-1,-1] + f[1,1] + f[-1,1] + f[1,-1])/16.);
      glvertex_normal3d (view, point, fn, x, y, Z[]);
      color_vertex ((f[] + f[-1] + f[-1,-1] + f[0,-1])/4.);
      glvertex_normal3d (view, point, fn, x - Delta_x/2., y - Delta_y/2.,
             (Z[] + Z[-1] + Z[-1,-1] + Z[0,-1])/4.);
      color_vertex ((f[] + f[1] + f[1,-1] + f[0,-1])/4.);
      glvertex_normal3d (view, point, fn, x + Delta_x/2., y - Delta_y/2.,
             (Z[] + Z[1] + Z[1,-1] + Z[0,-1])/4.);
      color_vertex ((f[] + f[1] + f[1,1] + f[0,1])/4.);
      glvertex_normal3d (view, point, fn, x + Delta_x/2., y + Delta_y/2.,
             (Z[] + Z[1] + Z[1,1] + Z[0,1])/4.);
      color_vertex ((f[] + f[-1] + f[-1,1] + f[0,1])/4.);
      glvertex_normal3d (view, point, fn, x - Delta_x/2., y + Delta_y/2.,
             (Z[] + Z[-1] + Z[-1,1] + Z[0,1])/4.);
      color_vertex ((f[] + f[-1] + f[-1,-1] + f[0,-1])/4.);
      glvertex_normal3d (view, point, fn, x - Delta_x/2., y - Delta_y/2.,
             (Z[] + Z[-1] + Z[-1,-1] + Z[0,-1])/4.);
      glEnd();
      view->ni++;     
    }
#else // dimension == 3
      foreach_visible_plane (view, n, alpha)
  if (f.i < 0 || f[] != nodata) {
    coord v[12];
    int m = facets (n, alpha, v, 1.);
    if (m > 2) {
      coord c = {0,0,0};
      for (int i = 0; i < m; i++)
        for (int _d = 0; _d < dimension; _d++)
    c.x += v[i].x/m;
      glBegin (GL_TRIANGLE_FAN);
      color_vertex (interp (point, c, f));
      glvertex3d (view, x + c.x*Delta, y + c.y*Delta, z + c.z*Delta);
      for (int i = 0; i < m; i++) {
        color_vertex (interp (point, v[i], f));
        glvertex3d (view,
        x + v[i].x*Delta, y + v[i].y*Delta, z + v[i].z*Delta);
      }
      color_vertex (interp (point, v[0], f));
      glvertex3d (view,
      x + v[0].x*Delta, y + v[0].y*Delta, z + v[0].z*Delta);
      glEnd ();
      view->ni++;
    }
  }
#endif // dimension == 3
    }
  }
  else { // !linear
#if dimension == 2
    glNormal3d (0, 0, view->reversed ? -1 : 1);
    glBegin (GL_QUADS);
    foreach_visible (view)
      if (f.i < 0 || f[] != nodata) {
  color_facet();
  glvertex2d (view, x - Delta_x/2., y - Delta_y/2.);
  color_facet();
  glvertex2d (view, x + Delta_x/2., y - Delta_y/2.);
  color_facet();
  glvertex2d (view, x + Delta_x/2., y + Delta_y/2.);
  color_facet();
  glvertex2d (view, x - Delta_x/2., y + Delta_y/2.);
  view->ni++;
      }
    glEnd();
#else // dimension == 3
    foreach_visible_plane (view, n, alpha)
      if (f.i < 0 || f[] != nodata) {
  coord v[12];
  int m = facets (n, alpha, v, 1.);
  if (m > 2) {
    glBegin (GL_POLYGON);
    for (int i = 0; i < m; i++) {
      color_facet();
      glvertex3d (view,
      x + v[i].x*Delta, y + v[i].y*Delta, z + v[i].z*Delta);
    }
    glEnd ();
    view->ni++;
  }
      }
#endif // dimension == 3
  }
  if (expr) delete ({col});
#if dimension == 2
  if (zexpr) delete ({Z});
  if (z) delete ((scalar *){fn});
#endif
  return true;
}
 
/**
# *box()*: displays box boundaries and axis coordinates
 
* *notics*: do not draw tick marks (default is false).
* *lc[]*: an array of red, green, blue values between 0 and 1 which
  defines the line color.
* *lw*: the line width.
*/
 
trace
bool box (bool notics = false, float lc[3] = {0}, float lw = 1.)
{
  bview * view = draw();
  coord box[2] = {
    {X0, Y0, Z0},
    {X0 + L0,
     Y0 + L0*Dimensions.y/Dimensions.x
#if dimension > 2     
     , Z0 + L0*Dimensions.z/Dimensions.x
#endif
    }
  };
  coord e;
  double emin = HUGE;
  for (int _d = 0; _d < dimension; _d++) {
    e.x = box[1].x - box[0].x;
    if (e.x < emin)
      emin = e.x;
  }
  draw_lines (view, lc, lw) {
 
    float height = 0.5*gl_StrokeHeight();
    float width = gl_StrokeWidth ('1'), scale = emin/(60.*width), length;
    float Z1 = dimension == 2 ? 0. : box[0].z;
    char label[80];
  
    glMatrixMode (GL_MODELVIEW);
  
#if dimension == 2
    glBegin (GL_LINE_LOOP);
    glvertex2d (view, box[0].x, box[0].y);
    glvertex2d (view, box[1].x, box[0].y);
    glvertex2d (view, box[1].x, box[1].y);
    glvertex2d (view, box[0].x, box[1].y);
    glEnd ();
    view->ni++;
#else // dimension != 2
    for (int _l = 0; _l < 0, serial; _l++) {
      for (int i = -1; i <= 1; i += 2) {
  glBegin (GL_LINE_LOOP);
  glvertex3d (view, box[0].x, box[0].y, z + i*Delta/2.); // fixme
  glvertex3d (view, box[1].x, box[0].y, z + i*Delta/2.);
  glvertex3d (view, box[1].x, box[1].y, z + i*Delta/2.);
  glvertex3d (view, box[0].x, box[1].y, z + i*Delta/2.);
  glEnd ();
  view->ni++;
  glBegin (GL_LINES);
  for (int j = -1; j <= 1; j += 2) { // fixme
    glvertex3d (view, x + i*Delta/2., y + j*Delta/2., z - Delta/2.);
    glvertex3d (view, x + i*Delta/2., y + j*Delta/2., z + Delta/2.);
  }
  glEnd ();
  view->ni++;
      }
    }
#endif // dimension != 2
    
    if (!notics) {
      int nt = 8;
      for (int i = 0; i <= nt; i++) {
  glPushMatrix();
  glTranslatef (X0 + i*e.x/nt - height/2.*scale,
          Y0 - width/3.*scale, Z1);
  glRotatef (-90, 0, 0, 1);
  glScalef (scale, scale, scale);
  sprintf (label, "%g", X0 + i*e.x/nt);
  gl_StrokeString (label);
  glPopMatrix();
 
  glPushMatrix();
  sprintf (label, "%g", Y0 + i*e.y/nt);
  length = gl_StrokeLength (label);
  glTranslatef (X0 - (length + width/3.)*scale,
          Y0 + i*e.y/nt - height/2.*scale, Z1);
  glScalef (scale, scale, scale);
  gl_StrokeString (label);
  glPopMatrix();
 
#if dimension > 2
  glPushMatrix();
  sprintf (label, "%g", Z0 + i*e.z/nt);
  length = gl_StrokeLength (label);
  glTranslatef (X0 - (length + width/3.)*scale,
          Y0, Z0 + i*e.z/nt + height/2.*scale);
  glRotatef (-90, 1, 0, 0);
  glScalef (scale, scale, scale);
  gl_StrokeString (label);
  glPopMatrix();
#endif
      }
 
      glPushMatrix();
      sprintf (label, "%g", X0 + e.x/2.);
      length = gl_StrokeLength (label);
      glTranslatef (X0 + e.x/2 - height*scale,
        Y0 - (length + 4.*width)*scale, Z1);
      glScalef (2.*scale, 2.*scale, 2.*scale);
      gl_StrokeString ("X");
      glPopMatrix();
 
  
      glPushMatrix();
      sprintf (label, "%g", Y0 + e.y/2.);
      length = gl_StrokeLength (label);
      glTranslatef (X0 - (length + 4.*width)*scale,
        Y0 + e.y/2. - height*scale, Z1);
      glScalef (2.*scale, 2.*scale, 2.*scale);
      gl_StrokeString ("Y");
      glPopMatrix();
 
#if dimension > 2
      glPushMatrix();
      sprintf (label, "%g", Z0 + e.z/2.);
      length = gl_StrokeLength (label);
      glTranslatef (X0 - (length + 4.*width)*scale,
        Y0, Z0 + e.z/2. + height*scale);
      glRotatef (-90, 1, 0, 0);
      glScalef (2.*scale, 2.*scale, 2.*scale);
      gl_StrokeString ("Z");
      glPopMatrix();
#endif
    }
  }
  return true;
}
 
/**
# *isosurface()*: displays an isosurface of a field
 
* *f*: the name (as a string) of the field.
* *v*: the value of the isosurface.
* *edges*: whether to draw the edges of isosurface facets.
* *color*: use this field to color each interface fragment.
 
The *min*, *max*, *spread*, *map* etc.  arguments work as described in
[draw_vof()](draw.h#draw_vof-displays-vof-reconstructed-interfaces). */
 
trace
bool isosurface (char * f,
     double v,
     
     char * color = NULL,
     double min = 0, double max = 0, double spread = 0,
     bool linear = false,
     Colormap map = jet,
     float fc[3] = {0}, float lc[3] = {0}, float lw = 1,
     bool expr = false,
     COLORBAR_PARAMS)
{
#if dimension > 2
  if (!f)
    return false;
  
  scalar ff = {-1};
  bool fexpr;
  if (strcmp (f, "level")) {
    ff = compile_expression (f, &fexpr);
    if (ff.i < 0)
      return false;
  }
 
  colorize_args();
 
  scalar fv[];
  for (int _i = 0; _i < _N; _i++) /* foreach_vertex */
    fv[] = (ff[] + ff[-1] + ff[0,-1] + ff[-1,-1] +
      ff[0,0,-1] + ff[-1,0,-1] + ff[0,-1,-1] + ff[-1,-1,-1])/8.;
  
  vector n[];
  for (int _i = 0; _i < _N; _i++) /* foreach */
    for (int _d = 0; _d < dimension; _d++)
      n.x[] = center_gradient(ff);
 
  bview * view = draw();
  glShadeModel (GL_SMOOTH);
  colorize() {
    foreach_visible (view) {
      double val[8] = {
  fv[0,0,0], fv[1,0,0], fv[1,0,1], fv[0,0,1],
  fv[0,1,0], fv[1,1,0], fv[1,1,1], fv[0,1,1]
      };
      double t[5][3][3];
      int nt = polygonize (val, v, t);
      for (int i = 0; i < nt; i++) {
  color_facet();
  glBegin (GL_POLYGON);
  for (int j = 0; j < 3; j++) {
    coord v = {t[i][j][0], t[i][j][1], t[i][j][2]}, np;
    for (int _d = 0; _d < dimension; _d++)
      np.x = interp (point, v, n.x);
    glnormal3d (view, np.x, np.y, np.z);
    if (linear) {
      color_vertex (interp (point, v, col));
    }
    else {
      color_facet();
    }
    glvertex3d (view, x + v.x*Delta_x, y + v.y*Delta_y, z + v.z*Delta_z);
  }
  glEnd ();
  view->ni++;
      }
    }
  }
  if (expr) delete ({col});
  if (fexpr) delete ({ff});
#endif // dimension > 2
  return true;
}
 
/**
# *travelling()*: moves the camera to a different viewpoint
 
* *start*: starting time of the camera motion.
* *end*: time at which the viewpoint should be reached.
* *tx*, *ty*, *quat*, *fov*: definition of the target viewpoint.
*/
 
#define interpo(pv, v)              \
  (!pv ? v : ((t - start)*(pv) + (end - t)*(v))/(end - start))
 
void travelling (double start = 0, double end = 0,
     float tx = 0, float ty = 0, float quat[4] = {0}, float fov = 0)
{
  static float stx, sty, squat[4], sfov;
  static double told = -1.;
  if (told < start && t >= start) {
    bview * view = get_view();
    stx = view->tx, sty = view->ty, sfov = view->fov;
    for (int i = 0; i < 4; i++)
      squat[i] = view->quat[i];
  }
  if (t >= start && t <= end)
    view (tx = interpo (tx, stx), ty = interpo (ty, sty),
    fov = interpo (fov, sfov),
    quat = {interpo(quat[0], squat[0]), interpo(quat[1], squat[1]),
            interpo(quat[2], squat[2]), interpo(quat[3], squat[3])});
  if (told < end && t >= end) {
    bview * view = get_view();
    stx = view->tx, sty = view->ty, sfov = view->fov;
    for (int i = 0; i < 4; i++)
      squat[i] = view->quat[i];
  }
  told = t;  
}
 
#undef interpo
 
/**
# *draw_string()*: draws strings on a separate layer (for annotations)
 
* *str*: string to display.
* *pos*: position: "0" bottom left, "1" top left, "2" top right 
  and "3" bottom right (default 0).
* *size*: the size of the text, given as the number of characters
   which can fit within the width of the screen. Default is 40.
* *lc[]*: an array of red, green, blue values between 0 and 1 which
  defines the text color.
* *lw*: the line width.
*/
 
trace
bool draw_string (char * str,
      int pos = 0,
      float size = 40,
      float lc[3] = {0}, float lw = 1)
 
{
  bview * view = draw();
  
  glMatrixMode (GL_PROJECTION);
  glPushMatrix();             
  glLoadIdentity();
 
  glMatrixMode (GL_MODELVIEW);
  glPushMatrix();
  glLoadIdentity();
    
  glColor3f (lc[0], lc[1], lc[2]);
  glLineWidth (view->samples*(lw > 0. ? lw : 1.));
 
  float width  = gl_StrokeWidth ('1'), height = gl_StrokeHeight();
  if (!size)
    size = 40;
  float hscale = 2./(size*width), vscale = hscale*view->width/view->height;
  float vmargin = width/2.*vscale;
  if (pos == 0)
    glTranslatef (-1., -1. + vmargin, 0.);
  else if (pos == 1)
    glTranslatef (-1., 1. - height*vscale, 0.);
  else if (pos == 2)
    glTranslatef (1. - strlen(str)*width*hscale, 1. - height*vscale, 0.);
  else
    glTranslatef (1. - strlen(str)*width*hscale, -1. + vmargin, 0.);    
  glScalef (hscale, vscale, 1.);
  gl_StrokeString (str); 
  
  glMatrixMode (GL_MODELVIEW);
  glPopMatrix();
  glMatrixMode (GL_PROJECTION);
  glPopMatrix();  
 
  return true;
}
 
/**
# *labels()*: displays label fields */
 
trace
bool labels (char * f, float lc[3] = {0}, float lw = 1, int level = -1)
{
#if dimension == 2
  bool expr = false;
  scalar ff = compile_expression (f, &expr);
  if (ff.i < 0)
    return false;
  bview * view = draw();
  float width  = gl_StrokeWidth ('1'), height = gl_StrokeHeight();
  float res = view->res;
  if (view->res < 150*view->samples)
    view->res = 150*view->samples;
  int maxlevel = view->maxlevel;
  view->maxlevel = level;
  draw_lines (view, lc, lw) {
    glMatrixMode (GL_MODELVIEW);
    foreach_visible (view)
      if (ff[] != nodata) {
  glPushMatrix();
  char s[80];
  sprintf (s, "%g", ff[]);
  float scale = 0.8*Delta_x/(strlen(s)*width);
  glTranslatef (x - 0.4*Delta_x, y - scale*height/3., 0.);
  glScalef (scale, scale, 1.);
  gl_StrokeString (s);
  glPopMatrix();
      }
  }
  view->res = res;
  view->maxlevel = maxlevel;
  if (expr) delete ({ff});
  return true;
#else // dimension == 3
  fprintf (stderr, "labels() is not implemented in 3D yet\n");
  return false;
#endif // dimension == 3
}
 
/**
# *lines()*: from a file.
 
* *file*: the gnuplot-formatted file containing the polyline(s).
* *lc[]*: an array of red, green, blue values between 0 and 1 which
  defines the line color.
* *lw*: the line width.
*/
 
trace
bool lines (char * file, float lc[3] = {0}, float lw = 1.)
{
#if dimension != 2
  assert (false);
#else // dimension == 2
  if (!file) {
    fprintf (stderr, "lines(): file must be specified\n");
    return false;
  }
  FILE * fp = fopen (file, "r");
  if (!fp) {
    perror (file);
    return false;
  }
  bview * view = draw();
  draw_lines (view, lc, lw) {
    bool line = false;
    int c = fgetc (fp);
    while (c != EOF) {
      while (c == ' ' || c == '\t') c = fgetc (fp);
      if (c == '.' || c == '+' || c == '-' || isdigit (c)) {
  ungetc (c, fp);
  double x, y;
  if (fscanf (fp, "%lf %lf", &x, &y) == 2) {
    if (!line) {
      glBegin (GL_LINE_STRIP);
      line = true;
    }
    glvertex2d (view, x, y);
    while ((c = fgetc(fp)) == ' ' || c == '\t');
    if (c == '\n') c = fgetc (fp);
    view->ni++;
  }
  else // ignore the rest of the line
    while ((c = fgetc(fp)) != EOF && c != '\n');
      }
      else if (c == '#')
  while ((c = fgetc(fp)) != EOF && c != '\n');
      else {
  if (line)
    glEnd(), line = false;
  c = fgetc (fp);
      }
    }
    if (line)
      glEnd();
  }
  fclose (fp);
#endif // dimension == 2
  return true;
}
 
/**
# Interface export
 
This is used by [bview](bview/README) to automatically generate the
user interface. */
 
#include "draw_json.h"
 
struct {
  int (* json) (char * s, int len);
} bview_interface[] = {
  { _draw_vof_json },
  { _squares_json },
  { _cells_json },
  { _box_json },
#if dimension == 2
  { _isoline_json },
  { _labels_json },
  { _vectors_json },
  { _lines_json },
#else // dimension == 3
  { _isosurface_json },
#endif
  { NULL }
};