// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 0
// Copyright 2003 by Tor Olav Kristensen
// Email: t o r . o l a v . k [_A_T_] g m a i l . c o m
// http://subcube.com
// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 0

#version 3.5;

global_settings {
  ambient_light 1
}

#include "functions.inc" // For f_r() and f_torus
#include "textures.inc"
#include "glass.inc"

#declare Clamp_Fn = function(N, Min, Max) { max(Min, min(N, Max)) }

#declare Nil = 1E-9;

#declare Acos_Fn = function(Cos) { acos(Clamp_Fn(Cos, 0 + Nil, 1 - Nil)) }

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 0
// Macro to build "Angle Of Incidence" functions

#macro AOI_Function(Fn, pV)

  #local pVx = pV.x;
  #local pVy = pV.y;
  #local pVz = pV.z;

  #local H = 1E-4;
  #local H2 = H*2;

  #local Temp0Fn =
    function(vNx, vNy, vNz, vVx, vVy, vVz, Div) {
      abs(vNx*vVx + vNy*vVy + vNz*vVz)/Div
    }

  #local Temp1Fn =
    function(vNx, vNy, vNz, vVx, vVy, vVz) {
      Temp0Fn(
        vNx, vNy, vNz,
        vVx, vVy, vVz,
        f_r(vNx, vNy, vNz)*f_r(vVx, vVy, vVz)
      )
    }

  function {
    Temp1Fn(
      (Fn(x + H, y, z) - Fn(x - H, y, z))/H2,
      (Fn(x, y + H, z) - Fn(x, y - H, z))/H2,
      (Fn(x, y, z + H) - Fn(x, y, z - H))/H2,
      pVx - x,
      pVy - y,
      pVz - z
    )
  }

#end // macro AOI_Function

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 0

#macro PisoFunction(Pd, Ps, StdDev, Fn, pI, pR)

  #local pIx = pI.x;
  #local pIy = pI.y;
  #local pIz = pI.z;
  #local pRx = pR.x;
  #local pRy = pR.y;
  #local pRz = pR.z;

  #local Nil = 1E-9;
  #local SqSD = StdDev*StdDev;
  #local FourPi = 4*pi;
  #local H = 1E-4;
  #local H2 = H*2;

  #local ClampFn = function(N) { max(0, min(N, 1)) }
  #local DotFn = function(x, y, z, xx, yy, zz) { x*xx + y*yy + z*zz }

  #local PisoFn =
    function(CosAi, CosAr, Ah) {
      ClampFn(
        Pd/pi +
//        select(CosAr*CosAi, 0, 1)*
        Ps/(Nil + sqrt(abs(CosAi*CosAr)))*exp(-pow(tan(Ah), 2)/SqSD)/(FourPi*SqSD)
      )
    }

  #local AhFn =
    function(vNxn, vNyn, vNzn, vHx, vHy, vHz) {
      acos(
        DotFn(
          vNxn, vNyn, vNzn,
          vHx, vHy, vHz
        )/(Nil + f_r(vHx, vHy, vHz))
      )
    }

  #local TempFn0 =
    function(vNxn, vNyn, vNzn, vIxn, vIyn, vIzn, vRxn, vRyn, vRzn) {
      PisoFn(
        DotFn(
          vNxn, vNyn, vNzn,
          vIxn, vIyn, vIzn
        ),
        DotFn(
          vNxn, vNyn, vNzn,
          vRxn, vRyn, vRzn
        ),
        AhFn(
          vNxn, vNyn, vNzn,
          vIxn + vRxn, vIyn + vRyn, vIzn + vRzn
        )
      ) 
    }

  #local TempFn1 =
    function(vNx, vNy, vNz, lN, vIx, vIy, vIz, lI, vRx, vRy, vRz, lR) {
      TempFn0(
        vNx/lN, vNy/lN, vNz/lN,
        vIx/lI, vIy/lI, vIz/lI,
        vRx/lR, vRy/lR, vRz/lR
      )
    }

  #local TempFn2 =
    function(vNx, vNy, vNz, vIx, vIy, vIz, vRx, vRy, vRz) {
      TempFn1(
        vNx, vNy, vNz, Nil + f_r(vNx, vNy, vNz),
        vIx, vIy, vIz, Nil + f_r(vIx, vIy, vIz),
        vRx, vRy, vRz, Nil + f_r(vRx, vRy, vRz)
      ) 
    }
  
  function {
    TempFn2(
      (Fn(x + H, y, z) - Fn(x - H, y, z))/H2,
      (Fn(x, y + H, z) - Fn(x, y - H, z))/H2,
      (Fn(x, y, z + H) - Fn(x, y, z - H))/H2,

      pIx - x,
      pIy - y,
      pIz - z,

      pRx - x,
      pRy - y,
      pRz - z
    )
  }

