// A simple fluid solver.. perhaps extended with some other stuff (hopefully).
// Martin Michelsen
// Caltech CS 101c, spring 2010

#include <stdlib.h>
#include <stdio.h>
#include <string.h>

#include <GL/glew.h>

#if defined (__APPLE__) || defined(MACOSX)
#include <GLUT/glut.h>
#else
#include <GL/glut.h>
#endif

#include <cuda_gl_interop.h>

#include "lab8_kernel.cu"

////////////////////////////////////////////////////////////////////////////////
// first, some error handling functions

inline void _safecall(cudaError err, const char *file, const unsigned int line)
{
    if (err != cudaSuccess)
    {
        fprintf(stderr, "at %s/%u: CUDA runtime error: %d/%s\n",
                file, line, err, cudaGetErrorString(err));
        exit(-1);
    }
}
inline void _checkmsg(const char *msg, const char *file, const unsigned int line)
{
    cudaError_t err = cudaGetLastError();
    if (err != cudaSuccess)
    {
        fprintf(stderr, "at %s/%u: %s: %d/%s\n",
                file, line, msg, err, cudaGetErrorString(err));
        exit(-1);
    }
}
#define safecall(a) _safecall(a, __FILE__, __LINE__)
#define checkmsg(a) _checkmsg(a, __FILE__, __LINE__)

#define iDivUp(a, b) ((a % b != 0) ? (a / b + 1) : (a / b))

////////////////////////////////////////////////////////////////////////////////
// simulation setup parameters

// these are default parameters. they can be changed (TODO) via the command line.
long fieldDimension = 128; // dimension of simulation grid (it's a cube)
long windowWidth = 512; // dimensions of viewing window
long windowHeight = 512;
long numParticles = 0; // set this using the cmd line...
float dt = 0.001f; // step in time (try different values)

// view settings (from lab 6)
float3 viewRotation = make_float3(20.0f, 11.0f, 2.0f);//0.5, 0.5, 0.0);
float3 viewTranslation = make_float3(0.0, 0.0, -8.0f);
int renderFlags = DEFAULT_RENDER_FLAGS;
float transparencyFactor = 17.0f;

// some program parameters
int benchmark = 0;
char videoPrefix[64] = "";
long frameNum = 0;

// pointers to texture arrays
cudaArray *dVelocity = 0;
cudaArray *dPressure = 0;

// pointers to buffers for new computed data (after each iteration)
float4 *dNewVelocity = NULL;
float *dNewPressure = NULL;
float4 *dParticles = NULL;

// OpenGL pixel buffer
GLuint pbo = 0;

// memcpy params for moving new data into textures for next iteration
cudaExtent volumeSize;
cudaMemcpy3DParms dVelocityParams = {0};
cudaMemcpy3DParms dPressureParams = {0};

////////////////////////////////////////////////////////////////////////////////
// update & display functions

// compute a fluid step
void updateFluid()
{
    // grid dimensions for field computations
    dim3 blockSize(fieldDimension, 1, 1);
    dim3 gridSize = dim3(fieldDimension, fieldDimension, 1);

    cudaEvent_t start, stop;
    if (benchmark)
    {
        // start timing...
        cudaEventCreate(&start);
        cudaEventCreate(&stop);
        cudaEventRecord(start, 0);
    }

    // transport the velocity field
    VelocityTransportK<<< gridSize, blockSize >>>(dNewVelocity, dt);
    checkmsg("VelocityTransportK failed");
    safecall(cudaMemcpy3D(&dVelocityParams));

    // recalculate pressures
    PressureRecalculateK<<< gridSize, blockSize >>>(dNewPressure);
    checkmsg("PressureRecalculateK failed");
    safecall(cudaMemcpy3D(&dPressureParams));

    // project the velocity onto its divergence-free component
    VelocityProjectK<<< gridSize, blockSize >>>(dNewVelocity);
    checkmsg("VelocityProjectK failed");
    safecall(cudaMemcpy3D(&dVelocityParams));

    // transport particles (if there are any particles)
    if (numParticles)
    {
        dim3 particleBlockSize(128, 1, 1);
        dim3 particleGridSize(numParticles / particleBlockSize.x, 1, 1);
        ParticlesTransportK<<< particleGridSize, particleBlockSize >>>(dParticles, dt, ((float)rand())/RAND_MAX);
        checkmsg("ParticlesTransportK failed");
    }

    if (benchmark)
    {
        // done timing...
        cudaEventRecord(stop, 0);
        cudaEventSynchronize(stop);
        float elapsedTime;
        cudaEventElapsedTime(&elapsedTime, start, stop);
        printf("Recompute took %f ms\n", elapsedTime);
        cudaEventDestroy(start);
        cudaEventDestroy(stop);
    }
}

