!-----------------------------------------------------
! 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)
!-----------------------------------------------------
SUBROUTINE Mandelb(data, ext, Nx, Ny, max_iterations)
  IMPLICIT NONE ! Don't use any implicit names of variables!
  ! Function arguments
  REAL*8, intent(out) :: data(Nx,Ny)
  REAL*8, intent(in) :: ext(4)
  INTEGER, intent(in) :: max_iterations
  INTEGER, intent(in) :: Nx, Ny
  !f2py integer optional, intent(in)          :: max_iterations=1000
  !  !f2py integer intent(hide), depend(data) :: Nx=shape(data,0)
  !  !f2py integer intent(hide), depend(data) :: Ny=shape(data,1)
  ! Local variables
  INTEGER    :: i, j, itt
  COMPLEX*16 :: z0, z
  data(:,:) = 0.0
  DO i=1,Nx
     DO j=1,Ny
        z0 = dcmplx( ext(1) + (ext(2)-ext(1))*(i-1)/(Nx-1.), ext(3) + (ext(4)-ext(3))*(j-1)/(Ny-1.) )
        z=0
        DO itt=1,max_iterations
           IF (abs(z)>2.) THEN
              data(i,j) = 1./itt            ! result is number of iterations
              EXIT
           ENDIF
           z = z**2 + z0             ! f(z) = z**2+z0 -> z
        ENDDO
     ENDDO
  ENDDO
  RETURN
END SUBROUTINE Mandelb