#ifdef __APPLE__
#include <GLUT/glut.h>
#else
#include <GL/glut.h>
#endif

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

#define WINDOW_WIDTH    70
#define WINDOW_HEIGHT   70

#define LATTICE_WIDTH   50
#define LATTICE_HEIGHT  50
#define LATTICE_DEPTH   5

#define VECTOR_SPACING  2

#define POINT_SIZE      9
#define SPACING         1

#define Q               19
#define DIM             3

double omega = 2.0;

double lattice[LATTICE_DEPTH][LATTICE_HEIGHT][LATTICE_WIDTH][Q];
double tmp_lattice[LATTICE_DEPTH][LATTICE_HEIGHT][LATTICE_WIDTH][Q];
int e[Q][DIM] = {{0,0,0}, {1, 0, 0}, {-1, 0, 0}, {0, 1, 0}, {0, -1, 0},
                          {0, 0, 1}, {0, 0, -1}, 
                          {1, 1, 0}, {1, -1, 0}, {-1, 1, 0}, {-1, -1, 0}, 
                          {0, 1, 1}, {0, 1, -1}, {0, -1, 1}, {0, -1, -1}, 
                          {1, 0, 1}, {1, 0, -1}, {-1, 0, 1}, {-1, 0, -1}};
double weights[Q] = {1.0/3.0, 
                     1.0/18.0, 1.0/18.0, 1.0/18.0, 1.0/18.0, 1.0/18.0, 1.0/18.0,
                     1.0/18.0, 1.0/18.0, 1.0/18.0,
                     1.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0,
                     1.0/36.0, 1.0/36.0, 1.0/36.0};


void resolve_collisions( ) {
    int x,y,z, i,j;
    int rho;
    double u[DIM];
    double u_squared;
    double u_dot_e[Q];
    double f_eq;

    for(z = 0; z < LATTICE_DEPTH; z++) {
        for(x = 0; x < LATTICE_WIDTH; x++) {
            for(y = 0; y < LATTICE_HEIGHT; y++) {
                
                // Clear Variables
                rho = 0;
                u_squared = 0;
                u[0] = u[1] = u[2] = 0;
                for(i = 0; i < Q; i++) {
                    u_dot_e[i] = 0;
                }

                for(i = 0; i < Q; i++) {
                    rho += lattice[z][y][x][i];
                    for(j = 0; j < DIM; j++) {
                        u[j] += e[i][j] * lattice[z][y][x][i];
                    }
                }

                for(j = 0; j < DIM; j++) {
                    u[j] /= ((rho) + 0.1);
                    u_squared += u[j] * u[j];
                }

                for(i = 0; i < Q; i++) {
                    for(j = 0; j < DIM; j++) {
                        u_dot_e[i] += e[i][j] * u[j];
                    }
                }

                for(i = 0; i < Q; i++) {
                    f_eq = weights[i] * rho * (1 - 1.5 * u_squared + 3 * u_dot_e[i] + 
                                            4.5 * u_dot_e[i] * u_dot_e[i]);
                    tmp_lattice[z][y][x][i] = (1.0 - omega) * lattice[z][y][x][i] +
                                                omega * f_eq;
                }
            }
        }
    }
}

void propagate( ) {
    int x, y, z, i;
    int e_x, e_y, e_z;
    int new_x, new_y, new_z;

    for(z = 0; z < LATTICE_DEPTH; z++) {
        for(x = 0; x < LATTICE_WIDTH; x++) {
            for(y = 0; y < LATTICE_HEIGHT; y++) {
                for(i = 0; i < Q; i++) {
                    e_x = e[i][0];
                    e_y = e[i][1];
                    e_z   = e[i][2];
                    new_x = (x + e_x + LATTICE_WIDTH)  % LATTICE_WIDTH;
                    new_y = (y + e_y + LATTICE_HEIGHT) % LATTICE_HEIGHT;
                    new_z = (z + e_z + LATTICE_DEPTH)  % LATTICE_DEPTH;
                    lattice[new_z][new_y][new_x][i] = tmp_lattice[z][y][x][i];
                }
            }
        }
    }
}

void next( ) {
    resolve_collisions( );
    propagate( );
}

void display( ) {
    int x, y, z, i, j;
    int v_x, v_y, v_z;
    double slope[DIM];

    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
    glColor3f(1.0, 1.0, 1.0);
    glBegin(GL_LINE);
        for(z = 0; z < LATTICE_DEPTH; z += VECTOR_SPACING) {
            for(x = 0; x < LATTICE_WIDTH; x += VECTOR_SPACING) {
                for(y = 0; y < LATTICE_HEIGHT; y += VECTOR_SPACING) {
                    for(j = 0; j < DIM; j++) {
                        slope[j] = 0.0;
                    }
                    for(v_z = z; v_z < z + VECTOR_SPACING; v_z++) {
                        for(v_x = x; v_x < x + VECTOR_SPACING; v_x++) {
                            for(v_y = y; v_y < y + VECTOR_SPACING; v_y++) {
                                for(i = 0; i < Q; i++) {
                                    for(j = 0; j < DIM; j++) {
                                        slope[j] += e[i][j] * lattice[v_z][v_y][v_x][i];
                                    }
                                }
                            }
                        }
                    }

                    for(j = 0; j < DIM; j++) {
                        slope[j] /= (Q * VECTOR_SPACING * VECTOR_SPACING);
                    }
                    glVertex3f(x, y, z);
                    glVertex3f(x + slope[0], y + slope[1], z + slope[2]);
                }
            }
        }
    glEnd( );
/*    glutSolidSphere(0.1, 100, 80);*/
    glutSwapBuffers( );
}

void idle( ) {
    next( );
    glutPostRedisplay();
}

void init_vars( )
{
    int x,y,z, i;

    for(x = 0; x < LATTICE_WIDTH; x++) {
        for(y = 0; y < LATTICE_HEIGHT; y++) {
            for(z = 0; z < LATTICE_DEPTH; z++) {
                for(i = 0; i < Q; i++) {
                    lattice[z][y][x][i] = 0;
                    tmp_lattice[z][y][x][i] = 0;
                }
            }
        }
    }

    for(x = LATTICE_HEIGHT/3; x < 3*LATTICE_WIDTH/4; x++) {
        for(y = LATTICE_WIDTH/3; y < 3*LATTICE_HEIGHT/4; y++) {
            for(z = 0; z < LATTICE_DEPTH; z++) {
                for(i = 0; i < Q/2; i++) {
                    lattice[z][y][x][i] = 1;
                    tmp_lattice[z][y][x][i] = 1;
                }
            }
        }
    }
}

void init_gl()
{
    glClearColor(0.0, 0.0, 0.0, 0.0);

    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
/*    glFrustum(0, 2*WINDOW_WIDTH, 0, 2*WINDOW_HEIGHT, 80, -80.0);*/
    glFrustum(-WINDOW_WIDTH, 2*WINDOW_WIDTH, 
              -WINDOW_HEIGHT, 2*WINDOW_HEIGHT, 
              LATTICE_DEPTH, -LATTICE_DEPTH);

    glMatrixMode(GL_MODELVIEW);

    glShadeModel(GL_SMOOTH);
    glPointSize(POINT_SIZE);
    glLineWidth(POINT_SIZE);
}

int main(int argc, char** argv) {
    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
    glutInitWindowSize(WINDOW_WIDTH, WINDOW_HEIGHT);
    glutCreateWindow("D3Q19");
    glEnable(GL_DEPTH_TEST);
    init_vars();
    init_gl();
    glutDisplayFunc(display);
    glutIdleFunc(idle);
    glutMainLoop();
    return 0;
}