#!/usr/pkg/bin/python
"""minimalist program to create Mandelbrot set

   this creates a PPM file of the full set centered at -0.5+0j
   based on mandelgen.py, but this program is released to the 
   public domain by the author, jcomeau at unternet dot net
"""
import sys, os

current = {
 'x': -0.5, # center of this image, x (real) and y (imaginary)
 'y': 0.0,
 'width': 2048,  # in pixels
 'height': 1536,
 'ratio': 0,  # to be set in init()
 'zoom': 1.0,  # a zoom ratio of 1 shows every point approaching the set
 'output': sys.stdout,  # stdout means the output must be redirected to a file
 'iterations': None,  # to be set in init()
 'colors': None,  # set in init()
 'maxiterations': 0,  # brighter colors for fewer iterations; the Set is black
 'maxvalue': 4,  # x^2 + y^2
 'constants': None,  # to be set in init()
 'values': None,  # to be set in init()
}

def mandelmin():
 """main routine

    creates PPM image of Mandelbrot set at -0.5+0j
 """
 init()
 gen_by_pixel()
 ppmformat()

def gen_by_pixel():
 """this routine does the work, calculating z->z^2+z0 for each pixel
    either until it escapes (a^2+b^2>=2), or proves that it won't, 
    or the maximum iterations we decided upon have been exhausted
 """
 for index in range(len(current['constants'])):
  z, fifo = current['constants'][index], []
  while True:
   z = complex_add(complex_multiply(z, z), current['constants'][index])
   if z in fifo:  # value is oscillating, will never escape
    current['iterations'][index] = current['maxiterations']
    break
   fifo.append(z)
   try:
    if ((z[0] ** 2) + (z[1] ** 2)) < current['maxvalue']:
     current['iterations'][index] += 1
   except:
    break
   if current['iterations'][index] == current['maxiterations']:
    break

def complex_multiply(i, j):
 """multiply two complex numbers the old-fashioned way
    (before complex-number-aware languages like Python came along)"""
 [a, b], [c, d] = i, j
 return [(a * c) - (b * d), (b * c) + (a * d)]

def complex_add(i, j):
 """add two complex numbers the old-fashioned way
    (before complex-number-aware languages like Python came along)"""
 [a, b], [c, d] = i, j
 return [a + c, b + d]

def heightvalue(pixelheight):
 """determine the imaginary part of the number represented by a given pixel"""
 heightspan = current['top'] - current['bottom']
 return current['top'] - ((heightspan * pixelheight) / current['height'])

def widthvalue(pixelwidth):
 """determine the real part of the number represented by a given pixel"""
 widthspan = current['right'] - current['left']
 return current['left'] + ((widthspan * pixelwidth) / current['width'])

def init():
 """initialize any values not yet set in concrete"""
 if current['colors'] == None:
  current['colors'] = rgbvalues()
  # maximum iterations must be one less than the number of colors; if there
  # are just two colors, black and white, one iteration is sufficient
  current['maxiterations'] = len(current['colors']) - 1
  current['ratio'] = 640.0 / (current['width'] * 100)
  current['iterations'] = [0] * (current['width'] * current['height'])
  current['constants'] = list(current['iterations'])
  current['ratio'] = float(current['ratio']) / current['zoom']
  current['bottom'] = current['y'] - \
   ((current['height'] * current['ratio']) / 2)
  current['top'] = current['bottom'] + (current['height'] * current['ratio'])
  current['left'] = current['x'] - ((current['width'] * current['ratio']) / 2)
  current['right'] = current['left'] + (current['width'] * current['ratio'])
  for n in range(len(current['constants'])):
   x, y = n % current['width'], int(n / current['width'])
   current['constants'][n] = [widthvalue(x), heightvalue(y)]

def binary(binary_string):
 """convert a binary string to a number"""
 number = 0
 for character in binary_string:
  number <<= 1
  number |= (character == "1")
 return number

def rgbvalues():
 """define start color as cyan: r=0 g>=b

    decrement for each iteration such that 0 is black
    to match the colorForth program, g must be 6 bits and b 5, with the low
    bits of each cleared, and two subtracted each time
    also, found out that the g and b values must be shifted into the
    high-order bits of their respective bytes, to match colorForth results
 """
 cyan = binary("11111011110")
 divisor = binary("100000")
 values = [(0, (color / divisor) << 2,
  (color % divisor) << 3) for color in range(cyan, -2, -2)]
 return values

def ppmformat():
 """write out the pixel informated generated as a Portable Pixel Map"""
 current['output'].write('%s %d %d %d\n' % ('P6',
  current['width'], current['height'], 255))
 pixels = rgbvalues()
 for index in range(len(current['iterations'])):
  r, g, b = pixels[current['iterations'][index]]
  current['output'].write('%c%c%c' % (r, g, b))

if __name__ == '__main__':  # meaning this program was executed directly
 mandelmin()