        PROGRAM DIFFUS
C
C       --> Kompiliere mit "g77 diffusion.f -o diffusion"
C       Dabei ist /usr/bin/g77 ein f77-Compiler von "GNU".
C
C       Erstellt eine Histogramm fuer die Diffusion. Ein Teilchen wird
C       an der Stelle X = 0 mit einer Brownschen Bewegung, erzeugt durch
C       Zufallszahlen, in Bewegung "gesetzt". Dabei ist
C              Xneu = Xalt +- 1
C       und nur die Vorzeichen werden "ausgewuerfelt".
C       Die Diffusionzeit wird durch Nstep festgelegt und
C       dieser Prozess Ndiff wiederholt.
C
C       Fuer das Histogramm muesste man folgendes bekommen:
C          Hist(x)  proportional zu  exp[-x**2/(2*Nstep)]

        IMPLICIT NONE
        INTEGER NMAX, NSTEP, NDIFF
        PARAMETER(NMAX=40)
        INTEGER HISTO(-NMAX:NMAX), X, DI, N, NS
        REAL DX, RAND, SIG

C       ---Initialisierung--->
        DATA HISTO/NMAX*0,0,NMAX*0/

        WRITE(*,*) ' Wieviele Zeitschritte? [z.B 64]'
        READ (*,*) NSTEP
        WRITE(*,*) ' Wieviele Wiederholungen? [z.B 8192]'
        READ (*,*) NDIFF

C       ---eigentliche Rechnung--->
        DO N =  1, NDIFF
         X = 0
C        Einzelne Brownsche Bewegung, Beginn bei X = 0
         DO NS = 1, NSTEP
          DX = RAND() + 0.5 
          DI = INT(DX)
          DI = DI + DI - 1
          X = X + DI
         END DO
C        Ende der einzelne Brownsche Bewegung:
C        Position des Partikels nach Nstep bei X
          IF(IABS(X).LE.NMAX) THEN
           HISTO(X) = HISTO(X) + 1
          END IF
        END DO

C       Resultate ---------->
C       Eine mit "Gnuplot" lesbare Datei fort.11 wird erzeugt.
C       Informationen, wie man damit umgeht:
        WRITE(*,*) ' Histogramm --> fort.11\n'
        WRITE(*,*) ' "gnuplot" eingeben und dann folgendes schreiben:'
        WRITE(*,*) '>plot "fort.11" \n'
        WRITE(*,*) ' und eine Funktion zum Fitten definieren:'
        WRITE(*,*) '>f(x) = a*exp(-x**2/b)\n'
        WRITE(*,*) '>fit [-25:25] f(x) "fort.11" via a, b'
        WRITE(*,*) ' Der Fit braucht man Grenzen!\n  Fit anschauen:'
        WRITE(*,*) '>plot [-30:30] "fort.11", f(x) with lines'
        WRITE(*,*) ' Es sollte b = ', 2*NSTEP, ' sein.'
        WRITE(11,'(A39,I6)')'# Histogramm fuer Diffusion, Zeit:',NSTEP
C       Resultate werden nur geschrieben, wenn Histo != 0 ist.
        DO N = -NMAX, NMAX
         IF(HISTO(N).GT.0) THEN
          SIG = HISTO(N)
          SIG = SQRT(SIG)
          WRITE(11,'(2I5,F8.2)') N, HISTO(N), SIG
         END IF
        END DO

        STOP
        END
