#!/usr/bin/enb python
# -*- coding: utf-8 -*-
#
# Sample script to generate a PNM image
# of the Newton's fractal associated with
# x^\alpha - 1 polynomial.
#
 
def GetNewtonConvergenceSpeed(z, maxit = 255, eps = 10e-11, alpha = 3):
    """
    This function returns, given a complex number z,
    an integer between 0 and maxit representing the
    number of iteration to approximate a root of
    the polynomial x^\alpha - 1 = f(x) so that
    |f(x)| < eps.
    """
 
    print " => Iterating on %f + i%f" % (z.real, z.imag)
 
    fz = pow(z, alpha) - 1
    iterations = 0
 
    while (abs(fz) > eps and iterations < maxit):
        t = pow(z, alpha - 1)
        fz = t * z - 1
        fpz = alpha * t
 
        # Questo per prevenire la divisione per zero
        # (tanto la nostra derivata si annulla solo lì)
        if (fpz == 0):
            return maxit
 
        z = z - fz / fpz
        iterations += 1
 
    return iterations
 
def Newton(size = 200):
    """Compute the newton's method on a net of points
    in the set { x + iy \in C | |x| < 2 and |y| < 2 }
    and fill the image matrix with the value returned
    by the GetNewtonConvergenceSpeed ()
    """
 
    # Delta è l'ampiezza dell'intervallo fratto il numero di
    # divisioni - 1
    delta = 4 / float(size - 1)
    print " => delta = %f" % delta
 
    x_values = range(0 , size)
    y_values = x_values
 
    image = []
 
    for x in x_values:
        image.append([])
        for y in y_values:
            image[x_values.index(x)].append(
                GetNewtonConvergenceSpeed(complex(-2 + x * delta,
                                                  -2 + y * delta)))
 
    f = open("prova.pnm", 'w')
    f.write("P3\n")
    f.write(str(size) + " " + str(size) + "\n")
    f.write("255\n")
 
 
    for x in x_values:
        for y in y_values:
            value = image[x][y]
            f.write( str(3 * value) + " " +
                     str(5 * value) + " " +
                     str(9 * value) + "\n")
 
 
if __name__ == "__main__":
 
    Newton (1024)