;Translated from: JAVA REFERENCE IMPLEMENTATION OF IMPROVED NOISE - COPYRIGHT 2002 KEN PERLIN.
;http://mrl.nyu.edu/~perlin/noise/

function fade,t
  return, t * t * t * (t * (t * 6 - 15) + 10)
end
function lerp,t,a,b
  return, a + t * (b - a)
end
function grad,hash,x,y,z
  ;Convert lo 4 bits of hash code into 12 gradient directions
  h = hash and 15;
  u=x*0
  v=x*0
  hlt8=h lt 8

  whlt8=where(hlt8,n_whlt8,complement=wnhlt8,ncomplement=n_wnhlt8)
  if n_whlt8 gt 0 then u[whlt8]=x[whlt8]
  if n_wnhlt8 gt 0 then u[wnhlt8]=y[wnhlt8]

  part1=h lt 4
  part2=logical_or(h eq 12, h eq 14)

  wy=where(part1,n_wy)
  wx=where(logical_and(~part1,part2),n_wx)
  wz=where(logical_and(~part1,~part2),n_wz)
  if n_wy gt 0 then v[wy]=y[wy]
  if n_wx gt 0 then v[wx]=x[wx]
  if n_wz gt 0 then v[wz]=z[wz]
  acc=x*0;
  and1=(h and 1) eq 0
  and2=(h and 2) eq 0
  wand1=where(and1,n_wand1,complement=wnand1,ncomplement=n_wnand1)
  wand2=where(and2,n_wand2,complement=wnand2,ncomplement=n_wnand2)
  if n_wand1 gt 0 then acc[wand1]+=u[wand1]
  if n_wnand1 gt 0 then acc[wnand1]-=u[wnand1]
  if n_wand2 gt 0 then acc[wand2]+=v[wand2]
  if n_wnand2 gt 0 then acc[wnand2]-=u[wnand2]
  return,acc
end

function noise,x_,y_,z_
  common noise_static,p
  if n_elements(p) eq 0 then begin
    p = [ $
      151,160,137, 91, 90, 15,131, 13,201, 95, 96, 53,194,233,  7,225, $
      140, 36,103, 30, 69,142,  8, 99, 37,240, 21, 10, 23,190,  6,148, $
      247,120,234, 75,  0, 26,197, 62, 94,252,219,203,117, 35, 11, 32, $
       57,177, 33, 88,237,149, 56, 87,174, 20,125,136,171,168, 68,175, $
       74,165, 71,134,139, 48, 27,166, 77,146,158,231, 83,111,229,122, $
       60,211,133,230,220,105, 92, 41, 55, 46,245, 40,244,102,143, 54, $
       65, 25, 63,161,  1,216, 80, 73,209, 76,132,187,208, 89, 18,169, $
      200,196,135,130,116,188,159, 86,164,100,109,198,173,186,  3, 64, $
       52,217,226,250,124,123,  5,202, 38,147,118,126,255, 82, 85,212, $
      207,206, 59,227, 47, 16, 58, 17,182,189, 28, 42,223,183,170,213, $
      119,248,152,  2, 44,154,163, 70,221,153,101,155,167, 43,172,  9, $
      129, 22, 39,253, 19, 98,108,110, 79,113,224,232,178,185,112,104, $
      218,246, 97,228,251, 34,242,193,238,210,144, 12,191,179,162,241, $
       81, 51,145,235,249, 14,239,107, 49,192,214, 31,181,199,106,157, $
      184, 84,204,176,115,121, 50, 45,127,  4,150,254,138,236,205, 93, $
      222,114, 67, 29, 24, 72,243,141,128,195, 78, 66,215, 61,156,180  $
    ]
    p=[p,p]
  end
  XX = floor(x_) and 255                  ;// FIND UNIT CUBE THAT
  YY = floor(y_) and 255    ;              // CONTAINS POINT.
  ZZ = floor(z_) and 255;
  x=x_-floor(x_);                     ;           // FIND RELATIVE X,Y,Z
  y=y_-floor(y_);                      ;          // OF POINT IN CUBE.
  z=z_-floor(z_);
  u = fade(x);                       ;         // COMPUTE FADE CURVES
  v = fade(y);                        ;        // FOR EACH OF X,Y,Z.
  w = fade(z);
  A = p[XX  ]+YY      ;// HASH COORDINATES OF
  AA = p[A]+ZZ        ;// THE 8 CUBE CORNERS,
  AB = p[A+1]+ZZ
  B = p[XX+1]+YY
  BA = p[B]+ZZ
  BB = p[B+1]+ZZ;

      return,lerp(w, lerp(v, lerp(u, grad(p[AA  ], x  , y  , z   ),  $;// AND ADD
                                     grad(p[BA  ], x-1, y  , z   )), $;// BLENDED
                             lerp(u, grad(p[AB  ], x  , y-1, z   ),  $;// RESULTS
                                     grad(p[BB  ], x-1, y-1, z   ))),$;// FROM  8
                     lerp(v, lerp(u, grad(p[AA+1], x  , y  , z-1 ),  $;// CORNERS
                                     grad(p[BA+1], x-1, y  , z-1 )), $;// OF CUBE
                             lerp(u, grad(p[AB+1], x  , y-1, z-1 ),  $
                                     grad(p[BB+1], x-1, y-1, z-1 )))) ;
end

function mnoise,x,y,z,lambda=lambda,omega=omega,octaves=octaves
  ;These are all POV-Ray defaults
  if n_elements(lambda) eq 0 then lambda=2
  if n_elements(omega) eq 0 then omega=0.5
  if n_elements(octaves) eq 0 then octaves=6
  value=x*0
  o=1.0
  l=1.0
  for i=0,octaves-1 do begin
    value+=o*noise(l*x,l*y,l*z)
    o*=omega
    l*=lambda
  end
  return,value
end

function noisegrid,xrange,n_x,yrange,n_y,zrange,n_z,lambda=lambda,omega=omega,octaves=octaves
  x=linterp(0,xrange[0],n_x,xrange[1],findgen(n_x))
  y=linterp(0,yrange[0],n_y,yrange[1],findgen(n_y))
  z=linterp(0,zrange[0],n_z,zrange[1],findgen(n_z))

  x=rebin(reform(x,n_x,1,1,/overwrite),n_x,n_y,n_z)
  y=rebin(reform(y,1,n_y,1,/overwrite),n_x,n_y,n_z)
  z=rebin(reform(z,1,1,n_z,/overwrite),n_x,n_y,n_z)
  return,mnoise(x,y,z,lambda=lambda,omega=omega,octaves=octaves)
end