// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10
// Copyright 2001-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 ======= 10

#version 3.5;

#include "functions.inc" // For f_sphere(), f_torus and f_noise3d()
#include "rad_def.inc"

global_settings { radiosity { Rad_Settings(Radiosity_Normal, off, off) } }

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

#declare SquashFn = function(T, S) { select(S, -1, 1)*(1/(1 + exp(-abs(S)*T)) - 1) }


#macro TransformFunction(Fn, Transform)

  #local TT = function { transform { Transform inverse } }

  function { Fn(TT(x, y, z).x, TT(x, y, z).y, TT(x, y, z).z) }

#end // macro TransformFunction


#declare HyperTorusX_Fn =
  function(x, y, z, Rmaj, Rmin, A, B) {
    pow(pow(x, B) + pow(pow(pow(y, A) + pow(z, A), 1/A) - Rmaj, B), 1/B) - Rmin
  }

#declare NoiseFn = TransformFunction(function { 0.5 - f_noise3d(x, y, z) }, transform { scale 1.5 })

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

#declare StrgthI = 2;
#declare StrgthO = 6;

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

#declare RmajA = 11.9;
#declare RminA = 1.7;
#declare RadA = 3.2;

#declare pCtrSphereA1 = -(RmajA - 1)*y;
#declare pCtrSphereA2 =  (RmajA - 1)*y;

#declare vScaleA = <2.5, 1.0, 2.5>;


#declare TransA = transform { scale vScaleA }

#declare TorusFnA0 = function { HyperTorusX_Fn(x, y, z, RmajA, RminA, StrgthI, StrgthO) }
#declare TorusFnA1 = TransformFunction(TorusFnA0, transform { rotate  60*y })
#declare TorusFnA2 = TransformFunction(TorusFnA0, transform { rotate 120*y })

#declare SphereFnA = function { f_sphere(x, y, z, RadA) } 
#declare SphereFnA1 = TransformFunction(SphereFnA, transform { TransA translate pCtrSphereA1 })
#declare SphereFnA2 = TransformFunction(SphereFnA, transform { TransA translate pCtrSphereA2 })

#declare RmaxA = RmajA + RminA;
#declare pMinA = <-RmaxA, (pCtrSphereA1 - RadA*vScaleA).y, -RmaxA>;
#declare pMaxA = < RmaxA, (pCtrSphereA2 + RadA*vScaleA).y,  RmaxA>;

#declare IsoA =
  isosurface {
    function {
      SquashFn(TorusFnA0(x, y, z), 8)
     +SquashFn(TorusFnA1(x, y, z), 8)
     +SquashFn(TorusFnA2(x, y, z), 8)
     +SquashFn(SphereFnA1(x, y, z), 2)
     +SquashFn(SphereFnA2(x, y, z), 2)
     +0.1*NoiseFn(x, y, z)
    }
    threshold -0.5
    max_gradient 2
    contained_by { box { pMinA, pMaxA } }
    transform { TransA inverse }
    rotate -15*y
  }  

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

#declare RmajRing = 10;
#declare RminRing = 0.6;

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

#declare RmajB = RmajRing;
#declare RminB = RminA;
#declare RadB = 3.7;

#declare pCtrSphereB1 = -RmajB*y;
#declare pCtrSphereB2 =  RmajB*y;

#declare vScaleB = <1.0, 2.0, 1.0>;

#declare TransB = transform { scale vScaleB }

#declare TorusFnB = function { HyperTorusX_Fn(x, y, z, RmajB, RminB, StrgthI, StrgthO) }
#declare TorusFnB0 = TransformFunction(TorusFnB, transform { scale <0.6, 1, 1> rotate -60*y })
#declare TorusFnB1 = TransformFunction(TorusFnB, transform { scale <0.6, 1, 1> })
#declare TorusFnB2 = TransformFunction(TorusFnB, transform { scale <0.6, 1, 1> rotate  60*y })

#declare SphereFnB = function { f_sphere(x, y, z, RadB) } 
#declare SphereFnB1 = TransformFunction(SphereFnB, transform { TransB translate pCtrSphereB1 })
#declare SphereFnB2 = TransformFunction(SphereFnB, transform { TransB translate pCtrSphereB2 })

#declare RingFn = TransformFunction(function { f_torus(x, y, z, RmajRing, RminRing + 0.1) }, TransB)

#declare RmaxB = RmajB + RminB;
#declare pMinB = <-RmaxB, (pCtrSphereB1 - RadB*vScaleB).y, -RmaxB>;
#declare pMaxB = < RmaxB, (pCtrSphereB2 + RadB*vScaleB).y,  RmaxB>;

#declare IsoB =
  isosurface {
    function {
      SquashFn(TorusFnB0(x, y, z), 5)
     +SquashFn(TorusFnB1(x, y, z), 5)
     +SquashFn(TorusFnB2(x, y, z), 5)
     +SquashFn(SphereFnB1(x, y, z), 8)
     +SquashFn(SphereFnB2(x, y, z), 8)
     +SquashFn(RingFn(x, y, z), -3)
     +0.1*NoiseFn(x, y, z)
    }
    threshold -0.5
    max_gradient 2
    contained_by { box { pMinB, pMaxB } }
    transform { TransB inverse }
    rotate  15*y
  }                                     

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

#declare RmajC = RmajRing;
#declare RminC = RminRing;

#declare RmaxC = RmajC + RminC;
#declare pMinC = -<RmaxC, RminC, RmaxC>;
#declare pMaxC =  <RmaxC, RminC, RmaxC>;

#declare IsoC =
  isosurface {
    function { f_torus(x, y, z, RmajC, RminC) + 0.04*NoiseFn(x, y, z) }
    threshold 0
    max_gradient 2
    contained_by { box { pMinC, pMaxC } }
  }                                     

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

union {
  object { IsoA scale 2.5 }
  object { IsoB scale 2.0 }
  object { IsoC scale 2.0 }
  texture {
    pigment { color rgb <1, 1, 1> }
    finish {
      ambient color rgb <0, 0, 0>
      diffuse 1
    }
  }
}

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

sky_sphere {
  pigment {
    gradient y
    color_map {
      [ 0.0 color rgb <0.8, 0.8, 1.6> ]
      [ 0.5 color rgb <0.0, 0.0, 0.0> ]
      [ 1.0 color rgb <0.8, 0.8, 1.6> ]
    }
    scale <1, 2, 1>
    translate -y
  }
}

light_source {
  <-3, -5, -2>*20
  color rgb <0.95, 0.78, 0.42>
  area_light
  10*x, 10*z,
  4, 4
  jitter
  orient
}

camera {
  location vnormalize(< 3, 5, -14>)*90
  look_at <0, 4, 0>
}

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