Fortraneries/Fraktalism/map_henon.f

150 lines
3.6 KiB
FortranFixed
Raw Permalink Normal View History

2022-04-07 04:22:32 +11:00
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----------------------------------------------------------------