     double sigma_msmt[5] = {2.2e-4, 3.3e-4, 5.2e-4, 3.7e-4, 2.2e-4};  // Standard deviation for the measurement errors
 
     std::cout << "Generating some data..." << std::endl;
-
-	for (int i=0; i<ndata; i++) {
-        // Loop over the data indices
-        int diag_index = 0;
-		for (int j=0; j<p; j++) {
-            // Loop over the chi indices to generate true value of the characteristics
-			h_true_chi[i + ndata * j] = h_theta[j] + sqrt(exp(h_theta[p + j])) * snorm(rng);
-            // Initialize the covariance matrix of the chi proposal distribution to be 0.01 * identity(p)
-            h_prop_cholfact[i + ndata * diag_index] = 0.01;
-            diag_index += j + 2;
-        }
-        for (int k=0; k<mfeat; k++) {
-            // Loop over the feature indices to generate measurements, given the characteristics
-            double bbody_numer = 2.0 * hplanck * nu[k] * nu[k] * nu[k] / (clight * clight);
-            double temperature = exp(h_true_chi[ndata * 2 + i]); // Grab the temperature for this data point
-            double bbody_denom = exp(hplanck * nu[k] / (kboltz * temperature)) - 1.0;
-            double bbody = bbody_numer / bbody_denom;
-            // Grab the SED normalization and power-law index
-            double normalization = exp(h_true_chi[i]);
-            double beta = h_true_chi[ndata + i];
-            // Compute the model SED
-            double nu_over_nu0 = nu[k] / nu0;
-            double SED_ik = normalization * pow(nu_over_nu0, beta) * bbody;
-            // Generate measured SEDs
-			h_meas_unc[k * ndata + i] = sigma_msmt[k];
-            h_meas[k * ndata + i] = SED_ik + sigma_msmt[k] * snorm(rng);
-		}
-	}
 
     // Initialize the chis to their true values for now
     thrust::copy(h_true_chi.begin(), h_true_chi.end(), h_chi.begin());

static const double kboltz = 1.380658e-16;
static __constant__ double c_2h_over_csqr = 1.4745002e-47;  // 2 * hplanck / clight^2, used in the Planck function
static __constant__ double c_h_over_k = 4.7992157e-11;  // hplanck / kboltz

// Function to return the model SED at a input frequency, given the characteristics chi = (norm, beta, temp)
__device__

#include <iostream>
#include <math.h>
#include <vector>

// Boost includes
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/normal_distribution.hpp>
#include <boost/random/uniform_real_distribution.hpp>

//#include <thrust/sort.h>
//#include <thrust/copy.h>
