#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# Questo script genera il frattale di newton
# associato al polinomio x^alpha - 1; Può essere
# divertente cambiare alpha ed eventualemnte
# anche il numero di pixel dell'immagine.
#
# Ovviamente, essendo scritto in python, è piuttosto
# lento, ma comunque crea le immagini in tempi
# ragionevoli.
 
def GetNewtonConvergenceSpeed(z, maxit = 255, eps = 10e-11, alpha = 3):
    """
    Ritorna, dato un numero complesso z ed un eps > 0, il numero
    di iterazioni necessarie (del metodo di Newton) perché
    |f(x_k)| < eps, oppure maxit se dopo maxit iterazioni
    non si è ancora raggiunta la precisione desiderata.
    """
 
    fz = pow(z, alpha) - 1
    iterations = 0
 
    # Finchè non abbiamo la necessaria precisione continuiamo
    # ad iterare. Usciamo solamente se sforiamo il limite che
    # ci siamo dati.
    while (abs(fz) > eps and iterations < maxit):
 
        # Usiamo una var temporanea per non dover
        # ricalcolare le parti di comuni di funzione
        # e derivata.
        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 in z=0,
        # e perdere un punto non è poi tanto grave)
        if (fpz == 0):
            return maxit
 
        # Iterazione del metodo di newton
        z = z - fz / fpz
        iterations += 1
 
    return iterations
 
def Newton(size = 200, filename="newton.pnm"):
    """
    Calcola la velocità di convergenza su una rete di punti
    size x size sul quadrato di [-2,2]^2 attorno all'origine
    del piano complesso. La salva poi in un'immagine dal nome
    newton.pnm
    """
 
    # delta è l'ampiezza dell'intervallo fratto il numero di
    # divisioni - 1
    delta = 4 / float(size - 1)
 
    # Creiamo gli indici su cui iterare
    x_range = xrange(0 , size)
    y_range = x_range
 
 
    # Apriamo il file dove salveremo il nostro lavoro
    # e scriviamo l'intestazione del file PNM.
    f = open(filename, 'w')
    f.write("P3\n")
    f.write(str(size) + " " + str(size) + "\n")
    f.write("255\n")
 
    # Scriviamo la matrice su file
    for x in x_range:
        for y in y_range:
            # Otteniamo un valore "rappresentativo" della velocità di convergenza
            value = GetNewtonConvergenceSpeed(complex(-2 + x * delta,
                                                      -2 + y * delta))
            # Scriviamo i valori RGB moltiplicandoli per dei coefficienti in modo
            # da avere più blu.
            f.write( " ".join( [str(3 * value), str(5 * value), str(9*value)] )+ "\n")
    f.close ()
 
 
if __name__ == "__main__":
 
    # Creiamo un'immagine 800px x 800px.
    Newton (800, "newton.pnm")