// save frame to a numbered file
void saveFrame(void* deviceData)
{
    char filename[64];
    sprintf(filename, "%s%d.ppm", videoPrefix, frameNum);

    int size = windowWidth * windowHeight * sizeof(GLubyte) * 4;

    FILE* f = fopen(filename, "wt");
    if (!f) printf("failed to save frame: couldn\'t open file %s\n", filename);
    else {
        long* frameData = (long*)malloc(size);
        if (!frameData)
        {
            printf("failed to save frame: couldn\'t allocate frame data in main memory\n");
            fclose(f);
            return;
        }
        safecall(cudaMemcpy(frameData, deviceData, size, cudaMemcpyDeviceToHost));
        fprintf(f, "P6 %d %d 255\n", windowWidth, windowHeight);
        long x, y;
        for (y = windowHeight - 1; y >= 0; y--)
        {
            for (x = 0; x < windowWidth; x++)
            {
                fwrite(&frameData[y * windowWidth + x], 1, 3, f);
            }
        }
        fclose(f);
        free(frameData);
        frameNum++;
    }
}

// display results using OpenGL (called by GLUT)
void display()
{
    // do a time-step in the fluid...
    updateFluid();

    // use OpenGL to build view matrix
    GLfloat modelView[16];
    glMatrixMode(GL_MODELVIEW);
    glPushMatrix();
    glLoadIdentity();
    glRotatef(-viewRotation.x, 1.0f, 0.0f, 0.0f);
    glRotatef(-viewRotation.y, 0.0f, 1.0f, 0.0f);
    glTranslatef(-viewTranslation.x, -viewTranslation.y, -viewTranslation.z);
    glGetFloatv(GL_MODELVIEW_MATRIX, modelView);
    glPopMatrix();

    // transpose matrix to conform with OpenGL's notation
    GLfloat invViewMatrix[12];
    invViewMatrix[0] = modelView[0];
    invViewMatrix[1] = modelView[4];
    invViewMatrix[2] = modelView[8];
    invViewMatrix[3] = modelView[12];
    invViewMatrix[4] = modelView[1];
    invViewMatrix[5] = modelView[5];
    invViewMatrix[6] = modelView[9];
    invViewMatrix[7] = modelView[13];
    invViewMatrix[8] = modelView[2];
    invViewMatrix[9] = modelView[6];
    invViewMatrix[10] = modelView[10];
    invViewMatrix[11] = modelView[14];

    // set necessary constants in hardware
    safecall(cudaMemcpyToSymbol(dInvViewMatrix, invViewMatrix, sizeof(float4) * 3));
    safecall(cudaMemcpyToSymbol(dViewMatrix, modelView, sizeof(float4) * 3));

    // map PBO to get CUDA device pointer
    unsigned long* renderOutput;
    safecall(cudaGLMapBufferObject((void**)&renderOutput, pbo));

    // call CUDA kernel, writing results to PBO
    if (renderFlags & RENDER_FOG)
    {
        dim3 blockSize(32, 8, 1);
        dim3 gridSize = dim3(iDivUp(windowWidth, blockSize.x), iDivUp(windowHeight, blockSize.y), 1);
        RenderParameterFogK<<< gridSize, blockSize >>>(renderOutput, windowWidth, windowHeight);
        checkmsg("RenderParameterFogK failed");
    } else safecall(cudaMemset(renderOutput, 0, windowWidth * windowHeight * sizeof(long)));
    if (numParticles)
    {
        dim3 particleBlockSize(64, 1, 1);
        dim3 particleGridSize(numParticles / particleBlockSize.x, 1, 1);
        RenderParticlesK<<< particleGridSize, particleBlockSize >>>(renderOutput, dParticles, windowWidth, windowHeight);
        checkmsg("RenderParticlesK failed");
    }
    if (videoPrefix[0]) saveFrame(renderOutput);
    safecall(cudaGLUnmapBufferObject(pbo));

    // display results
    glClear(GL_COLOR_BUFFER_BIT);

    // draw image from PBO
    glDisable(GL_DEPTH_TEST);
    glRasterPos2i(0, 0);
    glBindBufferARB(GL_PIXEL_UNPACK_BUFFER_ARB, pbo);
    glDrawPixels(windowWidth, windowHeight, GL_RGBA, GL_UNSIGNED_BYTE, 0);
    glBindBufferARB(GL_PIXEL_UNPACK_BUFFER_ARB, 0);

    glutSwapBuffers();
    glutReportErrors();

    glutPostRedisplay();
}

