#!/usr/bin/env python
# v1.1 29 December 2007
# Numerical Python program for making a julia set in a PPM file.
# uses numerical arrays of Numeric, MA, or numpy
# Here is what the program does.  For every (x,y) point within
# a rectangular grid, the point is treated as complex number z.
# The following very simple equation is interated for up to
# itermax times:
# z=z^2 +c
# Depending on whether that iteration converges or diverges,
# a color is assigned to the (x,y) point.  By diverge ("breakout"), it is
# meant that |z|^2 becomes greater than the number breakout
# with itermax iterations.
import sys,getopt 
from math import sqrt
from os import system 
###############
#choose one of Numeric, MA or numpy:
from Numeric import *
#from MA import * #masked arrays will be *slower*
#from numpy import *
###############
#default arguments:
filename="fromnumeric.ppm"
h=400
w=600
magnify=10
itermax=100
breakout=64.
cr=-.7492
ci=.1
x0=.09950
y0=-.00062
displayppm=False
#change defaults from commandline
(options,arguments)=getopt.getopt(sys.argv[1:],'h:w:b:m:i:d',["cr=","ci=","x0","y0="])
for op,ar in options:
	if op=="--cr": cr=float(ar)
	if op=="--ci": ci=float(ar)
	if op=="--x0": ci=float(ar)
	if op=="--y0": ci=float(ar)
	if op=="-m": magnify=float(ar)
	if op=="-i": itermax=int(ar)
	if op=="-h": h=int(ar)
	if op=="-w": w=int(ar)
	if op=="-b": breakout=float(ar)
	if op=="-d": displayppm=True
rb=sqrt(1.*breakout)
x=resize(arange(1,w+1.),(h,w))
y=transpose(resize(arange(h+1.,1.,-1.),(w,h)))
x=4.*(x/w-.5)/magnify+x0
y=(4.*h/w)*(y/h-.5)/magnify+y0 
zm=zeros((h,w),'f')
nbr=ones((h,w),'f') 
iter=0
xl=array(x) # makes copy, same as x.copy()
yl=array(y)
while iter<itermax:
	xl=where(nbr,x,xl)
	yl=where(nbr,y,yl)
	x=xl*xl - yl*yl + cr
	y=2*xl*yl + ci
	zm=maximum(x*x+y*y,zm)
	nbr=less(zm,breakout)
#	x=masked_where(logical_not(nbr),x,copy=0) #use with MA
#	y=masked_where(logical_not(nbr),y,copy=0) #use with MA
	iter+=1
	print "iter=",iter,"  breakouts=",1.-1.*sum(nbr.flat)/size(nbr)
f0=zeros((h,w),'f')
red=where(nbr,f0,255.*(rb+xl)/(2*rb))
green=where(nbr,255.*zm/breakout,f0) 
blue=where(nbr,255.*zm/breakout,.5*255.*(rb+yl)/(2*rb))
rgb=zeros((h,w,3),'b') # b means byte 
header="P6 %d %d %d\n" % (w,h,255)
rgb[:,:,0]=red.astype('b')
rgb[:,:,1]=green.astype('b')
rgb[:,:,2]=blue.astype('b')
ppm=open(filename,'wb') # b means binary here, needed for Windows, irrelevant in Unix/Linux
ppm.write(header+rgb.tostring())
ppm.close()
print "wrote ",w,"by",h, "ppm file: ",filename 
if displayppm: system("display "+filename)