68 lines
1.5 KiB
Fortran
68 lines
1.5 KiB
Fortran
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
|
|
|