text: /**
  # Volumetric visualization
  
  A proof of a inefficient concept.
  
  ![It works a bit, perhaps a bit of shading helps to see
   depth.](plume/mov.mp4)
  */
  #include "grid/multigrid3D.h"
  #include "navier-stokes/centered.h"
  #include "tracer.h"
  #include "integrator2.h"
  
  scalar s[], * tracers = {s};
  
  double sigma = 0.1; // size
  double att = 0.1;   // Attenuation coefficient
  double prf = 4e3;    // Light intensity of some sort
  
  trace
  double get_c (ray r, scalar Sl, scalar Sc, scalar s) {
    double c = 0;
    foreach_ray_cell_intersection(r) {
      double l = 0;
      foreach_dimension()
        l += sq(_a[1].x - _a[0].x);
      l = l > 0 ? sqrt(l) : 0;
      c += l*
        interpolate(s, (_a[0].x +_a[1].x)/2, (_a[0].y +_a[1].y)/2, (_a[0].z +_a[1].z)/2)*
        interpolate(Sl, (_a[0].x +_a[1].x)/2, (_a[0].y +_a[1].y)/2, (_a[0].z +_a[1].z)/2)*
        interpolate(Sc, (_a[0].x +_a[1].x)/2, (_a[0].y +_a[1].y)/2, (_a[0].z +_a[1].z)/2);
    }
    return fabs(c);
  }
  
  trace
  void make_a_frame (scalar s, FILE * fp) {
    scalar Sl[], Sc[];
    coord ildir = {1, 1, 3};
    normalize (&ildir);
    integrate_dn (Sl, s, ildir);
    foreach()
      Sl[] = exp(-Sl[]/att);
    boundary ({Sl});
    coord icdir = {0.01, 0.01, 1};
    normalize (&icdir);
    integrate_dn (Sc, s, icdir);
    foreach()
      Sc[] = exp(-Sc[]/att);
    boundary ({Sc});
    fprintf (fp, "P6\n%d %d\n%d\n", cam.nx, cam.ny, 255); // .ppm header
    foreach_ray() {
      unsigned char px[3] = {10.*(ii + jj)/cam.nx + 20,
  			   25.*(ii + jj)/cam.nx + 120,
  			   50.*jj/cam.nx + 120};
      double c = get_c (_r, Sl, Sc, s);
      double occ = fabs(BC_int(_r.O, _r.dir, s));
      for (int i = 0; i < 3; i++) 
        px[i] =  min(max(px[i]*exp(-occ/att) + prf*c, 1), 255);
      fwrite (px, 1, 3, fp);
    }
    fflush (fp);
  }
  
  FILE * fp;
  int main() {
    X0 = Y0 = Z0 = -0.5;
    N = 64;
    fp = popen ("ppm2mp4 -r 10 mov.mp4", "w");
    run();
  }
  
  event init (t = 0) {
    fraction (s, sq(sigma) - (sq(x) + sq(y + 0.25) + sq(z)));
    foreach() 
      u.y[] = 0.4*s[];
    boundary ({s, u});
  }
  
  event frames (t += 0.2) {
    make_a_frame (s, fp);
  }
  
  event stop (t = 5) {
    pclose (fp);
    return 1;
  }
  