mandelbrot.f90

! 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