archival purpose only, redo from start

Fraktalism/map_henon.f

@@ -0,0 +1,149 @@
+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----------------------------------------------------------------