	model{ 
		theta[1] ~ dunif(theta[2], theta.max)
		theta[2] ~ dunif(theta[3], theta[1])
		theta[3] ~ dunif(theta[9], theta[2])
		theta[4] ~ dunif(theta[9], theta.max)
		theta[5] ~ dunif(theta[7], theta.max)
		theta[6] ~ dunif(theta[7], theta.max)
		theta[7] ~ dunif(theta[9], theta7max)
		theta7max <- min(theta[5], theta[6])
		theta[8] ~ dunif(theta[9], theta.max)
		theta[9] ~ dunif(theta[10], theta9max)
		theta9max <-min(min(theta[3], theta[4]),  min(theta[7], theta[8]))
		theta[10] ~ dunif(theta[11], theta[9])
		theta[11] ~ dunif(0 ,theta[10])
		
		bound[1] <- ranked(theta[1:8],  8)
		bound[2] <- ranked(theta[1:8], 1)
		bound[3] <- ranked(theta[9:11], 3)
		bound[4] <- ranked(theta[9:11], 1)
		
		for (j in 1 : 5){
			theta[j + 11] ~ dunif(0, theta.max)
			within[j, 1] <- 1 - step(bound[1] - theta[j + 11])
			for (k in 2 : 4){
				within[j, k] <- step(bound[k - 1] - theta[j + 11]) 
					- step(bound[k] - theta[j + 11])
			}
			within[j, 5] <- step(bound[4] - theta[j + 11])
		}


		for (i in 1:nDate){
			X[i] ~ dnorm(mu[i], tau[i])
			tau[i] <- 1/pow(sigma[i],2)
			mu[i] <- interp.lin(theta[i], calBP[], C14BP[])

	# monitor the following variable to smooth density of theta
			theta.smooth[i] <- 10 * round(theta[i] / 10)
		}
	}
