model {
	
	#### Priors and constraints ----
	
	## the recapture parameters: time periods 1:40
	for( t in 1:nyears ){
		for( a in 1:3 ){                                      
			eps.p1[a, t] ~ dnorm(0, tau.p1[a])
			logit(p[a, t]) <- mu.p1[a] + eps.p1[a, t]         
			q[a, t] <- 1 - p[a, t]                              
		} #a
	} #t
    
	mu.p1[1] <- 0
	tau.p1[1] <- 1
	for( a in 2:3 ){
		mean.p1[a] ~ dunif(0, 1)                               # (BPA: 186)
		mu.p1[a] <- logit(mean.p1[a])
		sigma.p1[a] ~ dt(0, 1, 2)T(0, )
		tau.p1[a] <- pow(sigma.p1[a], -2)
	} #a
  
	## the recapture parameters: time periods 41:...
	for( t in 41:(n.occ-1) ){
		for( a in 1:3 ){                                      
			eps.p2[a, t] ~ dnorm(0, tau.p2[a])
			logit(p[a, t]) <- mu.p2[a] + eps.p2[a, t]         
			q[a, t] <- 1 - p[a, t]                              
		} #a
	} #t
    
	mu.p2[1] <- 0
	tau.p2[1] <- 1
	for( a in 2:3 ){
		mean.p2[a] ~ dunif(0, 1)                               # (BPA: 186)
		mu.p2[a] <- logit(mean.p2[a])
		sigma.p2[a] ~ dt(0, 1, 1)T(0, )
		tau.p2[a] <- pow(sigma.p2[a], -2)
	} #a
  
	## the survival parameters
	# the means
	for( a in 1:3 ){                                  
		for( p in 1:2 ){
			mean.phi[a, p] ~ dunif(0, 1)
			mu.phi[a, p] <- logit(mean.phi[a, p])
		} #p
	} #a 

	# the annual values
	for( t in 1:(n.occ-1) ){                                
		for( a in 1:3 ){              
			logit(month.phi[a, t]) <- mu.phi[a, season[t]] + eps.phi[a, t]
		} #a
		
		# convert to the relevant time periods
		phi[1, t] <- pow(month.phi[1, t], jv.days[t])
		phi[2, t] <- pow(month.phi[2, t], ad.days[t])
		phi[3, t] <- pow(month.phi[3, t], ad.days[t])
		
		# the errors are sampled jointly
		eps.phi[1:3, t] ~ dmnorm.vcov(c(0,0,0), Sigma2[season[t], 1:3, 1:3])

	} #t
	
	for( p in 1:2 ){
		# Priors for the variances on the diagonal
		for( i in 1:3 ){
			sigma[p, i] ~ dunif(0, 5)
			Sigma2[p, i, i] <- sigma[p, i] * sigma[p, i]
			# there hapen to be three (off-diagonal) rho's
			rho[p, i] ~ dunif(-1, 1)   
		} #i
		
		# Specify the off diagonals individually
		Sigma2[p, 1, 2] <- sigma[p, 1] * sigma[p, 2] * rho[p, 1]
		Sigma2[p, 1, 3] <- sigma[p, 1] * sigma[p, 3] * rho[p, 2]
		Sigma2[p, 2, 3] <- sigma[p, 2] * sigma[p, 3] * rho[p, 3]
		Sigma2[p, 2, 1] <- Sigma2[p, 1, 2]
		Sigma2[p, 3, 1] <- Sigma2[p, 1, 3]
		Sigma2[p, 3, 2] <- Sigma2[p, 2, 3]
	} #p
	
	#### Define the multinomial likelihood ----
	
	for (t in 1:(n.occ-1)){
		marr.1[t, 1:n.occ] ~ dmulti(pr.1[t, ], rel.1[t])
		marr.2[t, 1:n.occ] ~ dmulti(pr.2[t, ], rel.2[t])
		marr.3[t, 1:n.occ] ~ dmulti(pr.3[t, ], rel.3[t])
		marr.4[t, 1:n.occ] ~ dmulti(pr.4[t, ], rel.4[t])
	} #t
  
	#### Define the cell probabilities of the m-arrays ----
	## Juveniles - always enter in autumn
	for( t in 1:(n.occ-1) )
		{ pr.1[t, t] <- phi[1, t] * p[2, t]  }                   
	for( t in 2:(n.occ-1) ) 
		{pr.1[(t-1), t] <- phi[1, (t-1)] * phi[2, t] * q[2, (t-1)] * p[2, t] }
	for( t in 3:(n.occ-1) ) 
		{pr.1[(t-2), t] <- phi[1, (t-2)] * prod(phi[2, (t-1):t]) * prod(q[2, (t-2):(t-1)]) * p[2, t] }
	for( t in 4:(n.occ-1) ) 
		{pr.1[(t-3), t] <- phi[1, (t-3)] * prod(phi[2, (t-2):t]) * prod(q[2, (t-3):(t-1)]) * p[2, t] }
	for( t in 5:(n.occ-1) ) 
		{pr.1[(t-4), t] <- phi[1, (t-4)] * prod(phi[2, (t-3):t]) * prod(q[2, (t-4):(t-1)]) * p[3, t] }
	for( t in 1:(n.occ-6) ){             # rows            
		for( u in (t+5):(n.occ-1) ){       # columns
			pr.1[t, u] <- phi[1, t] * prod(phi[2, (t+1):(t+4)]) * prod(phi[3, (t+5):u]) * prod(q[2, t:(t+3)]) * prod(q[3, (t+4):(u-1)]) * p[3, u] 
		} #u
	} #t
	
	## 2nd years - may enter in spring or autumn, so number of imm periods varies
	for( t in 1:(n.occ-1) ) 
		{ pr.2[t, t] <- phi[2, t] * p[2, t] }
	for( t in 2:(n.occ-1) ) 
		{ pr.2[(t-1), t] <- prod(phi[2, (t-1):t]) * q[2, (t-1)] * p[2, t] }
	for( t in 3:(n.occ-1) ) 
		{ pr.2[(t-2), t] <- prod(phi[2, (t-2):t]) * prod(q[2, (t-2):(t-1)]) * p[idx[t], t] }
	for( t in 4:(n.occ-1) ) 
		{ pr.2[(t-3), t] <- prod(phi[2, (t-3):(t-1)]) * phi[idx[(t-1)],t] * prod(q[2, (t-3):(t-2)]) * q[idx[(t-1)],t] * p[3, t] }
	for( t in 1:(n.occ-5) ){
		for( u in (t+4):(n.occ-1) ){   
			pr.2[t, u] <- prod(phi[2, t:(t+2)]) * phi[idx[(t+3)],(t+3)] * prod(phi[3, (t+4):u]) * prod(q[2, t:(t+1)]) * q[idx[(t+2)],(t+2)] * prod(q[3, (t+3):(u-1)]) * p[3, u] 
		} #u
	} #t
	
	## Third-years
	for( t in 1:(n.occ-1) ) 
		{ pr.3[t, t] <- phi[2, t] * p[idx[t], t] }
	for( t in 2:(n.occ-1) ) 
		{ pr.3[(t-1), t] <- phi[2, (t-1)] * phi[idx[(t-1)], t] * q[idx[(t-1)], (t-1)] * p[3, t] }
	for( t in 1:(n.occ-3) ){
		for( u in (t+2):(n.occ-1) ){   
			pr.3[t, u] <- phi[2, t] * phi[idx[t], (t+1)] * prod(phi[3, (t+2):u]) * q[idx[t], t] * prod(q[3, (t+1):(u-1)]) * p[3, u]
		} #u
	} #t
	
	## Adults
	for( t in 1:(n.occ-1) ) 
		{ pr.4[t, t] <- phi[3, t] * p[3, t] }
	for( t in 1:(n.occ-2) ){                                
		for( u in (t+1):(n.occ-1) ){   
			pr.4[t, u] <- prod(phi[3, t:u]) * prod(q[3, t:(u-1)]) * p[3, u] 
		} #u
	} #t


	#### Finishing off the matrices ----

	for (t in 1:(n.occ-1)){
		## Last column: probability of non-recapture
		pr.1[t, n.occ] <- 1 - sum(pr.1[t, 1:(n.occ-1)])       
		pr.2[t, n.occ] <- 1 - sum(pr.2[t, 1:(n.occ-1)])
		pr.3[t, n.occ] <- 1 - sum(pr.3[t, 1:(n.occ-1)])
		pr.4[t, n.occ] <- 1 - sum(pr.4[t, 1:(n.occ-1)])

		## the cells below the diagonals 
		for ( u in 1:(t-1) ){                                 
		  pr.1[t, u] <- 0                                     
		  pr.2[t, u] <- 0
		  pr.3[t, u] <- 0
		  pr.4[t, u] <- 0
		} #u
	} #t
  
}