#include "macrodef.h"
#include "bssn_rhs.h"
#include "fmisc.h"
#include "ricci_gamma.h"
#include "share_func.h"
#include "tool.h"
#include <vector>

#ifdef fortran1
#define f_constraint_bssn constraint_bssn
#define f_z4c_rhs_point z4c_rhs_point
#endif
#ifdef fortran2
#define f_constraint_bssn CONSTRAINT_BSSN
#define f_z4c_rhs_point Z4C_RHS_POINT
#endif
#ifdef fortran3
#define f_constraint_bssn constraint_bssn_
#define f_z4c_rhs_point z4c_rhs_point_
#endif

extern "C" void f_constraint_bssn(int *, double *, double *, double *,
                                  double *, double *,
                                  double *, double *, double *, double *, double *, double *,
                                  double *, double *, double *, double *, double *, double *,
                                  double *, double *, double *,
                                  double *, double *, double *, double *, double *, double *, double *, double *, double *, double *,
                                  double *, double *, double *, double *, double *, double *,
                                  double *, double *, double *, double *, double *, double *,
                                  double *, double *, double *, double *, double *, double *,
                                  double *, double *, double *, double *, double *, double *, double *, double *,
                                  double *, double *, double *,
                                  int &);

extern "C" void f_z4c_rhs_point(
    double &A11,
    double &A12,
    double &A13,
    double &A22,
    double &A23,
    double &A33,
    double &alpha,
    double &B1,
    double &B2,
    double &B3,
    double &beta1,
    double &beta2,
    double &beta3,
    double &chi,
    double &chiDivFloor,
    double &da1,
    double &dA111,
    double &dA112,
    double &dA113,
    double &dA122,
    double &dA123,
    double &dA133,
    double &da2,
    double &dA211,
    double &dA212,
    double &dA213,
    double &dA222,
    double &dA223,
    double &dA233,
    double &da3,
    double &dA311,
    double &dA312,
    double &dA313,
    double &dA322,
    double &dA323,
    double &dA333,
    double &db11,
    double &dB11,
    double &db12,
    double &dB12,
    double &db13,
    double &dB13,
    double &db21,
    double &dB21,
    double &db22,
    double &dB22,
    double &db23,
    double &dB23,
    double &db31,
    double &dB31,
    double &db32,
    double &dB32,
    double &db33,
    double &dB33,
    double &dchi1,
    double &dchi2,
    double &dchi3,
    double &dda11,
    double &dda12,
    double &dda13,
    double &dda22,
    double &dda23,
    double &dda33,
    double &ddb111,
    double &ddb112,
    double &ddb113,
    double &ddb121,
    double &ddb122,
    double &ddb123,
    double &ddb131,
    double &ddb132,
    double &ddb133,
    double &ddb221,
    double &ddb222,
    double &ddb223,
    double &ddb231,
    double &ddb232,
    double &ddb233,
    double &ddb331,
    double &ddb332,
    double &ddb333,
    double &ddchi11,
    double &ddchi12,
    double &ddchi13,
    double &ddchi22,
    double &ddchi23,
    double &ddchi33,
    double &deldelg1111,
    double &deldelg1112,
    double &deldelg1113,
    double &deldelg1122,
    double &deldelg1123,
    double &deldelg1133,
    double &deldelg1211,
    double &deldelg1212,
    double &deldelg1213,
    double &deldelg1222,
    double &deldelg1223,
    double &deldelg1233,
    double &deldelg1311,
    double &deldelg1312,
    double &deldelg1313,
    double &deldelg1322,
    double &deldelg1323,
    double &deldelg1333,
    double &deldelg2211,
    double &deldelg2212,
    double &deldelg2213,
    double &deldelg2222,
    double &deldelg2223,
    double &deldelg2233,
    double &deldelg2311,
    double &deldelg2312,
    double &deldelg2313,
    double &deldelg2322,
    double &deldelg2323,
    double &deldelg2333,
    double &deldelg3311,
    double &deldelg3312,
    double &deldelg3313,
    double &deldelg3322,
    double &deldelg3323,
    double &deldelg3333,
    double &delG11,
    double &delg111,
    double &delg112,
    double &delg113,
    double &delG12,
    double &delg122,
    double &delg123,
    double &delG13,
    double &delg133,
    double &delG21,
    double &delg211,
    double &delg212,
    double &delg213,
    double &delG22,
    double &delg222,
    double &delg223,
    double &delG23,
    double &delg233,
    double &delG31,
    double &delg311,
    double &delg312,
    double &delg313,
    double &delG32,
    double &delg322,
    double &delg323,
    double &delG33,
    double &delg333,
    double &dKhat1,
    double &dKhat2,
    double &dKhat3,
    double &dTheta1,
    double &dTheta2,
    double &dTheta3,
    double &G1,
    double &g11,
    double &g12,
    double &g13,
    double &G2,
    double &g22,
    double &g23,
    double &G3,
    double &g33,
    double &kappa1,
    double &kappa2,
    double &Khat,
    double &rA11,
    double &rA12,
    double &rA13,
    double &rA22,
    double &rA23,
    double &rA33,
    double &rchi,
    double &rG1,
    double &rg11,
    double &rg12,
    double &rg13,
    double &rG2,
    double &rg22,
    double &rg23,
    double &rG3,
    double &rg33,
    double &rKhat,
    double &rTheta,
    double &Theta);

