#include <fcntl.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <unistd.h>

#define eps 1e-3

int readfile(char *filename, void *buf, int count, int size) {
  int inputfd = open(filename, O_RDONLY);
  if (inputfd == -1) {
    printf("File open error!\n");
    return 1;
  }

  int expected = count * size;
  int result = read(inputfd, buf, expected);
  close(inputfd);

  if (result != expected) {
    printf("File read error: expected %d bytes but read %d bytes!\n", expected,
           result);
    return 1;
  }
  return 0;
}

int writefile(char *filename, void *buf, int count, int size) {
  int outputfd = open(filename, O_RDWR | O_CREAT, 0666);
  if (outputfd == -1) {
    printf("File create error!\n");
    return 1;
  }

  int expected = count * size;
  int result = write(outputfd, buf, expected);
  close(outputfd);

  if (result != expected) {
    printf("File write error: expected %d bytes but wrote %d bytes!\n",
           expected, result);
    return 1;
  }
  return 0;
}

int main(int ac, char *av[]) {
  int N;
  char *inputfile;
  int timesteps;
  double deltaT;

  if (ac != 6) {
    printf("Using default arguments '3000 ../../nbody_input.gal 100 1e-5 0'\n");
    // printf("Correct Format: %s N inputfile timesteps deltaT graphics\n",
    // av[0]);
    N = 3000;
    inputfile = "../../nbody_input.gal";
    timesteps = 100;
    deltaT = 1e-5;
  } else {
    N = atoi(av[1]);
    inputfile = av[2];
    timesteps = atoi(av[3]);
    deltaT = atof(av[4]);
  }

  const double G = 100.0 / N;

  clock_t start;
  start = clock();

  // Phase 1: File to Datastucture

  typedef struct {
    double x;
    double y;
    double m;
    double vx;
    double vy;
    double pad[1];
  } part;

  typedef struct {
    double ax;
    double ay;
  } accel;

  part *parts = (part *)malloc(N * sizeof(part));
  accel *accels = (accel *)malloc(N * sizeof(accel));

  printf("Reading input...\n");
  int err = readfile(inputfile, parts, N, sizeof(part));
  if (err != 0) {
    return err;
  }

  printf("Time for File to DS(s): %lf\n",
         (double)(clock() - start) / CLOCKS_PER_SEC);
  start = clock();

  // Phase 2: Processing
  int instant = 0;
  while (instant < timesteps) {
    memset(accels, 0, N * sizeof(accel));
    for (int row = 0; row < N; row++) {
      double mi = parts[row].m;
      double x = parts[row].x;
      double y = parts[row].y;
      double ax = 0.0;
      double ay = 0.0;

      for (int col = row + 1; col < N; col++) {
        double dx = x - parts[col].x;
        double dy = y - parts[col].y;
        double mj = parts[col].m;
        double rij = sqrt(dx * dx + dy * dy);
        double rije = rij + eps;
        double radius = rije * rije * rije;
        double force = G / radius;

        // Optimization: Newton's Third Law
        double fx = force * dx;
        double fy = force * dy;
        ax -= fx * mj;
        ay -= fy * mj;
        accels[col].ax += fx * mi;
        accels[col].ay += fy * mi;
      }

      accels[row].ax += ax;
      accels[row].ay += ay;

      parts[row].vx += deltaT * accels[row].ax;
      parts[row].vy += deltaT * accels[row].ay;
      parts[row].x += deltaT * parts[row].vx;
      parts[row].y += deltaT * parts[row].vy;
    }

    instant++;
  }

  printf("Time for Processing(s): %lf\n",
         (double)(clock() - start) / CLOCKS_PER_SEC);
  start = clock();

  // Phase 3: Datastucture to File
  err = writefile("results.gal", parts, N, sizeof(part));

  free(parts);
  free(accels);

  printf("Time for DS to File(s): %lf\n",
         (double)(clock() - start) / CLOCKS_PER_SEC);

  return 0;
}