Lorenz attractor divergence animation

#!/usr/bin/python3
#lorenz system demo
#larsupilami73

import os,sys
from scipy import *
from scipy import integrate
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

#parameters
sigma = 10.
rho = 28.
beta = 8/3.

tbegin = 0.
tend   = 200.
tstep = 0.01
T = linspace(tbegin,tend, int((tend-tbegin)/tstep))

#system
def dxdt(X, t=0):
    x = X[0]
    y = X[1]
    z = X[2]
    return array([sigma * (y - x),x * (rho - z ) - y, x*y - beta*z])


#initial run to find consistent initial values
X0=[0.5,0.5,0.5]
X, infodict = integrate.odeint(dxdt, X0, T, full_output=1)
#infodict['message'] # >>> 'Integration successful.' if argument full_output=True
x,y,z = X.T
#select closeby initial values to step over the attractor
numpoints = 20
initial_distance = 0.2
xp=[x[-1]+random.normal(scale=initial_distance) for dummy in range(numpoints)]
yp=[y[-1]+random.normal(scale=initial_distance) for dummy in range(numpoints)]
zp=[z[-1]+random.normal(scale=initial_distance) for dummy in range(numpoints)]

savetxt('out.dat',vstack((x,y,z,T)).T)
#sys.exit(0)

numframes = 200
try:
    os.mkdir('frames')
except:
    pass
try:
    os.system('rm -f ./frames/*.png')   #remove previous run frames
except:
    pass
fig = plt.figure(1, figsize=(8, 6)) #make the figure once
for k in range(numframes):
    for l in range(numpoints):
        T = [0,0.01,0.02,0.03,0.04] #range(10)* tstep
        X0=[xp[l],yp[l],zp[l]]
        X, infodict = integrate.odeint(dxdt, X0, T, full_output=1)
        #infodict['message'] # >>> 'Integration successful.' if argument full_output=True
        x1,y1,z1 = X.T

        #new initial values
        xp[l] = x1[-1]
        yp[l] = y1[-1]
        zp[l] = z1[-1]

    #plot it
    plt.cla()
    plt.clf()

    #plt.margins((0,0.02))
    #plt.rc('mathtext', fontset = 'cm') #make sure mathtext uses nice font.
    #plt.rc('mathtext', rm='serif')

    ax = Axes3D(fig, elev=-156, azim=110+k/2.)

    fig.set_facecolor('black')
    ax.set_facecolor('black')
    ax.grid(False)
    ax.w_xaxis.pane.fill = False
    ax.w_yaxis.pane.fill = False
    ax.w_zaxis.pane.fill = False

    #ax.set_title('rho=%f'%rho)
    no_points = len(x)
    cm = plt.get_cmap('binary_r')
    #this worked for previous matplotlib versions
    #then no c=c needed in ax.plot below
    #ax.set_color_cycle([cm(1.*i/(no_points-1)) for i in range(no_points-1)])
    for i in range(no_points-1):
        c = cm(int(255*i/(no_points-1)))
        bar = ax.plot(x[i:i+2],y[i:i+2],z[i:i+2], c=c,linewidth=0.6, alpha=0.5)

    ax.scatter(xp,yp,zp, c=range(numpoints), cmap=plt.cm.tab20, edgecolor='white', s=30.0,alpha=1.0,marker='o')

    ax.set_xlim(-26,26)
    ax.set_ylim(-20,20)
    ax.set_zlim(5,45)
    ax.set_xlabel('x', color='white')
    ax.set_ylabel('y',color='white')
    ax.set_zlabel('z',color='white')

    plt.savefig('./frames/frame%05d.png' %k)
    ax.clear()
    print('Frame number: %d' %k)


#now assemble all the pictures in a movie
#print('converting to gif!')
#os.system('convert -delay 40 -loop 0 ./frames/*.png fb_tot.gif')
#or to mp4
print('converting to mp4!')
os.system("ffmpeg -y -r 40 -i ./frames/frame%05d.png -c:v libx264 -vf fps=20 lorenz.mp4")