PROGRAM RandomWalk

*     RandomWalk auf einem quadratischen Gitter
****
 
      REAL R(97)
      INTEGER HOP,IBIT
      EXTERNAL RAN
 
      PARAMETER (NBIG=100000)

      REAL XS(NBIG),YS(NBIG)
      

  

C     Erzeugen von Zufallszahlen 0,1,2,3
      IBIT(X) = INT(4.*X)


      OPEN (10,FILE="RandomWalk.dat",
     &            FORM="FORMATTED",STATUS="UNKNOWN")
      OPEN (15,FILE="RandomWalk_t.dat",
     &            FORM="FORMATTED",STATUS="UNKNOWN")


      A=3.0

      PRINT *, 
     &     'Eingabe Anzahl der Zufallsschritte , Anfangszustand ISEED'
      PRINT *, 
     &     ' (2 natuerliche Zahlen):'
      READ *,M,ISEED

C     Startpunkte
      X=0.*A
      Y=0.*A

      WRITE (*,*) 'Nach Vollendung des Programmes kan man sich die'
      WRITE (*,*) 'berechnete Trajektorie mit'
      WRITE (*,*) 'xmgrace RandomWalk.dat anschauen'
      WRITE (*,*) '(auf Gross- und Kleinschreibung achten)'

      
      WRITE (10,*) 0,X,Y

C     RAN Zufallszahl zwischen 0 und 1
      DO 100, I=1,M
          HOP=IBIT(RAN(ISEED,IX1,IX2,IX3,R))
          IF (HOP.EQ.0) THEN
             X=X+1 
          ELSEIF (HOP.EQ.1) THEN
             X=X-1 
          ELSEIF (HOP.EQ.2) THEN
             Y=Y+1 
          ELSEIF (HOP.EQ.3) THEN
             Y=Y-1 
          ELSE
             PRINT *, 'Kein Sprung'
          ENDIF
           
          XS(I)=X*A       
          YS(I)=Y*A       
          WRITE (10,*) XS(I),YS(I)
          WRITE (15,*) I,XS(I),YS(I)

 100  CONTINUE
      


      END       


*
*
      REAL FUNCTION RAN(IDUM,IX1,IX2,IX3,R)

*     Returns a random deviate between 0.0 and 1.0. Set IDUM to
*     any negative value to initialize or reinitialize the sequence.
*     Numerical Recipes p.196

      REAL R(97)
      INTEGER  IX1,IX2,IX3,IC1,IC2,IC3,IA1,IA2,IA3,IFF
      PARAMETER (M1=259200,IA1=7141,IC1=54733,RM1=1./M1)
      PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=1./M2)
      PARAMETER (M3=243000,IA3=4561,IC3=51349)

      DATA IFF /0/

      IF (IDUM.LT.0.OR.IFF.EQ.0) THEN
         IFF=1
         IX1=MOD(IC1-IDUM,M1)
         IX1=MOD(IA1*IX1+IC1,M1)
         IX2=MOD(IX1,M2)
         IX1=MOD(IA1*IX1+IC1,M1)
         IX3=MOD(IX1,M3)
         DO 11, J=1,97
           IX1=MOD(IA1*IX1+IC1,M1)
           IX2=MOD(IA2*IX2+IC2,M2)
           R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
 11      CONTINUE
         IDUM=1

       ENDIF

       IX1=MOD(IA1*IX1+IC1,M1)
       IX2=MOD(IA2*IX2+IC2,M2)
       IX3=MOD(IA3*IX3+IC3,M3)
       
       J=1+(97*IX3)/M3
       
*       IF (J.GT.97.OR.J.LT.1) PAUSE

       RAN=R(J)  


       R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1

       RETURN
       END
*
*