(re-)starting work on Henon mapping

Fraktalism/.gitignore

@@ -4,6 +4,7 @@ pickover
 lorentz
 voxelize
 evolvopick
+henon
 
 WS/*.dat
 WS/*.txt

Fraktalism/Makefile

@@ -18,6 +18,9 @@ OBJS  = mods/spitpgm.o mods/points3d.o fraktals.o
 
 # ---------------------------------------------
 
+henon:		henon.f90 Makefile $(OBJS)
+	gfortran $(GFOPT) $<  $(OBJS)    -o $@
+
 julia:		julia.f90 Makefile $(OBJS)
 	gfortran $(GFOPT) $<  $(OBJS)    -o $@
 

Fraktalism/henon.f90

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