/* ewald_c.c
 * MPIなし・最小依存のEwald和（実空間＋逆格子＋自己項）を全サイト同時計算。
 * 返すのは：
 *   - mp[i] : サイト i のポテンシャル φ_i [1/m]（ewaldf_BM の MP と同じ次元）
 *   - F[i]  : サイト i の力ベクトル [N]
 *
 * ビルド例:
 *  Linux/macOS:
 *    g++ -O3 -fPIC -shared -o libewald_c.so ewald_c.c -lm
 *  Windows (MinGW):
 *    g++ -O3 -shared -o ewald_c.dll ewald_c.c -static-libgcc -static-libstdc++
 *
 * 使い方は同ディレクトリの ewald_fast_ctypes.py を参照。
 */

#include <stddef.h>
#include <stdint.h>
#include <string.h>

#ifdef __cplusplus
  #include <cmath>
  using std::exp;
  using std::sqrt;
  using std::erfc;
#else
  #include <math.h>
#endif

/* ---- 物理定数（ewaldf_BM と一致） ---- */
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

static const double PI   = M_PI;
static const double E_CH = 1.602176634e-19;      /* C */
static const double EPS0 = 8.854418782e-12;      /* C^2 N^-1 m^-2 */
static const double KE   = (E_CH*E_CH)/(4.0*M_PI*EPS0); /* J·m·C^-2 */
static const double ANG2M = 1e-10;

/* 行列 × ベクトル（転置） out = A^T * v, A は 3x3（行メジャ） */
static inline void matvec_T(const double *A, const double v[3], double out[3]) {
    /* A[r*3+c] */
    out[0] = A[0]*v[0] + A[3]*v[1] + A[6]*v[2];
    out[1] = A[1]*v[0] + A[4]*v[1] + A[7]*v[2];
    out[2] = A[2]*v[0] + A[5]*v[1] + A[8]*v[2];
}

static inline double dot3(const double a[3], const double b[3]) {
    return a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
}

static inline double norm3(const double v[3]) {
    return sqrt(dot3(v,v));
}

/* 2π*(h*x + k*y + l*z) */
static inline double phase_2pi(int h, int k, int l, const double frac[3]) {
    return 2.0*M_PI*( (double)h*frac[0] + (double)k*frac[1] + (double)l*frac[2] );
}

/* v^T * M * v （Mは3x3, 行メジャ） */
static inline double quad3(const double M[9], const double v[3]) {
    double t0 = M[0]*v[0] + M[1]*v[1] + M[2]*v[2];
    double t1 = M[3]*v[0] + M[4]*v[1] + M[5]*v[2];
    double t2 = M[6]*v[0] + M[7]*v[1] + M[8]*v[2];
    return v[0]*t0 + v[1]*t1 + v[2]*t2;
}

