#COARSE FRACTION PUMP SIMULATION
	
	a=c(.9,.1); b=.12; #
	index=1:50
	B=.8; C=.25
	x=a[1]* B*exp(-C*index) /sum(exp(-C*index) )
		#initialization: this gets worked through in the iteration
  	x=c(-sum(x),x)
	X=cbind(rep(.9,200),rep(.1,200))
	X[2:52,]=cbind(X[2:52,1]-x,X[2:52,2]+x)
	X=rep(list(X),1000);

   for( i in 2:1000) {
	  x=X[[i-1]][1,1]* B*exp(-C*index) /sum(exp(-C*index) )  
		#this moves fines down in a neg exponential snakes and ladders
	  x=c(x,rep(0,198-length(x)))	 
	  X[[i]][,1]= c(a[1],X[[i-1]][1,1]-sum(x),X[[i-1]][2:199,1]+x)  
	  X[[i]][,2]= 1-X[[i]][,1]	
		#this percolates coarse fraction back up for a material balance
	  }
	X=X[501:1000]
	
	y=apply( array(X[[500]][,2],dim=c(10,20)),2,mean);
		#this is an average of 10 years
	ts.plot(y,xlim=c(0,10),ylim=c(0,.45)) 
#Run your code, then make this figure.

x11();par(mfrow=c(3,3))
for(i in c(1,5,10,20,50,100,150,200,500)){
plot(X[[i]][,2], main=i)
}