public
Last active

Mandelbrot with NumPy and Fortran

  • Download Gist
mandelbrot.f90
FORTRAN
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
program Mandelbrot
use types, only: dp
use constants, only: I
use utils, only: savetxt, linspace, meshgrid
implicit none
 
integer, parameter :: ITERATIONS = 100
integer, parameter :: DENSITY = 1000
real(dp) :: x_min, x_max, y_min, y_max
real(dp), dimension(DENSITY, DENSITY) :: x, y
complex(dp), dimension(DENSITY, DENSITY) :: c, z
integer, dimension(DENSITY, DENSITY) :: fractal
integer :: n
x_min = -2.68_dp
x_max = 1.32_dp
y_min = -1.5_dp
y_max = 1.5_dp
 
call meshgrid(linspace(x_min, x_max, DENSITY), &
linspace(y_min, y_max, DENSITY), x, y)
c = x + I*y
z = c
fractal = 255
 
do n = 1, ITERATIONS
print "('Iteration ', i0)", n
where (abs(z) <= 10) z = z**2 + c
where (fractal == 255 .and. abs(z) > 10) fractal = 254 * (n-1) / ITERATIONS
end do
 
print *, "Saving..."
call savetxt("fractal.dat", log(real(fractal, dp)))
call savetxt("coord.dat", reshape([x_min, x_max, y_min, y_max], [4, 1]))
end program
mandelbrot.py
Python
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
import numpy as np
 
ITERATIONS = 100
DENSITY = 1000
x_min, x_max = -2.68, 1.32
y_min, y_max = -1.5, 1.5
 
x, y = np.meshgrid(np.linspace(x_min, x_max, DENSITY),
np.linspace(y_min, y_max, DENSITY))
c = x + 1j*y
z = c.copy()
fractal = np.zeros(z.shape, dtype=np.uint8) + 255
 
for n in range(ITERATIONS):
print "Iteration %d" % n
mask = abs(z) <= 10
z[mask] *= z[mask]
z[mask] += c[mask]
fractal[(fractal == 255) & (-mask)] = 254. * n / ITERATIONS
 
print "Saving..."
np.savetxt("fractal.dat", np.log(fractal))
np.savetxt("coord.dat", [x_min, x_max, y_min, y_max])

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.