Lect3 <- function()
{
 R0Val <- 1.43563216170687
 Steepness <- 0.550149
 surv <- exp(-0.4)
 AgeM <- 8
 MaxY <- 150
 SigmaR <- 0.6

 TheCatch <- read.csv("C:\\courses\\fish558_09\\Lectures\\Lec3Data.Csv",header=F)
 Nyear <- length(TheCatch[,2])
 
 Wght <- 0.0055*(230.3*(1-exp(-0.046*(0:7+0.825))))^3.084
 Fecu <- rep(0,AgeM)
 for (Ia in 1:AgeM)
  if (Ia < 4+1) Fecu[Ia] <- 0 else Fecu[Ia] <- Wght[Ia]
 Sel <- 1.0/(1.0+exp(-(0:7-3)*2))

 Nsim <- 1
 CatchSeq <- seq(from=0,to=200,by=1)
 Resu <- matrix(0,ncol=5,nrow=length(CatchSeq))
 for (II in 1:length(CatchSeq))
  {
   Above <- 0
   Below <- 0
   AveBio <- 0
   AveDepl <- 0
   for (Isim in 1:Nsim)
    {
     Eps <- c(rep(0,Nyear+20))
     RecMult <- exp(Eps)
     Outs <- ModelPrj(R0Val,Steepness,surv,AgeM,Wght,Fecu,Sel,TheCatch,Nyear,RecMult,20,CatchSeq[II])
     AveBio <- AveBio + Outs$Bfinal/Nsim
     AveDepl <- AveDepl + Outs$Depl*100/Nsim
     if (Outs$Depl > 0.4) Above <- Above + 100/Nsim
     if (Outs$Depl < 0.1) Below <- Below + 100/Nsim
    } 
   Resu[II,1] <- CatchSeq[II] 
   Resu[II,2] <- Above
   Resu[II,3] <- Below
   Resu[II,4] <- AveBio
   Resu[II,5] <- AveDepl
   cat(CatchSeq[II],Above,Below,AveBio,"\n") 
  }
 
 par(mfrow=c(2,2))
 plot(Resu[,1],Resu[,2],xlab="Catch limit",ylab="Probability of B > 0.4B0",ylim=c(0,105),type='l',pch=16,yaxs="i",lwd=2)
 plot(Resu[,1],Resu[,3],xlab="Catch limit",ylab="Probability of B < 0.1B0",ylim=c(0,105),type='l',pch=16,yaxs="i",lwd=2)
 plot(Resu[,1],Resu[,4],xlab="Catch limit",ylab="Average biomass",type='l',pch=16,yaxs="i",lwd=2)
 plot(Resu[,1],Resu[,5],xlab="Catch limit",ylab="Average depletion",type='l',ylim=c(0,105),pch=16,yaxs="i",lwd=2)

 Nsim <- 100
 CatchSeq <- seq(from=0,to=200,by=5)
 Resu <- matrix(0,ncol=5,nrow=length(CatchSeq))
 for (II in 1:length(CatchSeq))
  {
   Above <- 0
   Below <- 0
   AveBio <- 0
   AveDepl <- 0
   for (Isim in 1:Nsim)
    {
     Eps <- c(rep(0,Nyear),rnorm(20,0,SigmaR)-SigmaR^2/2)
     RecMult <- exp(Eps)
     Outs <- ModelPrj(R0Val,Steepness,surv,AgeM,Wght,Fecu,Sel,TheCatch,Nyear,RecMult,20,CatchSeq[II])
     AveBio <- AveBio + Outs$Bfinal/Nsim
     AveDepl <- AveDepl + Outs$Depl*100/Nsim
     if (Outs$Depl > 0.4) Above <- Above + 100/Nsim
     if (Outs$Depl < 0.1) Below <- Below + 100/Nsim
    } 
   Resu[II,1] <- CatchSeq[II] 
   Resu[II,2] <- Above
   Resu[II,3] <- Below
   Resu[II,4] <- AveBio
   Resu[II,5] <- AveDepl
   cat(CatchSeq[II],Above,Below,AveBio,"\n") 
  }
 
 par(mfrow=c(2,2))
 plot(Resu[,1],Resu[,2],xlab="Catch limit",ylab="Probability of B > 0.4B0",ylim=c(0,105),type='b',pch=16,yaxs="i")
 plot(Resu[,1],Resu[,3],xlab="Catch limit",ylab="Probability of B < 0.1B0",ylim=c(0,105),type='b',pch=16,yaxs="i")
 plot(Resu[,1],Resu[,4],xlab="Catch limit",ylab="Average biomass",type='b',pch=16,yaxs="i")
 plot(Resu[,1],Resu[,5],xlab="Catch limit",ylab="Average depletion",type='b',ylim=c(0,105),pch=16,yaxs="i")

 }
 
 # ===========================================================================================

ModelPrj<-function(R0Val,Steepness,surv,AgeM,Wght,Fecu,Sel,TheCatch,Nyear,RecMult,FutProj,Catch)
{

 MaxY <-  Nyear + FutProj

 # Compute Alpha and Beta
 Alpha <- (1-Steepness)/(4*Steepness*R0Val)
 Beta <- (5*Steepness-1)/(4*Steepness*R0Val)

 N <- matrix(0,ncol=AgeM,nrow=MaxY+1)
 SSB <- rep(0,nrow=MaxY+1)
 ExpBio <- rep(0,nrow=MaxY)
 ExpRate <- rep(0,nrow=MaxY)
 
 # Set Initial age-structue
 N[1,1] <- R0Val
 for (age in 2:AgeM) N[1,age] <- N[1,age-1]*surv
 N[1,AgeM] <- N[1,AgeM]/(1-surv)
 SSB[1] <- 0
 for (Ia in 1:AgeM) SSB[1] <- SSB[1] + N[1,Ia]*Fecu[Ia]

 # Project forward 
 for (Year in 1:MaxY)
  {
   ExpBio[Year] <- 0
   for (Ia in 1:AgeM) ExpBio[Year] <- ExpBio[Year] + N[Year,Ia]*Wght[Ia]*Sel[Ia]
   
   # Exploitation rate   
   if (Year <=Nyear)
    ExpRate[Year] <- TheCatch[Year,2]/ExpBio[Year]
   else
    ExpRate[Year] <- Catch/ExpBio[Year]
   if (ExpRate[Year] > 0.99) ExpRate[Year] <- 0.99 
       
   # Update age-classes
   N[Year+1,AgeM] <- N[Year,AgeM-1]*surv*(1-Sel[AgeM-1]*ExpRate[Year])+ N[Year,AgeM]*surv*(1-Sel[AgeM]*ExpRate[Year])
   for (Ia in 2:(AgeM-1))
    N[Year+1,Ia] <- N[Year,Ia-1]*surv*(1-Sel[Ia-1]*ExpRate[Year])
   SSB[Year+1] <- 0
   for (Ia in 1:AgeM) SSB[Year+1] <- SSB[Year+1] + N[Year+1,Ia]*Fecu[Ia]
   N[Year+1,1] <- SSB[Year+1]/SSB[1]/(Alpha+Beta*SSB[Year+1]/SSB[1])*RecMult[Year]
   
  }
 Bfinal <- SSB[(MaxY+1)]
 
 Outs <- NULL
 Outs$Bfinal <- Bfinal
 Outs$Depl <- Bfinal/SSB[1]
 return(Outs)
  
}

# =========================================================================================

Lect3()