(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/common.sh

@@ -1,4 +1,4 @@
 
 
-POVOPT=" -d +q9 +a +W1024 +H768 -v +WT4"
+POVOPT=" -d +q9 -a +W1920 +H1080 -v +WT4"
 

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
+

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----------------------------------------------------------------

Fraktalism/mkvoxvidz.sh

@@ -12,34 +12,36 @@ fi
 
 ./voxelize > WS/voxels.dat
 ./vox2inc.awk < WS/voxels.dat > WS/voxels.inc
+grep NbrVox WS/voxels.inc
 
 TMPNG="/dev/shm/voxvidz.png"
 NBIMG=360
+
 for idx in $( seq   0 $(( NBIMG - 1)) )
 do
 
 	dst=$(printf "frames/voxel/%05d.png" $idx)
-	delta=$( echo "scale=3 ; ${idx}/360" | bc -l)
-	echo "Renderbox work on " $dst "        tu"  $delta
+	delta=$( echo "scale=6 ; ${idx}/360" | bc -l)
+	echo "Renderbox work on "$dst"     delta = "  $delta
 
 	povray		-ishowvoxels.pov -K$idx ${POVOPT}	\
 			-O${TMPNG} 2> WS/toto.err
 	if [ $? -ne 0 ] ; then
-		echo "ERROR ERROR ERROR ERROR"
+		echo "ERROR ERROR ERROR ERROR ERROR ERROR"
 		tail -15 WS/toto.err
 		sleep 20
 		#exit 1
 	fi
 
-	titre='>>> Voxelisation - tTh - avril 2022'
+	titre='Voxelisation - tTh - avril 2022'
 	numbers=$(tail -1 WS/camvox.log |		\
-		  awk '{printf ">>>   %6.3f   %6.3f   %6.3f", $1, $2, $3}')
-	echo $numbers
+		  awk '{printf "  %6.3f   %6.3f   %6.3f", $1, $2, $3}')
+	echo "$numbers"
 
 	convert		${TMPNG}			\
-			-fill  Cyan			\
+			-fill  Gray70			\
 			-kerning 1			\
-			-pointsize 16			\
+			-pointsize 24			\
 			-font AvantGarde-Book		\
 			-gravity South-West		\
 			-annotate +20+32 "$titre"	\
@@ -51,4 +53,4 @@ do
 
 done
 
-./encode.sh frames/voxel/ voxel-1.mp4
+./encode.sh frames/voxel/ voxel-2.mp4

Fraktalism/showvoxels.pov

@@ -17,7 +17,7 @@ global_settings {
 #declare VOXEL = object
 {
 // sphere { 0, 1.18 }
-#local D = 1.38;
+#local D = 1.0111;
 box { <-D, -D, -D>, <D, D, D> }
 }
 
@@ -34,18 +34,18 @@ object {
 
 //----------------------------------------------------------------
 // light_source { <-12,  45, -25> color Gray70 }
-light_source { <-52,  5, -38>   color Yellow*0.55 	}
+light_source { <-52,  5, -38>   color Yellow*0.45 	}
 light_source { < 59,  5,  48>   color Gray20		}
-light_source { <  3, -59,  8>   color Red*0.75 		}
-light_source { <  8,  61,  3>   color Green*0.85 	}
+light_source { <  3,  59,  8>   color Red*0.65 		}
+light_source { <  8,  -61,  3>   color Green*0.75 	}
 
 #declare NormClock = (clock/360);
 
-#declare ECAM  =  88 - 23*NormClock;
+#declare ECAM  =  88 - (73*NormClock);
 #declare XCAM  =  ECAM * sin(radians(clock));
 #declare YCAM  =  -12;
 #declare ZCAM  =  ECAM * cos(radians(clock));
-#declare ACAM  =  65 + (27*NormClock);
+#declare ACAM  =  65 + (48*NormClock);
 
 #if (0 = clock)
   #fopen CL "WS/camvox.log" write

Fraktalism/vox2inc.awk

@@ -9,7 +9,7 @@ BEGIN	{
 	print "union {"
 	}
 
-$4 > 145 {
+$4 > 500 {
 	count = $4
 	value = $5
 	if (count > maxcount) 

Fraktalism/voxelize.f90

@@ -7,23 +7,25 @@
 program voxelize
   use fraktals
 
+  implicit none
+
   integer, parameter                          :: DIM = 100
   integer, dimension(:,:,:), allocatable      :: cube
   type(t_point3d), dimension(:), allocatable  :: points
   integer                                     :: errcode, foo
   integer                                     :: ix, iy, iz
-  double precision, dimension(4)              :: coefs
+  integer                                     :: nbr_points
+  double precision, dimension(4)              :: coefs, KB
   double precision                            :: dval
 
-  ! XXX foo = (DIM*DIM*DIM) / (1024)
-  ! XXX PRINT *, "memory request for cube (in Kwords) ", foo
+  write(0, *) "--- start of voxelize"
 
   allocate (cube(DIM,DIM,DIM), stat=errcode)
   if (0 .NE. errcode) then
     STOP " : NO ENOUGH MEMORY FOR CUBE"
   endif
 
-  nbr_points = 3100000
+  nbr_points = 2800000
   allocate(points(nbr_points), stat=errcode)
   if (0 .NE. errcode) then
     STOP " : NO ENOUGH MEMORY FOR POINTS"
@@ -31,7 +33,10 @@ program voxelize
 
   coefs(1) =  2.23     ;     coefs(2) =  0.42
   coefs(3) = -0.64     ;     coefs(4) = -2.42
-  call compute_pickover(points, coefs)
+
+  KB(1) =  1.51        ;     KB(2) = -1.89
+  KB(3) =  1.69        ;     KB(4) =  0.79
+  call compute_pickover(points, KB)
 
   call clear_cube(cube)
 
@@ -46,7 +51,7 @@ program voxelize
   enddo
 
   dval = DBLE(MAXVAL(cube))
-  write(0, *) "--- cube maximum = ", dval
+  write(0, *) "--- maximum of the cube= ", dval
 
   do foo=1, nbr_points
     call fcoor2icoor(points(foo)%x, ix)
@@ -57,6 +62,9 @@ program voxelize
              DBLE(cube(ix,iy,iz)) / dval
   enddo
 
+  dval = DBLE(MAXVAL(cube))
+  write(0, *) "--- end of voxelize"
+
 !-----------------------------------------------------
 contains
 !-----------------------------------------------------
@@ -76,7 +84,7 @@ subroutine fcoor2icoor(in, out)
   if (outvalue .GE. DIM)    out = DIM-1
 
 end subroutine
-!-----------------------------------------------------
+!------------------------------------------------------------
 subroutine clear_cube(cube)
   type(integer), dimension(:,:,:), intent(out)  :: cube
 
@@ -91,7 +99,7 @@ subroutine clear_cube(cube)
   enddo
 
 end subroutine
-!-----------------------------------------------------
+!------------------------------------------------------------
 subroutine spit_cube_as_union(fname, datas, notused)
   character(*)                                 :: fname
   type(integer), dimension(:,:,:), intent(in)  :: datas
@@ -106,6 +114,7 @@ subroutine spit_cube_as_union(fname, datas, notused)
     STOP ' : SPITUNION : FAIL OPEN OUTPUT FILE'
   endif
 
+
   close(fd)
 
 end subroutine