from matplotlib import pylab
from mpl_toolkits.mplot3d import axes3d

# time increment
dt = 0.001

a = 0.2
b = 0.2
c = 5.7

def f_dx(x, y, z):
    global a, b, c
    return -y - z

def f_dy(x, y, z):
    global a, b, c
    return x + a*y

def f_dz(x,y,z):
    global a, b, c
    return b + z*(x - c)

	
#	foutrh order Runge-Kutta routine	
def rk4(x, y, z):
    global dt
    k1x = dt * f_dx(x, y, z)
    k1y = dt * f_dy(x, y, z)
    k1z = dt * f_dz(x, y, z)

    k2x = dt * f_dx(x + k1x/2.0, y + k1y/2.0, z + k1z/2.0)
    k2y = dt * f_dy(x + k1x/2.0, y + k1y/2.0, z + k1z/2.0)
    k2z = dt * f_dz(x + k1x/2.0, y + k1y/2.0, z + k1z/2.0)

    k3x = dt * f_dx(x + k2x/2.0, y + k2y/2.0, z + k2z/2.0)
    k3y = dt * f_dy(x + k2x/2.0, y + k2y/2.0, z + k2z/2.0)
    k3z = dt * f_dz(x + k2x/2.0, y + k2y/2.0, z + k2z/2.0)

    k4x = dt * f_dx(x + k3x, y + k3y, z + k3z)
    k4y = dt * f_dy(x + k3x, y + k3y, z + k3z)
    k4z = dt * f_dz(x + k3x, y + k3y, z + k3z)

    x = x + k1x/6.0 + k2x/3.0 + k3x/3.0 +k4x/6.0
    y = y + k1y/6.0 + k2y/3.0 + k3y/3.0 +k4y/6.0
    z = z + k1z/6.0 + k2z/3.0 + k3z/3.0 +k4z/6.0

    return [x, y, z]
	
# list of all points to be plotted
all_x = []
all_y = []
all_z = []

i = 0

while i < 10:
	# set initial values
    x = -9 + (i - 5) / 2.0
    y = 0 + (i - 5) / 2.0
    z = 0 + (i - 5) / 2.0
	
	# list of points for one orbit
    list_x = [x]
    list_y = [y]
    list_z = [z]

    t = 1
    while t < 20000:
        [x, y, z] = rk4(list_x[t-1], list_y[t-1], list_z[t-1])
        list_x.append(x)
        list_y.append(y)
        list_z.append(z)
        t = t + 1
    all_x = all_x + list_x
    all_y = all_y + list_y
    all_z = all_z + list_z
    i = i + 1

# draw graph using [all_x, all_y, all_z]
fig = pylab.figure()
ax = axes3d.Axes3D(fig)
ax.scatter3D(all_x, all_y, all_z, s = 0.1)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
pylab.show()