!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  PROGRAM LogMap
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  (c) 2005, Christian B. Lang, Mathematische Methoden in der Physik
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Die Folgenglieder S_201 bis S_400 der nichtlinearen
!  Abbildung  a_n = alpha * a_(n-1) * (1-a_(n-1))
!  werden durch Punkte fuer verschiedene Werte des Parameters
!  alpha dargestellt. Die sich ergebende Struktur wurde von
!  M. Feigenbaum genau untersucht.
!  Der gewuenschte Bereich von alpha kann angegeben werden,
!  um so detailliertere Aufloesung der beobachteten Struktur
!  zu erlauben.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

PROGRAM LogMap

   implicit none
   real :: a1,a2,alpha,adel,an
   integer::i,j,n
   character (len=20)::str1,str2
   character (len=256)::gnudat,gnucmd
   
   print *,'Ich berechne die Folge a_n=alpha*a_(n-1)*(1-a_(n-1))'
   print *,'(fuer a_0=0.5) und stelle die Werte der Folgenglieder'
   print *,'a_201 bis a_400 durch Punkte fuer gegebene Werte des'
   print *,'Parameter alpha dar.'
   print *,'Bitte geben Sie den gewuenschten Wertebereich'
   print *,'von alpha ( >0, <4 ) an (Beispiel: 2.8, 3.99) :'
   print *,'alpha_min = '
   read(*,*) a1
   print *,'alpha_max = '
   read(*,*) a2
   
!  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='logmap_plot.dat'
   gnucmd='logmap_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,f10.6,a,f10.6,a)') 'set title "Feigenbaum-Plot fuer alpha von ',a1,' bis ',a2,'"'
   write(11,'(a)') 'set xlabel "alpha"'
   write(11,'(a)') 'set ylabel "x"'
   write(11,'(a)') 'plot "'//trim(gnudat)//'" using 1:2 with points pt 6'
   write(11,'(a)') 'pause -1'
   write(11,'(a)') 'quit'
   close(11)
   
   adel = (a2 - a1) / 400
   alpha = a1 - adel
   do n = 1,401
     alpha = alpha + adel ;
     an = 0.5 ;
!    Zuerst 200 Voriterationen
     do i = 1,200
        an = alpha * an * (1 - an);
     end do
!    Nun a_201 bis a_400
     do i = 201,400
        an = alpha * an * (1 - an);
        write(10,*)alpha, an
     end do
   end do
   
!  Wenn das Unterprogramm "system" nicht verfuegbar ist, so kommentieren
!  Sie die folgende Zeile aus. Dann muessen Sie die Grafik selbst mittels
!  "gnuplot logmap_plot.cmd" aufreufen.

   call system(trim('gnuplot '//trim(gnucmd)))
   print *,'Ende des Programms'

END PROGRAM logmap

