#Assignment 1 for Florida

A		<-	15
age		<-	1:A

ro		<-	100
kap		<-	12
m		<-	0.23
winf		<-	1
k		<-	0.153
adot		<-	4.5
sd.adot	<-	1.5
ahat		<-	5.0
sd.ahat	<-	1.25
fe		<-	0.0

##	-------------------------------------------------------------------
equil		<-	function(fe=0)
{
	wa		<-	winf*(1-exp(-k*age))^3
	ma		<-	plogis(age,adot,sd.adot)
	va		<-	plogis(age,ahat,sd.ahat)
	iota		<-	exp(-m)^(age-1); iota[A]=iota[A]/(1-exp(-m))
	
	iota.f	<-	iota
	for(i in 2:A) iota.f[i]=iota.f[i-1]*exp(-m-fe*va[i-1])
	iota.f[A]=iota.f[A]/(1-exp(-m-fe*va[A]))
	
	phiE		<-	sum(iota*wa*ma)
	phie		<-	sum(iota.f*wa*ma)
	phib		<-	sum(iota.f*wa)
	phiq		<-	sum(iota.f*wa*va*(1-exp(-m-fe*va))/(m+fe*va))
	
	re		<-	max(0,ro*(kap-phiE/phie)/(kap-1))
	ye		<-	fe*re*phiq
	ypr		<-	fe*phiq
	be		<-	re*phib
	spr		<-	phie/phiE
	
	return(list(re=re,ye=ye,be=be,spr=spr,ypr=ypr))
}

fn.ye	<- function(fe) equil(fe)$ye
fn.re <- function(fe) equil(fe)$re
fn.be <- function(fe) equil(fe)$be
fn.spr <- function(fe) equil(fe)$spr
fn.ypr <- function(fe) equil(fe)$ypr

get.fmsy	<-	function(fe)
{
	#A much more efficient algorithm could be used
	#ie. Newton Raphson method to find dYe/dfe=0
	ye<-sapply(fe,fn.ye)
	return(fe[which(ye==max(ye))])
}


fe	<- seq(0,1.5,by=0.01)
ye	<- sapply(fe,fn.ye)
be	<- sapply(fe,fn.be)
re	<- sapply(fe,fn.re)
spr	<- sapply(fe,fn.spr)
ypr	<- sapply(fe,fn.ypr)

im	<- which(ye==max(ye))
print(paste("MSY =",round(ye[im],2),"Fmsy =",fe[im],"Fmax = ",fe[min(which(ye[-1]==0))]))
print(paste("Bmsy/Bo =",round(be[im]/be[1],2),"Rmsy/Ro =",round(re[im]/re[1],2)))
print(paste("SPR at MSY =",round(spr[im],2)))

set.seed(1233)		
ftry	<- seq(0.1,0.5,by=0.01)
md=rnorm(200,m,0.1*m)
fn=function(md){m<<-md; get.fmsy(ftry)}
fmsy=sapply(md,fn)


##	-------------------------------------------------------------------
#SET UP GRAPHICS DEVICE FOR PLOTTING (page-up page-down to scroll through figs)
graphics.off()
if(exists(".SavedPlots",where=1)==T){rm(.SavedPlots,pos=1)}
#windows(record=T); 
par(las=1,ps=12,cex.main=14/12)
par(mfcol=c(2,2))
plot(fe,ye,type="l",xlab="Fishing rate (fe)",ylab="Yield")
plot(fe,be,type="l",xlab="Fishing rate (fe)",ylab="Biomass")
plot(fe,re,type="l",xlab="Fishing rate (fe)",ylab="Age-1 recruits")
plot(fe,spr,type="l",xlab="Fishing rate (fe)",ylab="SPR")

par(mfcol=c(1,1),pty="s",cex.lab=1.5,cex.axis=1.5)
plot(density(md,adjust=1.65),xlab="M",ylab="Prior density",main="",xlim=c(0.9*min(md,fmsy),1.1*max(md,fmsy)))
lines(density(fmsy,adjust=1.65),lty=2,col=2,xlab="Fmsy",ylab="Prior density",main="")
legend("topright",c("M","Fmsy"),lty=c(1,2),col=c(1,2),bty="n")
#dev.copy2eps(file="ImpliedFmsy.eps")