library(MASS)
library(scales)
library(ggplot2)
library(gridExtra)

####Amur Estuary S2
A = matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1799, ##fecundity divided by 4 because they reproduce ~quadrennially.
              0.00053, 0, 0, 0, 0, 0, 0, 0, 0, 0,
              0, 0.67, 0, 0, 0, 0, 0, 0, 0, 0,
              0, 0, 0.67, 0, 0, 0, 0, 0, 0, 0,
              0, 0, 0, 0.67, 0, 0, 0, 0, 0, 0,
              0, 0, 0, 0, 0.67, 0, 0, 0, 0, 0,
              0, 0, 0, 0, 0, 0.67, 0, 0, 0, 0,
              0, 0, 0, 0, 0, 0, 0.67, 0, 0, 0,
              0, 0, 0, 0, 0, 0, 0, 0.67, 0, 0,
              0, 0, 0, 0, 0, 0, 0, 0, 0.67, 0.57), nr = 10, byrow = TRUE)

# initial population vector N  
N0 = matrix(c(51919140, 11088, 11088, 11088, 11088, 11088, 11088, 11088, 11088, 28860), ncol = 1)

N1 = A %*% N0 #%% is a matrix multiplier command

#simulation length
years = 30
N.projected = matrix(0, nrow = nrow(A), ncol = years+1)
N.projected[, 1] = N0

#matrix multiplication
for (i in 1:years)
{
  N.projected[1,i] <- N.projected[1,i] + 1250000 #adding y0 fish; represents females only of 2x10^6 total
  N.projected[, i + 1] = A %*% N.projected[,i]
}

#dataframe for plots
df.Estuary = as.data.frame(t(N.projected))
df.Estuary[,11] = c(2014:2044)
colnames(df.Estuary) = c("Fry", "Y1", "Y2", "Y3", "Y4", "Y5", "Y6", "Y7", "Y8", "Mature", "Year")
df.Estuary$not.mature = rowSums(df.Estuary[,1:9]) #abundance of immature fish
df.Estuary$total = rowSums(df.Estuary[,1:10]) #total abundance
df.Estuary$one.up = rowSums(df.Estuary[,2:10]) #abundance of immature fish


####Amur Lower S2
A = matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1799, ##fecundity divided by 4 because they reproduce ~quadrennially.
              0.00053, 0, 0, 0, 0, 0, 0, 0, 0, 0,
              0, 0.67, 0, 0, 0, 0, 0, 0, 0, 0,
              0, 0, 0.67, 0, 0, 0, 0, 0, 0, 0,
              0, 0, 0, 0.67, 0, 0, 0, 0, 0, 0,
              0, 0, 0, 0, 0.67, 0, 0, 0, 0, 0,
              0, 0, 0, 0, 0, 0.67, 0, 0, 0, 0,
              0, 0, 0, 0, 0, 0, 0.67, 0, 0, 0,
              0, 0, 0, 0, 0, 0, 0, 0.67, 0, 0,
              0, 0, 0, 0, 0, 0, 0, 0, 0.67, 0.57), nr = 10, byrow = TRUE)

# initial population vector N  
N0 = matrix(c(764575, 1483, 1483, 1483, 1483, 1483, 1483, 1483, 1483, 425), ncol = 1)

N1 = A %*% N0 #%% is a matrix multiplier command

#Simulation length
years = 30
N.projected = matrix(0, nrow = nrow(A), ncol = years+1)
N.projected[, 1] = N0

#matrix multiplication
for (i in 1:years)
{
  N.projected[1,i] <- N.projected[1,i] + 1250000 #adding y0 fish; represents females only of 2x10^6 total
  N.projected[, i + 1] = A %*% N.projected[,i]
}

#df for plots
df.Lower = as.data.frame(t(N.projected))
df.Lower[,11] = c(2020:2050)
colnames(df.Lower) = c("Fry", "Y1", "Y2", "Y3", "Y4", "Y5", "Y6", "Y7", "Y8", "Mature", "Year")
df.Lower$not.mature = rowSums(df.Lower[,1:9])
df.Lower$total = rowSums(df.Lower[,1:10])
df.Lower$one.up = rowSums(df.Lower[,2:10])


