library(stockassessment)
load("run/model.RData")

# not used
FC<-list()

#for partial mortality F_ibc
oldwd<-setwd("data")
ibf<-read.ices("whg47d_ibf.dat")

#ref points
Fmsy<-0.41  
Fpa<-0.44
Fsq<-0.085

Bpa<-195927
Blim<-139948

# catch TAC for area 4 in  2025
TAC<-111861

# proportion caught in area 4
propc4<-0.8023


# use recent 3 years of fishing selectivity
selYears<-max(fit$data$years) +(-3:-1)  

# to calculate partial Fibc########################
# mean recent 3 years of proportion ibc
ibc_mean<-apply(ibf[length(ibf[,1])+(-2:0),],2,mean) 

#estimated total F at age
summary_F<-summary(fit)
f_est<-summary_F[,7]

#total selectivity
fsel <- faytable(fit)/f_est

# mean recent 3 years of F, scaled t final year
catch_mean<-apply(fsel[length(fsel[,1])+(-3:-1),],2,mean)*f_est[length(f_est)-1]  
IBC_sel<-ibc_mean*catch_mean
Fibc<-mean(IBC_sel[3:6]) #age 2-5



set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, ave.years=max(fit$data$years) +(-3:-1), rec.years=2002:max(fit$data$years),overwriteSelYears=selYears,fval=c(Fsq,Fsq,Fsq,Fsq),splitLD =FALSE, label="FSQ, all years")


set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, ave.years=max(fit$data$years) +(-3:-1), rec.years=2002:max(fit$data$years),overwriteSelYears=selYears, fval=c(Fsq,0.000001,0.000001, 0.000001),splitLD =FALSE, label="SQ then zero")

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, ave.years=max(fit$data$years) +(-3:-1), rec.years=2002:max(fit$data$years),overwriteSelYears=selYears, fval=c(Fsq,Fmsy,Fmsy, Fmsy),splitLD =FALSE, label="SQ then Fmsy")

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, ave.years=max(fit$data$years) +(-3:-1), rec.years=2002:max(fit$data$years),overwriteSelYears=selYears, fval=c(Fsq,Fpa,Fpa, Fpa),splitLD =FALSE, label="SQ then Fpa")


set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, ave.years=max(fit$data$years) +(-3:-1), rec.years=2002:max(fit$data$years),overwriteSelYears=selYears, fscale=c(NA,0.75,0.75,0.75),splitLD =FALSE, fval=c(Fsq,NA,NA, NA), label="SQ then 0.75*SQ")

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, ave.years=max(fit$data$years) +(-3:-1), rec.years=2002:max(fit$data$years),overwriteSelYears=selYears, fscale=c(NA,1.25,1.25,1.25), fval=c(Fsq,NA,NA, NA),splitLD =FALSE, label="SQ then 1.25*SQ")


set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, ave.years=max(fit$data$years) +(-3:-1), rec.years=2002:max(fit$data$years),overwriteSelYears=selYears, fval=c(Fsq,Fibc,Fibc,Fibc), customSel=IBC_sel,splitLD =FALSE, label="SQ then no HC fishery") # use partial mortality for IBC of SQ


set.seed(12345)
FC[[length(FC)+1]]<-forecast(fit, ave.years=max(fit$data$years) +(-3:-1),rec.years=2002:max(fit$data$years),overwriteSelYears=selYears, fval=c(Fsq,NA,NA,NA), nextssb=c(NA,Blim,Blim,Blim),splitLD =FALSE, label="SQ then Blim in year 3")

set.seed(12345)
FC[[length(FC)+1]]<-forecast(fit,ave.years=max(fit$data$years) +(-3:-1),rec.years=2002:max(fit$data$years),overwriteSelYears=selYears, fval=c(Fsq,NA,NA,NA), nextssb=c(NA,Bpa,Bpa,Bpa),splitLD =FALSE, label="SQ then Bpa in year 3")


setwd("..")
save(FC, file="run/forecast.RData")


