In [1]:

    

# Double pendulum formula translated from the C code at
# http://www.physics.usyd.edu.au/~wheat/dpend_html/solve_dpend.c

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation


    

In [2]:

    

import pickle
arquivo = open("../xyz", "r")

X = pickle.loads(arquivo.read())
print X.shape

joelho = X[0,:,:].T
trocanter = X[1,:,:].T
tibia = X[2,:,:].T

x_min = np.min(np.array([joelho[:, 0], trocanter[:, 0], tibia[:, 0]]))
x_max = np.max(np.array([joelho[:, 0], trocanter[:, 0], tibia[:, 0]]))

y_min = np.min(np.array([joelho[:, 1], trocanter[:, 1], tibia[:, 1]]))
y_max = np.max(np.array([joelho[:, 1], trocanter[:, 1], tibia[:, 1]]))

z_min = np.min(np.array([joelho[:, 2], trocanter[:, 2], tibia[:, 2]]))
z_max = np.max(np.array([joelho[:, 2], trocanter[:, 2], tibia[:, 2]]))

print joelho.shape


    

    




(3, 3, 66)
(66, 3)

In [3]:

    

fig = plt.figure()
ax = fig.add_subplot(111, autoscale_on=False, xlim=(x_min, x_max), ylim=(z_min, z_max))
ax.grid()

line, = ax.plot([], [], 'o-', lw=2)
time_template = 'time = %.1fs'
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)


    

In [4]:

    

def init():
    line.set_data([], [])
    time_text.set_text('')
    return line, time_text

def animate(i):
    thisx = [trocanter[i, 0], joelho[i, 0], tibia[i, 0]]
    thisy = [trocanter[i, 2], joelho[i, 1], tibia[i, 2]]

    line.set_data(thisx, thisy)
    
    return line, time_text

ani = animation.FuncAnimation(fig, animate, np.arange(0, 65),
    interval=25, blit=False, init_func=init)

#ani.save('double_pendulum.mp4', fps=15)
plt.show()


    

In [ ]: