rm(list=ls())

N=100 # number of particles
T=50 # final time
theta=2.898 # parameter of interaction
dt=0.00001 # time step
eps=0.0001 # a smoothing parameter : X/|X|^2 replaced by X/(|X|^2+eps)
deltaplot=100 # plots the configuration each "deltaplot" time steps

#########
I=floor(T/dt)
ect=sqrt(dt)
dtdrift=dt*theta/N
nextiplot=0

######### Initial condition
tt=seq(1,N)
X=0.1*rnorm(N)+3*(tt%%2==0)-1.5 
Y=0.1*rnorm(N)+3*(tt%%4==1)+3*(tt%%4==2)-1.5

plot(X,Y,main=paste("N = ",N,", theta = ",theta,", time = ",0*dt),xlab="", ylab="",xlim=c(-3,3), ylim=c(-3,3))
Sys.sleep(1)

for (i in seq(2,I))
{
X=X+ect*rnorm(N)
Y=Y+ect*rnorm(N)
A=cbind(X)%*%rbind(rep(1,N))-cbind(rep(1,N))%*%rbind(X)
B=cbind(Y)%*%rbind(rep(1,N))-cbind(rep(1,N))%*%rbind(Y)
C=A*A+B*B+eps
X=X-rowSums(A/C)*dtdrift
Y=Y-rowSums(B/C)*dtdrift
if (i>nextiplot) 
{
plot(X,Y,main=paste("N = ",N,", theta = ",theta,", time = ",i*dt),xlab="", ylab="",xlim=c(-3,3), ylim=c(-3,3))
nextiplot=nextiplot+deltaplot
}
}