The first experiments are working well, although they are certainly a bit untidy looking – different projections and more consistent behaviours of pendulums are needed.
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>