/*-------simple particle system tutorial

by Jonathan Rafael Ghiglia 
12th January 2002 

Desc: a basic particle flow. Only still. 
This is not meant to be the best particle flow algorithm. Nevertheless it's the simplest to understand. 

You can freely use this code. 
Report bugs to: jrg_pov@hotmail.com

--------------------------------------*/ 


#declare Version = 31; // -31 MegaPov; -35 POV-Ray 3.5 (It crashes with the present beta!)

#if (Version=31)
#version unofficial MegaPov 0.7;
#else #include "transforms.inc"
#end

#macro Blob (Threshold,Radius,Alpha)
#local Strength_Temp = (Alpha+1)*Threshold;
#local Radius_Temp = Radius/ sqrt ( 1 - sqrt (1/(Alpha+1)));
Radius_Temp
strength Strength_Temp
#end  

#macro round_cube (_x,_y,_z,_r)
merge {
        box {-1,1 scale <_x, _y-_r, _z-_r>}
        box {-1,1 scale <_x-_r, _y-_r, _z>}
        box {-1,1 scale <_x-_r, _y, _z-_r>}
        cylinder {<0,-_y+_r,0>, <0,_y-_r,0>, _r translate <-_x+_r,0,-_z+_r>} 
        cylinder {<0,-_y+_r,0>, <0,_y-_r,0>, _r translate <_x-_r,0,-_z+_r>}
        cylinder {<0,-_y+_r,0>, <0,_y-_r,0>, _r translate <-_x+_r,0,_z-_r>}
        cylinder {<0,-_y+_r,0>, <0,_y-_r,0>, _r translate <_x-_r,0,_z-_r>}
        cylinder {<-_x+_r,0,0>, <_x-_r,0,0>, _r translate <0,_y-_r,-_z+_r>}
        cylinder {<-_x+_r,0,0>, <_x-_r,0,0>, _r translate <0,-_y+_r,-_z+_r>}
        cylinder {<-_x+_r,0,0>, <_x-_r,0,0>, _r translate <0,_y-_r,_z-_r>}
        cylinder {<-_x+_r,0,0>, <_x-_r,0,0>, _r translate <0,-_y+_r,_z-_r>}
        cylinder {<0,0,-_z+_r>, <0,0,_z-_r>, _r translate <-_x+_r,_y-_r,0>}
        cylinder {<0,0,-_z+_r>, <0,0,_z-_r>, _r translate <_x-_r,_y-_r,0>}
        cylinder {<0,0,-_z+_r>, <0,0,_z-_r>, _r translate <-_x+_r,-_y+_r,0>}
        cylinder {<0,0,-_z+_r>, <0,0,_z-_r>, _r translate <_x-_r,-_y+_r,0>}
        sphere {0,1 scale _r translate <-_x+_r,_y-_r,_z-_r>}
        sphere {0,1 scale _r translate <_x-_r,_y-_r,_z-_r>}
        sphere {0,1 scale _r translate <-_x+_r,-_y+_r,_z-_r>}
        sphere {0,1 scale _r translate <_x-_r,-_y+_r,_z-_r>}
        sphere {0,1 scale _r translate <_x-_r,_y-_r,-_z+_r>}
        sphere {0,1 scale _r translate <-_x+_r,_y-_r,-_z+_r>}
        sphere {0,1 scale _r translate <_x-_r,-_y+_r,-_z+_r>}
        sphere {0,1 scale _r translate <-_x+_r,-_y+_r,-_z+_r>}
        }
#end        

#declare Environment=
union {
        plane {y,0}
union {        
object {round_cube (3,7,1,0.2)
        translate -10*z}
        cylinder {0,7*y,1.3 translate 10*z}
pigment {rgb x*.2}
finish {ambient 0 diffuse .8 specular .5 roughness .01 reflection .2}
}
}


#declare Time = 5;  
#declare Iterations = 100;
#declare Particle_Count = 4000; // must be > Iterations
#declare Elasticity = 0.35; 
#declare Acceleration = <0,-9.81,0>; 
#declare Particle_Radius = 0.2;   
#declare RS = seed (0);
#declare Random_Pos = 0.05; // max random amount added to position
#declare Friction = 0.98;   // 1= no friction
#declare dt = Time/ Iterations;   
#declare File_Option = 0; // -0 do not save; -1 save file; -2 load file
#declare File_Name = "par_tut_arr.inc"

#declare Particles_Array = array [Particle_Count][2]

#if (File_Option !=2)


#macro Start_Point ()  // this macro defines the starting position. Must return a vector.
vtransform(vtransform(0.8*x*rand(RS),transform {rotate 360*rand(RS)*z}),transform {translate <0,10,0.8>})
#end 

#macro Start_Vel ()  // this macro defines the starting velocity. Must return a vector.
vtransform(vtransform ((-15+5*rand(RS))*z,transform {rotate 10*rand(RS)*x}),transform {rotate 360*rand(RS)*z})
#end  

