/* nbody.c by Jack Carrozzo */

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

#define NUMPARTS   1000 //how many particles
#define GRID_X     640
#define GRID_Y     480

#define MASS       1 //plus 10 zeros
#define CENTERMASS 100000 //plus 10 zeros
#define TIMESTEP   0.05 //seconds per step //0.05,1 is good
#define IMAGESTEP  1 //make bmp every n steps

#define rad(d)     (3.1415926 * d) / 180 //radians from degress
#define DEBUG      1 //0 off, 1 on

struct PartType
{
    int               x; //x coord for body
    int               y; //y coord
    unsigned long int m;
    double            vi;//i-hat of velocity vector
    double            vj;//j-hat "                "
};

struct ForceType
{
    double i;
    double j;
};

struct PartType Particles[NUMPARTS];
struct ForceType ForceSum;
struct ForceType MoveForce[NUMPARTS]; //where forces are stored until move()
int ImageGrid[GRID_X][GRID_Y];
int Steps,ERR;
//double G=0.667; //has 10 zeros taken off
double G=0.667;
int ImgSteps=0;

void InitParticles( void )
{
    int i,j,r=100;

    srand(time()*getpid());
    Particles[0].x = (int)(GRID_X / 2);
    Particles[0].y = (int)(GRID_Y / 2);
    Particles[0].m = CENTERMASS;
    Particles[0].vi = 0;
    Particles[0].vj = 0;

    /*Particles[1].x = (Particles[0].x);
      Particles[1].y = (Particles[0].y + r);
      Particles[1].m = MASS;
    //Particles[1].vi = (sqrt((G * CENTERMASS) / (double)(r)) / 10);
    //Particles[1].vi = (sqrt((G * CENTERMASS) / (double)(r)));
    Particles[1].vi = 5; //seems to work
    printf("STARTING Vi is %f.\n",Particles[1].vi);
    Particles[1].vj = 0; */

    for (i=1;i<NUMPARTS;i++) //starts at 1, not 0
    {
        Particles[i].x = ((double)rand()/((double)pow(2,31))) * GRID_X;
        Particles[i].y = ((double)rand()/((double)pow(2,31))) * GRID_Y;
        Particles[i].m = MASS;
        Particles[i].vi = 0;
        Particles[i].vj = 0;
    }

    for (i=0;i<GRID_X;i++)
    {
        for (j=0;j<GRID_Y;j++)
        {
            ImageGrid[i][j]=0;
        }
    }
}

void ComputeForces( void )
{
    int i,j,fxsign,fysign; //signs: 0 +, 1 -
    double r,F,theta,dx,dy,fx,fy;

    for (i=0;i<NUMPARTS;i++)
    {
        ForceSum.i=0;
        ForceSum.j=0;
        for (j=0;j<NUMPARTS;j++)
        {
            if (j!=i)
            {
                fxsign=0;
                fysign=0;
                if ((Particles[j].x - Particles[i].x) < 0) fxsign++;
                if ((Particles[j].y - Particles[i].y) < 0) fysign++;
                dx = abs(Particles[i].x - Particles[j].x);
                dy = abs(Particles[i].y - Particles[j].y);
                theta = atan(dy/dx); //theta in radians
                r = sqrt((dx * dx) + (dy * dy)); //distance
                F = ((G * (double)(Particles[i].m) * (double)(Particles[j].m)) / (r * r));
                if (DEBUG) printf("F= %e, G= %e m1= %d, m2= %d, r= %f\n",F,G,Particles[i].m,Particles[j].m,r);
                fx = cos(theta) * F;
                fy = sin(theta) * F;
                if (DEBUG) if (i==2) printf("compareing to %d, r - %f, force x %f, force y %f, theta %f\n",j,r,fx,fy,theta);
                if (fxsign==0)
                    ForceSum.i += fx;
                else
                    ForceSum.i -= fx;

                if (fysign==0)
                    ForceSum.j += fy;
                else
                    ForceSum.j -= fy;
            }
        }
        MoveForce[i].i = ForceSum.i;
        MoveForce[i].j = ForceSum.j;
    }
}

void MoveParticles( void )
{
    int i;
    double A; //acceleration

    for (i=0;i<NUMPARTS;i++)
    {
        //F=MA, A=F/M
        A = (MoveForce[i].i / (Particles[i].m));
        Particles[i].vi += (A * TIMESTEP);
        A = (MoveForce[i].j / (Particles[i].m));
        //A *= 25; //don't know why, but need this
        Particles[i].vj += (A * TIMESTEP);

        Particles[i].x += (int)(Particles[i].vi);
        Particles[i].y += (int)(Particles[i].vj);

        if (DEBUG) if (i==2) printf("part 2, y %f, f %f, %f.\n",Particles[i].vj,MoveForce[i].i,MoveForce[i].j);
    }
    Steps++;
}

