nfish<-5000		#5000 fish...
ndays<-50		#...over 50 days
sizes<-growth<-matrix(0,nrow=ndays,ncol=nfish)	
alive<-array(1,dim=nfish)
sigmag<-0.04		#Standard deviation for each day's growth rate

sizes[1,1:nfish]<-12	#Each fish starts at 12mm
growth[1,1:nfish]<-0.2	#Mean for each day's growth rate

for(i in 1:(ndays-1)){
	meetpred<-runif(nfish)<0.2		#Does fish meet a predator?
	pdie<--0.33+0.15*(90/sizes[i,])		#Probability a fish is consumed if it met a predator
	eaten<-runif(nfish)<pdie		#Fish is consumed
	m<-unique(sort((meetpred*eaten)*(1:nfish)))[-1]	#Stores index of consumed fish
	alive[m]<-0
	sizes[i+1,]<-sizes[i,]+growth[i,]	#New size
	growth[i+1,]<-growth[i,] + sigmag*rnorm(nfish)	#New growth rate
	growth[growth<0]<-0
	growth[growth>0.8]<-0.8
}

m<-unique(sort(alive*(1:nfish)))[-1]
growthrate<-(sizes[ndays,m]-sizes[1,m])/ndays
hist(growthrate)