// Computes the movement of a sphere influenced by gravity and wind.
// Sigmund Kyrre Ås <as@stud.ntnu.no> 2001

// 1 unit = 1 meter
// Note that bouyancy and hydrodynamic mass is left out
// This is a good estimate as long as rof << ros

#version unofficial MegaPov 0.7;
// (global variables)
$Time =2      ; // Total time [s]
$nu   =1.5e-5  ; // Kinematical viscosity [m^2/s]
$rof  =1.22    ; // Density of fluid [kg/m^3]
$ros  =1e3     ; // Density of sphere [kg/m^3]
$d    =.5    ; // Diameter of sphere [m]
$fps  =20      ; // number of frames per second
$V0   =<0,1,1>*10;  // Initial velocity [m/s]
$X0   =<0,0,0> ; // Initial position
$Vf   =30*<1,-.7,-2>;  // fluid speed

#macro Skrivutvektor(V)  #render concat("<",str(V.x,2,2),",",str(V.y,2,2),",",str(V.z,2,2),">\n")#end

// arrow showing wind direction
#if (vlength(Vf)!=0)
   union {
      $tn=vnormalize(Vf)*2;
      cone {0 .5 tn 0}
      cylinder {-tn 0 .2}
      pigment {rgb 1} finish {specular .3 roughness .2}
      translate <5,2,4>
   }
#end

$dt=1/fps; // Timestep [s]
$g = 9.81;
$ro= rof/ros;
$C = 0.75*ro/d;
$T=clock*Time;
#render concat("Clock    ",str(T,2,3),"\n")
#render concat("Stepsize ",str(dt,2,3),"\n")

disc{ 0,y,25 
   pigment {checker rgb .6, rgb .3 scale 5} 
   finish {diffuse .7 brilliance 1.1 specular .3 roughness .2} 
   no_shadow
}
light_source { 50 rgb .4} light_source { 100*y rgb 1}
camera {location <5,5,-5> look_at <3,1,4>}   

#include "num.inc" 
// use initial values on first run, else use previously computed value.
#ifndef (X) $X=X0; #else $X=F_a[0]; #end
#ifndef (V) $V=V0; #else $V=F_a[1]; #end
$F_a = rk4kule(X,V,dt)

#render "Pos now: " Skrivutvektor(X)
sphere {X, d/2 pigment {rgb <.5,0,.6>} finish {specular .5 roughness .1}}