#macro Update_Particle (Particle)
#declare Particles_Array [Particle][0] = Particles_Array [Particle][0] + Particles_Array [Particle][1]*dt +0.5*Acceleration*dt*dt +<-Random_Pos+2*Random_Pos*rand(RS),-Random_Pos+2*Random_Pos*rand(RS),-Random_Pos+2*Random_Pos*rand(RS)>;
#declare Particles_Array [Particle][1] = Particles_Array [Particle][1] + Acceleration*dt;
#end

#macro Rebound (Particle,Normal)
#declare Particles_Array [Particle][1] = (Particles_Array [Particle][1] - vdot (Particles_Array [Particle][1],vnormalize (Normal))*vnormalize (Normal)*(1+Elasticity/Friction))*Friction;
#end  

//-----main loop

#declare Present_Iteration = 1;
#while (Present_Iteration <= Iterations)

#declare Present_Particle = 0;
#while (Present_Particle < floor (Particle_Count/Iterations)*Present_Iteration)

//---- create new particle

#ifndef (Particles_Array [Present_Particle][0])
#declare Particles_Array [Present_Particle][0] = Start_Point ();
#declare Particles_Array [Present_Particle][1] = Start_Vel ();

#else

#declare Old_Position = Particles_Array [Present_Particle][0];
Update_Particle (Present_Particle)

//---- check for collisions

#declare Normal = <0,0,0>;
#declare Vel_Dir = Particles_Array [Present_Particle][0]-Old_Position; 
#if (vlength(trace (Environment,Old_Position,Vel_Dir,Normal)-Old_Position)-vlength (Vel_Dir)<Particle_Radius)
#declare Particles_Array [Present_Particle][0] = Particles_Array [Present_Particle][0]-(Particle_Radius-(vlength(trace (Environment,Old_Position,Vel_Dir,Normal)-Old_Position)-vlength (Vel_Dir)))*vnormalize (Vel_Dir);
#if (vdot (Vel_Dir,Normal)<0)
Rebound (Present_Particle, Normal)
#end
#end

#end

#declare Present_Particle = Present_Particle +1 ;
#end   

#debug concat ("Iterations: ", str (Present_Iteration,0,0), ". Going to ", str (Iterations,0,0), ".\n")
#declare Present_Iteration = Present_Iteration +1;
#end   

#if (File_Option = 1 )
#fopen arr_file File_Name write  
#write (arr_file,"//------Particles Data------\n") 
#local i=0;
#while (i<Particle_Count)
#ifdef (Particles_Array[i][0])
#write (arr_file, "#declare Particles_Array[",i,"][0] = ",Particles_Array[i][0],";\n") 
#end
#local i=i+1;
#end
#fclose arr_file
#end

#else

#include File_Name
#end
 


  

blob {
        threshold 0.5
        #declare j=0;
        #while (j<Particle_Count)
        
        sphere {0,Blob (0.5,Particle_Radius,.2)
                        translate Particles_Array[j][0]   
                        }

        #declare j=j+1;
        #end 
                        pigment {rgb 1}
                        finish {diffuse 0.75 specular 0.7 roughness 0.008 
                            } 
} 

global_settings {
    assumed_gamma 1
    ambient_light 0
}

//--- set

union { 
object {Environment
    pigment {rgb 0}
    finish {ambient 0
        diffuse 0.8  specular 0.5 roughness 0.01 reflection 0.2
    } 
} 
      union {  
         difference {
                superellipsoid {<1,0.25> scale <1.5,1.5,0.8>}
                cylinder {-2*z,2*z,1}
                translate 0.8*z
                translate 10*y} 
        difference {
                cylinder {0.5*z,7*z,1.3}
                cylinder {-z,8*z,1}
                translate 10*y}
        intersection {
                torus {3,1.3 rotate 90*z}
                plane {-y,0}
                plane {-z,0}
                translate <0,10-3,7>} 
pigment {rgb .2*x}
finish {
        ambient 0
        diffuse .8
        specular .5
        roughness .01
        reflection .2
}
      }
}





//----lights

light_source {
    <0,140,0>
    rgb 1
    area_light 30*x,30*z,5,5 adaptive 0 jitter
    fade_power 2
    fade_distance 100
}

light_source {
    <30,80,-99>
    rgb 1 
    area_light 10*x,10*y,5,5 adaptive 0 jitter circular orient
    fade_power 2
    fade_distance 120
}

light_source {
    <0,50,99>
    rgb 0.8
    area_light 10*x,10*y,5,5 adaptive 0 jitter circular orient
    fade_power 2
    fade_distance 100
} 

light_source {
        <0,25,-40>
        rgb 0.35
        rotate -100*y
        shadowless}

camera {
    location <0,25,-40>
    right image_width/image_height*x
    look_at 5*y 
    angle 40
    rotate -100*y
}