MODULE quadrato
   IMPLICIT NONE
!
!  La matrice "figura" rappresenta una tabella quadrata di punti
!
   INTEGER, PARAMETER :: nd = 64 
   INTEGER, DIMENSION(nd,nd) :: figura
END MODULE quadrato

PROGRAM frattale
   USE quadrato
   IMPLICIT NONE
   INTEGER :: esponente, passi, n, i, j, r1, c1, ios
!
!  Si azzerano tutti i punti di figura
!
   figura = 0

!
!  Ad ogni passo si dividera' il quadrato in 4 parti uguali
!
   PRINT*,"Inserisci il numero di passi"
   READ*,passi
   PRINT*," passi = ",passi
!
   n = nd
   r1 = 1
   c1 = 1
!
!  All'inizio l'intero quadrato viene passato a  quadrato_frattale()
!
   CALL quadrato_frattale(r1,c1,n,passi)

   IF ( nd < 80 ) THEN
      DO i = 1, nd
         WRITE(*,'(80i1)') figura(i,:)
      END DO
   ELSE
!     Se il quadrato e' grande, viene prodotta un'immagine pbm, 
!     visualizzabile ad esempio con Paint Shop o Photo Shop 
!     oppure convertibile in altri formati con Image Magick 
!     (http://www.imagemagick.org/)
!
      OPEN(11,FILE="frattale.pbm",IOSTAT=ios)
      IF ( ios /= 0 ) THEN
         PRINT*,"Errore accesso scrittura del file  frattale.pbm"
         STOP
      ENDIF
      WRITE(11,"(A)") "P1"   ! Il codice "magico" per immagini PBM
      WRITE(11,"(2(I6,1X))") nd, nd   ! Dimensioni dell'immagine
      DO i = 1, nd
         write(11,*) figura(i,:)
      END DO
      CLOSE(11)
   ENDIF
!
   STOP
END PROGRAM frattale


RECURSIVE SUBROUTINE quadrato_frattale(r0,c0,n,passi)
   USE quadrato
   IMPLICIT NONE
   INTEGER :: r0, c0, n, passi
   INTEGER :: r1, c1
  
!   PRINT*," in figura n, passi = ",n, n/2+1, passi
!   IF ( n < 2 ) RETURN


!
!  Il quadrato viene diviso in 4 parti, di cui solo 
!  le 3 parti in basso a sinistra vengono passate a 
!  quadrato_disegna()
!
   r1 = r0; c1 = c0
   CALL quadrato_disegna(r1,c1,n/2,passi-1)
!      PRINT*,"figura(",r1,c1,")",n/2,passi-1
      
   r1 = r0 + n/2; c1 = c0
   CALL quadrato_disegna(r1,c1,n/2,passi-1)
!      print*,"figura(",r1,c1,")",n/2,passi-1
      
   r1 = r0 + n/2; c1 = c0 + n/2
   CALL quadrato_disegna(r1,c1,n/2,passi-1)
!   print*,"figura(",r1,c1,")",n/2,passi-1
   
   RETURN
END SUBROUTINE quadrato_frattale

RECURSIVE SUBROUTINE quadrato_disegna(r1,c1,n,passi)
   USE quadrato
   IMPLICIT NONE
   INTEGER :: r1, c1, n, passi
   INTEGER :: i, j
   
!
!  Il quadrato viene "disegnato" effettivamente solo se
!  sufficientemente piccolo, altrimenti viene passato alla
!  quadrato_frattale()
!
!   print*," in disegna r1, c1, n, n/2, passi = ",r1, c1, n, n/2+1, passi
   IF ( n <= 1 .or. passi <= 1 ) THEN
      DO i = 1, n
         DO j = 1, i
	    figura(r1+i-1,c1+j-1) = 1
	 END DO
      END DO
   ELSE
      CALL quadrato_frattale(r1,c1,n,passi)
   END IF

   RETURN
END SUBROUTINE quadrato_disegna