Involute of cubic curve

#!/usr/bin/env python

import sys
import os
from math import *

from pyx import canvas, path, deco, trafo, style, text, color, deformer

rgb = color.rgb
rgbfromhexstring = color.rgbfromhexstring

red, green, blue, yellow = (
    rgbfromhexstring("#d00000"),
    rgbfromhexstring("#006000"),
    rgb.blue, rgb(0.75, 0.75, 0))


def axis(c, x0, y0, x1, y1, extra):
    w = x1-x0
    mx = 0.01*w

    h = y1-y0
    my = 0.01*h

    tick = 1.
    while w/tick<10.:
        tick *= 0.1
    x = floor(x0)-1.
    assert x<x0
    assert x<x1
    while x < x1:
        x += tick
        if x<x0:
            continue
        count = int(round(x/tick))
        m = 0.5*my if (count%5) else my
        c.stroke(path.line(x, -m, x, +m), [style.linewidth.thin]+extra)

    c.stroke(path.line(x0-tick, 0, x1+2*tick, 0),
        extra+[deco.earrow(size=0.4)]) #[style.linewidth.thin]+extra)


def plot(c, f, df, x0, x1):

    n = 100
    dx = 1.*(x1-x0)/n
    x = x0
    y = y0 = y1 = f(x)

    ps = [path.moveto(x, y)]
    for i in range(n):

        x2 = x+dx
        y2 = f(x2)

        y0 = min(y0, y2)
        y1 = max(y1, y2)

        r = dx
        ps.append(path.lineto(x2, y2))

        x, y = x2, y2

    c.stroke(path.path(*ps), [red]+[style.linewidth.THick]+extra)

    axis(c, x0, y0, x1, y1, extra)



def main():

    os.chdir("frames")

    def f(x):
        return x**3

    def df(x):
        return 3*x**2

    global extra
    sx, sy = 5., 5.
    extra = [trafo.scale(sx=sx, sy=sy)]

    pts = []

    delta = 1.7 # path length

    N = 200
    dx = 1.0 / N
    x0 = -0.7
    i = 0

    r = 1.2
    clip = path.rect(-r*sx, -sy, 2*r*sx, 2*sy)

    while 1:

        c = canvas.canvas([canvas.clip(clip)])

        #x0 = -0.9 + i*dx
        y0 = f(x0)
        #print "x0=%.3f y0=%.3f" % (x0, y0)
        s = (1 + df(x0)**2)**0.5
        dy = f(x0+dx)-y0
        #print dx, dy, delta,
        delta -= (dx**2 + dy**2)**0.5
        print i, "%.4f"%delta

        x1 = x0 + delta/s

        y1 = y0 + (x1-x0)*df(x0)
        c.stroke(path.line(x0, y0, x1, y1), [style.linewidth.Thick]+extra)

        plot(c, f, df, -1., 1.)

        r = 0.02
        c.fill(path.circle(x1, y1, r), extra)

        pts.append((x1, y1))

        if x1 > 1.:
            break

        x0 += dx

        ps = [path.moveto(pts[0][0], pts[0][1])] + [path.lineto(p[0], p[1]) for p in pts[1:]]
        if len(pts)>1:
            c.stroke(path.path(*ps), [blue, style.linewidth.Thick]+extra)

        name = "pic-%.3d.pdf"%i
        assert name.endswith(".pdf")
        if i%10 == 0:
            c.writePDFfile(name)
            os.system("pdftocairo -png %s"%(name))

        i += 1

    print


if __name__ == "__main__":

    main()