Function for drawing the Feigenbaum diagram

feigenbaum.py

import random

from PIL import Image


def feigenbaum_diagram(scale=2000, r_ax=None, x_ax=None, name="feigenbaum_diagram.png"):
    if not x_ax:
        x_ax = [0.0, 1.0]
    if not r_ax:
        r_ax = [0.0, 1.0]

    xdiff = abs(r_ax[1] - r_ax[0])
    ydiff = abs(x_ax[1] - x_ax[0])

    if xdiff > ydiff:
        box_ratio = ydiff / xdiff
        width = scale
        height = round(scale * box_ratio)
    else:
        box_ratio = xdiff / ydiff
        height = scale
        width = round(scale * box_ratio)

    img = Image.new("RGB", (width + 1, height + 1), "white")

    start = random.uniform(0, 1)

    print(width, height)

    for i in range(width):
        old = start
        r = r_ax[0] + i * (r_ax[1] - r_ax[0]) / width
        for iteration in range(200):
            new = 4 * r * old * (1 - old)
            old = new
            if iteration > 100:
                x = old * height
                # flip the diagram
                if x < height / 2:
                    x = round(x + (height - 2 * x))
                else:
                    x = round(x - (height - 2 * (height - x)))

                if 0 <= i <= width and 0 <= x <= height:
                    img.putpixel((i, x), (0, 0, 0))

    img.save(name)


feigenbaum_diagram(2000, [0.65, 1.0], [0, 1], name="feigenbaum.png")
feigenbaum_diagram(2000, [0.85, 0.90], [0.8, 0.9], name="feigenbaum_zoom.png")