!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  PROGRAM Monte
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  (c) 2005, Christian B. Lang, Mathematische Methoden in der Physik
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Monte-Carlo Integration. In diesem Programm werden zufaellig Punkte
!  in einer quadratischen Flaeche (1x1) placiert. Die Zahl der Punkte,
!  die ins Innere einer Kreisflaeche (r=0.5) treffen, werden gezaehlt.
!  Das Verhaeltnis Treffer/Gesamtzahl konvergiert gegen das Verhaeltnis
!  Kreisflaeche/Quadratflaeche (=Pi/4=0.78539816...), also die Kreisflaeche!
!  In den mit * auskommentierten Zeilen ist beschrieben, wie man  mit
!  diesem Programm auch eine Massendichteverteilung der Kreisscheibe
!  einbauen kann und so die Masse der Scheibe als Ergebnis erhalten
!  (exakter Wert= Pi*ln(1.25) = 0.70102614...).   
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


PROGRAM Monte

   implicit none
   
   real :: x,y,Versuche,Treffer,Schaetzwert,rho
   integer :: samples,i
   character (len=256)::gnucmd,gnudat1,gnudat2
   
   print *,'Wieviele Punkte soll ich berechnen?'
   read(*,*) samples
   
   gnudat1='monte_plot1.dat'
   gnudat2='monte_plot2.dat'
   gnucmd='monte_plot.cmd'
   
   print *,'Ich erzeuge drei Files:'
   print *,trim(gnudat1)
   print *,trim(gnudat2)
   print *,trim(gnucmd)
   print *,'Falls diese Files schon existieren, erfolgt eine Fehlermeldung.'
   open(unit=9 ,file=trim(gnudat1),status='new')
   open(unit=10,file=trim(gnudat2),status='new')
   open(unit=11,file=trim(gnucmd),status='new')
   write(11,'(a)') 'set title "Monte Carlo"'
   write(11,'(a)') 'set xlabel "x"'
   write(11,'(a)') 'set ylabel "y"'
   write(11,'(a)') 'set size ratio 1'
   write(11,'(a)') 'plot "'//trim(gnudat1)//'" using 1:2 with points pt 4 , "'&
   &                      //trim(gnudat2)//'" using 1:2 with points pt 1'
   write(11,'(a)') 'pause -1'
   write(11,'(a)') 'quit'
   close(11)
   
   Versuche = 0
   Treffer = 0
   
   do i=1,samples
     Versuche  = Versuche + 1
!    Zufallskoordinaten im Quadrat:
     call random_number(x)
     x  = x - 0.5
     call random_number(y)
     y  = y - 0.5
   
!    Treffer im Kreis mit r =0.5 ?

!    Im Falle einer Massendichte: Ersetze letzen Befehl durch:
!    if(x*x + y*y < 0.25) Treffer  = Treffer + rho(x, y)
   
     if(x*x + y*y < 0.25) then
       Treffer  = Treffer + 1
       write(9,*)x,y
     else
      write(10,*)x,y
     endif
   
   end do
   
   Schaetzwert  = real(Treffer)/ real(Versuche)
   
   close(9)
   close(10)
   
   write(*,*) 'Schaetzwert bei ',Versuche,' Versuchen  = ',Schaetzwert
!  Wenn das Unterprogramm "system" nicht verfuegbar ist, so kommentieren
!  Sie die folgende Zeile aus. Dann muessen Sie die Grafik selbst mittels
!  "gnuplot monte_plot.cmd" aufreufen.
   
   call system(trim('gnuplot '//trim(gnucmd)))
   
   print *,'Berechnung und Darstellung sind abgeschlossen.'

END PROGRAM Monte

REAL FUNCTION rho(x, y)

   implicit none
   real :: x,y

!  Bei Bedarf: Massendichte der Kreisscheibe. z.B. :

  rho = 1/(1 + x*x + y*y)
   
END FUNCTION rho