#end // macro PisoFunction

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 0

#macro PhongFunction(Fn, pI, pV, Ks, Ls, Shininess)

  #local pIx = pI.x;
  #local pIy = pI.y;
  #local pIz = pI.z;
  #local pVx = pV.x;
  #local pVy = pV.y;
  #local pVz = pV.z;

  #local Nil = 1E-9;
  #local H = 1E-4;
  #local H2 = H*2;

  #local ClampFn = function(N) { max(0, min(N, 1)) }
  #local DotFn = function(x, y, z, xx, yy, zz) { x*xx + y*yy + z*zz }

  #local SpecularFn0 =
    function(Phi) {
      select(Phi, 0, 1)*ClampFn(Ks*Ls*pow(abs(Phi), Shininess), 0, 1)
    }

  #local SpecularFn1 =
    function(S, vNxn, vNyn, vNzn, vIxn, vIyn, vIzn, vVxn, vVyn, vVzn) {
      SpecularFn0(
        DotFn(
          vVxn, vVyn, vVzn,
          S*vNxn - vIxn, S*vNyn - vIyn, S*vNzn - vIzn
        )
      ) 
    }

  #local SpecularFn2 =
    function(vNxn, vNyn, vNzn, vIxn, vIyn, vIzn, vVxn, vVyn, vVzn) {
      SpecularFn1(
        DotFn(
          vIxn, vIyn, vIzn,
          vNxn, vNyn, vNzn
        )*2,
        vNxn, vNyn, vNzn,
        vIxn, vIyn, vIzn,
        vVxn, vVyn, vVzn
      ) 
    }

  #local SpecularFn3 =
    function(vNx, vNy, vNz, lN, vIx, vIy, vIz, lI, vVx, vVy, vVz, lV) {
      SpecularFn2(
        vNx/lN, vNy/lN, vNz/lN,
        vIx/lI, vIy/lI, vIz/lI,
        vVx/lV, vVy/lV, vVz/lV
      ) 
    }

  #local SpecularFn4 =
    function(vNx, vNy, vNz, vIx, vIy, vIz, vVx, vVy, vVz) {
      SpecularFn3(
        vNx, vNy, vNz, Nil + f_r(vNx, vNy, vNz),
        vIx, vIy, vIz, Nil + f_r(vIx, vIy, vIz),
        vVx, vVy, vVz, Nil + f_r(vVx, vVy, vVz)
      ) 
    }

  function {
    SpecularFn4(
      (Fn(x + H, y, z) - Fn(x - H, y, z))/H2, // vNx
      (Fn(x, y + H, z) - Fn(x, y - H, z))/H2, // vNy
      (Fn(x, y, z + H) - Fn(x, y, z - H))/H2, // vNz

      pIx - x, // vIx
      pIy - y, // vIy
      pIz - z, // vIz

      pVx - x, // vVx
      pVy - y, // vVy
      pVz - z  // vVz
    )
  }

#end // macro PhongFunction

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 0