////////////////////////////////////////////////////////////////////////////////
// initialization functions

// initializes the contents of the velocity and pressure fields
void initFieldContents()
{
    // get some local storage to fill with random crap
    float4* initialVelocityData = (float4*)malloc(sizeof(float4) * fieldDimension * fieldDimension * fieldDimension);
    float* initialPressureData = (float*)malloc(sizeof(float) * fieldDimension * fieldDimension * fieldDimension);
    if (!initialVelocityData || !initialPressureData)
    {
        printf("can\' initialize field contents: out of memory\n");
        exit(-1);
    }

    // velocities are all initially 0
    memset(initialVelocityData, 0, sizeof(float4) * fieldDimension * fieldDimension * fieldDimension);

    // initialize pressure data
    float dx, dy, dz, intensity;
    int w, x, y, z;
    int centersX[4] = {rand() % fieldDimension,
                       rand() % fieldDimension,
                       rand() % fieldDimension,
                       rand() % fieldDimension};
    int centersY[4] = {rand() % fieldDimension,
                       rand() % fieldDimension,
                       rand() % fieldDimension,
                       rand() % fieldDimension};
    int centersZ[4] = {rand() % fieldDimension,
                       rand() % fieldDimension,
                       rand() % fieldDimension,
                       rand() % fieldDimension};
    for (w = 0; w < 4; w++)
    {
        for (x = 0; x < fieldDimension; x++)
        {
            dx = x - centersX[w];
            for (y = 0; y < fieldDimension; y++)
            {
                dy = y - centersY[w];
                for (z = 0; z < fieldDimension; z++)
                {
                    dz = z - centersZ[w];
                    intensity = 10.0f * exp(-((dx * dx) + (dy * dy) + (dz * dz)) / (fieldDimension * 4));
                    initialPressureData[z * fieldDimension * fieldDimension + y * fieldDimension + x] +=
                        intensity;
                }
            }
        }
    }

    // copy this data to the card
    cudaMemcpy3DParms params;
    memset(&params, 0, sizeof(cudaMemcpy3DParms));
    params.srcPtr   = make_cudaPitchedPtr(initialPressureData, fieldDimension * sizeof(float), fieldDimension, fieldDimension);
    params.dstArray = dPressure;
    params.extent   = volumeSize;
    params.kind     = cudaMemcpyHostToDevice;
    safecall(cudaMemcpy3D(&params));

    // and copy that data to the card
    params.srcPtr   = make_cudaPitchedPtr(initialVelocityData, fieldDimension * sizeof(float4), fieldDimension, fieldDimension);
    params.dstArray = dVelocity;
    safecall(cudaMemcpy3D(&params));

    // initialize particles, if present
    if (numParticles)
    {
        // get local storage
        float4* initialParticleData = (float4*)malloc(numParticles * sizeof(float4));

        // make random data (particles start in the cube of dimension 0..1 on all sides)
        for (x = 0; x < numParticles; x++)
        {
            initialParticleData[x].x = ((float)rand())/RAND_MAX;
            initialParticleData[x].y = ((float)rand())/RAND_MAX;
            initialParticleData[x].z = ((float)rand())/RAND_MAX;
            initialParticleData[x].w = 1.0f;
        }

        // copy data to graphics card and clean up
        safecall(cudaMemcpy(dParticles, initialParticleData, numParticles * sizeof(float4), cudaMemcpyHostToDevice));
        free(initialParticleData);
    }

    free(initialVelocityData);
    free(initialPressureData);
}

