from pylab import *
from scipy import weave

# The C++ code
code="""
     using namespace std; // for using cout when debugging
     for (int i=0; i<Nx; i++){
         for (int j=0; j<Ny; j++){
             complex<double> z0( ext(0)+(ext(1)-ext(0))*i/(Nx-1.), ext(2)+(ext(3)-ext(2))*j/(Ny-1.) );
             complex<double> z=0.0;
             for (int itt=0; itt<max_iterations; itt++){
                 if (norm(z)>4.) { data(i,j)=1./itt; break; }
                 z = z*z + z0;     // if |z|>2 the point is not part of mandelbrot set
             }
         }
     }
     """

Nx = 400
Ny = 400
data = zeros((Nx,Ny), 'd')
max_iterations=1000
ext=array([-1.8,-1.72,-0.05,0.05]) # The range of the mandelbrot plot [x0,x1,y0,y1]

weave.inline(code, ['data', 'Nx', 'Ny', 'max_iterations', 'ext'],
             type_converters=weave.converters.blitz, compiler = 'gcc')
        
# Using python's pylab, we display pixels to the screen!
data = data.transpose()
imshow(data, interpolation='bilinear', cmap=cm.hot, origin='lower', extent=ext, aspect=1.)
colorbar()
show()