First experiments with virtual Harmonograph

The first experiments are working well, although they are certainly a bit untidy looking – different projections and more consistent behaviours of pendulums are needed.

Individual harmonograph

Image with individual pendulum routes and the harmonograph beneath.

Script taking the csv files of several pendulums into one harmonograph plot, using theory from previous post:

 import numpy
 import matplotlib
 import pandas
 import matplotlib.pyplot as plt
 from mpl_toolkits.mplot3d import Axes3D</pre>
course1 = pandas.read_csv('C:\Users\david\Desktop\spherical_pendulum_0.98m.csv')
course2 = pandas.read_csv('C:\Users\david\Desktop\spherical_pendulum_0.99m.csv')
course3 = pandas.read_csv('C:\Users\david\Desktop\spherical_pendulum_1.00m.csv')

#check each of course1 course2 and course3 have the same number of timesteps
if (course1.shape[0]==course2.shape[0]) and (course2.shape[0]==course3.shape[0]):
print("all the same shape, good news!")
if (course1.shape[1]==course2.shape[1]) and (course2.shape[1]==course3.shape[1]):
print("same number of columns, good news!")
if list(course1.columns.values) == list(course2.columns.values) and list(course2.columns.values) == list(course3.columns.values):
print("same column titles, good news!")

numberOfTimeSteps = course1.shape[0]
timeSteps = numpy.arange(numberOfTimeSteps)
columnTitles = ['t', 'x1', 'y1', 'z1', 'x2', 'y2', 'z2', 'x3', 'y3', 'z3', 'x', 'y', 'z']
constants = {'lambda1':0.8, 'lambda2':0.8, 'lambda3':0.8, 'r1':0.8, 'r2':0.8, 'r3':0.8, 'xsupport1':0.0,'ysupport1':0.0,'zsupport1':0.0,'xsupport2':1.0,'ysupport2':0.0,'zsupport2':0.0,'xsupport3':0.5,'ysupport3':0.86,'zsupport3':0.0,}
harmonographCourse = pandas.DataFrame(numpy.zeros([numberOfTimeSteps,13],dtype=long), index = timeSteps, columns = columnTitles)

harmonographCourse['t'] = course1['t']
harmonographCourse['x1'] = constants['lambda1']*course1['x']+constants['xsupport1']
harmonographCourse['y1'] = constants['lambda1']*course1['y']+constants['ysupport1']
harmonographCourse['z1'] = constants['lambda1']*course1['z']+constants['zsupport1']
harmonographCourse['x2'] = constants['lambda2']*course2['x']+constants['xsupport2']
harmonographCourse['y2'] = constants['lambda2']*course2['y']+constants['ysupport2']
harmonographCourse['z2'] = constants['lambda2']*course2['z']+constants['zsupport2']
harmonographCourse['x3'] = constants['lambda3']*course3['x']+constants['xsupport3']
harmonographCourse['y3'] = constants['lambda3']*course3['y']+constants['ysupport3']
harmonographCourse['z3'] = constants['lambda3']*course3['z']+constants['zsupport3']

xsupport1 = numpy.array([constants['xsupport1'],constants['ysupport1'],constants['zsupport1']])

for timeStep in timeSteps:
x1 = numpy.array([harmonographCourse.loc[timeStep,'x1'],harmonographCourse.loc[timeStep,'y1'],harmonographCourse.loc[timeStep,'z1']])
x2 = numpy.array([harmonographCourse.loc[timeStep,'x2'],harmonographCourse.loc[timeStep,'y2'],harmonographCourse.loc[timeStep,'z2']])
x3 = numpy.array([harmonographCourse.loc[timeStep,'x3'],harmonographCourse.loc[timeStep,'y3'],harmonographCourse.loc[timeStep,'z3']])
vector_a = x2 - x1
vector_c = numpy.cross((x2-x1),(x3-x1))
vector_b = numpy.cross((x2-x1),vector_c)
norm_vector_a = vector_a / numpy.sqrt(numpy.sum(numpy.square(vector_a)))
norm_vector_b = vector_b / numpy.sqrt(numpy.sum(numpy.square(vector_b)))
norm_vector_c = vector_c / numpy.sqrt(numpy.sum(numpy.square(vector_c)))
d = numpy.dot((x2-x1),norm_vector_a)
i = numpy.dot((x3-x1),norm_vector_a)
j = numpy.dot((x3-x1),norm_vector_b)
a = (numpy.square(constants['r1']) - numpy.square(constants['r2']) + numpy.square(d)) / (2*d)
b = (numpy.square(constants['r1']) - numpy.square(constants['r3']) + numpy.square(i) + numpy.square(j))/(2*j) - i/j*a
c = -1*numpy.sqrt(numpy.square(constants['r1'])-numpy.square(a)-numpy.square(b))
x = xsupport1 + x1 + a*norm_vector_a + b*norm_vector_b + c*norm_vector_c
harmonographCourse.loc[timeStep,'x'] = x[0]
harmonographCourse.loc[timeStep,'y'] = x[1]
harmonographCourse.loc[timeStep,'z'] = x[2]
print harmonographCourse

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(harmonographCourse['x'], harmonographCourse['y'], harmonographCourse['z'], label = 'harmonograph')
ax.plot(harmonographCourse['x1'], harmonographCourse['y1'], harmonographCourse['z1'], label = '0.98m pendulum')
ax.plot(harmonographCourse['x2'], harmonographCourse['y2'], harmonographCourse['z2'], label = '0.99m pendulum')
ax.plot(harmonographCourse['x3'], harmonographCourse['y3'], harmonographCourse['z3'], label = '1.00m pendulum')
ax.set_xlabel('x value')
ax.set_ylabel('y value')
ax.set_zlabel('z value')
plt.axis('equal')
plt.show()
<pre>

Leave a Reply