#macro ShadeFunction(Fn, pI, Kd, Ld, AttFn)

  #local pIx = pI.x;
  #local pIy = pI.y;
  #local pIz = pI.z;

  #local Nil = 1E-9;
  #local H = 1E-4;
  #local H2 = H*2;

  #local ClampFn = function(N) { max(0, min(N, 1)) }
  #local DotFn = function(x, y, z, xx, yy, zz) { x*xx + y*yy + z*zz }

  #local DiffuseFn0 =
    function(Theta, lI) {
      ClampFn(Kd*Ld*Theta*AttFn(lI))
    }

  #local DiffuseFn1 =
    function(vNxn, vNyn, vNzn, vIxn, vIyn, vIzn, lI) {
      DiffuseFn0(
        DotFn(
          vIxn, vIyn, vIzn,
          vNxn, vNyn, vNzn
        ),
        lI
      ) 
    }

  #local DiffuseFn2 =
    function(vNx, vNy, vNz, lN, vIx, vIy, vIz, lI) {
      DiffuseFn1(
        vNx/lN, vNy/lN, vNz/lN,
        vIx/lI, vIy/lI, vIz/lI,
        lI
      ) 
    }

  #local DiffuseFn3 =
    function(vNx, vNy, vNz, vIx, vIy, vIz) {
      DiffuseFn2(
        vNx, vNy, vNz, Nil + f_r(vNx, vNy, vNz),
        vIx, vIy, vIz, Nil + f_r(vIx, vIy, vIz)
      ) 
    }

  function {
    DiffuseFn3(
      (Fn(x + H, y, z) - Fn(x - H, y, z))/H2, // vNx
      (Fn(x, y + H, z) - Fn(x, y - H, z))/H2, // vNy
      (Fn(x, y, z + H) - Fn(x, y, z - H))/H2, // vNz

      pIx - x, // vIx
      pIy - y, // vIy
      pIz - z  // vIz
    )
  }

#end // macro ShadeFunction


#macro PhongFunction(ShapeFn, pCamera, pLight, AttFn, Iambient, Idiffuse, Ispecular, Shininess)

  #local pLx = pLight.x;
  #local pLy = pLight.y;
  #local pLz = pLight.z;
  #local pCx = pCamera.x;
  #local pCy = pCamera.y;
  #local pCz = pCamera.z;

  #local Nil = 1E-9;
  #local H = 1E-4;
  #local H2 = H*2;

  #local ClampFn = function(N) { max(0, min(N, 1)) }
  #local DotFn = function(x, y, z, xx, yy, zz) { x*xx + y*yy + z*zz }

  #local PhongFn0 =
    function(lI, Theta, Phi) {
      ClampFn(
        AttFn(lI)*
        (
          Iambient
         +select(Theta, 0, Idiffuse)*Theta
         +select(Phi, 0, Ispecular)*pow(abs(Phi), Shininess)
        )
      )
    }

  #local PhongFn1 =
    function(lI, vLxn, vLyn, vLzn, vNxn, vNyn, vNzn, vCxn, vCyn, vCzn, S) {
      PhongFn0(
        lI,
        DotFn(
          vLxn, vLyn, vLzn,
          vNxn, vNyn, vNzn
        ),
        DotFn(
          vCxn, vCyn, vCzn,
          S*vNxn - vLxn, S*vNyn - vLyn, S*vNzn - vLzn
        )
      ) 
    }

  #local PhongFn2 =
    function(lI, vLxn, vLyn, vLzn, vNxn, vNyn, vNzn, vCxn, vCyn, vCzn) {
      PhongFn1(
        lI,
        vLxn, vLyn, vLzn,
        vNxn, vNyn, vNzn,
        vCxn, vCyn, vCzn,
        DotFn(
          vLxn, vLyn, vLzn,
          vNxn, vNyn, vNzn
        )*2
      ) 
    }

  #local PhongFn3 =
    function(vLx, vLy, vLz, lI, vNx, vNy, vNz, lN, vCx, vCy, vCz, lV) {
      PhongFn2(
        lI,
        vLx/lI, vLy/lI, vLz/lI,
        vNx/lN, vNy/lN, vNz/lN,
        vCx/lV, vCy/lV, vCz/lV
      ) 
    }

  #local PhongFn4 =
    function(vLx, vLy, vLz, vNx, vNy, vNz, vCx, vCy, vCz) {
      PhongFn3(
        vLx, vLy, vLz, Nil + f_r(vLx, vLy, vLz),
        vNx, vNy, vNz, Nil + f_r(vNx, vNy, vNz),
        vCx, vCy, vCz, Nil + f_r(vCx, vCy, vCz)
      ) 
    }

  function {
    PhongFn4(
      pLx - x, // vLx
      pLy - y, // vLy
      pLz - z, // vLz
      (ShapeFn(x + H, y, z) - ShapeFn(x - H, y, z))/H2, // vNx
      (ShapeFn(x, y + H, z) - ShapeFn(x, y - H, z))/H2, // vNy
      (ShapeFn(x, y, z + H) - ShapeFn(x, y, z - H))/H2, // vNz
      pCx - x, // vCx
      pCy - y, // vCy
      pCz - z  // vCz
    )
  }