static inline void z4c_contract_gamma(
    const double gxx, const double gxy, const double gxz,
    const double gyy, const double gyz, const double gzz,
    const double gxxx, const double gxyx, const double gxzx,
    const double gyyx, const double gyzx, const double gzzx,
    const double gxxy, const double gxyy, const double gxzy,
    const double gyyy, const double gyzy, const double gzzy,
    const double gxxz, const double gxyz, const double gxzz,
    const double gyyz, const double gyzz, const double gzzz,
    double &Gamxa, double &Gamya, double &Gamza)
{
    double det = gxx * gyy * gzz + gxy * gyz * gxz + gxz * gxy * gyz -
                 gxz * gyy * gxz - gxy * gxy * gzz - gxx * gyz * gyz;
    const double gupxx = (gyy * gzz - gyz * gyz) / det;
    const double gupxy = -(gxy * gzz - gyz * gxz) / det;
    const double gupxz = (gxy * gyz - gyy * gxz) / det;
    const double gupyy = (gxx * gzz - gxz * gxz) / det;
    const double gupyz = -(gxx * gyz - gxy * gxz) / det;
    const double gupzz = (gxx * gyy - gxy * gxy) / det;

    const double Gamxxx = 0.5 * (gupxx * gxxx + gupxy * (2.0 * gxyx - gxxy) + gupxz * (2.0 * gxzx - gxxz));
    const double Gamyxx = 0.5 * (gupxy * gxxx + gupyy * (2.0 * gxyx - gxxy) + gupyz * (2.0 * gxzx - gxxz));
    const double Gamzxx = 0.5 * (gupxz * gxxx + gupyz * (2.0 * gxyx - gxxy) + gupzz * (2.0 * gxzx - gxxz));

    const double Gamxyy = 0.5 * (gupxx * (2.0 * gxyy - gyyx) + gupxy * gyyy + gupxz * (2.0 * gyzy - gyyz));
    const double Gamyyy = 0.5 * (gupxy * (2.0 * gxyy - gyyx) + gupyy * gyyy + gupyz * (2.0 * gyzy - gyyz));
    const double Gamzyy = 0.5 * (gupxz * (2.0 * gxyy - gyyx) + gupyz * gyyy + gupzz * (2.0 * gyzy - gyyz));

    const double Gamxzz = 0.5 * (gupxx * (2.0 * gxzz - gzzx) + gupxy * (2.0 * gyzz - gzzy) + gupxz * gzzz);
    const double Gamyzz = 0.5 * (gupxy * (2.0 * gxzz - gzzx) + gupyy * (2.0 * gyzz - gzzy) + gupyz * gzzz);
    const double Gamzzz = 0.5 * (gupxz * (2.0 * gxzz - gzzx) + gupyz * (2.0 * gyzz - gzzy) + gupzz * gzzz);

    const double Gamxxy = 0.5 * (gupxx * gxxy + gupxy * gyyx + gupxz * (gxzy + gyzx - gxyz));
    const double Gamyxy = 0.5 * (gupxy * gxxy + gupyy * gyyx + gupyz * (gxzy + gyzx - gxyz));
    const double Gamzxy = 0.5 * (gupxz * gxxy + gupyz * gyyx + gupzz * (gxzy + gyzx - gxyz));

    const double Gamxxz = 0.5 * (gupxx * gxxz + gupxy * (gxyz + gyzx - gxzy) + gupxz * gzzx);
    const double Gamyxz = 0.5 * (gupxy * gxxz + gupyy * (gxyz + gyzx - gxzy) + gupyz * gzzx);
    const double Gamzxz = 0.5 * (gupxz * gxxz + gupyz * (gxyz + gyzx - gxzy) + gupzz * gzzx);

    const double Gamxyz = 0.5 * (gupxx * (gxyz + gxzy - gyzx) + gupxy * gyyz + gupxz * gzzy);
    const double Gamyyz = 0.5 * (gupxy * (gxyz + gxzy - gyzx) + gupyy * gyyz + gupyz * gzzy);
    const double Gamzyz = 0.5 * (gupxz * (gxyz + gxzy - gyzx) + gupyz * gyyz + gupzz * gzzy);

    Gamxa = gupxx * Gamxxx + gupyy * Gamxyy + gupzz * Gamxzz +
            2.0 * (gupxy * Gamxxy + gupxz * Gamxxz + gupyz * Gamxyz);
    Gamya = gupxx * Gamyxx + gupyy * Gamyyy + gupzz * Gamyzz +
            2.0 * (gupxy * Gamyxy + gupxz * Gamyxz + gupyz * Gamyyz);
    Gamza = gupxx * Gamzxx + gupyy * Gamzyy + gupzz * Gamzzz +
            2.0 * (gupxy * Gamzxy + gupxz * Gamzxz + gupyz * Gamzyz);
}

