■

なつかしいマンデルブロ集合をば。

PROGRAM Mandel
USE plotter
IMPLICIT NONE
INTEGER, PARAMETER :: kd = SELECTED_REAL_KIND(15)
INTEGER, PARAMETER :: m = 1000, nwinx = 800 , nwiny = 600
INTEGER :: i, j, imax, jmax, maxiter, icount
REAL (kd) :: xmin, xmax, ymin, ymax 
REAL (kd) :: xmin1, xmax1, ymin1, ymax1 
REAL (kd) :: x, y, a, b, a2, b2, dx, dy
INTEGER  :: icol(0:m + 1)
xmin =  -1.0
xmax =   2.2
ymin =  -1.2
ymax =   1.2
maxiter = 100
!
dx = xmax - xmin
dy = ymax - ymin
IF (dx <= 0.0_kd .OR. dy <= 0.0_kd .OR. maxiter <= 0 .OR. maxiter > M) STOP 'input error'
IF (dx * nwinx > dy * nwiny) THEN
 imax = nwinx
 jmax = INT(nwinx * dy / dx + 0.5_kd)
ELSE
 imax = INT(nwiny * dx / dy + 0.5_kd)
 jmax = INT(nwiny)
END IF
dx = dx / REAL(imax, kd)
dy = dy / REAL(jmax, kd)
icol(0) = 0 ! black
j = irgb(255, 255, 255)
icol(0) = 0
DO i = maxiter, 1, -1
 icol(i) = j 
 IF (j > 1) j = j - irgb(255, 255, 255) / maxiter
END DO
!
CALL gr_on('Mandelbrot', nwinx, nwiny)
CALL gr_show()
DO i = 0, imax
 x = xmin + i * dx
 DO j = 0, jmax
  y = ymax - j * dy
  a = x
  b = y
  a2 = a * a
  b2 = b * b
  icount = maxiter
  DO WHILE (a2 + b2 <= 4.0_kd .AND. icount > 0) 
   b = 2.0_kd * a * b - y
   a = a2 - b2 - x
   a2 = a * a
   b2 = b * b
   icount = icount - 1
  END DO
  CALL gr_dot(i, j, icol(icount))
 END DO
 CALL gr_show()
 CALL sleep(5)
END DO
CALL gr_off()
STOP
END PROGRAM Mandel