 
     // Cuda grid launch
     dim3 nThreads(256);
-    dim3 nBlocks((n + nThreads.x-1) / nThreads.x);
+    dim3 nBlocks((ndata + nThreads.x-1) / nThreads.x);
     printf("nBlocks: %d\n", nBlocks.x);  // no more than 64k blocks!
     if (nBlocks.x > 65535)
     {

// Constants used for the model SED
static __constant__ double c_nu0 = 2.3e11;  // reference frequency, nu0 = 230 GHz
static const double nu0 = 2.3e11;
static const int m = 5;  // the number of measured features per data point
static __constant__ int c_m = 5;
    // observational frequencies, corresponding to the Herschel PACS and SPIRES bands
static const double nu[m] = {4.286e12, 1.765e12, 1.200e12, 8.571e11, 6.000e11};
static __constant__ double c_nu[m] = {4.286e12, 1.765e12, 1.200e12, 8.571e11, 6.000e11};
static const double hplanck = 6.6260755e-27;  // Planck's constant, in cgs
static const double clight = 2.997925e10;
static const double kboltz = 1.380658e-16;

        double* p_prop_cholfact = thrust::raw_pointer_cast(&d_prop_cholfact[0]);
        double* p_logdens = thrust::raw_pointer_cast(&d_logdens[0]);
        int current_iter = i + 1;
        update_chi<<<nBlocks,nThreads>>>(p_chi, p_meas, p_meas_unc, ndata, p_logdens, p_prop_cholfact, devStates, current_iter);

        // Generate new theta in parallel with GPU calculation above
        thrust::host_vector<double> proposed_theta(dim_theta);
