program mandel
      implicit none

      integer max, res
      parameter (max=1000, res=255)

      integer i,j,r,start,stop,red,green,blue,n
      character*30 dummy
      character*30 filename

      integer colour(max,max)
      real  zr(max,max), zi(max,max), cr(max,max)
      real  ci(max,max), zrs(max,max), zis(max,max)

      read (*,*) dummy, n
      read (*,*) dummy, filename

      print *, "calculating the mandelbrot set"
      print *, "NxN GRID, with N:", n
      print *, "\n"

C      initialise zr and zi using forall

      do j=1,n
        do i=1,n
	 cr(i,j) = real(i-1)/real(n-1)
	 ci(i,j) = real(j-1)/real(n-1)
	end do
      end do
      	
      
C      forall (i=1:n, j=1:n) cr(i,j) = real(i-1)/real(n-1)
C      forall (i=1:n, j=1:n) ci(i,j) = real(j-1)/real(n-1)

      
C     initialise other arrays

      do j=1,n
        do i=1,n
	zr(i,j) = cr(i,j)
	zi(i,j) = ci(i,j) 
	end do
      end do
            
C      zr = cr
C      zi = ci


C      zrs = zr*zr
C      zis = zi*zi
C      colour = 0 
      
      do j=1,n
        do i=1,n
	zrs(i,j) = zr(i,j)*zr(i,j)
        zis(i,j)= zi(i,j)*zi(i,j)
	colour(i,j) = 0 
	end do
      end do
      
    

C     main loop : resolution number of iterations

C      do r = 0, resolution

       do j=1,n
        do i=1,n 
         r=-1

         do while ( (zrs(i,j) + zis(i,j)) .le. 4.0 .and. r.lt.res)
	
	   zrs(i,j) = zr(i,j)*zr(i,j)
	   zis(i,j) = zi(i,j)*zi(i,j)

	   zi(i,j) = 2.0*zr(i,j)*zi(i,j) + ci(i,j)
	   zr(i,j) = zrs(i,j) - zis(i,j) + cr(i,j)

           r=r+1
	   colour(i,j) = r	 
         end do
C	 print *, j,i, "\t", colour(i,j)
        end do 
       print *, "The ", j, "th row has completed"
       call sleep(1) 
       end do 
C     end do

      print *,"The imagemap is being written to: ", filename


      open(unit=10, file=filename)
      write(10, fmt='(''P3'',/,i3,2x,i3,/,i3)') n, n, res
      
C    some pretty colour map, maybe
      
      do i = 1, n
	 do j = 1, n
	    if (colour(i,j) == 255) then
	       red=0; green=0; blue=0
	    else if( colour(i,j) >= 0 .and. colour(i,j) <= 10 ) then
	       red = colour(i,j) * 25; green = 0; blue = 0
	    else if( colour(i,j) >= 11 .and. colour(i,j) <= 15 ) then
	       red = 255; green = 0; blue = 0
	    else
	       red = res; green = res; blue = res; 
	    end if
	    write(10,*) red, green, blue
	 end do
      end do
     
      close(unit=10)

      print *," Mandelbrot is successfully completed"
      end