"seasonal.logistic"<-
function(r=0.2,a=0.5,K=50,omega=2*pi/25,nsteps=1000,grid=10000,initial.X=10,iseed=123)
{
	#
	# Default values come from Renshaw, Modelling Biological Populations in Space and Time, p. 227
	#
	#
	# Set the random number seed for iseed in [0,1023].
	set.seed(iseed)
	#
	#  
	# 
	X <- rep(0, nsteps)
	sojourn <- X
	W <- X
	X[1] <- initial.X
	W[1] <- sojourn[1]
	#
	#
	# Create a grid of time points covering over 99.99% of the range of Exponential(1)
	#
	tt<-(0:grid)*10/grid
	#
	#
	# Note that this is time = 0.
	#
	for(i in 2:nsteps) {
		# Model has a non-seasonal birth rate r*X and seasonal death rate (r*X^2/K)*(1+a*cos(omega*t)),
		# so it exhibits logistic growth with seasonal fluctuations.  Need to integrate the total
		# birth+death rate and apply the inverse of this integral to a standard exponential.
		#
		t0<-W[i-1]
		ff<-a*(sin(omega*(tt+t0))-sin(omega*t0))
		B<-r*X[i-1]
		tmp<-r*X[i-1]^2/K
		Lambda<-(B+tmp)*tt+tmp*ff
		#
		# Draw standard exponential to be transformed into sojourn time.
		#
		S.star<-rexp(1,1)		
		#
		# Evaluate inverse of Lambda at S.star.
		#
		index<-grid+2-length(tt[Lambda>S.star])
		sojourn[i]<-tt[index]
		birth<-B
		death<-tmp*(1+a*cos(omega*tt[index]))
		#
		# In case the sojourn time was even longer than anticipated...
		#
		if(is.na(sojourn[i])){
			sojourn[i]<-10;
			birth<-B[grid+1];
			death<-tmp[grid+1]*(1+a*cos(omega*tt[grid+1]))
						}
		#
		# Note that the transition probabilities depend on the instantaneous rates,
		# not the integrated rates.
		#
		rates<-c(birth,death)	
		#
		# Make sure all rates are non-negative.
		#
		rates[rates < 0] <- 0
		#
		# Check to see if the population is extinct. 
		#
		overall.rate <- sum(rates)
		if(overall.rate <= 0 || X[i-1] == 0) {
			# population is extinct
			sojourn[i] <- 0
			W[i] <- W[i - 1]
		}
		#end if extinct
		if(overall.rate > 0) {
			W[i] <- W[i - 1] + sojourn[i]
			#
			# Determine the probabilities of each state transition and 
			# draw a random state transition. 
			#
			st <- sample(1:2, size = 1, prob = rates/overall.rate)
			if(st == 1) {
				# birth
				X[i] <- X[i - 1] + 1
			}
			if(st == 2) {
				# death
				X[i] <- X[i - 1] - 1
			}
		}
	}
	# end for loop on i
	plot(c(0,W),c(0,X), type = "n", xlab = "Time", ylab = "Population Size")
	lines(W,X,type="s")
	grid<-(0:1000)*max(W)/1000
	top<-max(X)/5
	f<-1+a*cos(omega*grid)
	f<-(f-min(f))/max(f)
	#lines(grid,top*f,col=8)
	ncycles<-(max(W)*omega/(2*pi))%/%1
	abline(v=(0:ncycles)*2*pi/omega,col=8)
		return()
}