/* declare some usefull constants
 * regarding dodcahedral symmetry
*/
#declare g = (1+sqrt(5))/2;
#declare V1 = (3*g+1)*0.02;
#declare V2 = (4*g+3)*0.02;

/* parameter controling how sharp the edges will be*/
#declare sharpness = 26;                           

/* parameter for declaring the type of polyhedra:

 Possible threshold values:
  a 0.5 threshold will create a regular dodcahedron {5, 3}
  a 1.5 threshold will create a small stellated dodecahedron {5/2, 5}
  a 2.5 will create a great dodecahedron {5, 5 / 2}
  a 3.5 will create a great stellated dodecahedron {5/2, 3}  

*/

#declare myThreshold = 1.5;

/************************************************************************************************************** 
 * Declare a sigmoid function (
 * f : R -> {y| 0 < y < 1}
 *
 *                         1
 * sigmoid(s) =  -----------------------
 *               1 + exp(-sharpness * s)
 *
 * A function that almost equals 0 in negative areas of s 
 * and alomst equals 1 in positive areas.
 * It slopes smoothly up from 0 to 1 around (s=0) 
 * (the higher 'sharpness' is, the steeper the slope is)
 *
 * Its values never reach zero (0) nor one (1), but are always between them.
 * Its values mainly distribute
 * around 1 (in positive s values) and 0 (in negetive s values)
 

  S:------0+++++++
  +
  1        ,~~~~~~
  +       |         
  + _____,
  0
  -
        
*/ 
#declare sigmoid = function(s) { 1 / (1 + exp(-sharpness * s))}                                         
/**************************************************************************/


// number of planes
#declare n = 12;

/* declare 12 plane functions (one for each face of the dodcahedron)

 (The dodcahedron is the intersection of these plane or -inside- all of the planes)

 each function is a 3d function whose values are
 negative "iside" the plane,
 positive "outside" the plane
 and zero on its surface
*/
 
#declare planes = array[n]{
        function{(-V1*x-V2*z - 1)},
        function{(-V1*y-V2*x - 1)},
        function{(-V1*z-V2*y - 1)},

        function{(-V1*x+V2*z - 1)},
        function{(-V1*y+V2*x - 1)},
        function{(-V1*z+V2*y - 1)},

        function{(+V1*x-V2*z - 1)},
        function{(+V1*y-V2*x - 1)},
        function{(+V1*z-V2*y - 1)},

        function{(+V1*x+V2*z - 1)},
        function{(+V1*y+V2*x - 1)},
        function{(+V1*z+V2*y - 1)}
}                         


/************************************************************************************************************** 

 dodca : R3 -> {y| 0 < y < 12}
 dodca(x,y,z) = sigma{i: 0 -> 11}(sigmoid(planes[i]))
 
 
 create a function accumulating sigmoids of values of all functions.
 accumulating sigmoids creates a 3d function in space matching each point
 with the number of planes the point is "outside" of.
 point 0, for instance is inside all planes so it will recieve 
 a value almost equal to zero.                
 (notice: it does not equal "n" but is a bit smaller, since the sigmoid never really equals 1)
 
 function's values spread between 0 and 12, are never whole integers but are 
 distributed mainly around integers.
 */
#declare dodca = function {
                #declare j = 0;                                             
                #while (j < n)                                                                                
                              sigmoid(planes[j](x, y, z)) +
                              #declare j = j + 1;
                #end                             
                0    
}                                                                                                              

/************************************************************************************************************** 
 create an isosurface using the "dodca" function
 since the dodca function count the number of planes a
 given point is "outside".
 setting the threshold on 1.5
 divides 3d space between all points that are -outside- 1 or less than one plane
 and all points that are -outside- two ore more planes       
 
 in other words,
 
 the space is divided into two regions, one containg 
 points that are either -inside- all but one plane or are -inside- all plane,
 the other region containing all the rest (points that are inside less planes).
 
 By that we are creating a small stellated dodecahedron {5/2, 5}
   
*/

isosurface {
        function {
                 dodca(x,y,z)
        }           
        threshold myThreshold
        max_gradient 2.4
        contained_by {sphere {0, 10}}

        pigment{color rgb<1,1,1>}
        no_shadow
        
}
        
                                                       
/**************************************************************************************************************                                                       
 lights, camera, action !!
*/
light_source {<17,19,12> color <1,1,1>}
camera {
        location <14,21,16>
        look_at 0
        }