150 lines
3.6 KiB
Fortran
150 lines
3.6 KiB
Fortran
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----------------------------------------------------------------
|