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>