#end // macro PhongFunction

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 0

#declare pCameraLocation = <0.0, 3.4, -5.1>;
#declare pCameraLookAt = <0.0, -2.6, 0.0>;

#declare pLightLocation =
//  pCameraLocation;
//  <3, 3, -2>*50;
  <6/2, 2, -2*2>;

camera {
  location pCameraLocation
  look_at pCameraLookAt
  angle 100
}

light_source { pLightLocation color rgb <1, 1, 1> shadowless }

sky_sphere {
  pigment {
    gradient y
    turbulence 0.1
    color_map {
      [ 0.0 color rgb <1.5, 1.5, 1.5>/2 ]
      [ 0.4 color rgb <0.0, 0.0, 0.5>/2 ]
      [ 0.6 color rgb <0.0, 0.0, 0.5>/2 ]
      [ 1.0 color rgb <1.5, 1.5, 1.5>/2 ]
    }
  }
}

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 0
// Isosurface with "Angle Of Incidence" controlled texture

#declare Ambient = 0.1;
#declare Diffuse = 0.5;
#declare Specular = 0.4;
#declare Shininess = 50;
#declare Attenuation_Fn = function(D) { 1/(1 + 0.00*D*D) }

#declare Rmaj = 4;
#declare Rmin = 1;

#declare Torus_Fn = function { f_torus(x, y, z, Rmaj, Rmin) }
#declare Sphere_Fn = function { f_sphere(x, y, z, Rmaj) }

#declare GlassMaterial =
  material {
    texture {
      AOI_Function(Torus_Fn, pCameraLocation)
      texture_map {
        [ 0 Orange_Glass ]
        [ 1 Vicks_Bottle_Glass ]
      }
    }
    interior { I_Glass }
  }


#declare Ambient_ = 1;
#declare Diffuse_ = 0;

#declare TexA = 
  texture {
    pigment { color rgb <0, 0, 0> }
    finish { ambient Ambient_ diffuse Diffuse_ }
  }

#declare TexB = 
  texture {
    pigment { color rgb <1, 1, 1> }
    finish { ambient Ambient_ diffuse Diffuse_ }
  }

#declare RedTex = 
  texture {
    pigment { color red 1 }
    finish { ambient Ambient_ diffuse Diffuse_ }
  }

#declare GreenTex = 
  texture {
    pigment { color green 1 }
    finish { ambient Ambient_ diffuse Diffuse_ }
  }

#declare BlueTex = 
  texture {
    pigment { color blue 1 }
    finish { ambient Ambient_ diffuse Diffuse_ }
  }

#declare YellowTex = 
  texture {
    pigment { color rgb <1, 1, 0> }
    finish { ambient Ambient_ diffuse Diffuse_ }
  }

#declare OrangeTex = 
  texture {
    pigment { color rgb <0.8, 0.7, 0> }
    finish { ambient Ambient_ diffuse Diffuse_ }
  }

#declare Simple_Map =
  texture_map {
    [ 0 TexA ]
    [ 1 TexB ]
  }

#declare AOI_Map =
  texture_map {
    [ 0.00 RedTex ]
    [ 0.50 TexA ]
    [ 0.90 TexB ]
    [ 1.00 GreenTex ]
  }