void DataOut( void )
{
    int i,j;
    char filename[100];

    if (DEBUG)
    {
        printf("\n   [+] Time: %d [+]\n",Steps);
        for (i=0;i<NUMPARTS;i++)
        {
            printf("----------------------\n");
            printf("Particle Num     : %d\n",i);
            printf("Particle Position: (%d,%d)\n",Particles[i].x,Particles[i].y);
            printf("Particle Velocity: (%f,%f)\n",Particles[i].vi,Particles[i].vj);
            printf("----------------------\n");
        }
    }

    if ((Steps % IMAGESTEP)==0)
    {
        CreateBitmap(GRID_X,GRID_Y);

        for (i=0;i<NUMPARTS;i++)
        {
            for (j=0;j<NUMPARTS;j++)
            {
                SetPixel(i, j, 0, 0, 0); //blackout image
            }
        }

        for (i=0;i<NUMPARTS;i++)
        {
            if (Particles[i].x >= 0 && Particles[i].y >= 0 && Particles[i].x < GRID_X && Particles[i].y < GRID_Y)
            {
                if (i==0)
                {
                    SetPixel((Particles[i].x), (Particles[i].y), 255, 255, 255); //big dot for big piece
                    SetPixel((Particles[i].x + 1), (Particles[i].y), 255, 255, 255);
                    SetPixel((Particles[i].x + 1), (Particles[i].y - 1),  255, 255, 255);
                    SetPixel((Particles[i].x), (Particles[i].y -1), 255, 255, 255);

                    SetPixel((Particles[i].x - 1), (Particles[i].y + 1), 255, 255, 255);
                    SetPixel((Particles[i].x - 1), (Particles[i].y), 255, 255, 255);
                    SetPixel((Particles[i].x - 1), (Particles[i].y - 1), 255, 255, 255);
                    SetPixel((Particles[i].x - 1), (Particles[i].y - 2), 255, 255, 255);

                    SetPixel((Particles[i].x), (Particles[i].y + 1), 255, 255, 255, 255);
                    SetPixel((Particles[i].x + 1), (Particles[i].y + 1), 255, 255, 255);
                    SetPixel((Particles[i].x + 2), (Particles[i].y + 1), 255, 255, 255);

                    SetPixel((Particles[i].x + 2), (Particles[i].y), 255, 255, 255);
                    SetPixel((Particles[i].x + 2), (Particles[i].y - 1), 255, 255, 255);
                    SetPixel((Particles[i].x + 2), (Particles[i].y - 2), 255, 255, 255);

                    SetPixel((Particles[i].x), (Particles[i].y - 2), 255, 255, 255);
                    SetPixel((Particles[i].x + 1), (Particles[i].y - 2), 255, 255, 255);
                } else {
                    SetPixel((Particles[i].x), (Particles[i].y), 255, 255, 255);
                    SetPixel((Particles[i].x + 1), (Particles[i].y), 255, 255, 255);
                    SetPixel((Particles[i].x + 1), (Particles[i].y - 1), 255, 255, 255);
                    SetPixel((Particles[i].x), (Particles[i].y -1), 255, 255, 255);
                }
            }
        }

        sprintf(filename,"img/step%d.bmp",ImgSteps);
        ImgSteps++;

        WriteBitmapToFile(&filename);

        DeleteBitmap();
    }
}

//go to bottom to see main

// bitmap.c
// --------
// author: Levi Lansing
// email: levi_is@hotmail.com
// last updated: 10/15/03
//
// DESCRIPTION:
// - a simple set of bitmap opperations
// - outputs a standard windows bitmap file

//#include <stdio.h>

typedef int BOOL;
#define TRUE	1
#define FALSE	0

typedef unsigned short WORD;
typedef unsigned long DWORD;
typedef unsigned char BYTE;


//bitmap headers
typedef struct tag_BITMAPFILEHEADER {
    WORD     bfType;
    DWORD      bfSize;
    WORD     bfReserved1;
    WORD     bfReserved2;
    DWORD      bfOffBits;
} BITMAPFILEHEADER;

typedef struct tag_BITMAPINFOHEADER{
    DWORD      biSize;
    DWORD      biWidth;
    DWORD      biHeight;
    WORD     biPlanes;
    WORD     biBitCount;
    DWORD      biCompression;
    DWORD      biSizeImage;
    DWORD      biXPelsPerMeter;
    DWORD      biYPelsPerMeter;
    DWORD      biClrUsed;
    DWORD      biClrImportant;
    DWORD      bmiColors;
} BITMAPINFOHEADER;



//some global variables to keep track of the bitmap info
//don't mess with this stuff
DWORD b_height;
DWORD b_width;
DWORD b_widthbytes;
BYTE* b_pBits = NULL;