static int compute_rhs_z4c_cartesian(
    int *ex, double &T, double *X, double *Y, double *Z,
    double *chi_state, double *chi_constraints, double *trK,
    double *dxx, double *gxy, double *gxz, double *dyy, double *gyz, double *dzz,
    double *Axx, double *Axy, double *Axz, double *Ayy, double *Ayz, double *Azz,
    double *Gamx, double *Gamy, double *Gamz,
    double *Lap, double *betax, double *betay, double *betaz,
    double *dtSfx, double *dtSfy, double *dtSfz,
    double *TZ,
    double *chi_rhs, double *trK_rhs,
    double *gxx_rhs, double *gxy_rhs, double *gxz_rhs, double *gyy_rhs, double *gyz_rhs, double *gzz_rhs,
    double *Axx_rhs, double *Axy_rhs, double *Axz_rhs, double *Ayy_rhs, double *Ayz_rhs, double *Azz_rhs,
    double *Gamx_rhs, double *Gamy_rhs, double *Gamz_rhs,
    double *Lap_rhs, double *betax_rhs, double *betay_rhs, double *betaz_rhs,
    double *dtSfx_rhs, double *dtSfy_rhs, double *dtSfz_rhs,
    double *TZ_rhs,
    double *rho, double *Sx, double *Sy, double *Sz,
    double *Sxx, double *Sxy, double *Sxz, double *Syy, double *Syz, double *Szz,
    double *Gamxxx, double *Gamxxy, double *Gamxxz, double *Gamxyy, double *Gamxyz, double *Gamxzz,
    double *Gamyxx, double *Gamyxy, double *Gamyxz, double *Gamyyy, double *Gamyyz, double *Gamyzz,
    double *Gamzxx, double *Gamzxy, double *Gamzxz, double *Gamzyy, double *Gamzyz, double *Gamzzz,
    double *Rxx, double *Rxy, double *Rxz, double *Ryy, double *Ryz, double *Rzz,
    double *Hcon, double *Mxcon, double *Mycon, double *Mzcon, double *Gmxcon, double *Gmycon, double *Gmzcon,
    int &Symmetry, int &Lev, double &eps, int &co)
{
    (void)T;

    const int nx = ex[0];
    const int ny = ex[1];
    const int nz = ex[2];
    const int all = nx * ny * nz;

    double alpn1[all], chin1[all], gxx[all], gyy[all], gzz[all];
    double chix[all], chiy[all], chiz[all], chixx[all], chixy[all], chixz[all], chiyy[all], chiyz[all], chizz[all];
    double gxxx[all], gxyx[all], gxzx[all], gyyx[all], gyzx[all], gzzx[all];
    double gxxy[all], gxyy[all], gxzy[all], gyyy[all], gyzy[all], gzzy[all];
    double gxxz[all], gxyz[all], gxzz[all], gyyz[all], gyzz[all], gzzz[all];
    double gxxxx[all], gxxxy[all], gxxxz[all], gxxyy[all], gxxyz[all], gxxzz[all];
    double gxyxx[all], gxyxy[all], gxyxz[all], gxyyy[all], gxyyz[all], gxyzz[all];
    double gxzxx[all], gxzxy[all], gxzxz[all], gxzyy[all], gxzyz[all], gxzzz[all];
    double gyyxx[all], gyyxy[all], gyyxz[all], gyyyy[all], gyyyz[all], gyyzz[all];
    double gyzxx[all], gyzxy[all], gyzxz[all], gyzyy[all], gyzyz[all], gyzzz[all];
    double gzzxx[all], gzzxy[all], gzzxz[all], gzzyy[all], gzzyz[all], gzzzz[all];
    double Lapx[all], Lapy[all], Lapz[all], Lapxx[all], Lapxy[all], Lapxz[all], Lapyy[all], Lapyz[all], Lapzz[all];
    double betaxx[all], betaxy[all], betaxz[all], betayx[all], betayy[all], betayz[all], betazx[all], betazy[all], betazz[all];
    double dBxx[all], dBxy[all], dBxz[all], dByx[all], dByy[all], dByz[all], dBzx[all], dBzy[all], dBzz[all];
    double sfxxx[all], sfxxy[all], sfxxz[all], sfxyy[all], sfxyz[all], sfxzz[all];
    double sfyxx[all], sfyxy[all], sfyxz[all], sfyyy[all], sfyyz[all], sfyzz[all];
    double sfzxx[all], sfzxy[all], sfzxz[all], sfzyy[all], sfzyz[all], sfzzz[all];
    double Gamxx[all], Gamxy[all], Gamxz[all], Gamyx[all], Gamyy[all], Gamyz[all], Gamzx[all], Gamzy[all], Gamzz[all];
    double Kx[all], Ky[all], Kz[all], TZx[all], TZy[all], TZz[all];
    double Axxx[all], Axxy[all], Axxz[all], Axyx[all], Axyy[all], Axyz[all];
    double Axzx[all], Axzy[all], Axzz[all], Ayyx[all], Ayyy[all], Ayyz[all];
    double Ayzx[all], Ayzy[all], Ayzz[all], Azzx[all], Azzy[all], Azzz[all];

    const double SSS[3] = {1.0, 1.0, 1.0};
    const double AAS[3] = {-1.0, -1.0, 1.0};
    const double ASA[3] = {-1.0, 1.0, -1.0};
    const double SAA[3] = {1.0, -1.0, -1.0};
    const double ASS[3] = {-1.0, 1.0, 1.0};
    const double SAS[3] = {1.0, -1.0, 1.0};
    const double SSA[3] = {1.0, 1.0, -1.0};

    const double ONE = 1.0;
    const double TWO = 2.0;
    const double ZEO = 0.0;
    double chiDivfloor = 1.0e-5;

    double kappa1 = 2.0e-2;
    double kappa2 = 0.0;
    double FF = 0.75;
    double eta = 2.0;

    for (int idx = 0; idx < all; ++idx)
    {
        alpn1[idx] = Lap[idx] + ONE;
        chin1[idx] = chi_state[idx] + ONE;
        gxx[idx] = dxx[idx] + ONE;
        gyy[idx] = dyy[idx] + ONE;
        gzz[idx] = dzz[idx] + ONE;
    }

    fderivs(ex, betax, betaxx, betaxy, betaxz, X, Y, Z, -1.0, 1.0, 1.0, Symmetry, Lev);
    fderivs(ex, betay, betayx, betayy, betayz, X, Y, Z, 1.0, -1.0, 1.0, Symmetry, Lev);
    fderivs(ex, betaz, betazx, betazy, betazz, X, Y, Z, 1.0, 1.0, -1.0, Symmetry, Lev);
    fderivs(ex, dtSfx, dBxx, dBxy, dBxz, X, Y, Z, -1.0, 1.0, 1.0, Symmetry, Lev);
    fderivs(ex, dtSfy, dByx, dByy, dByz, X, Y, Z, 1.0, -1.0, 1.0, Symmetry, Lev);
    fderivs(ex, dtSfz, dBzx, dBzy, dBzz, X, Y, Z, 1.0, 1.0, -1.0, Symmetry, Lev);
    fderivs(ex, chi_state, chix, chiy, chiz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
    fderivs(ex, dxx, gxxx, gxxy, gxxz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
    fderivs(ex, gxy, gxyx, gxyy, gxyz, X, Y, Z, -1.0, -1.0, 1.0, Symmetry, Lev);
    fderivs(ex, gxz, gxzx, gxzy, gxzz, X, Y, Z, -1.0, 1.0, -1.0, Symmetry, Lev);
    fderivs(ex, dyy, gyyx, gyyy, gyyz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
    fderivs(ex, gyz, gyzx, gyzy, gyzz, X, Y, Z, 1.0, -1.0, -1.0, Symmetry, Lev);
    fderivs(ex, dzz, gzzx, gzzy, gzzz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);

    fdderivs(ex, dxx, gxxxx, gxxxy, gxxxz, gxxyy, gxxyz, gxxzz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
    fdderivs(ex, dyy, gyyxx, gyyxy, gyyxz, gyyyy, gyyyz, gyyzz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
    fdderivs(ex, dzz, gzzxx, gzzxy, gzzxz, gzzyy, gzzyz, gzzzz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
    fdderivs(ex, gxy, gxyxx, gxyxy, gxyxz, gxyyy, gxyyz, gxyzz, X, Y, Z, -1.0, -1.0, 1.0, Symmetry, Lev);
    fdderivs(ex, gxz, gxzxx, gxzxy, gxzxz, gxzyy, gxzyz, gxzzz, X, Y, Z, -1.0, 1.0, -1.0, Symmetry, Lev);
    fdderivs(ex, gyz, gyzxx, gyzxy, gyzxz, gyzyy, gyzyz, gyzzz, X, Y, Z, 1.0, -1.0, -1.0, Symmetry, Lev);

    fderivs(ex, Gamx, Gamxx, Gamxy, Gamxz, X, Y, Z, -1.0, 1.0, 1.0, Symmetry, Lev);
    fderivs(ex, Gamy, Gamyx, Gamyy, Gamyz, X, Y, Z, 1.0, -1.0, 1.0, Symmetry, Lev);
    fderivs(ex, Gamz, Gamzx, Gamzy, Gamzz, X, Y, Z, 1.0, 1.0, -1.0, Symmetry, Lev);

    fderivs(ex, Lap, Lapx, Lapy, Lapz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
    fderivs(ex, trK, Kx, Ky, Kz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
    fderivs(ex, TZ, TZx, TZy, TZz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);

    fdderivs(ex, betax, sfxxx, sfxxy, sfxxz, sfxyy, sfxyz, sfxzz, X, Y, Z, -1.0, 1.0, 1.0, Symmetry, Lev);
    fdderivs(ex, betay, sfyxx, sfyxy, sfyxz, sfyyy, sfyyz, sfyzz, X, Y, Z, 1.0, -1.0, 1.0, Symmetry, Lev);
    fdderivs(ex, betaz, sfzxx, sfzxy, sfzxz, sfzyy, sfzyz, sfzzz, X, Y, Z, 1.0, 1.0, -1.0, Symmetry, Lev);

    fdderivs(ex, chi_state, chixx, chixy, chixz, chiyy, chiyz, chizz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
    fdderivs(ex, Lap, Lapxx, Lapxy, Lapxz, Lapyy, Lapyz, Lapzz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);

    fderivs(ex, Axx, Axxx, Axxy, Axxz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
    fderivs(ex, Axy, Axyx, Axyy, Axyz, X, Y, Z, -1.0, -1.0, 1.0, Symmetry, Lev);
    fderivs(ex, Axz, Axzx, Axzy, Axzz, X, Y, Z, -1.0, 1.0, -1.0, Symmetry, Lev);
    fderivs(ex, Ayy, Ayyx, Ayyy, Ayyz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
    fderivs(ex, Ayz, Ayzx, Ayzy, Ayzz, X, Y, Z, 1.0, -1.0, -1.0, Symmetry, Lev);
    fderivs(ex, Azz, Azzx, Azzy, Azzz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);

    for (int idx = 0; idx < all; ++idx)
    {
        double point_kappa1 = 0.0;
        f_z4c_rhs_point(
            Axx[idx], Axy[idx], Axz[idx], Ayy[idx], Ayz[idx], Azz[idx],
            alpn1[idx], dtSfx[idx], dtSfy[idx], dtSfz[idx],
            betax[idx], betay[idx], betaz[idx],
            chin1[idx], chiDivfloor,
            Lapx[idx],
            Axxx[idx], Axyx[idx], Axzx[idx], Ayyx[idx], Ayzx[idx], Azzx[idx],
            Lapy[idx],
            Axxy[idx], Axyy[idx], Axzy[idx], Ayyy[idx], Ayzy[idx], Azzy[idx],
            Lapz[idx],
            Axxz[idx], Axyz[idx], Axzz[idx], Ayyz[idx], Ayzz[idx], Azzz[idx],
            betaxx[idx], dBxx[idx], betayx[idx], dByx[idx], betazx[idx], dBzx[idx],
            betaxy[idx], dBxy[idx], betayy[idx], dByy[idx], betazy[idx], dBzy[idx],
            betaxz[idx], dBxz[idx], betayz[idx], dByz[idx], betazz[idx], dBzz[idx],
            chix[idx], chiy[idx], chiz[idx],
            Lapxx[idx], Lapxy[idx], Lapxz[idx], Lapyy[idx], Lapyz[idx], Lapzz[idx],
            sfxxx[idx], sfyxx[idx], sfzxx[idx],
            sfxxy[idx], sfyxy[idx], sfzxy[idx],
            sfxxz[idx], sfyxz[idx], sfzxz[idx],
            sfxyy[idx], sfyyy[idx], sfzyy[idx],
            sfxyz[idx], sfyyz[idx], sfzyz[idx],
            sfxzz[idx], sfyzz[idx], sfzzz[idx],
            chixx[idx], chixy[idx], chixz[idx], chiyy[idx], chiyz[idx], chizz[idx],
            gxxxx[idx], gxyxx[idx], gxzxx[idx], gyyxx[idx], gyzxx[idx], gzzxx[idx],
            gxxxy[idx], gxyxy[idx], gxzxy[idx], gyyxy[idx], gyzxy[idx], gzzxy[idx],
            gxxxz[idx], gxyxz[idx], gxzxz[idx], gyyxz[idx], gyzxz[idx], gzzxz[idx],
            gxxyy[idx], gxyyy[idx], gxzyy[idx], gyyyy[idx], gyzyy[idx], gzzyy[idx],
            gxxyz[idx], gxyyz[idx], gxzyz[idx], gyyyz[idx], gyzyz[idx], gzzyz[idx],
            gxxzz[idx], gxyzz[idx], gxzzz[idx], gyyzz[idx], gyzzz[idx], gzzzz[idx],
            Gamxx[idx], gxxx[idx], gxyx[idx], gxzx[idx],
            Gamyx[idx], gyyx[idx], gyzx[idx],
            Gamzx[idx], gzzx[idx],
            Gamxy[idx], gxxy[idx], gxyy[idx], gxzy[idx],
            Gamyy[idx], gyyy[idx], gyzy[idx],
            Gamzy[idx], gzzy[idx],
            Gamxz[idx], gxxz[idx], gxyz[idx], gxzz[idx],
            Gamyz[idx], gyyz[idx], gyzz[idx],
            Gamzz[idx], gzzz[idx],
            Kx[idx], Ky[idx], Kz[idx],
            TZx[idx], TZy[idx], TZz[idx],
            Gamx[idx], gxx[idx], gxy[idx], gxz[idx],
            Gamy[idx], gyy[idx], gyz[idx],
            Gamz[idx], gzz[idx],
            point_kappa1, kappa2,
            trK[idx],
            Axx_rhs[idx], Axy_rhs[idx], Axz_rhs[idx], Ayy_rhs[idx], Ayz_rhs[idx], Azz_rhs[idx],
            chi_rhs[idx],
            Gamx_rhs[idx], gxx_rhs[idx], gxy_rhs[idx], gxz_rhs[idx],
            Gamy_rhs[idx], gyy_rhs[idx], gyz_rhs[idx],
            Gamz_rhs[idx], gzz_rhs[idx], trK_rhs[idx], TZ_rhs[idx], TZ[idx]);
    }

    for (int idx = 0; idx < all; ++idx)
        Lap_rhs[idx] = -TWO * alpn1[idx] * trK[idx];

#if (GAUGE == 0)
    for (int idx = 0; idx < all; ++idx)
    {
        betax_rhs[idx] = FF * dtSfx[idx];
        betay_rhs[idx] = FF * dtSfy[idx];
        betaz_rhs[idx] = FF * dtSfz[idx];
        dtSfx_rhs[idx] = Gamx_rhs[idx] - eta * dtSfx[idx];
        dtSfy_rhs[idx] = Gamy_rhs[idx] - eta * dtSfy[idx];
        dtSfz_rhs[idx] = Gamz_rhs[idx] - eta * dtSfz[idx];
    }
#elif (GAUGE == 1)
    for (int idx = 0; idx < all; ++idx)
    {
        betax_rhs[idx] = Gamx[idx] - eta * betax[idx];
        betay_rhs[idx] = Gamy[idx] - eta * betay[idx];
        betaz_rhs[idx] = Gamz[idx] - eta * betaz[idx];
        dtSfx_rhs[idx] = ZEO;
        dtSfy_rhs[idx] = ZEO;
        dtSfz_rhs[idx] = ZEO;
    }
#else
#error "z4c_rhs_c.C currently supports GAUGE == 0 or GAUGE == 1 for Z4C"
#endif

    lopsided(ex, X, Y, Z, gxx, gxx_rhs, betax, betay, betaz, Symmetry, SSS);
    lopsided(ex, X, Y, Z, gxy, gxy_rhs, betax, betay, betaz, Symmetry, AAS);
    lopsided(ex, X, Y, Z, gxz, gxz_rhs, betax, betay, betaz, Symmetry, ASA);
    lopsided(ex, X, Y, Z, gyy, gyy_rhs, betax, betay, betaz, Symmetry, SSS);
    lopsided(ex, X, Y, Z, gyz, gyz_rhs, betax, betay, betaz, Symmetry, SAA);
    lopsided(ex, X, Y, Z, gzz, gzz_rhs, betax, betay, betaz, Symmetry, SSS);

    lopsided(ex, X, Y, Z, Axx, Axx_rhs, betax, betay, betaz, Symmetry, SSS);
    lopsided(ex, X, Y, Z, Axy, Axy_rhs, betax, betay, betaz, Symmetry, AAS);
    lopsided(ex, X, Y, Z, Axz, Axz_rhs, betax, betay, betaz, Symmetry, ASA);
    lopsided(ex, X, Y, Z, Ayy, Ayy_rhs, betax, betay, betaz, Symmetry, SSS);
    lopsided(ex, X, Y, Z, Ayz, Ayz_rhs, betax, betay, betaz, Symmetry, SAA);
    lopsided(ex, X, Y, Z, Azz, Azz_rhs, betax, betay, betaz, Symmetry, SSS);

    lopsided(ex, X, Y, Z, chi_state, chi_rhs, betax, betay, betaz, Symmetry, SSS);
    lopsided(ex, X, Y, Z, trK, trK_rhs, betax, betay, betaz, Symmetry, SSS);

    lopsided(ex, X, Y, Z, Gamx, Gamx_rhs, betax, betay, betaz, Symmetry, ASS);
    lopsided(ex, X, Y, Z, Gamy, Gamy_rhs, betax, betay, betaz, Symmetry, SAS);
    lopsided(ex, X, Y, Z, Gamz, Gamz_rhs, betax, betay, betaz, Symmetry, SSA);

    lopsided(ex, X, Y, Z, Lap, Lap_rhs, betax, betay, betaz, Symmetry, SSS);
    lopsided(ex, X, Y, Z, betax, betax_rhs, betax, betay, betaz, Symmetry, ASS);
    lopsided(ex, X, Y, Z, betay, betay_rhs, betax, betay, betaz, Symmetry, SAS);
    lopsided(ex, X, Y, Z, betaz, betaz_rhs, betax, betay, betaz, Symmetry, SSA);
#if (GAUGE == 0)
    lopsided(ex, X, Y, Z, dtSfx, dtSfx_rhs, betax, betay, betaz, Symmetry, ASS);
    lopsided(ex, X, Y, Z, dtSfy, dtSfy_rhs, betax, betay, betaz, Symmetry, SAS);
    lopsided(ex, X, Y, Z, dtSfz, dtSfz_rhs, betax, betay, betaz, Symmetry, SSA);
#endif
    lopsided(ex, X, Y, Z, TZ, TZ_rhs, betax, betay, betaz, Symmetry, SSS);

    for (int idx = 0; idx < all; ++idx)
    {
        double Gamxa = 0.0, Gamya = 0.0, Gamza = 0.0;
        z4c_contract_gamma(
            gxx[idx], gxy[idx], gxz[idx], gyy[idx], gyz[idx], gzz[idx],
            gxxx[idx], gxyx[idx], gxzx[idx], gyyx[idx], gyzx[idx], gzzx[idx],
            gxxy[idx], gxyy[idx], gxzy[idx], gyyy[idx], gyzy[idx], gzzy[idx],
            gxxz[idx], gxyz[idx], gxzz[idx], gyyz[idx], gyzz[idx], gzzz[idx],
            Gamxa, Gamya, Gamza);

        TZ_rhs[idx] -= alpn1[idx] * (TWO + kappa2) * kappa1 * TZ[idx];
        trK_rhs[idx] += alpn1[idx] * kappa1 * (ONE - kappa2) * TZ[idx];
        Gamx_rhs[idx] -= TWO * alpn1[idx] * kappa1 * (Gamx[idx] - Gamxa);
        Gamy_rhs[idx] -= TWO * alpn1[idx] * kappa1 * (Gamy[idx] - Gamya);
        Gamz_rhs[idx] -= TWO * alpn1[idx] * kappa1 * (Gamz[idx] - Gamza);
    }

    if (eps > 0.0)
    {
        kodis(ex, X, Y, Z, chi_state, chi_rhs, SSS, Symmetry, eps);
        kodis(ex, X, Y, Z, trK, trK_rhs, SSS, Symmetry, eps);
        kodis(ex, X, Y, Z, gxx, gxx_rhs, SSS, Symmetry, eps);
        kodis(ex, X, Y, Z, gxy, gxy_rhs, AAS, Symmetry, eps);
        kodis(ex, X, Y, Z, gxz, gxz_rhs, ASA, Symmetry, eps);
        kodis(ex, X, Y, Z, gyy, gyy_rhs, SSS, Symmetry, eps);
        kodis(ex, X, Y, Z, gyz, gyz_rhs, SAA, Symmetry, eps);
        kodis(ex, X, Y, Z, gzz, gzz_rhs, SSS, Symmetry, eps);
        kodis(ex, X, Y, Z, Axx, Axx_rhs, SSS, Symmetry, eps);
        kodis(ex, X, Y, Z, Axy, Axy_rhs, AAS, Symmetry, eps);
        kodis(ex, X, Y, Z, Axz, Axz_rhs, ASA, Symmetry, eps);
        kodis(ex, X, Y, Z, Ayy, Ayy_rhs, SSS, Symmetry, eps);
        kodis(ex, X, Y, Z, Ayz, Ayz_rhs, SAA, Symmetry, eps);
        kodis(ex, X, Y, Z, Azz, Azz_rhs, SSS, Symmetry, eps);
        kodis(ex, X, Y, Z, Gamx, Gamx_rhs, ASS, Symmetry, eps);
        kodis(ex, X, Y, Z, Gamy, Gamy_rhs, SAS, Symmetry, eps);
        kodis(ex, X, Y, Z, Gamz, Gamz_rhs, SSA, Symmetry, eps);
        kodis(ex, X, Y, Z, Lap, Lap_rhs, SSS, Symmetry, eps);
        kodis(ex, X, Y, Z, betax, betax_rhs, ASS, Symmetry, eps);
        kodis(ex, X, Y, Z, betay, betay_rhs, SAS, Symmetry, eps);
        kodis(ex, X, Y, Z, betaz, betaz_rhs, SSA, Symmetry, eps);
#if (GAUGE == 0)
        kodis(ex, X, Y, Z, dtSfx, dtSfx_rhs, ASS, Symmetry, eps);
        kodis(ex, X, Y, Z, dtSfy, dtSfy_rhs, SAS, Symmetry, eps);
        kodis(ex, X, Y, Z, dtSfz, dtSfz_rhs, SSA, Symmetry, eps);
#endif
        kodis(ex, X, Y, Z, TZ, TZ_rhs, SSS, Symmetry, eps);
    }

    if (co == 0)
    {
#if (ABV == 0)
        f_ricci_gamma(ex, X, Y, Z,
                      chi_constraints,
                      dxx, gxy, gxz, dyy, gyz, dzz,
                      Gamx, Gamy, Gamz,
                      Gamxxx, Gamxxy, Gamxxz, Gamxyy, Gamxyz, Gamxzz,
                      Gamyxx, Gamyxy, Gamyxz, Gamyyy, Gamyyz, Gamyzz,
                      Gamzxx, Gamzxy, Gamzxz, Gamzyy, Gamzyz, Gamzzz,
                      Rxx, Rxy, Rxz, Ryy, Ryz, Rzz,
                      Symmetry);
#endif
        f_constraint_bssn(ex, X, Y, Z,
                          chi_constraints, trK,
                          dxx, gxy, gxz, dyy, gyz, dzz,
                          Axx, Axy, Axz, Ayy, Ayz, Azz,
                          Gamx, Gamy, Gamz,
                          Lap, betax, betay, betaz, rho, Sx, Sy, Sz,
                          Gamxxx, Gamxxy, Gamxxz, Gamxyy, Gamxyz, Gamxzz,
                          Gamyxx, Gamyxy, Gamyxz, Gamyyy, Gamyyz, Gamyzz,
                          Gamzxx, Gamzxy, Gamzxz, Gamzyy, Gamzyz, Gamzzz,
                          Rxx, Rxy, Rxz, Ryy, Ryz, Rzz,
                          Hcon, Mxcon, Mycon, Mzcon, Gmxcon, Gmycon, Gmzcon,
                          Symmetry);
    }

    return 0;
}

extern "C" int f_compute_rhs_Z4c(int *ex, double &T,
                                 double *X, double *Y, double *Z,
                                 double *chi, double *trK,
                                 double *dxx, double *gxy, double *gxz, double *dyy, double *gyz, double *dzz,
                                 double *Axx, double *Axy, double *Axz, double *Ayy, double *Ayz, double *Azz,
                                 double *Gamx, double *Gamy, double *Gamz,
                                 double *Lap, double *betax, double *betay, double *betaz,
                                 double *dtSfx, double *dtSfy, double *dtSfz,
                                 double *TZ,
                                 double *chi_rhs, double *trK_rhs,
                                 double *gxx_rhs, double *gxy_rhs, double *gxz_rhs, double *gyy_rhs, double *gyz_rhs, double *gzz_rhs,
                                 double *Axx_rhs, double *Axy_rhs, double *Axz_rhs, double *Ayy_rhs, double *Ayz_rhs, double *Azz_rhs,
                                 double *Gamx_rhs, double *Gamy_rhs, double *Gamz_rhs,
                                 double *Lap_rhs, double *betax_rhs, double *betay_rhs, double *betaz_rhs,
                                 double *dtSfx_rhs, double *dtSfy_rhs, double *dtSfz_rhs,
                                 double *TZ_rhs,
                                 double *rho, double *Sx, double *Sy, double *Sz,
                                 double *Sxx, double *Sxy, double *Sxz, double *Syy, double *Syz, double *Szz,
                                 double *Gamxxx, double *Gamxxy, double *Gamxxz, double *Gamxyy, double *Gamxyz, double *Gamxzz,
                                 double *Gamyxx, double *Gamyxy, double *Gamyxz, double *Gamyyy, double *Gamyyz, double *Gamyzz,
                                 double *Gamzxx, double *Gamzxy, double *Gamzxz, double *Gamzyy, double *Gamzyz, double *Gamzzz,
                                 double *Rxx, double *Rxy, double *Rxz, double *Ryy, double *Ryz, double *Rzz,
                                 double *Hcon, double *Mxcon, double *Mycon, double *Mzcon, double *Gmxcon, double *Gmycon, double *Gmzcon,
                                 int &Symmetry, int &Lev, double &eps, int &co)
{
    return compute_rhs_z4c_cartesian(
        ex, T, X, Y, Z,
        chi, chi, trK,
        dxx, gxy, gxz, dyy, gyz, dzz,
        Axx, Axy, Axz, Ayy, Ayz, Azz,
        Gamx, Gamy, Gamz,
        Lap, betax, betay, betaz,
        dtSfx, dtSfy, dtSfz,
        TZ,
        chi_rhs, trK_rhs,
        gxx_rhs, gxy_rhs, gxz_rhs, gyy_rhs, gyz_rhs, gzz_rhs,
        Axx_rhs, Axy_rhs, Axz_rhs, Ayy_rhs, Ayz_rhs, Azz_rhs,
        Gamx_rhs, Gamy_rhs, Gamz_rhs,
        Lap_rhs, betax_rhs, betay_rhs, betaz_rhs,
        dtSfx_rhs, dtSfy_rhs, dtSfz_rhs,
        TZ_rhs,
        rho, Sx, Sy, Sz,
        Sxx, Sxy, Sxz, Syy, Syz, Szz,
        Gamxxx, Gamxxy, Gamxxz, Gamxyy, Gamxyz, Gamxzz,
        Gamyxx, Gamyxy, Gamyxz, Gamyyy, Gamyyz, Gamyzz,
        Gamzxx, Gamzxy, Gamzxz, Gamzyy, Gamzyz, Gamzzz,
        Rxx, Rxy, Rxz, Ryy, Ryz, Rzz,
        Hcon, Mxcon, Mycon, Mzcon, Gmxcon, Gmycon, Gmzcon,
        Symmetry, Lev, eps, co);
}

extern "C" int f_compute_rhs_Z4cnot(int *ex, double &T,
                                    double *X, double *Y, double *Z,
                                    double *chi, double *trK,
                                    double *dxx, double *gxy, double *gxz, double *dyy, double *gyz, double *dzz,
                                    double *Axx, double *Axy, double *Axz, double *Ayy, double *Ayz, double *Azz,
                                    double *Gamx, double *Gamy, double *Gamz,
                                    double *Lap, double *betax, double *betay, double *betaz,
                                    double *dtSfx, double *dtSfy, double *dtSfz,
                                    double *TZ,
                                    double *chi_rhs, double *trK_rhs,
                                    double *gxx_rhs, double *gxy_rhs, double *gxz_rhs, double *gyy_rhs, double *gyz_rhs, double *gzz_rhs,
                                    double *Axx_rhs, double *Axy_rhs, double *Axz_rhs, double *Ayy_rhs, double *Ayz_rhs, double *Azz_rhs,
                                    double *Gamx_rhs, double *Gamy_rhs, double *Gamz_rhs,
                                    double *Lap_rhs, double *betax_rhs, double *betay_rhs, double *betaz_rhs,
                                    double *dtSfx_rhs, double *dtSfy_rhs, double *dtSfz_rhs,
                                    double *TZ_rhs,
                                    double *rho, double *Sx, double *Sy, double *Sz,
                                    double *Sxx, double *Sxy, double *Sxz, double *Syy, double *Syz, double *Szz,
                                    double *Gamxxx, double *Gamxxy, double *Gamxxz, double *Gamxyy, double *Gamxyz, double *Gamxzz,
                                    double *Gamyxx, double *Gamyxy, double *Gamyxz, double *Gamyyy, double *Gamyyz, double *Gamyzz,
                                    double *Gamzxx, double *Gamzxy, double *Gamzxz, double *Gamzyy, double *Gamzyz, double *Gamzzz,
                                    double *Rxx, double *Rxy, double *Rxz, double *Ryy, double *Ryz, double *Rzz,
                                    double *Hcon, double *Mxcon, double *Mycon, double *Mzcon, double *Gmxcon, double *Gmycon, double *Gmzcon,
                                    int &Symmetry, int &Lev, double &eps, int &co, double &chitiny)
{
    const int all = ex[0] * ex[1] * ex[2];
    std::vector<double> chi_clamped(chi, chi + all);
    f_lowerboundset(ex, chi_clamped.data(), chitiny);

    const int ret = compute_rhs_z4c_cartesian(
        ex, T, X, Y, Z,
        chi_clamped.data(), chi, trK,
        dxx, gxy, gxz, dyy, gyz, dzz,
        Axx, Axy, Axz, Ayy, Ayz, Azz,
        Gamx, Gamy, Gamz,
        Lap, betax, betay, betaz,
        dtSfx, dtSfy, dtSfz,
        TZ,
        chi_rhs, trK_rhs,
        gxx_rhs, gxy_rhs, gxz_rhs, gyy_rhs, gyz_rhs, gzz_rhs,
        Axx_rhs, Axy_rhs, Axz_rhs, Ayy_rhs, Ayz_rhs, Azz_rhs,
        Gamx_rhs, Gamy_rhs, Gamz_rhs,
        Lap_rhs, betax_rhs, betay_rhs, betaz_rhs,
        dtSfx_rhs, dtSfy_rhs, dtSfz_rhs,
        TZ_rhs,
        rho, Sx, Sy, Sz,
        Sxx, Sxy, Sxz, Syy, Syz, Szz,
        Gamxxx, Gamxxy, Gamxxz, Gamxyy, Gamxyz, Gamxzz,
        Gamyxx, Gamyxy, Gamyxz, Gamyyy, Gamyyz, Gamyzz,
        Gamzxx, Gamzxy, Gamzxz, Gamzyy, Gamzyz, Gamzzz,
        Rxx, Rxy, Rxz, Ryy, Ryz, Rzz,
        Hcon, Mxcon, Mycon, Mzcon, Gmxcon, Gmycon, Gmzcon,
        Symmetry, Lev, eps, co);

    if (ret != 0 || co != 0)
        return ret;

#if (ABV == 0)
    f_ricci_gamma(ex, X, Y, Z,
                  chi,
                  dxx, gxy, gxz, dyy, gyz, dzz,
                  Gamx, Gamy, Gamz,
                  Gamxxx, Gamxxy, Gamxxz, Gamxyy, Gamxyz, Gamxzz,
                  Gamyxx, Gamyxy, Gamyxz, Gamyyy, Gamyyz, Gamyzz,
                  Gamzxx, Gamzxy, Gamzxz, Gamzyy, Gamzyz, Gamzzz,
                  Rxx, Rxy, Rxz, Ryy, Ryz, Rzz,
                  Symmetry);
#endif
    f_constraint_bssn(ex, X, Y, Z,
                      chi, trK,
                      dxx, gxy, gxz, dyy, gyz, dzz,
                      Axx, Axy, Axz, Ayy, Ayz, Azz,
                      Gamx, Gamy, Gamz,
                      Lap, betax, betay, betaz, rho, Sx, Sy, Sz,
                      Gamxxx, Gamxxy, Gamxxz, Gamxyy, Gamxyz, Gamxzz,
                      Gamyxx, Gamyxy, Gamyxz, Gamyyy, Gamyyz, Gamyzz,
                      Gamzxx, Gamzxy, Gamzxz, Gamzyy, Gamzyz, Gamzzz,
                      Rxx, Rxy, Rxz, Ryy, Ryz, Rzz,
                      Hcon, Mxcon, Mycon, Mzcon, Gmxcon, Gmycon, Gmzcon,
                      Symmetry);
    return ret;
}
