			model {

			# Set up 'data' to define spatial dependence structure
			# =====================================
				for(i in 1 : N) { 
					m[i] <- 1/E[i]       # scaling factor for variance in each cell
				} 

			# The vector C[] required as input into the car.proper distribution is a vector
			# respresention of the weight matrix with elements Cij. The first J1 elements of the C[]
			# vector contain the weights for the J1 neighbours of area i=1; the (J1+1) to J2
			# elements of the C[] vector contain the weights for the J2 neighbours of area i=2;
			# etc. To set up this vector, we need to define a variable cumsum, which gives the
			# values of J1, J2, etc.; we then set up an index matrix pick[,] with N columns
			#  corresponding to the  i=1,...,N areas, and with the same number of rows as there are
			# elements in the C[] vector (i.e. sumNumNeigh). The elements 
			#C[ (cumsum[i]+1):cumsum[i+1] ] correspond to 
			# the set of weights Cij associated with area i, and so we set up ith column of the
			# matrix pick[,]to have a 1 in all the rows k for which 
			#cumsum[i] < k <= cumsum[i+1], and 0's elsewhere. 
			# For example, let N=4 and cumsum=c(0,3,5,6,8), so area i=1 has 3 neighbours, area
			# i=2 has 2 neighbours, area i=3 has 1 neighbour and area i=4 has 2 neighbours. The
			# the matrix pick[,] is:
			#                pick                  
			#             1, 0, 0, 0,                   
			#             1, 0, 0, 0,                   
			#             1, 0, 0, 0,                    
			#             0, 1, 0, 0,                   
			#             0, 1, 0, 0,                   
			#             0, 0, 1, 0,                   
			#             0, 0, 0, 1,                   
			#             0, 0, 0, 1,                   
			#
			# We can then use the inner product (inprod(,)) function in WinBUGS and the kth
			# row of pick to select which area corresponds to the kth element in the vector C[];
			# likewise, we can use inprod(,) 
			# and the ith column of pick to select the elements of C[] which correspond to area i.
			#
			# Note: this way of setting up the C vector is somewhat convoluted!!!! In future
			# versions, we hope the GeoBUGS adjacency matrix tool will be able to dump out the
			# relevant  vectors required. Alternatively, the C vector could be created using another
			# package (e.g. Splus) and read into WinBUGS as data.
			#
				cumsum[1] <- 0
				for(i in 2:(N+1)) {
					cumsum[i] <- sum(num[1:(i-1)])
				}	
				for(k in 1 : sumNumNeigh) { 
					for(i in 1:N) {
						pick[k,i] <- step(k - cumsum[i] - epsilon)  * step(cumsum[i+1] - k)   
						#  pick[k,i] = 1    if     cumsum[i] < k <= cumsum[i=1];  otherwise, pick[k,i] = 0
					}                                                       
					C[k] <- sqrt(E[adj[k]] / inprod(E[], pick[k,]))    # weight for each pair of neighbours
				}
				epsilon <- 0.0001

			# Model
			# =====

			# Likelihood
				for (i in 1 : N) {
					O[i]  ~ dpois(mu[i])
					log(mu[i]) <- log(E[i]) + S[i]
			# Area-specific relative risk 
					RR[i] <- exp(S[i])
					theta[i] <- alpha
				}

			# Proper CAR prior distribution for spatial random effects: 
				S[1:N] ~ car.proper(theta[], C[], adj[], num[], m[], prec, gamma)

			# Other priors:
				alpha  ~ dnorm(0, 0.0001)  
			# prior on precision
				prec  ~ dgamma(0.5, 0.0005)
				v <- 1/prec				# variance
				sigma <- sqrt(1 / prec)				# standard deviation

				gamma.min <- min.bound(C[], adj[], num[], m[])
				gamma.max <- max.bound(C[], adj[], num[], m[])
				gamma ~ dunif(gamma.min, gamma.max)

			}
