###########################################
# FIRST, WORK WITH CONVEX CURVES          #
# Scenario 1 - Stable population Lambda=1 #
###########################################
L=1 				    # Lambda
c=seq(from=0, to=5,by=.1)   # Current repro value
f=(c+.1)^(-1)		    # Future  repro value
v=c+f/L			    # Lifetime reproductive output

par(mfrow=c(3,2))
plot(c,f,type='l',
	xlab='Fecundity at age x',
	ylab='Residual Repro. Value')
text(2,2,bquote(lambda==.(L)))
lines(c,(max(v)-L*c),lty=2,col='blue')
abline(v=0,lty=2,col='green')
plot(c,v,type='l',
	xlab='Fecundity at age x',
	ylab='Lifetime Reproductive Output')
text(2,4,bquote(lambda==.(L)))
abline(h=max(v),lty=2,col='blue')
abline(v=0,lty=2,col='green')

# Scenario 2 - Lambda>1 = de-value 
#     residual reproductive value	
L=2
c=seq(from=0, to=5,by=.1)
f=(c+.1)^(-1)
v=c+f/L
plot(c,f,type='l',
	xlab='Fecundity at age x',
	ylab='Residual Repro. Value')
abline(v=0,lty=2,col='green')
abline(v=5,lty=2,col='green')
text(2,2,bquote(lambda==.(L)))
lines(c,(max(f)-L*c),lty=2,col='blue')
plot(c,v,type='l',
	xlab='Fecundity at age x',
	ylab='Lifetime Reproductive Output')
text(2,2,bquote(lambda==.(L)))
abline(h=max(v),lty=2,col='blue')
abline(v=0,lty=2,col='green')
abline(v=5,lty=2,col='green')

# Scenario 3 - Lambda<1 = increase 
#     residual reproductive value	
L=0.5
c=seq(from=0, to=5,by=.1)
f=(c+.1)^(-1)
v=c+f/L
plot(c,f,type='l',
	xlab='Fecundity at age x',
	ylab='Residual Repro. Value')
abline(v=0,lty=2,col='green')
text(2,2,bquote(lambda==.(L)))
lines(c,(max(f)-L*c),lty=2,col='blue')
plot(c,v,type='l',
	xlab='Fecundity at age x',
	ylab='Lifetime Reproductive Output')
text(2,7,bquote(lambda==.(L)))
abline(h=max(v),lty=2,col='blue')
abline(v=0,lty=2,col='green')

###########################################
# NEXT, WORK WITH CONCAVE CURVES          #
# Scenario 1 - Stable population Lambda=1 #
###########################################
windows()
L=1 				    # Lambda
c=seq(from=0, to=5,by=.1)   # Current repro value
f=10-.02*c^2-.03*c^3	    # Future  repro value
v=c+f/L			    # Lifetime reproductive output
par(mfrow=c(3,2))
plot(c,f,type='l',
	xlab='Fecundity at age x',
	ylab='Residual Repro. Value')
text(1,7,bquote(lambda==.(L)))
lines(c,(max(v)-L*c),lty=2,col='blue')
abline(v=3.1,lty=2,col='green')
plot(c,v,type='l',
	xlab='Fecundity at age x',
	ylab='Lifetime Reproductive Output')
text(1,10.5,bquote(lambda==.(L)))
abline(h=max(v),lty=2,col='blue')
abline(v=3.1,lty=2,col='green')

# Scenario 2 - Lambda>1 = de-value 
#     residual reproductive value	
L=2
c=seq(from=0, to=5,by=.1)
f=10-.02*c^2-.03*c^3	    
v=c+f/L
plot(c,f,type='l',
	xlab='Fecundity at age x',
	ylab='Residual Repro. Value')
text(1,7,bquote(lambda==.(L)))
lines(c,(15.861-L*c),lty=2,col='blue')
abline(v=4.5,lty=2,col='green')
plot(c,v,type='l',
	xlab='Fecundity at age x',
	ylab='Lifetime Reproductive Output')
text(1,5.5,bquote(lambda==.(L)))
abline(h=max(v),lty=2,col='blue')
abline(v=4.5,lty=2,col='green')

# Scenario 3 - Lambda<1 = increase 
#     residual reproductive value	
L=0.5
c=seq(from=0, to=5,by=.1)
f=10-.02*c^2-.03*c^3
v=c+f/L
plot(c,f,type='l',
	xlab='Fecundity at age x',
	ylab='Residual Repro. Value')
text(1,7,bquote(lambda==.(L)))
lines(c,(10.68-L*c),lty=2,col='blue')
abline(v=2.15,lty=2,col='green')
plot(c,v,type='l',
	xlab='Fecundity at age x',
	ylab='Lifetime Reproductive Output')
text(1,18,bquote(lambda==.(L)))
abline(h=max(v),lty=2,col='blue')
abline(v=2.15,lty=2,col='green')