source: branches/numpy/pypar-numeric/mandelbrot_example/mandelbrot.py @ 6883

Last change on this file since 6883 was 5779, checked in by steve, 16 years ago

Added the old version of pypar which works with Numeric. Necessary for parallel code until we move anuga to numpy (and then we can use pypar as distribute via sourceforge).

File size: 1.8 KB
Line 
1"""Fundamentals for computing the Mandelbrot set
2"""
3
4def calculate_point(c, kmax):
5    """Calculate one point of the Mandelbrot set by iterating the
6    governing equation
7   
8        z_{k+1} = z_k^2 + c,  k = 0, 1, ... kmax-1
9        z_0 = 0 + 0i       
10   
11    for as long as k < kmax  and |z_k| <= 2
12   
13    Inputs:
14      c:    A complex number for which the iteration is computed
15      kmax: The maximal number of iterations
16     
17    Output:
18      count: The value of k after the iteration has completed.
19    """
20
21    z = complex(0,0) #Create complex number with initial value
22    k = 0            #Initialise iteration counter
23
24    while k < kmax and abs(z) <= 2:
25        z = z*z + c
26        k = k+1
27
28    return k
29
30
31
32def calculate_region(real_min, real_max, imag_min, imag_max, kmax, M, N):
33    """Calculate the mandelbrot set in a given rectangular subset of the complex plane.
34    Inputs:
35       real_min: Left boundary
36       real_max: Right boundary
37       imag_min: Lower boundary
38       imag_max: Upper boundary       
39       kmax: Maximal iteration count
40       M: Number of points along the real axis
41       N: Number of points along the imaginary axis
42    Output:
43       Matrix A (MxN) containing iteration counts for each point   
44    """
45   
46    from Numeric import zeros
47    from mandel_ext import calculate_point #Faster C-implementation
48
49    #Compute resolution
50    real_step = (real_max-real_min)/
51    imag_step = (imag_max-imag_min)/N
52   
53   
54    A = zeros((M, N))   # Create M x N matrix
55   
56    #Compute Mandelbrot iteration for each point in the rectangular subset
57    #and store iteration counts in matrix A                   
58    for i in range(M):
59        for j in range(N):
60            c = complex(real_min + i*real_step, imag_min + j*imag_step)
61            A[i,j] = calculate_point(c, kmax)
62           
63    return A
Note: See TracBrowser for help on using the repository browser.