// Sigmund Kyrre Ås <as@stud.ntnu.no> 2001

#version unofficial MegaPov 0.7;
#macro qparab(Y,X)
   (Y[0]+X*(Y[1]+(Y[2]+ Y[3]*X)*X))
#end

// Computes the drag-coefficient of a sphere as a function
// of the Reynolds nunber Re. Curve fitted after fig. A-56
// in "Fluid Mechanics & Hydraulics", Schaum's Solved Problems
#macro CDkule(Re) 
   %a2=array[4]{0.658, 34.87, -14.05, 4.22}
   %a3=array[4]{2.41, -17.08, 43.72, -30.41}
   %a4=array[4]{11.932, -8.472, 2.031, -0.1584}
   %a5=array[4]{14.93, -7.332, 1.1905, -0.06335}
   #if (Re<=0) %CD = 0; #end
   #if (Re> 8.0e6) %CD=0.2; #end
   %x0 = log(Re)/log(10.0);
   #if (Re> 0 & Re<=0.5) %CD=24/Re; #end
   #if (Re> 0.5 & Re<=100) %CD=qparab(a2,1/Re); #end
   #if (Re> 100 & Re<=1.0e4) %CD=qparab(a3,1/x0); #end
   #if (Re> 1.0e4 & Re<=3.35e5) %CD=qparab(a4,x0); #end
   #if (Re> 3.35e5 & Re <=5.0e5)
      %x1 = log(Re/4.5e5)/log(10);
      %CD=91.08*x1^4 + 0.0764;
   #end
   #if (Re> 5.0e5 & Re<=8.0e6) %CD=qparab(a5,x0); #end
   CD
#end

// equation of motion
#macro kule(Vm)
   %Vr=V-Vf; // realative fluid speed
   %vr=vlength(Vr);
   %Re= vr*d/nu;
   %CD= CDkule(Re);
   %f=C*vr*CD;
   //<-f*Vr.x, -g-f*Vr.y, -f*Vr.z> //dVdt
   (-f*Vr-g*y)
#end

#macro rk4kule(Xm,Vm,h)
   // Input:
   // X  position <,,>
   // V  speed  <,,>
   // h  step size
   // Output: Array containing F and F' at the end of the integration step.  
   
   // compute next position
   %k1 = h*Vm;
   %k2 = h*(Vm+k1/2);
   %k3 = h*(Vm+k2/2);
   %k4 = h*(Vm+k3);
   %Xnew = Xm+(k1+2*k2+2*k3+k4)/6;
   
   // compute next speed
   %k1 = h*kule(Vm);
   %k2 = h*kule(Vm+k1/2);
   %k3 = h*kule(Vm+k2/2);
   %k4 = h*kule(Vm+k3);
   %Vnew = Vm+(k1+2*k2+2*k3+k4)/6;
   #render "Speed:   " Skrivutvektor(Vnew)
   #render "Pos  :   " Skrivutvektor(Xnew)
   array[2] {Xnew,Vnew}
#end