#include "rmannotes.sl"

displacement 
boiling(
	float     minfreq = 3,
                  maxfreq = 10, 
                  plenitude = 0.3,    /* range: [0,1] */
                  Km = 0.2,
                  boil = 0.2,
	          ftime = 0.5;
	)
{
  float   surface_mag,
          layer_mag;
  float   fuzz = 0.005;
  float   ss, tt;
  float   row, col;
  float   freq, mag;
  point   diff_N;
  float   r, theta, angle;
  float   d, d0, d1;
  float   cx, cy;
  float   noifreq = 5;
  float   noiscale = 0.1;
  float   bub;
  float   sample;


  surface_mag = 0.0;


  /* bubble layer */

  for ( freq=minfreq; freq<maxfreq; freq=freq+1 )
    {
      angle = PI * snoise(freq * 16.31456);
      rotate2d(s, t, angle, 0.5, 0.5, cx, cy);
      ss = repeat(cx, freq);
      tt = repeat(cy, freq);
      col = whichtile(cx, freq);
      row = whichtile(cy, freq);

      bub = 0.5 + snoise2(row * 578.2 + col * 13.2 + 0.053,
			  col * 325.7 + row * 5.7 + 0.051);


      if (abs(bub - ftime) < 0.5 * boil)
	{
	  ss = ss + noiscale * snoise(snoise2(s * noifreq, t * noifreq) + 912);
	  tt = tt + noiscale * snoise(snoise2(s * noifreq, t * noifreq) + 333);
	  cx = 0.5 + 0.1 * snoise2(row * 12.312, col * 21.773);
	  cy = 0.5 + 0.1 * snoise2(row * 28.398, col * 62.112);
	  d = distance((cx, cy, 0), (ss, tt, 0));

	  if (d > 0.3)
	    mag = 0.0;
	  else
	    mag = ((0.5 * boil - abs(bub - ftime)) / boil)
	      * 60 * (0.09 - d * d) 
		* ((maxfreq - freq) / maxfreq);

	  layer_mag = mag;
	  surface_mag = max(surface_mag, layer_mag);
	}
    }


  /* waves layer */

  topolar2d(s - 0.5, t - 0.5, r, theta);
  sample = r * 12 * PI - ftime * 8 * PI;
  layer_mag = 0.4 * sin(sample);

  surface_mag += layer_mag;



  diff_N = normalize(N) - normalize(Ng);
  P += Km * surface_mag * normalize(N);
  N = normalize(calculatenormal(P)) + diff_N;
}