# Integrator.py: This module contains the 3D 
#                fourth-order Runge-Kutta method.
#
# For a 3D ODE
#     dx/dt = f(x,y,z)
#     dy/dt = g(x,y,z)
#     dz/dt = h(x,y,z)

# 3D Fourth-Order Runge-Kutta Integrator
def RK3DIntegrator(alpha,wx,wy,wz,x,y,z,f,g,h,dt):
   k1x = dt * f(alpha,wx,wy,wz,x,y,z)
   k1y = dt * g(alpha,wx,wy,wz,x,y,z)
   k1z = dt * h(alpha,wx,wy,wz,x,y,z)
   k2x = dt * f(alpha,wx,wy,wz,x + k1x / 2.0,y + k1y / 2.0,z + k1z / 2.0)
   k2y = dt * g(alpha,wx,wy,wz,x + k1x / 2.0,y + k1y / 2.0,z + k1z / 2.0)
   k2z = dt * h(alpha,wx,wy,wz,x + k1x / 2.0,y + k1y / 2.0,z + k1z / 2.0)
   k3x = dt * f(alpha,wx,wy,wz,x + k2x / 2.0,y + k2y / 2.0,z + k2z / 2.0)
   k3y = dt * g(alpha,wx,wy,wz,x + k2x / 2.0,y + k2y / 2.0,z + k2z / 2.0)
   k3z = dt * h(alpha,wx,wy,wz,x + k2x / 2.0,y + k2y / 2.0,z + k2z / 2.0)
   k4x = dt * f(alpha,wx,wy,wz,x + k3x,y + k3y,z + k3z)
   k4y = dt * g(alpha,wx,wy,wz,x + k3x,y + k3y,z + k3z)
   k4z = dt * h(alpha,wx,wy,wz,x + k3x,y + k3y,z + k3z)
   x += ( k1x + 2.0 * k2x + 2.0 * k3x + k4x ) / 6.0
   y += ( k1y + 2.0 * k2y + 2.0 * k3y + k4y ) / 6.0
   z += ( k1z + 2.0 * k2z + 2.0 * k3z + k4z ) / 6.0
   return x,y,z