#-----------------------------------------------------------------||||||||||||--
# Name:          chaos.py
# Purpose:       chaotic and fractal utilities.
#
# Authors:       Christopher Ariza
#
# Copyright:     (c) 2005 Christopher Ariza
# License:       GPL
#-----------------------------------------------------------------||||||||||||--



import math



#-----------------------------------------------------------------||||||||||||--
def _fibonacciNumber(goldenUpper, goldenLower, i):
    """return i-th number in the Fibonacci series
    """ 
    i = i + 1
    n = ((math.pow(goldenUpper, i) - math.pow(goldenLower, i)) /
            math.sqrt(5.0))
    return math.floor(n)

def fibonacciSeries(j, k):
    """return a subset of the fibonacci series.
    Syntax:
    fibonacciSeries(n)
    Return Fiboacci terms from 0 to n-1.
    fibonacciSeries(j, k)
    Return Fibonacci terms from j to k-1
    """
    goldenUpper = (1.0 + math.sqrt(5.0)) / 2.0
    goldenLower = (1.0 - math.sqrt(5)) / 2.0
    return [_fibonacciNumber(goldenUpper, goldenLower, i) for i in range(j, k)] 
    
def fibonacciSuccessor(n):
    """given a number from the Fibonacci series, return its
    sucessor in the same series.
    """
    return math.floor((n + 1.0 + math.sqrt(5.0 * n ** 2)) / 2.0)

#-----------------------------------------------------------------||||||||||||--
def verhulst(p, x):
    """logistic map or verhulst equation, w/o state variables"""
    return p * x * (1 - x)


#-----------------------------------------------------------------||||||||||||--
# note: usful to get values for x, (x, y) alternating, and x, y, z, alternating

# x(n+1)=y(n)+1-ax(n)^2
# y(n+1)=bx(n)
# X and Y are the initial values for the system and should both be a number close to zero. A and B control the behaviour of the system. Both should  be a positive value. The system will be stable if the values for A and B    are set following this general guideline:
#A<=(2.0-B)
# or; A+B <= 2.0
#It is not possible to predict the maximum and minimum of the output     values. To 
     
# test values
# a, b
# 1.4, 0.3: chaotic, b/n -1 and 1.2
# .4, .3: period 2, b/n 0 and 1
# 0, .3: cons, b/n 0 and 1
# .1, .1: cons, b/n 0 and 1
# 1, 1: overflow error
# 2, 2: negative infinity
# .5, .5: periodic 2, b/n -1 and 1
# 1.5, .2: chaotic, b/n -1 and 1
# .2, 1.01: used mathworld example

# reccomended x,y start points for strange attractor
# x = 0.63135448, y = 0.18940634.

class Henon:

    def __init__(self, a=1.4, b=0.3, x=0.63135448, y=0.18940634):
        self.a = a
        self.b = b
        self.x = x
        self.y = y
        
    def __call__(self, a=None, b=None):
        """can provide dynamic a and b values
        values are not scaled or altered from original
        in case of overflow error, will distort to accomodate"""
        if a != None: self.a = float(a)
        if b != None: self.b = float(b)
        try:
            xTemp = 1.0 - (self.a * pow(self.x, 2.0)) + self.y
        except OverflowError:
            xTemp = 0.0
        yTemp = self.b * self.x
        self.x = xTemp
        self.y = yTemp
        return self.x, self.y
        
        

#-----------------------------------------------------------------||||||||||||--


# x'=sigma*(y-x)
# y'=Rx-y-xz
# z'=xy-Bz
# 
# X, Y, and Z set the initial conditions. SIGMA, R, and B control the  behaviour of the system. In general, higher values for R and B increase  the chaos, while the same is true for lower values for SIGMA:
# It is not possible to predict the maximum and minimum of the output  values.

# The Lorenz attractor is an attractor that arises in a simplified system of equations describing the two-dimensional flow   of fluid of uniform depth, with an imposed temperature difference DeltaT, under gravity g, with buoyancy a, thermal diffusivity k, and kinematic viscosity v.
# ... representing one of the earliest discoveries of the so-called butterfly effect.
# where X is proportional to convective intensity, Y to the temperature  difference between descending and ascending currents, and Z to the difference  in vertical temperature profile from linearity.
# where [sigma] is the Prandtl number, Ra is the Rayleigh number, Ra, is the critical Rayleigh number, and b is a geometric factor  (Tabor 1989, p. 206). Lorenz took b=8/3 and [sigma]=10 (2.666).

# stable values for positive r if sigma greater than b + 1
# r = 28 is chaotic, r 99.96 is periodic

# b is often 8/3 (2.6666666666666665): others 3/8 (.375)
# act examples:
# s, r, b
# 10 28 1.0
# 10 28 3/8
# 10 28 5.2
# 2  28 3/8
# 10 40 3/8 

class Lorenz:
    """this does not work: produces values that are too small"""
    def __init__(self, r=28, s=10.0, b=2.6666666666, 
                             x=1.0, y=1.0, z=1.0):   
        self.r = float(r)
        self.s = float(s)
        self.b = float(b)

        self.x = float(x)
        self.y = float(y)
        self.z = float(z)
    
        # dont know where this value of d is coming from
        # little uses .003
        # others use 1.0/128 (0.0078125)
        # mikelson uses .001
        # berg uses .01
        self.d = 0.01
        
    def __call__(self, r=None, s=None, b=None):
        if r != None: self.r = float(r)
        if s != None: self.s = float(s)
        if b != None: self.b = float(b)
        
        xTemp = self.s * (self.y - self.x)
        yTemp = (self.r * self.x) - self.y - (self.x * self.z)
        zTemp = (self.x * self.y) - (self.b * self.z)
        self.x = self.x + (self.d * xTemp)
        self.y = self.y + (self.d * yTemp)
        self.z = self.z + (self.d * zTemp)
        return self.x, self.y, self.z
    




#-----------------------------------------------------------------||||||||||||--


if __name__ == "__main__":
    a = Henon()
    b = Lorenz()
    for obj in [a, b]:
        print '\n', obj.__class__
        for i in range(500):
            out = obj()
            if len(out) == 2:
                x, y = out
                print str(x).ljust(20), str(y).ljust(20)
            elif len(out) == 3:      
                x, y, z = out
                print str(x).ljust(20), str(y).ljust(20), str(z).ljust(20)
        


# end
