/* 
 * Galaxy Collision in OpenGL by Andrey Mirtchovski (aam396@mail.usask.ca)
 */


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <GL/glut.h>
#include <math.h>



#define SIZE    0.05
#define BOUNDS    400
#define NUMGAL 2

enum {
    SCROLL_LEFT = 1,
    SCROLL_RIGHT,
    SCROLL_UP,
    SCROLL_DOWN
} type = SCROLL_RIGHT;


typedef struct _star {
    float pos[3];
    float vel[3];
    float acc[3];
    float mass;
} star;



float ep = 0.001;

float x1, x2, x3;

double d;

int n = 2;
int gal = 1000;

star **glx;
star *s;

int h, w;
int spin_x=0, spin_y=0, old_x = 0, old_y = 0;
int move_z = 0;
int spin_z = -1500;

float boundary[8][3] = {
        {-BOUNDS, -BOUNDS, -BOUNDS},
        {-BOUNDS,  BOUNDS, -BOUNDS},
        { BOUNDS, -BOUNDS, -BOUNDS},
        { BOUNDS,  BOUNDS, -BOUNDS},
        {-BOUNDS, -BOUNDS,  BOUNDS},
        {-BOUNDS,  BOUNDS,  BOUNDS},
        { BOUNDS, -BOUNDS,  BOUNDS},
        { BOUNDS,  BOUNDS,  BOUNDS},
};

void init();

void drawboundary() {

    glColor3ub(70, 70, 70);

    glBegin(GL_LINES);
    glVertex3fv(boundary[4]);
    glVertex3fv(boundary[6]);
    glVertex3fv(boundary[5]);
    glVertex3fv(boundary[7]);
    glVertex3fv(boundary[4]);
    glVertex3fv(boundary[5]);
    glVertex3fv(boundary[6]);
    glVertex3fv(boundary[7]);
    glVertex3fv(boundary[6]);
    glVertex3fv(boundary[2]);
    glVertex3fv(boundary[4]);
    glVertex3fv(boundary[0]);
    glVertex3fv(boundary[0]);
    glVertex3fv(boundary[2]);
    glVertex3fv(boundary[5]);
    glVertex3fv(boundary[1]);
    glVertex3fv(boundary[7]);
    glVertex3fv(boundary[3]);
    glVertex3fv(boundary[1]);
    glVertex3fv(boundary[3]);
    glVertex3fv(boundary[0]);
    glVertex3fv(boundary[1]);
    glVertex3fv(boundary[2]);
    glVertex3fv(boundary[3]);
    glEnd();
}


void movegals(void) {
    int i, j;
    float ep_t_mass;

    for(i = 0; i < n; i++) {
        for(j = 0; j< n; j++) {

            if(j == i)
                continue;
    
             x1 = s[i].pos[0] - s[j].pos[0];
            x2 = s[i].pos[1] - s[j].pos[1];
            x3 = s[i].pos[2] - s[j].pos[2];

            d = sqrt (x1*x1 + x2*x2 + x3*x3);
    
            if(d > 10)
                ep_t_mass = -ep*s[j].mass/(d*d);
            else
                ep_t_mass = 0;
            
            s[i].acc[0] = ep_t_mass * x1; 
            s[i].acc[1] = ep_t_mass * x2;
            s[i].acc[2] = ep_t_mass * x3;

            s[i].vel[0] += s[i].acc[0];
            s[i].vel[1] += s[i].acc[1];
            s[i].vel[2] += s[i].acc[2];

            s[i].pos[0] += s[i].vel[0];
            s[i].pos[1] += s[i].vel[1];
            s[i].pos[2] += s[i].vel[2];
                                  
        }
    }
}