/* 共有ライブラリ公開：Cリンケージ */
#ifdef __cplusplus
extern "C" {
#endif

/* 戻り値: 0=成功, <0=エラー
 * 引数（Row-major 前提）:
 *   N                : サイト数
 *   aij[9]           : 実格子ベクトル行列（Å）
 *   Raij[9]          : 逆格子ベクトル行列（Å^-1）
 *   gij[9]           : 実空間metric（Å^2）
 *   Rgij[9]          : 逆空間metric（Å^-2）
 *   volume_A3        : 体積（Å^3）
 *   frac[3N]         : 各サイトの分率座標
 *   q[N]             : 各サイトの電荷（単位: e）
 *   alpha_invA       : Ewald α（Å^-1）
 *   nrmax[3]         : 実空間反復の最大（各軸の±範囲）
 *   hgmax[3]         : 逆格子反復の最大（各軸の±範囲）
 *   rmin_A           : 近接スキップ閾値（Å）
 *   halfspace        : 1=G空間で半空間 (l>=0、l!=0 は2倍)、0=全空間
 *   out_mp[N]        : 出力 φ_i（1/m）
 *   out_F_N[3N]      : 出力 力（N）
 */
int ewald_all_sites(
    int N,
    const double *aij,
    const double *Raij,
    const double *gij,
    const double *Rgij,
    double volume_A3,
    const double *frac,
    const double *q,
    double alpha_invA,
    const int *nrmax,
    const int *hgmax,
    double rmin_A,
    int halfspace,
    double *out_mp,
    double *out_F_N
) {
    if (N<=0 || !aij || !Raij || !gij || !Rgij || !frac || !q || !nrmax || !hgmax || !out_mp || !out_F_N) {
        return -1;
    }

    /* 初期化 */
    for (int i=0;i<N;i++) {
        out_mp[i] = 0.0;
        out_F_N[3*i+0] = 0.0;
        out_F_N[3*i+1] = 0.0;
        out_F_N[3*i+2] = 0.0;
    }

    const int nrx = nrmax[0], nry = nrmax[1], nrz = nrmax[2];
    const int hgx = hgmax[0], hgy = hgmax[1], hgz = hgmax[2];

    const double vol_m3 = volume_A3 * (ANG2M*ANG2M*ANG2M);
    const double Krec   = 1.0/(M_PI*vol_m3);          /* 1/(πV) */
    const double Kexp   = (M_PI*M_PI)/(alpha_invA*alpha_invA);
    const double alpha_m = alpha_invA * 1e10;         /* [1/m] */

    /* ---------- 実空間 ---------- */
    for (int ix=-nrx; ix<=nrx; ++ix) {
        for (int iy=-nry; iy<=nry; ++iy) {
            for (int iz=-nrz; iz<=nrz; ++iz) {
                for (int i=0; i<N; ++i) {
                    const double *ri_frac = &frac[3*i];
                    const double qi = q[i];
                    for (int j=0; j<N; ++j) {
                        /* 同一点（i==j && T==0）はスキップ */
                        if (ix==0 && iy==0 && iz==0 && i==j) continue;

                        double rj_frac[3] = { frac[3*j+0] + ix,
                                              frac[3*j+1] + iy,
                                              frac[3*j+2] + iz };
                        double d_frac[3] = { ri_frac[0] - rj_frac[0],
                                             ri_frac[1] - rj_frac[1],
                                             ri_frac[2] - rj_frac[2] };
                        /* r_vec[Å] = aij^T * d_frac */
                        double rA[3]; matvec_T(aij, d_frac, rA);
                        double r_ang = norm3(rA);
                        if (r_ang < 1e-12 || r_ang < rmin_A) continue;

                        /* φ_i += q_j erfc(αr)/r  [1/m] */
                        out_mp[i] += q[j] * erfc(alpha_invA * r_ang) / (r_ang * ANG2M);

                        /* 力: F_i += qi*qj*Ke*[erfc(αr)/r^3 + (2α/√π)e^{-α^2r^2}/r^2] * r_vec[m] */
                        double rm[3] = { rA[0]*ANG2M, rA[1]*ANG2M, rA[2]*ANG2M };
                        double r_m = norm3(rm);
                        if (r_m > 0.0) {
                            double expterm = exp(-(alpha_m*r_m)*(alpha_m*r_m));
                            double G = erfc(alpha_m*r_m)/(r_m*r_m*r_m)
                                     + (2.0*alpha_m/sqrt(M_PI))*expterm/(r_m*r_m);
                            double coeff = qi * KE * q[j] * G;
                            out_F_N[3*i+0] += coeff * rm[0];
                            out_F_N[3*i+1] += coeff * rm[1];
                            out_F_N[3*i+2] += coeff * rm[2];
                        }
                    }
                }
            }
        }
    }

    /* ---------- 逆格子 ---------- */
    /* l 半空間（l>=0）: l!=0 のとき寄与を 2 倍。全空間にしたい場合は halfspace=0 */
    const int l_start = halfspace ? 0 : -hgz;
    for (int l = l_start; l<=hgz; ++l) {
        for (int k=-hgy; k<=hgy; ++k) {
            for (int h=-hgx; h<=hgx; ++h) {
                if (h==0 && k==0 && l==0) continue;
                /* |G|^2 = [h k l]^T Rgij [h k l] */
                double hkltmp[3] = { (double)h, (double)k, (double)l };
                double G2_Ainv2 = quad3(Rgij, hkltmp);
                if (G2_Ainv2 <= 0.0) continue;

                /* amp = exp(-π^2 G2/α^2) / (G2 * 1e20)  → [m^2] */
                double amp = exp(-Kexp * G2_Ainv2) / (G2_Ainv2 * 1e20);

                /* 構造因子（電荷重み） */
                double Csum = 0.0, Ssum = 0.0;
                for (int j=0; j<N; ++j) {
                    double th = phase_2pi(h,k,l, &frac[3*j]);
                    Csum += q[j]*cos(th);
                    Ssum += q[j]*sin(th);
                }

                /* φ_i, ∇φ_i を評価 */
                for (int i=0; i<N; ++i) {
                    double thi = phase_2pi(h,k,l, &frac[3*i]);
                    double cti = cos(thi), sti = sin(thi);
                    double fcal = cti*Csum + sti*Ssum;       /* 実部: Re(S)*cos + Im(S)*sin */
                    double gph  = Ssum*cti - Csum*sti;       /* 位相微分に出る組み合わせ */

                    double mult = (halfspace && l!=0) ? 2.0 : 1.0;
                    double coef = Krec * amp;

                    /* φ_i [1/m] */
                    out_mp[i] += mult * coef * fcal;

                    /* 勾配：coef*(2π)*gph * G_vec[m^-1] */
                    double GvA[3] = {
                        h*Raij[0] + k*Raij[3] + l*Raij[6],
                        h*Raij[1] + k*Raij[4] + l*Raij[7],
                        h*Raij[2] + k*Raij[5] + l*Raij[8]
                    }; /* Å^-1 */
                    double Gvm[3] = { GvA[0]*1e10, GvA[1]*1e10, GvA[2]*1e10 };
                    double gcoef = mult * coef * (2.0*M_PI) * gph;

                    /* F = - q_i * Ke * grad_phi */
                    double scale = - q[i] * KE * gcoef;
                    out_F_N[3*i+0] += scale * Gvm[0];
                    out_F_N[3*i+1] += scale * Gvm[1];
                    out_F_N[3*i+2] += scale * Gvm[2];
                }
            }
        }
    }

    /* ---------- 自己項（サイト独立） ----------
     * φ_self = 2 α / √π  [1/Å → 1/m へ変換済み上で適用]
     * ここでは φ_i から self_term を減算。
     */
    const double self_term = 2.0 * (alpha_invA*1e10) / sqrt(M_PI);
    for (int i=0;i<N;++i) out_mp[i] -= self_term;

    return 0;
}

#ifdef __cplusplus
} /* extern "C" */
#endif