ccccccccccccccccccccccccc     Program 7.2     cccccccccccccccccccccccccc
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                      c
c Please Note:                                                         c
c                                                                      c
c (1) This computer program is part of the book, "An Introduction to   c
c     Computational Physics," written by Tao Pang and published and    c
c     copyrighted by Cambridge University Press in 1997.               c
c                                                                      c
c (2) No warranties, express or implied, are made for this program.    c
c                                                                      c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      SUBROUTINE MXWL(N,M,T,V)
C
C This subroutine assigns velocities according to the Maxwell
C distribution. N is the total number of velocity components
C and M is the total number of degrees of freedom. T is the
C system temperature in the reduced units.
C
      DIMENSION V(N)
      COMMON /CSEED/ISEED
C
C Assign a Gaussian distribution to each velocity component
C
      DO 100   I = 1, N-1, 2
        CALL GRNF(V1,V2)
        V(I)   = V1
        V(I+1) = V2
  100 CONTINUE
C
C Scale the velocity to satisfy the partition theorem
C
      EK = 0.0
      DO 200   I = 1, N
        EK = EK+V(I)*V(I)
  200 CONTINUE
      VS = SQRT(EK/(M*T))
      DO 300   I = 1, N
        V(I) = V(I)/VS
  300 CONTINUE
      RETURN
      END
C
      SUBROUTINE GRNF(X,Y)
      COMMON /CSEED/ ISEED
        PI = 4.0*ATAN(1.0)
        R1 = -ALOG(1.0-RANF())
        R2 = 2.0*PI*RANF()
        R1 = SQRT(2.0*R1)
        X  = R1*COS(R2)
        Y  = R1*SIN(R2)
      RETURN
      END
C
      FUNCTION RANF()
      DATA IA/16807/,IC/2147483647/,IQ/127773/,IR/2836/
      COMMON /CSEED/ ISEED
        IH = ISEED/IQ
        IL = MOD(ISEED,IQ)
        IT = IA*IL-IR*IH
        IF(IT.GT.0) THEN
          ISEED = IT
        ELSE
          ISEED = IC+IT
        END IF
        RANF = ISEED/FLOAT(IC)
      RETURN
      END