//Pella-Tomlinson surplus production model 

data{
  int<lower=1> N;
  vector[N] LnCPUE;
  vector<lower=0>[N] Catch;
  real LnK_mu;
  real<lower=0> LnK_sd;
  real Lnr_mu;
  real<lower=0> Lnr_sd;
}

parameters{
  real LnK;
  real Lnr;
  real Lnq1;
  real<lower=0> m;
  real<lower=0> h;
  vector[N] mu_raw; //put here so we can use in transformed params
  real<lower=0> sigma;
  real<lower=0> tau;
  real<lower=0> phi;
}

transformed parameters{
 real<lower=0> K;
 real<lower=0> r;
 real<lower=0> Bnew;
 real<lower=0> Pold;
 real<lower=0> Pnew;
 vector<lower=0>[N] B; 
 vector[N] LnCPUE_hat;
 vector[N] mu;
 
 K = exp(LnK);
 r = exp(Lnr);
 mu = mu_raw*sigma; // implies mu ~ normal(0, sigma);

 // Population model. But putting it in transformed params we 
 // speed up computations
 B[1] = K * phi;
 
 for (t in 2:N){
      Pold = B[t-1]/K;
      Pnew = Pold+(r/(m-1))*Pold*(1-Pold^(m-1))-Catch[t-1]/K;
      Bnew = Pnew*K*exp(mu[t]);
  
   if (Bnew < 0.001){
      B[t] = 0.001;
   } else {
      B[t] = Bnew;
   }
 }
 
  LnCPUE_hat = Lnq1 + h*log(B);

}

model{ 
  // Priors
  LnK ~ normal(LnK_mu, LnK_sd);
  Lnr ~ normal(Lnr_mu, Lnr_sd);
  m ~ normal(2,0.5);
  h ~ uniform(0.7, 1.3);
  Lnq1 ~ uniform(log(0.7), log(1.3));
  sigma ~ uniform(0.001, 0.3);
  tau ~ uniform(0.001, 0.3);
  phi ~ uniform(0.5, 1.0);

  // process errors
  mu_raw ~ normal(0, 0.5);

  // sampling 
  LnCPUE ~ normal(LnCPUE_hat, tau);

}

generated quantities{
  real<lower=0> Hmsy;
  real<lower=0> Bmsy;
  real<lower=0> MSY;
  vector<lower=0>[N] Brel; 

  Hmsy =  (r/(m-1)) * (1 - (1/m));
  Bmsy = K * (m^(-1/((m-1))));
  MSY = Hmsy * Bmsy;
  Brel = B/K;
}