// nbody.c
// gcc -O3 -shared -o nbody.dll nbody.c
// gcc -O3 -fPIC -shared nbody.c -o libnbody.so
// gcc -O3 -fPIC -shared nbody.c -o libnbody.so -fopenm

// nbody.c
#include <math.h>
#include <stddef.h>

#ifdef _WIN32
#define API __declspec(dllexport)
#else
#define API
#endif

#ifndef SOFTEN2
#define SOFTEN2 1e-18   // r^2=0 の事故を避ける微小ソフトニング
#endif

// a[3*i+k] に加速度を書き込む (in-place)
API void accel_inplace(const double *m, const double *r, double *a, int N, double G1) {
    for (int i = 0; i < N; ++i) {
        double xi = r[3*i+0], yi = r[3*i+1], zi = r[3*i+2];
        double ax = 0.0, ay = 0.0, az = 0.0;
        for (int j = 0; j < N; ++j) {
            if (i == j) continue;
            double dx = r[3*j+0] - xi;
            double dy = r[3*j+1] - yi;
            double dz = r[3*j+2] - zi;
            double r2 = dx*dx + dy*dy + dz*dz + SOFTEN2;
            double invr  = 1.0 / sqrt(r2);
            double invr3 = invr * invr * invr;
            double s = G1 * m[j] * invr3;
            ax += s * dx; ay += s * dy; az += s * dz;
        }
        a[3*i+0] = ax; a[3*i+1] = ay; a[3*i+2] = az;
    }
}

// Euler を k ステップ前進
API void step_euler_k(const double *m, double *r, double *v, double *a,
                      double dt, int k, int N, double G1) {
    for (int step = 0; step < k; ++step) {
        // v(t+dt) = v + dt a
        for (int i=0; i<N; ++i) {
            v[3*i+0] += dt * a[3*i+0];
            v[3*i+1] += dt * a[3*i+1];
            v[3*i+2] += dt * a[3*i+2];
        }
        // r(t+dt) = r + dt v(t+dt)
        for (int i=0; i<N; ++i) {
            r[3*i+0] += dt * v[3*i+0];
            r[3*i+1] += dt * v[3*i+1];
            r[3*i+2] += dt * v[3*i+2];
        }
        // a(t+dt)
        accel_inplace(m, r, a, N, G1);
    }
}

// 位置Verlet（r_prev 付き）を k ステップ前進
API void step_verlet_k(const double *m, double *r_prev, double *r, double *v, double *a,
                       double dt, int k, int N, double G1) {
    for (int step = 0; step < k; ++step) {
        // r_new = 2 r - r_prev + dt^2 a
        for (int i=0; i<N; ++i) {
            double r_new_x = 2.0*r[3*i+0] - r_prev[3*i+0] + (dt*dt)*a[3*i+0];
            double r_new_y = 2.0*r[3*i+1] - r_prev[3*i+1] + (dt*dt)*a[3*i+1];
            double r_new_z = 2.0*r[3*i+2] - r_prev[3*i+2] + (dt*dt)*a[3*i+2];
            // v ≈ (r_new - r_prev)/(2 dt)
            v[3*i+0] = (r_new_x - r_prev[3*i+0]) / (2.0*dt);
            v[3*i+1] = (r_new_y - r_prev[3*i+1]) / (2.0*dt);
            v[3*i+2] = (r_new_z - r_prev[3*i+2]) / (2.0*dt);
            // rotate buffers: r_prev <- r, r <- r_new
            r_prev[3*i+0] = r[3*i+0]; r[3*i+0] = r_new_x;
            r_prev[3*i+1] = r[3*i+1]; r[3*i+1] = r_new_y;
            r_prev[3*i+2] = r[3*i+2]; r[3*i+2] = r_new_z;
        }
        accel_inplace(m, r, a, N, G1);
    }
}

// 速度Verletを k ステップ前進
API void step_vverlet_k(const double *m, double *r, double *v, double *a,
                        double dt, int k, int N, double G1) {
    for (int step = 0; step < k; ++step) {
        // v(t+dt/2)
        for (int i=0; i<N; ++i) {
            v[3*i+0] += 0.5 * dt * a[3*i+0];
            v[3*i+1] += 0.5 * dt * a[3*i+1];
            v[3*i+2] += 0.5 * dt * a[3*i+2];
        }
        // r(t+dt)
        for (int i=0; i<N; ++i) {
            r[3*i+0] += dt * v[3*i+0];
            r[3*i+1] += dt * v[3*i+1];
            r[3*i+2] += dt * v[3*i+2];
        }
        // a(t+dt)
        accel_inplace(m, r, a, N, G1);
        // v(t+dt)
        for (int i=0; i<N; ++i) {
            v[3*i+0] += 0.5 * dt * a[3*i+0];
            v[3*i+1] += 0.5 * dt * a[3*i+1];
            v[3*i+2] += 0.5 * dt * a[3*i+2];
        }
    }
}

API double kinetic_energy(const double *m, const double *v, int N) {
    double K = 0.0;
    for (int i=0; i<N; ++i) {
        double vx = v[3*i+0], vy = v[3*i+1], vz = v[3*i+2];
        K += 0.5 * m[i] * (vx*vx + vy*vy + vz*vz);
    }
    return K;
}

API double potential_energy(const double *m, const double *r, int N, double G1) {
    double U = 0.0;
    for (int i=0; i<N; ++i) {
        double xi = r[3*i+0], yi = r[3*i+1], zi = r[3*i+2];
        for (int j=i+1; j<N; ++j) {
            double dx = r[3*j+0]-xi, dy = r[3*j+1]-yi, dz = r[3*j+2]-zi;
            double dist = sqrt(dx*dx + dy*dy + dz*dz + SOFTEN2);
            U += G1 * m[i] * m[j] / dist;
        }
    }
    return U;
}

API void momentum(const double *m, const double *v, int N,
                  double *Px, double *Py, double *Pz, double *Pmsm) {
    double px = 0.0, py = 0.0, pz = 0.0;
    for (int i=0; i<N; ++i) {
        px += m[i] * v[3*i+0];
        py += m[i] * v[3*i+1];
        pz += m[i] * v[3*i+2];
    }
    *Px = px; *Py = py; *Pz = pz;
    double p2 = px*px + py*py + pz*pz;
    *Pmsm = sqrt(p2 / (3.0 * (N > 0 ? N : 1)));
}
