PROGRAM VICSEK
! T VICSEK, PRL, 75, 6, 1226-1229, (1995)
! TO RUN THE PROGRAM:
! gfortran vicsek.f95; ./a.out
IMPLICIT NONE
INTEGER, PARAMETER :: N = 300
REAL*8 PI,LENGTH,LX,LY,V0,VX,VY,ETA,RCUT
REAL*8 DX,DY,R,ARG,AVGTHETA,NOISE,SUMX,SUMY,AVGVEL
REAL*8, DIMENSION (N) :: X,Y,THETA,NUM,DEN
INTEGER I,J,TIME,TMAX,TIMESTEP
INTEGER,PARAMETER :: SEED = 99999999
CALL SRAND(SEED)
PI = 4.D0*ATAN(1.D0)
OPEN(UNIT=5,FILE='System_information.dat',STATUS='UNKNOWN')
! USER INPUTS.
LENGTH = 5.D0
LX = LENGTH/2.D0
LY = LX
V0 = 0.03D0
ETA = 0.1D0
RCUT = 1.0D0
TMAX = 100
TIMESTEP = 1
! INITIAL POSITION & VELOCITY OF THE PARTICLES.
DO I = 1,N,1
X(I) = 2.D0*(RAND()-0.5D0)*LX
Y(I) = 2.D0*(RAND()-0.5D0)*LY
THETA(I) = 2.D0*RAND()*PI
VX = V0*DCOS(THETA(I))
VY = V0*DSIN(THETA(I))
WRITE(10,*) X(I),Y(I),THETA(I),VX,VY
ENDDO ! I
! TIME ITERATION
AVGVEL = 0.D0
DO TIME = 1,TMAX,TIMESTEP
SUMX = 0.D0
SUMY = 0.D0
DO I = 1,N,1
! INCLUDE Ith PARTICLE IN ARCTAN(THETA)
NUM(I) = DSIN(THETA(I))
DEN(I) = DCOS(THETA(I))
DO J = 1,N,1
! EVALUATE ARCTAN(THETA) OVER THE CIRCLE OF RADIUS RCUT
! MINIMUM IMAGE CONVENTION INTRODUCED
DX = (X(I)-X(J)) - NINT((X(I)-X(J))/(2.D0*LX))*2.D0*LX
DY = (Y(I)-Y(J)) - NINT((Y(I)-Y(J))/(2.D0*LY))*2.D0*LY
R = DSQRT(DX**2 + DY**2)
IF (R <= RCUT) THEN
NUM(I) = NUM(I) + DSIN(THETA(J))
DEN(I) = DEN(I) + DCOS(THETA(J))
ENDIF
ENDDO !J
! EVALUATE THETA
ARG = ABS(NUM(I)/DEN(I))
IF (NUM(I) >= 0.D0 .AND. DEN(I) >= 0.D0) THEN
AVGTHETA = DATAN(ARG)
ELSEIF (NUM(I) >= 0.D0 .AND. DEN(I) <= 0.D0) THEN
AVGTHETA = PI - DATAN(ARG)
ELSEIF (NUM(I) <= 0.D0 .AND. DEN(I) <= 0.D0) THEN
AVGTHETA = PI + DATAN(ARG)
ELSEIF (NUM(I) <= 0.D0 .AND. DEN(I) >= 0.D0) THEN
AVGTHETA = 2.D0*PI - DATAN(ARG)
ENDIF
! NOISE OR TEMPERATURE
NOISE = ETA*(RAND()-0.5D0)
! EVALUATE FINAL THETA
THETA(I) = AVGTHETA + NOISE
! TIME UPDATE
VX = V0*DCOS(THETA(I))
VY = V0*DSIN(THETA(I))
X(I) = X(I) + VX*TIMESTEP
Y(I) = Y(I) + VY*TIMESTEP
! PERIODIC BOUNDARY CONDITION
X(I) = X(I) - 2.D0*(INT(X(I)/LX))*LX
Y(I) = Y(I) - 2.D0*(INT(Y(I)/LY))*LY
ENDDO ! I
DO I = 1,N,1
! EVALUATE AVERAGED NORMALISED VELOCITY
SUMX = SUMX + V0*DCOS(THETA(I))
SUMY = SUMY + V0*DSIN(THETA(I))
WRITE(TIME+100,*) X(I),Y(I),THETA(I),VX,VY
ENDDO
CLOSE (TIME+100)
AVGVEL = AVGVEL + DSQRT(SUMX*SUMX+SUMY*SUMY)/(V0*DFLOAT(N))
ENDDO ! TIME
! EVALUATE ORDER PARAMETER
AVGVEL = AVGVEL/DFLOAT(TMAX)
WRITE(5,*) "Length of box=", LENGTH, "Noise/Temperature (Eta)=", ETA, "r=", RCUT, "Density (Rho)=", LX*LY/DFLOAT(N)
WRITE(5,*) "Average Normalised Velocity (V_a)=", AVGVEL
WRITE(5,*) "Total Timestep=", TMAX
CLOSE(5)
END PROGRAM VICSEK