!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  PROGRAM Mandel;
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  (c) 2005, Christian B. Lang, Mathematische Methoden in der Physik
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Die Mandelbrotmenge wird berechnet. Dabei wird ein
!  Graphikfeld von 200 x 200 Punkten dargestellt.
!  Die maximale Zahl der Iterationen pro Punkt und der
!  Darstellungsbereich wird vorher eingegeben.
!  Empfohlene Werte  Max. Iterationen= 50
!                    Re(c) := -2   ... 1
!                    Im(c) := -1.5 ... 1.5
!  oder              Max. Iterationen= 200
!                    Re(c) := -.28   ... -.24
!                    Im(c) :=  .75 ...    .79     }
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

PROGRAM Mandel

   implicit none

!  Plotgrid in Pixel
   integer, parameter :: xnpix=200,ynpix=200
   
   real ::cra,crb,cia,cib,cr,ci,crdel,cidel
   integer ::sw,n,niter,x,y
   complex :: c,z
   character (len=1):: ans
   character (len=256)::gnucmd,gnudat
   

   write(*,*) 'Ich zeichne die Mandelbrotmenge.'
   write(*,*) 'Maximale Zahl der Iterationen (z.B. 40) : '
   read(*,*) niter
   
1  write(*,*) 'Soll ich auch die Iterationstiefe darstellen (j=0/n=1)?'
   read(*,*) sw
   if((sw.lt.0).or.(sw.gt.1))goto 1

   write(*,*) 'Bereich fuer Re(c) : untere Grenze (z.b. -2)   = '
   read(*,*) cra
   write(*,*) '                     obere  Grenze (z.b.  1)   = '
   read(*,*) crb
   write(*,*) 'Bereich fuer Im(c) : untere Grenze (z.b. -1.5) = '
   read(*,*) cia
   write(*,*) '                     obere  Grenze (z.b. 1.5)  = '
   read(*,*) cib
   
!  Ich verwende gnuplot, um die Daten zu plotten.
!  Wenn gnuplot nicht installiert ist, so koennen die Daten im File
!  folgen_plot.dat' mit einem anderen Plottingprogramm dargestellt
!  werden.
   
   gnudat='mandel_plot.dat'
   gnucmd='mandel_plot.cmd'
   
   print *,'Ich erzeuge zwei Files:'
   print *,trim(gnudat)
   print *,trim(gnucmd)
   print *,'Falls diese Files schon existieren, erfolgt eine Fehlermeldung.'
   open(unit=10,file=trim(gnudat),status='new')
   open(unit=11,file=trim(gnucmd),status='new')
   write(11,'(a)') 'set title "Mandelbrot Menge"'
   write(11,'(a)') 'set xlabel "Re"'
   write(11,'(a)') 'set ylabel "Im"'
   write(11,'(a,f5.1)') 'set size ratio ',abs((cia-cib)/(cra-crb))
   write(11,'(a)') 'plot "'//trim(gnudat)//'" using 1:2 with points pt 4'
   write(11,'(a)') 'pause -1'
   write(11,'(a)') 'quit'
   close(11)
   
!  Nun die Berechnung:
   crdel = (crb - cra) / (xnpix-1)
   cidel = (cib - cia) / (ynpix-1)
   cr = cra - crdel
   do x = 1,xnpix
      cr = cr + crdel
      ci = cia - cidel
      do y = 1,ynpix
         ci = ci + cidel
         call Iteration(cr, ci, niter, n)
         n = niter - n
         if(sw.eq.0)  n = mod(n,2)
         if(n.eq.0) write(10,*)cr,ci
      end do
   end do
   
   close(10)
   
!  Wenn das Unterprogramm "system" nicht verfuegbar ist, so kommentieren
!  Sie die folgende Zeile aus. Dann muessen Sie die Grafik selbst mittels
! "gnuplot mandel_plot.cmd" aufreufen.

   call system(trim('gnuplot '//trim(gnucmd)))

   print *,'Berechnung und Darstellung sind abgeschlossen.'

END PROGRAM Mandel

SUBROUTINE Iteration (cr,ci,niter,n);
!  Konstruiere Folge z_n bis abs(z_n)>4 oder maximale
!  Zahl der Iterationen erreicht wird. 
   implicit none
   
   real, intent(in) ::cr, ci
   complex ::z,c
   integer, intent(in)::niter
   integer, intent(out)::n
   integer::i

   c=cmplx(cr,ci)
   z = 0
   n=0
  
   do while((abs(z).lt.2).and.(n.lt.niter))
     n=n+1
     z=c+z*z
   end do
   
END SUBROUTINE Iteration


