# condensation3.R
#
# Written by Chris McCarthy, Mathematics, Borough of Manhattan Community College, New York NY
#
#######################################
# Function Definitions
#######################################
####################################### 3
moveP3 = function(p,sites, dt = .3, directionChange = .2,  #
maxX = 100, maxY = 100, edge = .2, tol = 2){
           theta = runif(1, 0, 2*pi);
           pout = p;
           if(runif(1,0,1) < directionChange){
              s = sqrt(p[3]^2 + p[4]^2); 
              pout[3] = cos(theta)*s;
              pout[4] = sin(theta)*s;
           }
           pout[1] = p[1] + pout[3]*dt;
           pout[2] = p[2] + pout[4]*dt;
           if( (pout[1] < edge) || (pout[1] > (maxX- edge)) ){
              pout[3] = -pout[3];
              pout[1] = p[1] + pout[3]*dt;
              }
           if( (pout[2] < edge) || (pout[2] > (maxY - edge) )  ){
             pout[4] = -pout[4];
             pout[2] = p[2] + pout[4]*dt;
             }
           if(min(apply(sites, 1, function(s){sum(abs(s- p[1:2] + p[5]))}))<tol)
             {pout[1:2] = p[1:2]; pout[3] = 0; pout[4]=0; pout[5] =1; } 
          return(pout);}
####################################### 3
its = 30;
numParticles = 20;
p1 = c(50, 50, 5, 0, 0) # x, y, vx, vy stuck
pm = p1
for(i in 1:numParticles){
  xcord = sample(1:99, 1);
  ycord = sample(1:99, 1);
  theta = runif(1, 0, 2*pi);
  xvel =  6*cos(theta);
  yvel =  6*sin(theta);
   st = 0;
  pm = rbind(pm,c(xcord, ycord, xvel, yvel, st))  
}
Xsite =  c(1,0);   
for(i in 2:99){ Xsite = rbind(Xsite, c(i,0));}
interval = .01;
nmax =its;
clist = 0; 
ilist = 0;
###############################################################
# To create png image files that you can make into a movie:
###############################################################
# Uncomment the 2 lines starting with setwd and png 
# to produce png's which we can turn into a movie. 
# Change C:\\Users\\CM\\Desktop\\... 
# to the directory where you want to save the png files.
# Also uncomment the last line of this script dev.off().
################################################################  
# setwd("C:\\Users\\CM\\Desktop\\sync\\Spring 2016\\MSEIP Summer 2016\\png2");
# png(filename = "Rplot%03d.png", width = 960, height = 480, units = "px")
################################################################  
for (i in 1:nmax){
    pm = t(apply(pm, 1,function(pin){moveP3(p = pin, sites = Xsite);}))
    ilist = c(ilist,i)
    clist = c(clist,sum(pm[,5]))
    par(mfrow = c(1,2), pch = 19);
    plot(pm[,1], pm[,2], col = "blue", xlim = c(0,100), ylim = c(0,100),
            xlab = "", ylab = ""); 
    plot( ilist, clist, 
          xlim = c(0,nmax), ylim = c(0,numParticles),
          xlab = "time", ylab = "particles removed",  
          col = "red", grid());  
  };
#######################################
# First order model y' = k(numParticles - y) solution: y(x) = numParticles(1- exp(-k*x)
x = ilist; y = clist;
df <- data.frame(x, y);
mCondensationFirstOrd <- nls(y ~ (numParticles*(1-exp(-k*x))), #algorithm = "port",
data = df,
start = list(k = .1),
trace = T);
#s <- seq(from = 0, to = 300, length = 50);
s = ilist;
lines(s, predict(mCondensationFirstOrd, list(x = s)), col = "black", lwd = 3 );
########################################
summary(mCondensationFirstOrd)
########################################
# dev.off() # to allow plots to screen, etc.