code <- nimbleCode({ 
  #### likelihood ####
  for(i in 1:M){
    # First occasion
    # State proces
    z[i,1]~dbern(gamma[1])
    mu1[i]<-z[i,1]*p[i,1]
    # Observation proces
    Y[i,1]~dbern(mu1[i])
    
    for (t in 2:Occasions) { #subsequent occasions 
      # State proces
      q[i,t-1]<-1-z[i, t-1]
      logit(phi[i, t-1])<-beta0+beta1*F.Obs[i]
      
      mu2[i,t]<-phi[i, t-1]*z[i,t-1]+gamma[t]*prod(q[i,1:(t-1)])
      z[i,t] ~ dbern(mu2[i,t])
      
      # Observed proces
      mu3[i,t]<-z[i,t]*p[i,t]
      Y[i,t]~dbern(mu3[i,t])
    }
  }
  
  # Claculate deriven populaiton parameters
  for(t in 1:Occasions){
    qgamma[t]<-1-gamma[t]
  }
  
  cprob[1]<-gamma[1]
  for(t in 2:Occasions){
    cprob[t]<-gamma[t]*prod(qgamma[1:(t-1)])
  }
  
  psi<-sum(cprob[1:Occasions]) #inclussion probability
  for(t in 1:Occasions){
    b[t]<-cprob[t]/psi #entry probability
  }
  
  for(i in 1:M){
    recruit[i,1]<-z[i,1]
    for(t  in 2:Occasions){
      recruit[i,t]<-(1-z[i,t-1])*z[i,t]
    }
  }
  
  for(t in 1:Occasions){
    N[t]<-sum(z[1:M, t]) #actual pop size
    B[t]<-sum(recruit[1:M, t]) #Number of entries
  }
  
  for(i in 1:M){
    Nind[i]<-sum(z[i, 1:Occasions])
    Nalive[i]<-1-equals(Nind[i], 0)
  }
  
  Nsuper<-sum(Nalive[1:M]) #size of super populaiton
  
  #### Priors and Constrains ####
  for (i in 1:M) {
    for(t in 1:Occasions){
      p[i,t]<-mean.p
    }
  }
  
  beta0~dflat() 
  beta1~dflat()
  
  mean.p~dunif(0, 1)
  
  for(t in 1:Occasions){
    gamma[t]~dunif(0, 1)
  }  
  
  })