 #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>

    }
};


int main(void)
{
    /* TODO: Should use the boost random number generator libraries for doing this on the host */

    // Random number generator and distribution needed to simulate some data
	// unsigned int seed = 98724732;
    unsigned int seed = 123455;
	static thrust::minstd_rand rng(seed);
	thrust::random::experimental::normal_distribution<double> snorm(0., 1.0);
    thrust::random::uniform_real_distribution<double> uniform(0.0, 1.0);

    /*
     Population level parameters: theta, where theta parameterizes the distribution of chi.

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

// Global random number generator and distributions for generating random numbers on the host. The random number generator used
// is the Mersenne Twister mt19937 from the BOOST library.

boost::random::mt19937 rng;
boost::random::normal_distribution<> snorm(0.0, 1.0); // Standard normal distribution
boost::random::uniform_real_distribution<> uniform(0.0, 1.0); // Uniform distribution from 0.0 to 1.0

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