+// Cuda Includes
 #include "cuda_runtime.h"
 #include "device_launch_parameters.h"
-
 #include <cuda.h>
 #include <curand.h>
 #include <curand_kernel.h>
+
+// Thrust includes
 #include <thrust/device_vector.h>
 #include <thrust/host_vector.h>
-
 #include <thrust/random/linear_congruential_engine.h>
 #include <thrust/random/normal_distribution.h>
 #include <thrust/random.h>
 #include <thrust/reduce.h>
 #include <thrust/transform_reduce.h>
 
+// Standard includes
 #include <stdio.h>
 #include <fstream>
 #include <stdlib.h>
+#include <math.h>
+#include <vector>
 
 //#include <thrust/sort.h>
 //#include <thrust/copy.h>

//#include <algorithm>
//#include <cstdlib>

#include <vector>
/*
 typedef thrust::device_vector<double> dvec;
 typedef std::vector<dvec> vdvec;
 
 struct wrapvec
 {
 vdvec v;
 
 wrapvec(int m, int n) : v(m)
 {
 for (int i=0; i<m; i++)	{
 v[i].reserve(n);
 }
 }
 
 double** ptrs()
 {
 thrust::host_vector<double*> h_ptr(v.size());
 for (unsigned int i=0; i<v.size(); i++)
 h_ptr[i] = (double*) thrust::raw_pointer_cast(&(v[i][0]));
 
 thrust::device_vector<double*> d_ptr = h_ptr;
 return (double**) thrust::raw_pointer_cast(&d_ptr[0]);
 }
 };
 */
#define CUDA_CALL(x) do { if((x) != cudaSuccess) { \
printf("Error at %s:%d\n",__FILE__,__LINE__); \
return EXIT_FAILURE;}} while(0)


// 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];

// Constants used for the model SED
static __constant__ double c_nu0 = 2.3e11;  // nu0 = 230 GHz
static const double nu0 = 2.3e11;
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;
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__
double modified_bbody(double nu, double normalization, double beta, double temperature)
{
    double bbody_numer = c_2h_over_csqr * nu * nu * nu;
    double bbody_denom = exp(c_h_over_k * nu / temperature) - 1.0;
    double bbody = bbody_numer / bbody_denom;
    double nu_over_nu0 = nu / c_nu0;
    double SED = normalization * pow(nu_over_nu0, beta) * bbody;
    return SED;
}


// Initialize the random number generator state
__global__
void setup_kernel(curandState *state)
{
    int id = threadIdx.x + blockIdx.x * blockDim.x;
    /* Each thread gets same seed, a different sequence
