WinBUGS Code for Bayesian hierarchical model

Question

I’m reading this paper (which aims to model the nanowires (NW) growth using a Bayesian hierarchical approach). In page 7, the author proposed a model to describe the growth of nanowires. I’m trying to write a WinBUGS code (Please see below), but I’m not sure whether it’s correct.

I quoted this from the paper as a summary:

The case studies will demonstrate the procedure of hierarchical modeling and estimation of the NW growth process under uncertainties. ... The data therein were collected over time (t =15 s, 30 s, 180 s, 900 s) under six growth conditions (T = 365 ◦C, 380 ◦C, 400 ◦C, 420 ◦C, 430 ◦C, 440 ◦C). The first four conditions are used for model building. We do not consider the two high-temperature conditions because there are no observations at times 180 s and 900 s. “

So we’re dealing with three variables, including Time (in seconds), Temperature (in C), and Length of nanowires (in cm). There are four levels of Temperature (365 ◦C, 380 ◦C, 400 ◦C, and 420 ◦C), and four time points (15 s, 30 s, 180 s, 900 s). I don’t have the exact values of length under each condition, so I generated some values by looking at Figure 6 in the paper.

Here is the WinBUGS code. It runs perfectly, but I’m not sure if it’s correct. I appreciate any input. Thanks!

model
{
    for(j in 1 : 4) {
        for(i in 1 : 4) {
            X[i,j] ~ dnorm(mu[i,j], tau)
            mu[i,j] <- (step(t0[j] - t[i]))  * ((a1[j]*(exp(-alpha2[j]/t[i])))) + (step(t[i] - t0[j])) * ((a3[j]*t[i]) + (a1[j]*(exp(-alpha2[j]/t0[j]))) -a3[j]*t0[j])
        }

        t0[j] ~ dunif(1,40)
        alpha2[j] ~ dunif(1,20)

        a1[j]~ dnorm(mua1[j], tau1)
        a3[j]~ dnorm(mua3[j], tau3)

        alpha1[j] <- exp(a11)
        alpha3[j] <- exp(a33)

        log(mua1[j]) <- log(alpha1[j]) - (230/(T[j]*1.3806503*pow(10, -23)))
        log(mua3[j]) <- log(alpha3[j]) - (230/(T[j]*1.3806503*pow(10, -23)))
    }

    a11~ dnorm(0,0.001)
    a33~ dnorm(0, 0.001)

    tau1 <- pow(sigma1, -2)
    sigma1 ~ dunif(0,10)

    tau3 <- pow(sigma3, -2)
    sigma3 ~ dunif(0,10)

    tau <- pow(sigma, -2)
    sigma ~ dunif(0,20)
}

list(X= structure(.Data= c(0.00011E+00, 0.00022E+00, 2.00000E-02, 3.50000E-02, 7.00000E-03, 7.00000E-03, 3.00000E-01, 3.00000E-01, 3.50000E-02, 7.00000E-02, 2.50000E-01, 5.00000E-01, 2.00000E-01, 2.25000E-01, 5.50000E-01, 2.30000E+00), .Dim=c(4, 4)), T=c(3.65000E+02, 3.80000E+02, 4.00000E+02, 4.20000E+02), t=c(1.50000E+01, 3.00000E+01, 1.80000E+02, 9.00000E+02))