program henon

  implicit none

  integer                      :: passe
  double precision             :: vx, vy

  integer                      :: w, h
  integer                      :: foo, bar
  double precision             :: px, py
  w = 2000   ;   h=1600

  write(0, *) "###### Mapping of Henon "

  do foo=1, 16
    px = dble(foo) / 16.0
    do bar=1,16
      py = dble(bar) / 16.0
      call compute_pixel_henon(px, py, 1700,    &
                           passe, dble(0.5), vx, vy)
      write(0, fmt=*) "passe ", passe, vx, vy
    enddo
  end do

!-----------------------------------------------------
contains
!-----------------------------------------------------

!-----------------------------------------------------
subroutine compute_pixel_henon(a, b, maxpasse, passe, limit, rx, ry)
  implicit none
  double precision, intent(in)         :: a, b, limit
  integer, intent(in)                  :: maxpasse
  integer, intent(out)                 :: passe
  double precision, intent(out)        :: rx, ry

  double precision                     :: x, y, x2, y2

  write(0, fmt="('compute pixel:', (2F8.3, I6, F8.3))")       &
                              a, b, maxpasse, limit

  x = 0.0
  y = 0.0

  do passe=1, maxpasse

    x2 = 1d0 + y - (a * x * x)
    y2 = b * x
    x = x2
    y = y2
    write(0, fmt="(i4, 2F8.3)") passe, x, y
    if (x .lt. -limit) exit
    if (x .gt.  limit) exit
    if (y .lt. -limit) exit
    if (y .gt.  limit) exit

  enddo

  rx = x
  ry = y

end subroutine

!-----------------------------------------------------

end program