c----------------------------------------------------------------
c           *****************************************
c           *   THIS FILE IS ARCHIVED FOR HISTORY   *
c           *      DON'T TRY TO USE IT AT HOME      *
c           *****************************************
c----------------------------------------------------------------

c       	mapping of the Henon Diagram
c       a escape time fractal from 'Oulala Software'
c       http://www.chez.com/oulala/g77/map_henon.html
c
c----------------------------------------------------------------
        program	map_henon
c
c       control the image dimension
c
      integer		w, h
      parameter	(w=2000, h=1600)

c       maximum number of iterations by pixel
c       around 5000 seems a good value
c
      integer  maxpasse
      parameter (maxpasse=17000)

      double precision amp, limite
      parameter  (amp=1.42d0, limite=4.0d9)

c	control the position of the view window
c
      double precision a, b
      parameter  (a=0.65d0, b=0.15d0)

      integer 	x, y, image, passe
      integer 	masque, flag
      integer 	valr, valg, valb
      integer 	in, out
      double precision dx, dy, rx, ry
      double precision minrx, minry, maxrx, maxry
      real 	temps

      integer IMG_CREATE
c	---------------------------------------------------------

      image = IMG_CREATE(w, h, 3)
      masque = IMG_CREATE(w, h, 3)

      in = 0
      out = 0
      minrx = -9d99
      minry = -9d99
      maxrx =  9d99
      maxry =  9d99

      do 5 x=0, w-1
          dx = (a-amp) + (amp*2/(w-1)*x)
          do 4 y=0, h-1
              dy = (b-amp) + (amp*2/(h-1)*y)

             call henon(dx, dy, maxpasse, passe, limite, rx, ry)

       if (passe .gt. maxpasse) then
           in = in + 1
           valr = 128
           valg = abs(int(rx*1d2))
           valb = abs(int(ry*1d2))
           flag = 0
       else
           out = out + 1
           if (passe .le. 2) then
        valr = 0
        valg = 80 * passe
        valb = 10
               flag = 255
           else
            valg = (passe*255)/maxpasse
            valr = mod(passe,255)
            valb = (255-((passe*255)/maxpasse))/3
               flag = 128
           endif
       endif

C	        write(6,101) x, y, dx, dy, passe, valr, valg, valb, rx, ry

          	call IMG_PLOT(image, x, y, valr, valg, valb)
          	call IMG_PLOT(masque, x, y, flag, flag, flag)

4	    continue
      write(6,107) x, w, in, out
5	continue

c
c	Here, you can change the name of the saved TGA file.
c	Just modify the second parameter...
c
      call IMG_TGA_SA(image, "map_henon.tga", 0)
      call IMG_TGA_SA(masque, "msq_henon.tga", 0)

      temps = second() / 60.0
      write(6,120) temps

101   format(2i4, 2(2x,f9.5), i7, 3(i4), 2(2x,f9.5))
107   format(i5, ' / ', i5, 5x, 2(i12) )
120   format('total time: ', f8.2, ' minutes')

      end
c----------------------------------------------------------------
c	this subroutine do the iterations
c   in pars:	a
c		b
c		maxpasse
c		limit
c   out pars:	passe
c		rx
c		ry
c

      subroutine henon(a, b, maxpasse, passe, limit, rx, ry)
      double precision a, b, limit
      integer maxpasse, passe
      double precision rx, ry
      
      double precision x, x2, y, y2

      x = 0
      y = 0

      do 5 passe=1, maxpasse

          x2 = 1d0 + y - a * x * x
          y2 = b * x
          x = x2
          y = y2

          if (x .lt. -limit) goto 8
          if (x .gt.  limit) goto 8
          if (y .lt. -limit) goto 8
          if (y .gt.  limit) goto 8

5     continue

8     continue

      rx = x
      ry = y

      end

c----------------------------------------------------------------