#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "utilities.h"
#include "particle.h"

double axilrod_teller(double i, double j, double k) {

    double r_ij, r_ik, r_kj;
    double rsq_ij, rsq_ik, rsq_kj;
    double mul, num;

    r_ij = fmax(fabs(i - j), 1e-10);
    r_ik = fmax(fabs(i - k), 1e-10);
    r_kj = fmax(fabs(k - j), 1e-10);

    rsq_ij = r_ij * r_ij;
    rsq_ik = r_ik * r_ik;
    rsq_kj = r_kj * r_kj;

    mul = r_ij * r_ik * r_kj;

    num = 3 * (-rsq_ij + rsq_ik + rsq_kj)
            * (rsq_ij - rsq_ik + rsq_kj)
            * (rsq_ij + rsq_ik - rsq_kj);

    return 1 / pow(mul, 3.0) + num / pow(mul, 5.0);
}

void compute_forces_dim(double i, double j, double k, double *i_f, double *j_f, double *k_f) {

    double h = 4.69041575982343e-08;
    double h_i, h_j, h_k;
    double d_i, d_j, d_k;
    double i_f1, i_f2, j_f1, j_f2, k_f1, k_f2;

    h_i = h * i;
    h_j = h * j;
    h_k = h * k;

    d_i = (i + h_i) - (i - h_i);
    d_j = (j + h_j) - (j - h_j);
    d_k = (k + h_k) - (j - h_k);

    i_f1 = axilrod_teller(i + h_i, j, k);
    i_f2 = axilrod_teller(i - h_i, j, k);
    j_f1 = axilrod_teller(i, j + h_j, k);
    j_f2 = axilrod_teller(i, j - h_j, k);
    k_f1 = axilrod_teller(i, j, k + h_k);
    k_f2 = axilrod_teller(i, j, k - h_k);

    *i_f = - (i_f1 - i_f2) / d_i;
    *j_f = - (j_f1 - j_f2) / d_j;
    *k_f = - (k_f1 - k_f2) / d_k;
}


void  compute_forces(struct particle *i, struct particle *j, struct particle *k) {

    compute_forces_dim(i->x, j->x, k->x, &i->fx, &j->fx, &k->fx);
    compute_forces_dim(i->y, j->y, k->y, &i->fy, &j->fy, &k->fy);
    compute_forces_dim(i->z, j->z, k->z, &i->fz, &j->fz, &k->fz);
}

void update(struct particle *p, double deltatime) {

    p->x += p->vx * deltatime + p->ax * deltatime * deltatime / 2;
    p->y += p->vy * deltatime + p->ay * deltatime * deltatime / 2;
    p->z += p->vz * deltatime + p->az * deltatime * deltatime / 2;

    p->ax += p->fx;
    p->ay += p->fy;
    p->az += p->fz;

    p->vx += (2 * p->ax - p->fx) * deltatime / 2;
    p->vy += (2 * p->ay - p->fy) * deltatime / 2;
    p->vz += (2 * p->az - p->fz) * deltatime / 2;

    p->fx = 0;
    p->fy = 0;
    p->fz = 0;
}

int main() {
    int allParticlesCount;
    struct particle *allParticles;

    allParticlesCount = parse_input("./input.txt", &allParticles);

    for (int i = 0; i < allParticlesCount; ++i) {
        for (int j = i + 1;  j < allParticlesCount; ++j) {
            for (int k = j + 1; k < allParticlesCount; ++k) {
                compute_forces(&allParticles[i], &allParticles[j], &allParticles[k]);
            }
        }
    }

    for (int i = 0; i < allParticlesCount; ++i) {
        update(&allParticles[i], 0.5);
    }

    write_output("./output.txt", allParticlesCount, allParticles);

    free(allParticles);
}