void move(void) {
    int i, j, k;
    float ep_t_mass;


    for(j = 0; j < n; j++) {
           for(i = 0; i < gal; i++) {
               for(k = 0; k < n; k++) {

                 x1 = glx[j][i].pos[0] - s[k].pos[0];
                x2 = glx[j][i].pos[1] - s[k].pos[1];
                x3 = glx[j][i].pos[2] - s[k].pos[2];
    
                d = sqrt (x1*x1 + x2*x2 + x3*x3);
        
                if(d > 10)
                    ep_t_mass = -ep*s[k].mass/(d*d);
                else 
                    ep_t_mass = 0;
                
                glx[j][i].acc[0] = ep_t_mass * x1; 
                glx[j][i].acc[1] = ep_t_mass * x2;
                glx[j][i].acc[2] = ep_t_mass * x3;

                glx[j][i].vel[0] += glx[j][i].acc[0];
                glx[j][i].vel[1] += glx[j][i].acc[1];
                glx[j][i].vel[2] += glx[j][i].acc[2];
    
                glx[j][i].pos[0] += glx[j][i].vel[0];
                glx[j][i].pos[1] += glx[j][i].vel[1];
                glx[j][i].pos[2] += glx[j][i].vel[2];
                                  
            }
        }
    }
}

void
reshape(int width, int height)
{

    h = height;
    w = width;

    glViewport(0, 0, width, height);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    //glOrtho(-width, width, -height, height, -5000, 5000);
}

void
display(void)
{
    int i, j;

    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity ();

    glPushMatrix();
        glRotatef(spin_x, 0, 1, 0);
        glRotatef(spin_y, 1, 0, 0);
    
        glMatrixMode(GL_PROJECTION);
        glLoadIdentity ();

        glPushMatrix();
            gluPerspective (50, w / h, 0.1, 10000);
            glTranslatef(0, 0, spin_z);

            drawboundary();

            /* calculate new star positions */
            move();
            /* calculate new galaxy positions */
            movegals();
            /* draw stars */
            glBegin(GL_POINTS);
            glColor3ub(0, 150, 150);
            for(i = 0; i < n; i++) 
                for(j = 0; j < gal; j++) 
                    glVertex3f(glx[i][j].pos[0], glx[i][j].pos[1], glx[i][j].pos[2]);
            glEnd();

        glPopMatrix();
    glPopMatrix();


//    glMatrixMode(GL_PROJECTION);
//    glLoadIdentity ();

    
    glutSwapBuffers();
}

void
idle(void)
{
    glutPostRedisplay();
}

void
bail(int code)
{
    free(glx);
    exit(code);
}

void
mouse(int button, int state, int x, int y)
{

    switch(button) {
        case 0:
            old_x = x - spin_x;
            old_y = y - spin_y;
            break;
        case 2:
            old_y = y - spin_z;
            move_z = (move_z ? 0 : 1);
    }
            

    glutPostRedisplay();

}

void 
motion(x, y) {

    if(!move_z) {
        spin_x = x - old_x;
        spin_y = y - old_y;
    } else {
        spin_z = y - old_y;
    }

    glutPostRedisplay();
}

void
initgal(void) {
    int i;

    free(s);
    free(glx);

    s = (star *)malloc(sizeof(star) * n);
    glx = (star **)malloc(sizeof(star *) * n);

    for(i = 0; i < n; i++)
        glx[i] = (star*)malloc(sizeof(star) * gal);
}

void
keyboard(unsigned char key, int x, int y)
{
    static int old_x = 50;
    static int old_y = 50;
    static int old_width = 512;
    static int old_height = 512;

    switch (key) {
        case 27:
                bail(0);
            break;
        case ' ':
                init();
                break;
        case 'w':
                glutPositionWindow(old_x, old_y);
                glutReshapeWindow(old_width, old_height);
            break;
        case 'f':
            if (glutGet(GLUT_WINDOW_WIDTH) < glutGet(GLUT_SCREEN_WIDTH)) {
                old_x = glutGet(GLUT_WINDOW_X);
                old_y = glutGet(GLUT_WINDOW_Y);
                old_width = glutGet(GLUT_WINDOW_WIDTH);
                old_height = glutGet(GLUT_WINDOW_HEIGHT);
                glutFullScreen();
            }
            break;
        case '1':
            n = 1;
            initgal();
            init();
            break;
        case '2':
            n = 2;
            initgal();
            init();
            break;
        case '3':
            n = 3;
            initgal();
            init();
            break;
        case '4':
            n = 4;
            initgal();
            init();
            break;
        case '5':
            n = 5;
            initgal();
            init();
            break;
        case '6':
            n = 6;
            initgal();
            init();
            break;
        case '7':
            n = 7;
            initgal();
            init();
            break;
        case '8':
            n = 8;
            initgal();
            init();
            break;
        case '9':
            n = 9;
            initgal();
            init();
            break;
    }
}