#declare S_Map =
  texture_map {
    [ 0.10 TexA ]
    [ 0.20 RedTex ]
    [ 0.30 TexA ]
    [ 0.40 TexB ]
    [ 0.50 RedTex ]
    [ 0.60 TexB ]
    [ 0.90 TexA ]
    [ 1.00 GreenTex ]
  }

#declare Shade_Map =
  texture_map {
    [ 0.00 TexA ]
    [ 0.10 RedTex ]
    [ 0.20 TexA ]
    [ 0.95 TexB ]
    [ 1.00 GreenTex ]
  }

#declare Ring_Map =
  texture_map {
    [ 0.0 RedTex ]
    [ 0.001 TexA ]
    [ 0.2 TexB ]
    [ 0.3 TexA ]
    [ 0.4 TexB ]
    [ 0.5 TexA ]
    [ 0.5 TexB ]
    [ 0.6 TexA ]
    [ 0.7 TexB ]
    [ 0.8 TexA ]
    [ 0.999 TexB ]
    [ 1.0 GreenTex ]
  }

#declare Metal_Map =
  texture_map {
//    [ 0.0 Bright_Bronze ]
    [ 0.2 Soft_Silver ]
//    [ 0.2 Bright_Bronze ]
    [ 0.25 OrangeTex ]
//    [ 0.4 GreenTex ]
//    [ 0.5 New_Penny ]
    [ 0.3 Bright_Bronze ]
    [ 0.7 RedTex ]
    [ 0.9 TexB ]
  }

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 0

#declare Shade_Fn =
  ShadeFunction(
    Torus_Fn,
    pLightLocation,
    1, // Kd
    1, // Ld
    function(D) { 1/(1 + 0.008*D*D) }
  )

#declare Phong_Fn =
  PhongFunction(
    Torus_Fn,
    pCameraLocation,
    pLightLocation,
    Attenuation_Fn,
    Ambient,
    Diffuse,
    Specular,
    Shininess
  )

#declare Piso_Fn =
  PisoFunction(
    0.4, // Pd
    0.6, // Ps
    0.6, // StdDev
    Torus_Fn,
    pLightLocation,
    pCameraLocation
  )

#declare AOI_Fn =
  AOI_Function(
    Torus_Fn,
    pCameraLocation
  )

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 0

#declare SimpleMaterial =
  material {
    texture {
//      function { Shade_Fn(x, y, z) }
      function { Phong_Fn(x, y, z) }
//      function { Piso_Fn(x, y, z) }
//      function { AOI_Fn(x, y, z) }

//      texture_map { Shade_Map }
//      texture_map { AOI_Map }
//      texture_map { S_Map }
//      texture_map { Simple_Map }
//      texture_map { Ring_Map }
      texture_map { Metal_Map }
    }
  }

#declare SimpleTexture =
  texture {
    pigment { color rgb <1, 1, 1> }
    finish { ambient Ambient diffuse Diffuse phong 10 }
  }

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 0

#declare pMin = -<Rmaj + Rmin, Rmin, Rmaj + Rmin>;
#declare pMax =  <Rmaj + Rmin, Rmin, Rmaj + Rmin>;

isosurface {
  function { Torus_Fn(x, y, z) }
  contained_by { box { pMin, pMax } }
//  contained_by { sphere { <0, 0, 0>, Rmaj } }
//  material { GlassMaterial }
  material { SimpleMaterial }
//  texture { SimpleTexture }
}

/*
torus {
  Rmaj, Rmin
  pigment {
    AOI_Function(Torus_Fn, pCameraLocation)
    color_map {
      [ 0.0 color rgb <1.5, 0.0, 0.0> ]
      [ 1.0 color rgb <0.0, 1.5, 0.0> ]
    }
  }
}
*/

/*
#declare SlopeX = function { pigment { slope x } }
#declare SlopeY = function { pigment { slope y } }
#declare SlopeZ = function { pigment { slope z } }

torus {
  Rmaj, Rmin
  pigment {
    slope x
    color_map {
      [ 0.0 color rgb <1.5, 0.0, 0.0> ]
      [ 1.0 color rgb <0.0, 1.5, 0.0> ]
    }
  }
}
*/

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 0