// initializes the memory and textures used by the CUDA kernels
void initCuda()
{
    // create 3D cudaArrays for the velocity and pressure textures
    volumeSize = make_cudaExtent(fieldDimension, fieldDimension, fieldDimension);
    cudaChannelFormatDesc channelDescVelocity = cudaCreateChannelDesc<float4>();
    cudaChannelFormatDesc channelDescPressure = cudaCreateChannelDesc<float>();
    safecall(cudaMalloc3DArray(&dVelocity, &channelDescVelocity, volumeSize));
    safecall(cudaMalloc3DArray(&dPressure, &channelDescPressure, volumeSize));

    // create writeable 3D data arrays in global memory for updating
    safecall(cudaMalloc((void**)&dNewVelocity, sizeof(float4) * fieldDimension * fieldDimension * fieldDimension));
    safecall(cudaMalloc((void**)&dNewPressure, sizeof(float) * fieldDimension * fieldDimension * fieldDimension));
    if (numParticles) safecall(cudaMalloc((void**)&dParticles, numParticles * sizeof(float4)));

    // initialize contents of fields
    initFieldContents();

    // set up copy params for future copying
    memset(&dVelocityParams, 0, sizeof(cudaMemcpy3DParms));
    memset(&dPressureParams, 0, sizeof(cudaMemcpy3DParms));
    dVelocityParams.srcPtr   = make_cudaPitchedPtr(dNewVelocity, fieldDimension * sizeof(float4), fieldDimension, fieldDimension);
    dVelocityParams.dstArray = dVelocity;
    dVelocityParams.extent   = volumeSize;
    dVelocityParams.kind     = cudaMemcpyDeviceToDevice;
    dPressureParams.srcPtr   = make_cudaPitchedPtr(dNewPressure, fieldDimension * sizeof(float), fieldDimension, fieldDimension);
    dPressureParams.dstArray = dPressure;
    dPressureParams.extent   = volumeSize;
    dPressureParams.kind     = cudaMemcpyDeviceToDevice;

    // set texture parameters
    velocityTex.normalized = true;
    velocityTex.filterMode = cudaFilterModeLinear;
    velocityTex.addressMode[0] = cudaAddressModeClamp;
    velocityTex.addressMode[1] = cudaAddressModeClamp;
    velocityTex.addressMode[2] = cudaAddressModeClamp;
    pressureTex.normalized = true;
    pressureTex.filterMode = cudaFilterModeLinear;
    pressureTex.addressMode[0] = cudaAddressModeClamp;
    pressureTex.addressMode[1] = cudaAddressModeClamp;
    pressureTex.addressMode[2] = cudaAddressModeClamp;

    // bind arrays to 3D textures
    safecall(cudaBindTextureToArray(velocityTex, dVelocity, channelDescVelocity));
    safecall(cudaBindTextureToArray(pressureTex, dPressure, channelDescPressure));
}

// deletes all the stuff we've used
void cleanup()
{
    // don't need to check if these are null since atexit is called after they've
    // all been allocated
    safecall(cudaFree(dNewVelocity));
    safecall(cudaFree(dNewPressure));
    safecall(cudaFreeArray(dVelocity));
    safecall(cudaFreeArray(dPressure));
    safecall(cudaGLUnregisterBufferObject(pbo));    
    glDeleteBuffersARB(1, &pbo);
    if (numParticles) safecall(cudaFree(dParticles));
}

