!-----------------------------------------------------
! Produces Mandelbrot plot in the range [-2:1]x[-1:1]
! It uses formula z = z*z + z0 iteratively until
! asb(z) gets bigger than 2 
! (deciding that z0 is not in mandelbrot)
! The value returned is 1/(#-iterations to escape)
!-----------------------------------------------------

REAL*8 FUNCTION Mandelb(z0, max_iterations)
  IMPLICIT NONE ! Don't use any implicit names of variables!
  ! Function arguments
  COMPLEX*16 :: z0   ! Need to declare all variables first
  INTEGER    :: max_iterations
  ! Local variables
  INTEGER    :: i
  COMPLEX*16 :: z
  z=0                 ! Implementation after declarations
  DO i=1,max_iterations
     IF (abs(z)>2.) THEN
        Mandelb = i            ! result is number of iterations
        RETURN
     ENDIF
     z = z**2 + z0             ! f(z) = z**2+z0 -> z
  ENDDO
  Mandelb = max_iterations*1e3 ! choose some very large number
  RETURN
END FUNCTION Mandelb

!---------- Main program -----------------!
PROGRAM mand     ! Main part of the program
  IMPLICIT NONE  ! Don't use any implicit names of variables!
  REAL*8    :: Mandelb
  INTEGER   :: Nx, Ny, i, j
  REAL*8    :: x, y
  COMPLEX*16:: z0
  INTEGER   :: max_iterations
  Nx = 400                ! the mesh in complex plane
  Ny = 400
  max_iterations = 1000
  
  DO i=1,Nx
     DO j=1,Ny
        x = -2 + 3.*(i-1)/(Nx-1.)
        y = -1 + 2.*(j-1)/(Ny-1.)
        z0 = dcmplx(x,y)
        print *, x, y, 1/Mandelb(z0,max_iterations) ! call to the function
     ENDDO
  ENDDO
END PROGRAM mand