 
 // Pointer to the current value of theta. This is stored in constant memory so that all the threads on the GPU
 // can access the same theta quickly.
-static __constant__ double c_theta[6];
+static __constant__ double c_theta[dim_theta];
 
 // Constants used for the model SED
 static __constant__ double c_nu0 = 2.3e11;  // nu0 = 230 GHz

    
    // Cholesky factor of proposal distribution for each data point, used by the Metropolis algorithm. The proposal
    // covariance matrix for each data point is Covar_i = PropChol_i * transpose(PropChol_i).
    dim_cholfactor = p * p - ((p - 1) * p) / 2  // cholesky factor is a lower-diagonal matrix
    thrust::host_vector<double> h_prop_cholfact(n * dim_chol_factor);
    thrust::fill(h_prop_cholfact.begin(), h_prop_cholfact.end(), 0.0);
    

    curand_init(1234, id, 0, &state[id]);
}

// Perform the MHA update on the n characteristics, done in parallel on the GPU.
__global__
void update_chi(double* chi, double* meas, int m, double* meas_unc, int n, double* logdens, double* prop_cholfact,
                curandState* state, int current_iter)
{
	int i = blockDim.x * blockIdx.x + threadIdx.x;
	if (i < n)
	{
        // Copy variables to local memory for efficiency
        curandState localState = state[i];
        dim_cholfactor = c_p * c_p - ((c_p - 1) * c_p) / 2;  // cholesky factor is a lower-diagonal matrix
        double* cholfactor[dim_cholfactor];
        for (int j=0; j<dim_cholfactor; j++) {
            cholfactor[j] = prop_cholfact[n * j + i];
        }

        // get the unit proposal
        double* snorm_deviate[c_p];
        for (int j=0; j<c_p; j++) {
            snorm_deviate[j] = curand_normal_double(&localState);
        }
        
        // propose a new chi value
        double* new_chi[c_p];
        cholfact_index = 0
        for (int j=0; j<c_p; j++) {
            double centered_proposal = 0.0;
            for (int k=0; k<(j+1); k++) {
                // transform the unit proposal to the centered proposal, drawn from a multivariate normal
                // or t-distribution.
                centered_proposal += cholfactor[n * cholfact_index + i] * snorm_deviate[j];
                cholfact_index++
            }
            new_chi[i] = chi[i] + centered_proposal;
        }

        // Compute the conditional log-posterior of the proposed parameter values for this data point
        double centered_meas = meas_i - new_chi;