// initialize some parameters we use for rendering (there used to be more here)
void initGLParameters()
{
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    glOrtho(0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
}

// initialize (or re-initialize) the pixel buffer
void initPixelBuffer()
{
    // if buffer already exists, get rid of it
    if (pbo)
    {
        safecall(cudaGLUnregisterBufferObject(pbo));
        glDeleteBuffersARB(1, &pbo);
    }

    // create pixel buffer object for display
    glGenBuffersARB(1, &pbo);
    glBindBufferARB(GL_PIXEL_UNPACK_BUFFER_ARB, pbo);
    glBufferDataARB(GL_PIXEL_UNPACK_BUFFER_ARB, windowWidth * windowHeight * sizeof(GLubyte) * 4, 0, GL_STREAM_DRAW_ARB);
    glBindBufferARB(GL_PIXEL_UNPACK_BUFFER_ARB, 0);
    safecall(cudaGLRegisterBufferObject(pbo));
}

// set constants on device
void setConsts()
{
    safecall(cudaMemcpyToSymbol(dRenderFlags, &renderFlags, sizeof(int)));
    safecall(cudaMemcpyToSymbol(dFieldDimension, &fieldDimension, sizeof(int)));
    safecall(cudaMemcpyToSymbol(dTransparencyFactor, &transparencyFactor, sizeof(float)));
}

////////////////////////////////////////////////////////////////////////////////
// user interaction functions

void idle()
{
    // do nothing
}

void keyboard(unsigned char key, int x, int y)
{
    switch (key)
    {
      case 27:
        // quit
        exit(0);
        break;

      case 'c':
        // change color of particles
        renderFlags = (renderFlags & ~RENDER_COLOR_MASK) | ((renderFlags + RENDER_COLOR_INCREMENT) & RENDER_COLOR_MASK);
        if ((renderFlags & RENDER_COLOR_MASK) == RENDER_COLOR_VELOCITY) printf("color function: velocity\n");
        if ((renderFlags & RENDER_COLOR_MASK) == RENDER_COLOR_VEL_SPECT) printf("color function: velocity spectrum\n");
        if ((renderFlags & RENDER_COLOR_MASK) == RENDER_COLOR_POS_SPECT) printf("color function: position spectrum\n");
        if ((renderFlags & RENDER_COLOR_MASK) == RENDER_COLOR_GREEN) printf("color function: green\n");
        break;

      case 'f':
        // enable/disable fog
        renderFlags ^= RENDER_FOG;
        printf("pressure and velocity fog: %s\n", (renderFlags & RENDER_FOG) ? "enabled" : "disabled");
        break;

      case 'i':
        // switch view projection
        renderFlags ^= RENDER_ISOMETRIC;
        printf("projection mode: %s\n", (renderFlags & RENDER_ISOMETRIC) ? "isometric" : "perspective");
        break;

      case 'l':
        // change where particles are "recreated"
        renderFlags ^= RENDER_RAND_LOC;
        printf("particles start: %s\n", (renderFlags & RENDER_RAND_LOC) ? "somewhere in the field" : "at one of the walls");
        break;
    }

    setConsts();
    glutPostRedisplay();
}

// click locations and button state (for tracking mouse motion)
int ox, oy;
int buttonState = 0;
int moveMode = 0;
// move modes: rotate, translate, zoom, force

void mouse(int button, int state, int x, int y)
{
    if (state == GLUT_DOWN) buttonState |= 1 << button;
    else if (state == GLUT_UP) buttonState = 0;

    ox = x;
    oy = y;
    glutPostRedisplay();
}

void motion(int x, int y)
{
    float dx, dy;
    dx = x - ox;
    dy = y - oy;

    if (buttonState == 3) {
        // left+middle = zoom
        viewTranslation.z += dy / 100.0;
    }
    else if (buttonState & 2) {
        // middle = translate
        viewTranslation.x += dx / 100.0;
        viewTranslation.y -= dy / 100.0;
    }
    else if (buttonState & 1) {
        // left = rotate
        viewRotation.x += dy / 5.0;
        viewRotation.y += dx / 5.0;
    }

    ox = x;
    oy = y;
    glutPostRedisplay();
}

void reshape(int x, int y)
{
    // update windowWidth, windowHeight
    windowWidth = x;
    windowHeight = y;

    // reinitialize display state with new size
    initPixelBuffer();

    glViewport(0, 0, windowWidth, windowHeight);
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
}

////////////////////////////////////////////////////////////////////////////////
// Program main

// This function returns the best GPU (with maximum GFLOPS)
// (taken from cutil_inline_runtime.h)
int cutGetMaxGflopsDeviceId()
{
    int current_device   = 0, sm_per_multiproc = 0;
    int max_compute_perf = 0, max_perf_device  = 0;
    int device_count     = 0, best_SM_arch     = 0;
    int arch_cores_sm[3] = {1, 8, 32};
    cudaDeviceProp deviceProp;

    cudaGetDeviceCount(&device_count);

    // Find the best major SM Architecture GPU device
    while (current_device < device_count)
    {
        cudaGetDeviceProperties(&deviceProp, current_device);
        if (deviceProp.major > 0 && deviceProp.major < 9999)
            best_SM_arch = max(best_SM_arch, deviceProp.major);
        current_device++;
    }

    // Find the best CUDA capable GPU device
    for (current_device = 0; current_device < device_count; current_device++)
    {
        cudaGetDeviceProperties(&deviceProp, current_device);
        if (deviceProp.major == 9999 && deviceProp.minor == 9999)
            sm_per_multiproc = 1;
        else if (deviceProp.major <= 2)
            sm_per_multiproc = arch_cores_sm[deviceProp.major];
        else
            sm_per_multiproc = arch_cores_sm[2];

        int compute_perf = deviceProp.multiProcessorCount * sm_per_multiproc * deviceProp.clockRate;
        if (compute_perf > max_compute_perf)
        {
            // If we find GPU with SM major > 2, search only these
            if (best_SM_arch > 2)
            {
                // If our device==dest_SM_arch, choose this, or else pass
                if (deviceProp.major == best_SM_arch)
                {
                    max_compute_perf = compute_perf;
                    max_perf_device = current_device;
                }
            } else {
                max_compute_perf = compute_perf;
                max_perf_device = current_device;
            }
        }
    }
    return max_perf_device;
}

// (taken from cutil_inline_runtime.h)
void cutilDeviceInit(int deviceID)
{
    int deviceCount;
    safecall(cudaGetDeviceCount(&deviceCount));
    if (deviceCount == 0)
    {
        fprintf(stderr, "error: no devices supporting CUDA\n");
        exit(-1);
    }

    if (deviceID < 0) deviceID = 0;
    if (deviceID > deviceID - 1) deviceID = deviceCount - 1;
    cudaDeviceProp deviceProp;
    safecall(cudaGetDeviceProperties(&deviceProp, deviceID));
    if (deviceProp.major < 1)
    {
        fprintf(stderr, "error: GPU device does not support CUDA\n");
        exit(-1);
    }
    printf("using CUDA device %d: %s\n", deviceID, deviceProp.name);
    safecall(cudaSetDevice(deviceID));
}

int main(int argc, char* argv[])
{
    int deviceID = -1;

    // parse arguments
    for (int x = 1; x < argc; x++)
    {
        if (argv[x][0] != '-') continue;
        switch (argv[x][1])
        {
          case 'b':
          case 'B':
            benchmark = 1;
            break;
          case 's':
          case 'S':
            fieldDimension = atoi(&argv[x][2]);
            break;
          case 'p':
          case 'P':
            numParticles = atoi(&argv[x][2]);
            renderFlags &= ~RENDER_FOG; // turn off fog if particles are on
            break;
          case 'd':
          case 'D':
            deviceID = atoi(&argv[x][2]);
            break;
          case 'h':
          case 'H':
            windowHeight = atoi(&argv[x][2]);
            break;
          case 'w':
          case 'W':
            windowWidth = atoi(&argv[x][2]);
            break;
          case 'r':
          case 'R':
            sscanf(&argv[x][2], "%X", &renderFlags);
            break;
          case 't':
          case 'T':
            sscanf(&argv[x][2], "%f", &dt);
            break;
          case 'v':
          case 'V':
            strcpy(videoPrefix, &argv[x][2]);
            break;
          case 'z':
          case 'Z':
            sscanf(&argv[x][2], "%f", &transparencyFactor);
            break;
        }
    }

    // use command-line specified CUDA device, otherwise use device with highest Gflops/s
    if (deviceID != -1) cutilDeviceInit(deviceID);
    else {
        deviceID = cutGetMaxGflopsDeviceId();
        cudaSetDevice(deviceID);
    }

    struct cudaDeviceProp prop;
    cudaGetDeviceProperties(&prop, deviceID);
    if (!strncmp("Tesla", prop.name, 5))
    {
        printf("This sample needs a card capable of OpenGL and display.\n");
        printf("Please choose a different device with the -device=x argument.\n");
        exit(0);
    }

    srand(time(NULL));

    // initialize GLUT callback functions
    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE);
    glutInitWindowSize(windowWidth, windowHeight);
    glutCreateWindow("fluids");
    glutDisplayFunc(display);
    glutKeyboardFunc(keyboard);
    glutMouseFunc(mouse);
    glutMotionFunc(motion);
    glutReshapeFunc(reshape);
    glutIdleFunc(idle);

    glewInit();
    if (!glewIsSupported("GL_VERSION_2_0 GL_ARB_pixel_buffer_object"))
    {
        fprintf(stderr, "Required OpenGL extensions missing.");
        exit(-1);
    }

    // set up all the necessary data
    initCuda();
    initGLParameters();
    initPixelBuffer();
    setConsts();

    // give away our keys before we start drinking
    atexit(cleanup);

    // woo, let's go!
    glutMainLoop();

    cudaThreadExit();
    return 0;
}