void init(void) {
    int i, j;
    float dy, dx, dist;

    for(i = 0; i < n; i++) {
        s[i].pos[0] =  (rand()%2 ? -1 : 1)*rand()%w/2;
        s[i].pos[1] =  (rand()%2 ? -1 : 1)*rand()%h/2;
        s[i].pos[2] =  (rand()%2 ? -1 : 1)*rand()%w/2;
        s[i].vel[0] = 0;
        s[i].vel[1] = 0;
        s[i].vel[2] = 0;
        s[i].acc[0] = 0;
        s[i].acc[1] = 0;
        s[i].acc[2] = 0;
        s[i].mass = 5000;
    }



    for(j = 0; j < n; j++) {
        for(i = 0; i < gal; i++) {
            
            glx[j][i].pos[0] = s[j].pos[0] + (rand()%2 ? -1 : 1) * rand()%w/4;
            glx[j][i].pos[1] = s[j].pos[1] + (rand()%2 ? -1 : 1) * rand()%h/4;
            glx[j][i].pos[2] = s[j].pos[2] + (rand()%2 ? -1 : 1) * rand()%w;

            dist = sqrt( pow(glx[j][i].pos[0] - s[j].pos[0], 2) +
                        pow(glx[j][i].pos[1] - s[j].pos[1], 2) +
                        pow(glx[j][i].pos[2] - s[j].pos[2], 2));
    
            if(dist > w/6 || dist > h/6) {
                i--;
                continue;
            }

            glx[j][i].pos[2] = s[j].pos[2] + (rand()%2 ? -1 : 1) * (rand()%w*50)/(dist*dist);
        
            dx = glx[j][i].pos[0] - s[j].pos[0];
            dy = glx[j][i].pos[1] - s[j].pos[1];

            dist = sqrt(dx*dx + dy*dy);
    
            glx[j][i].vel[0] = -(dy*1.6 + s[j].vel[0])/dist;
            glx[j][i].vel[1] =  (dx*1.6 + s[j].vel[1])/dist;
            glx[j][i].vel[2] =  0;
    
            glx[j][i].acc[0] = glx[j][i].acc[1] = glx[j][i].acc[2] = 0;
            glx[j][i].mass = 1;
        }    
    }
}

int
main(int argc, char** argv)
{
    int i;

    srand(time());

    glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
    glutInitWindowPosition(50, 50);
    glutInitWindowSize(512, 512);
    glutInit(&argc, argv);

    glutCreateWindow("Galaxies");
    glutDisplayFunc(display);
    glutReshapeFunc(reshape);

    glutKeyboardFunc(keyboard);
    glutMouseFunc(mouse);
    glutMotionFunc(motion);

    glEnable (GL_DEPTH_TEST);

    if (argc == 3) {
        sscanf(argv[1], "%d", &gal);
        sscanf(argv[2], "%d", &n);
    } else if(argc == 2) {
        if (strcmp(argv[1], "-h") == 0) {
            fprintf(stderr, "%s [stars] [galaxies]\n", argv[0]);
            exit(0);
        }
        sscanf(argv[1], "%d", &gal);
    }


    initgal();


    h = w = 512;

    init();

    glutIdleFunc(idle);
    glutMainLoop();
    return 0;
}