####Middle Amur Scenario 2
A = matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1799, ##fecundity divided by 4 because they reproduce ~quadrennially.
              0.00053, 0, 0, 0, 0, 0, 0, 0, 0, 0,
              0, 0.67, 0, 0, 0, 0, 0, 0, 0, 0,
              0, 0, 0.67, 0, 0, 0, 0, 0, 0, 0,
              0, 0, 0, 0.67, 0, 0, 0, 0, 0, 0,
              0, 0, 0, 0, 0.67, 0, 0, 0, 0, 0,
              0, 0, 0, 0, 0, 0.67, 0, 0, 0, 0,
              0, 0, 0, 0, 0, 0, 0.67, 0, 0, 0,
              0, 0, 0, 0, 0, 0, 0, 0.67, 0, 0,
              0, 0, 0, 0, 0, 0, 0, 0, 0.67, 0.57), nr = 10, byrow = TRUE)

# initial population vector N  
N0 = matrix(c(152915, 297, 297, 297, 297, 297, 297, 297, 297, 85), ncol = 1)

N1 = A %*% N0 #%% is a matrix multiplier command

#Scenario length
years = 30
N.projected = matrix(0, nrow = nrow(A), ncol = years+1)
N.projected[, 1] = N0

#matrix multiplication
for (i in 1:years)
{
  N.projected[1,i] <- N.projected[1,i] + 1250000 #adding y0 fish; represents females only of 2x10^6 total
  N.projected[, i + 1] <- A %*% N.projected[,i]
}

#df for plots
df.Middle = as.data.frame(t(N.projected))
df.Middle[,11] = c(2020:2050)
colnames(df.Middle) = c("Fry", "Y1", "Y2", "Y3", "Y4", "Y5", "Y6", "Y7", "Y8", "Mature", "Year")
df.Middle$not.mature = rowSums(df.Middle[,1:9])
df.Middle$total = rowSums(df.Middle[,1:10])
df.Middle$one.up = rowSums(df.Middle[,2:10])


####Upper Amur Scenario 2
A = matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1799, ##fecundity divided by 4 because they reproduce ~quadrennially.
             0.00053, 0, 0, 0, 0, 0, 0, 0, 0, 0,
             0, 0.67, 0, 0, 0, 0, 0, 0, 0, 0,
             0, 0, 0.67, 0, 0, 0, 0, 0, 0, 0,
             0, 0, 0, 0.67, 0, 0, 0, 0, 0, 0,
             0, 0, 0, 0, 0.67, 0, 0, 0, 0, 0,
             0, 0, 0, 0, 0, 0.67, 0, 0, 0, 0,
             0, 0, 0, 0, 0, 0, 0.67, 0, 0, 0,
             0, 0, 0, 0, 0, 0, 0, 0.67, 0, 0,
             0, 0, 0, 0, 0, 0, 0, 0, 0.67, 0.57), nr = 10, byrow = TRUE)

# initial population vector N  
N0 <- matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0), ncol = 1)

N1 = A %*% N0 #%% is a matrix multiplier command
N1

#Scenario length
years = 30
N.projected = matrix(0, nrow = nrow(A), ncol = years+1)
N.projected[, 1] = N0

#matrix multiplication
for (i in 1:years)
{
  N.projected[1,i] <- N.projected[1,i] + 1250000 #adding y0 fish; represents females only of 2x10^6 total
  N.projected[, i + 1] <- A %*% N.projected[,i]
}

#df for plots
df.Upper = as.data.frame(t(N.projected))
df.Upper[,11] = c(2020:2050)
colnames(df.Upper) = c("Fry", "Y1", "Y2", "Y3", "Y4", "Y5", "Y6", "Y7", "Y8", "Mature", "Year")
df.Upper$not.mature = rowSums(df.Upper[,1:9]) #abundance of immature fish
df.Upper$total = rowSums(df.Upper[,1:10]) #total abundance
df.Upper$one.up = rowSums(df.Upper[,2:10]) #fish at least 1 y old