model
{
    ###Priors and constraints for first latent state process###
    #Transition probabilities among states (logit scale)
    for (h in 1:3) {
        a2[h] ~ dnorm(0, 1)
        a3[h] ~ dnorm(0, 1)
        b2[h] ~ dnorm(0.5, 1)
        b3[h] ~ dnorm(-0.2, 100)
    }
    #Distribution of distances from the closest boat (mean and SD)
    phi[1] ~ dnorm(3, 1)  T(phi[3], )
    phi[2] <- phi[1]
    phi[3] ~ dnorm(1, 1)  T(0, )
    sigma[1] ~ dunif(0, 1)
    sigma[2] <- sigma[1]
    sigma[3] ~ dunif(0, sigma[1])
    for (h in 1:3) {
        tau[h] <- 1/sigma[h]/sigma[h] #BUGS uses precision instead of SD
    }

    ###Priors and constraints for second latent state process###
    #Bearing concentration parameter
    rho[1] ~ dunif(0.5, 1)
    rho[2] ~ dunif(0, rho[1])
    #Step length distribution
    logbeta[1] ~ dunif(-1, log.maxbeta)
    logbeta[2] ~ dunif(-1, log.maxbeta)
    logalpha[1] ~ dunif(-10, log.maxalpha)
    logalpha[2] ~ dunif(-10, loga[1])
    for (k in 1:2) {
        alpha[k] <- exp(logalpha[k])
        beta[k] <- exp(logbeta[k])
        w[k] <- 1/pow(alpha[k], beta[k]) #BUGS uses a different parameterisation for the Weibull distribution
    }
    #Transition probabilities among states
    for (k in 1:2) {
        delta[k, 1:2] ~ ddirch(phiprior[1:2])
    }
    Pin <- (-Pi)

    ###Model###
    #Loop over trips
    for (i in 1:no.trips) {

        #Prior of trip-specific concentration parameter for state 1
        epsi[i] ~ dunif(0, 1)

        #States at first time step
        s[1, i] <- 1
        r[1, i] <- 1

        #Counters to monitor state convergence
        s.cnt[1, 1, i] <- equals(s[1, i], 1)
        s.cnt[2, 1, i] <- equals(s[1, i], 2)
        s.cnt[3, 1, i] <- equals(s[1, i], 3)
        r.cnt[1, 1, i] <- equals(r[1, i], 1)
        r.cnt[2, 1, i] <- equals(r[1, i], 2)

        #Loop over time steps
        for (t in 2:(steps[i] - 1)) {

            #Elapsed proportion of the trip
            p[t, i] <- t/steps[i]

            #Linear predictors for transition probabilities
            lp[t, i, 1] <- a2[s[t - 1, i]] + b2[s[t - 1, i]] * p[t, i]
            lp[t, i, 2] <- a3[s[t - 1, i]] + b3[s[t - 1, i]] * d[t, i]

            #Transition probabilities (multinomial logit parameterisation)
            gamma[t, i, 2] <- (exp(lp[t, i, 1])/(1 + exp(lp[t, i, 1]) + exp(lp[t, i, 2])))
            gamma[t, i, 3] <- (exp(lp[t, i, 2])/(1 + exp(lp[t, i, 1]) + exp(lp[t, i, 2])))
            gamma[t, i, 1] <- 1 - gamma[t, i, 2] - gamma[t, i, 3]

            #Bearing to sea (mean and variability is trip specific)
            mu[t, i, 1] ~ dunif(Pin, Pi)
            #"Ones trick" (OpenBUGS manual) to draw mean bearing under state 1 from wrapped Cauchy distribution
            onesa[t, i] ~ dbern(wCa[t, i])
            wCa[t, i] <- (1/(2 * Pi) * (1 - epsi[i] * epsi[i])/(1 + epsi[i] * epsi[i] - 2 * epsi[i] * cos(mu[t, i, 1] - upsi[i])))/1000

            #Bearing to colony
            deltay[t, i] <- ny[1, i] - ny[t, i]
            deltax[t, i] <- nx[1, i] - nx[t, i]
            mu[t, i, 2] <- 2 * arctan(deltay[t, i]/(sqrt(pow(deltax[t, i], 2) + pow(deltay[t, i], 2)) + deltax[t, i]))

            #Bearing to closest boat
            deltayB[t, i] <- by[t, i] - ny[t, i]
            deltaxB[t, i] <- bx[t, i] - nx[t, i]
            mu[t, i, 3] <- 2 * arctan(deltayB[t, i]/(sqrt(pow(deltaxB[t, i], 2) + pow(deltayB[t, i], 2)) + deltaxB[t, i]))

            #Draw states
            s[t, i] ~ dcat(gamma[t, i, 1:3])        #First latent state process
            r[t, i] ~ dcat(delta[r[t - 1, i], 1:2]) #Second latent state process

            #Mean bearing (given state s)
            mu.f[t, i] <- mu[t, i, s[t, i]]

            #Distance from closest boat (given state s)
            d[t, i] ~ dlnorm(phi[s[t, i]], tau[s[t, i]])

            #Step length (given state r)
            x[t, i] ~ dweib(beta[r[t, i]], w[r[t, i]])

            #"Ones trick" (OpenBUGS manual) to draw observed bearing from wrapped Cauchy distribution
            ones[t, i] ~ dbern(wC[t, i])
            wC[t, i] <- (1/(2 * Pi) * (1 - rho[r[t, i]] * rho[r[t, i]])/(1 + rho[r[t, i]] * rho[r[t, i]] - 2 * rho[r[t, i]] * cos(theta[t, i] - mu.f[t, i])))/1000

            #Counters to monitor state convergence
            s.cnt[1, t, i] <- equals(s[t, i], 1)
            s.cnt[2, t, i] <- equals(s[t, i], 2)
            s.cnt[3, t, i] <- equals(s[t, i], 3)
            r.cnt[1, t, i] <- equals(r[t, i], 1)
            r.cnt[2, t, i] <- equals(r[t, i], 2)
        }

        #Counters to monitor state convergence
        sf.cnt[1, i] <- sum(s.cnt[1, 1:(steps[i] - 1), i])
        sf.cnt[2, i] <- sum(s.cnt[2, 1:(steps[i] - 1), i])
        sf.cnt[3, i] <- sum(s.cnt[3, 1:(steps[i] - 1), i])
        rf.cnt[1, i] <- sum(r.cnt[1, 1:(steps[i] - 1), i])
        rf.cnt[2, i] <- sum(r.cnt[2, 1:(steps[i] - 1), i])
    }

    #Proportion of states (to monitor state convergence)
    lambda1[1] <- sum(sf.cnt[1, ])/(sum(steps[]) - no.trips)
    lambda1[2] <- sum(sf.cnt[2, ])/(sum(steps[]) - no.trips)
    lambda1[3] <- sum(sf.cnt[3, ])/(sum(steps[]) - no.trips)
    lambda2[1] <- sum(rf.cnt[1, ])/(sum(steps[]) - no.trips)
    lambda2[2] <- sum(rf.cnt[2, ])/(sum(steps[]) - no.trips)
}