fibonacci-spiral-python-matplotlib

fibonacci-spiral.py

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
import math

# Ported from Matlab code provided in "The “Other” Fibonacci Spiral and Binet Spirals" - Cye H. Waldman ([email protected])
# http://old.nationalcurvebank.org//waldman7/WaldmanOtherFibonacciSpiralandBinetSpirals.pdf

def HyperbolicGeneralBinet(n: np.ndarray, a: int, b: int, c: int) -> (np.ndarray, np.ndarray): # (chb, shb)
    """
    generalized hyperbolic Binet function
    c=-1 is the Fibonacci function
    f(0)=0, f(1)=1, f(n+1)=a*f(n)+b*f(n-1) for n>=2
    c=-1 is the Lucas function
    f(0)=2, f(1)=a, f(n+1)=a*f(n)+b*f(n-1) for n>=2
    """
    
    alpha = (a + math.sqrt(a**2 + 4 * b)) / 2
    beta = (a - math.sqrt(a**2 + 4 * b)) / 2
    
    chb = (alpha**n + abs(beta)**n)
    shb = (alpha**n - abs(beta)**n)
    
    if c == -1:
        chb = chb / (alpha + (c * beta))
        shb = shb / (alpha + (c * beta))
        
    return chb, shb
    
def GeneralBinetFn(n: np.ndarray, a: int , b: int, c: int) -> np.ndarray: # (y)
    """
    % generalized Binet function
    % c=-1 is the Fibonacci function
    % f(0)=0, f(1)=1, f(n+1)=a*f(n)+b*f(n-1) for n>=2
    % c=-1 is the Lucas function
    % f(0)=2, f(1)=a, f(n+1)=a*f(n)+b*f(n-1) for n>=2
    """
    
    alpha = (a + math.sqrt(a^2 + 4 * b)) / 2
    beta = (a - math.sqrt(a^2 + 4 * b)) / 2
    
    y = (alpha**n + c * beta**n)
    
    if c == -1:
        y = y / (alpha + c * beta)
        
    return y
        
def ShofarSurfaces() -> None:
    """
    b=1 and a=1,2,3 produces gold, silver, bronze shofars, respectively
    """
    
    a = 1
    b = 1
    
    # c = -1 --> Fibonacci-type
    # c = 1  --> Lucas-type
    c = -1
    
    # 1 --> enabled bounding curve plot
    # 0 --> disable bounding curve plot
    plotcurve = 0
    
    xpts = 501
    xmax = 5
    x = np.linspace(-xmax, xmax, xpts)
    cfs, sfs = HyperbolicGeneralBinet(x, a, b, c)
    ctr = (cfs + sfs) / 2
    r = (cfs - ctr) * 0.99
    X = np.meshgrid(x)
    t = np.linspace(0, 2 * math.pi, xpts)
    [R, T] = np.meshgrid(r, t)
    CTR = np.meshgrid(ctr)
    
    # in physical space this should be y-coordinate
    Z = np.add(np.multiply(R, np.sin(T)), CTR)
    
    # in physical space this should be z-coordinate
    Y = R * np.cos(T)
    
    #figure
    fig = plt.figure(figsize=(9, 6))
    ax = fig.gca(projection='3d')
    
    #h = surf(X, Y, Z);
    # consider -X when the bugle is backwards
    surf = ax.plot_surface(X, Y, Z, cmap='plasma',
                           linewidth=0, antialiased=False)

    if c == 1:
        ax.set_title(f'Generalized Binet-Lucas: [a b] = [{a} {b}]') #', 'FontSize', 14)
    else:
        ax.set_title(f'Generalized Binet-Fibonacci: [a b] = [{a} {b}]') #','FontSize',14)

    if plotcurve:
        #hold on
        z = GeneralBinetFn(x, a, b, c)
        ax.plot3D(x, np.imag(z), np.real(z)) #, 'r', 'LineWidth', 2)
        Zs = np.zeros(x.shape[0])
        ax.plot3D(x, Zs, cfs) #, 'y', 'LineWidth', 2)
        ax.plot3D(x, Zs, sfs) #, 'y', 'LineWidth', 2)

    plt.show()

ShofarSurfaces()