public
Created

  • 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 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94
! Condensed modules with only the necessary functions
module types
implicit none
private
public dp
integer, parameter :: dp=kind(0.d0)
end module
module constants
use types, only: dp
implicit none
private
public I
complex(dp), parameter :: I = (0, 1)
end module
module utils
use types, only: dp
implicit none
private
public savetxt, linspace, meshgrid
contains
subroutine savetxt(filename, d)
character(len=*), intent(in) :: filename
real(dp), intent(in) :: d(:, :)
integer :: s, i
open(newunit=s, file=filename, status="replace")
do i = 1, size(d, 1)
write(s, *) d(i, :)
end do
close(s)
end subroutine
function linspace(a, b, n) result(s)
real(dp), intent(in) :: a, b
integer, intent(in) :: n
real(dp) :: s(n)
integer :: i
do i = 1, n
s(i) = a + (b-a) * (i-1) / (n-1)
end do
end function
subroutine meshgrid(x, y, x2, y2)
real(dp), intent(in) :: x(:), y(:)
real(dp), intent(out) :: x2(:, :), y2(:, :)
x2 = spread(x, 1, size(y))
y2 = spread(y, 2, size(x))
end subroutine
end module
 
 
 
! Main program
 
 
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
real(dp) :: t1, t2, t3
call cpu_time(t1)
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 cpu_time(t2)
call savetxt("fractal.dat", log(real(fractal, dp)))
call savetxt("coord.dat", reshape([x_min, x_max, y_min, y_max], [4, 1]))
call cpu_time(t3)
print *, "Times:"
print *, "Calculation:", t2-t1
print *, "Saving:", t3-t2
print *, "Total:", t3-t1
end program

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.