Simulate an arbitrary number of linked oscillators (harmonograph), drawing output to .png file.

Harmonograph.py

import numpy as np
import math
import matplotlib.pyplot as plt
import webbrowser

def makePlot(xvals,yvals,fname="harmonograph.png",fopen=False,w=50,h=50,width=0.1):
    maxs=[np.max(np.abs(xvals)),np.max(np.abs(yvals))]
    farthest=max(maxs)

    drw=plt.figure(1)
    drw.set_figheight(50)
    drw.set_figwidth(50)
    rect=drw.patch
    rect.set_facecolor('white')
    ax=drw.add_subplot(1,1,1)
    ax.plot(xvals/farthest,yvals/farthest,linewidth=width,color='black')
    ax.axis('off')
    drw.savefig(fname,facecolor=drw.get_facecolor(),edgecolor="none")
    if fopen:
        webbrowser.open(fname)

def osc(A=1,f=1,p=0,d=0):
    def dynamics(t):
        return(A*math.sin(t*f+p)*math.exp(-1*d*t))
    return(dynamics)

def makeDosc(pdicts):
    dynList=[osc(pd["A"],pd["f"],pd["p"],pd["d"]) for pd in pdicts]
    def dosc(t):
        return(sum(dyn(t) for dyn in dynList))
    return(dosc)

def randomParams(Amin=0.5,Amax=1.5,fmin=0.5,fmax=1.5,pmin=0,pmax=math.pi,dmin=0,dmax=0.01):
    pdict=dict({"A":random.uniform(Amin,Amax),"f":random.uniform(fmin,fmax),"p":random.uniform(pmin,pmax),"d":random.uniform(dmin,dmax)})
    return(pdict)

xdicts=[
    dict({"A":1,"f":2,"p":1,"d":0.005}),
    dict({"A":0.25,"f":4,"p":1,"d":0.004}),
    dict({"A":0.5,"f":2,"p":1,"d":0.003}),
    dict({"A":1,"f":4,"p":1,"d":0.002}),
    dict({"A":0.6,"f":2,"p":1,"d":0})
    ]
ydicts=[
    dict({"A":1,"f":1,"p":1,"d":0.005}),
    dict({"A":0.25,"f":2,"p":1,"d":0.004}),
    dict({"A":0.5,"f":1,"p":1,"d":0.003}),
    dict({"A":1,"f":2,"p":1,"d":0.002}),
    dict({"A":0.25,"f":1,"p":1,"d":0})
    ]

x=makeDosc(xdicts)
y=makeDosc(ydicts)

tmax=10000             
tvals=np.arange(0,tmax,0.005)

xvals=[x(t) for t in tvals]
yvals=[y(t) for t in tvals]

makePlot(xvals,yvals,fname="harmonograph.png",fopen=True,w=2000,h=2000,width=0.3)