# Project sizes of population growing deterministically at per capita finite rate R = 1 + r per time step
# By R. Gomulkiewicz (8 Sep 2015)
# Note: everything to the right of # on a line is not interpreted in R, i.e., is regarded as a "comment"

initial.size <- 10 #starting population size
steps <- 50 
r <- -0.01 #intrinsic rate of growth or decline

# This is a function that computes the population size at the next time step assuming each 
# individual has R = 1+r descendants
next.size <- function(n, r) {n*(1+r)}
		

# create a list to hold the population sizes
# "populate" with a bunch of zeros
n <- matrix(0,steps+1,1)

#store the initial population size in the first position of the list (zeros everywhere else)
n[1] <- initial.size

#this is a "loop" that projects the population sizes for "steps"
for(i in 1:steps){ n[i+1] = next.size(n[i],r)}

#plot the results
plot(0:steps,n)

################################################
#Compare projections for 3 different values of r
#create a (steps+1)x3 table to hold the results
nr <- matrix(0,steps+1,3)

#store the initial population size in the first position (row) of each column
nr[1,] <- initial.size

#loop that projects the population sizes for r = 0.01
# results will be stored in the 1st column
for(i in 1:steps){ nr[i+1,1] = next.size(nr[i,1],0.01)}

#loop that projects the population sizes for r = 0.0
# results will be stored in the 2nd column
for(i in 1:steps){ nr[i+1,2] = next.size(nr[i,2],0.00)}

#loop that projects the population sizes for r = -0.01
# results will be stored in the 3rd column
for(i in 1:steps){ nr[i+1,3] = next.size(nr[i,3],-0.01)}


#plot all three trajectories together
matplot(0:steps,nr)

#line plot of the same
matplot(0:steps,nr,type="l")

#point & line plot of the same, with x and y axis labels
matplot(0:steps,nr,type="b",xlab="step",ylab="population size")


######################################################
#plot relative frequency of "clone 1" versus "clone 3"
freq <- nr[,1]/(nr[,1]+nr[,2])
plot(0:steps,freq,type="l",xlab="step",ylab="fraction of clone 1")