// CALL THIS FIRST!
// creates a bitmap and allocates the memory for it
void CreateBitmap(int width,int height)
{
    if(b_pBits)
        free(b_pBits);

    b_width = width;
    b_height = height;
    b_widthbytes = b_width*3;
    if(b_width%4)
        b_widthbytes += 4-(b_width%4);
    b_pBits = (BYTE*)malloc(b_widthbytes*height);
}

//frees the allocated memory for the bitmap
void DeleteBitmap()
{
    if(b_pBits)
        free(b_pBits);
    b_pBits = NULL;
}

//sets the color of a single pixel in the bitmap
// colors red,green,and blue accept values between 0 (dark) and 255 (light) only!
void SetPixel(int x, int y, BYTE red, BYTE green, BYTE blue)
{
    if(!b_pBits)
    {
        printf("Error: you forgot to create the bitmap first!\n");
        return;
    }
    if(x>=b_width || y>=b_height || x<0 || y<0)
    {
        printf("pixel (%d,%d) is out of bounds- bitmap is only %dx%d!\n",x,y,b_width,b_height);
        return;
    }

    BYTE* pBits = b_pBits + y*b_widthbytes + x*3;
    pBits[0] = (BYTE) blue;
    pBits[1] = (BYTE) green;
    pBits[2] = (BYTE) red;
}

BOOL WriteBitmapToFile(char* filename)
{
    if(!b_pBits)
    {
        printf("Error: you forgot to create the bitmap first!\n");
        return;
    }

    BITMAPFILEHEADER bmfh;
    BITMAPINFOHEADER bmi;
    int i;

    FILE* fp = fopen(filename,"wb");
    if(!fp)
        return FALSE;

    bmfh.bfType = 0x4D42;
    bmfh.bfSize = 40+14+b_widthbytes*b_height;
    bmfh.bfReserved1 = 0;
    bmfh.bfReserved2 = 0;
    bmfh.bfOffBits = 14+40;

    //write the file header
    fwrite(&bmfh.bfType,2,1,fp);
    fwrite(&bmfh.bfSize,4,1,fp);
    fwrite(&bmfh.bfReserved1,2,1,fp);
    fwrite(&bmfh.bfReserved2,2,1,fp);
    fwrite(&bmfh.bfOffBits,4,1,fp);

    bmi.biSize = 40;
    bmi.biBitCount = 24;
    bmi.biClrImportant = 0;
    bmi.biClrUsed = 0;
    bmi.biCompression = 0;
    bmi.biHeight = b_height;
    bmi.biPlanes = 1;
    bmi.biSizeImage = b_width*b_height*3;
    bmi.biWidth = b_width;
    bmi.biXPelsPerMeter = 3780;
    bmi.biYPelsPerMeter = 3780;
    bmi.bmiColors = 0;

    //write the bitmap header
    fwrite(&bmi.biSize,40,1,fp);

    //bitmap files are upside down! flip it over!
    BYTE* pBits = b_pBits + (b_height-1)*b_widthbytes;
    for(i=0; i<b_height; i++)
    {
        fwrite((void*)pBits,b_widthbytes,1,fp);
        pBits = pBits - b_widthbytes;
    }
    fclose(fp);

    return TRUE;
}

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

    printf("Simple NBODY solver by Jack Carrozzo.\n");

    if (argc!=2)
    {
        printf("Usage: ./nbody <frames>\n");
        exit(1);
    }

    frames = atoi(argv[1]);

    printf("Will go for %d iterations, writing %d frames.\n",(int)(frames*IMAGESTEP),frames);

    InitParticles();

    if (DEBUG)
    {
        printf("\n   [+] Time: %d [+]\n",Steps);
        for (i=0;i<NUMPARTS;i++)
        {
            printf("----------------------\n");
            printf("Particle Num     : %d\n",i);
            printf("Particle Position: (%d,%d)\n",Particles[i].x,Particles[i].y);
            printf("Particle Velocity: (%f,%f)\n",Particles[i].vi,Particles[i].vj);
            printf("----------------------\n");
        }
    }
    DataOut();

    for (i=0;i<(frames*IMAGESTEP);i++)
    {
        ComputeForces();
        MoveParticles();
        DataOut();
    }

    if (DEBUG)
    {
        printf("\n   [+] Time: %d [+]\n",Steps);
        for (i=0;i<NUMPARTS;i++)
        {
            printf("----------------------\n");
            printf("Particle Num     : %d\n",i);
            printf("Particle Position: (%d,%d)\n",Particles[i].x,Particles[i].y);
            printf("Particle Velocity: (%f,%f)\n",Particles[i].vi,Particles[i].vj);
            printf("----------------------\n");
        }
    }
    printf("Done.\n");

    return 0;
}