#include <math.h>

#include "cctk.h"

void apply_dissipation(const cGH *gh, const double *var, double *rhs, double dx[3],
                       int order, double *epsdis)
{
    double inv_dx[3] = { 1.0 / dx[0], 1.0 / dx[1], 1.0 / dx[2] };

    int ni = gh->cctk_lsh[0];
    int nj = gh->cctk_lsh[1];
    int nk = gh->cctk_lsh[2];

    double *y;
    int y_idx = -1;

    y = CCTK_VarDataPtr(gh, 0, "grid::y");
    for (int j = 0; j < nj; j++) {
        if (fabs(y[CCTK_GFINDEX3D(gh, 0, j, 0)]) < 1e-8) {
            y_idx = j;
            break;
        }
    }

    switch (order) {
    case 1:
#pragma omp parallel for
        for (int k = 1; k < nk - 1; k++)
            for (int i = 1; i < nk - 1; i++) {
                int idx   = CCTK_GFINDEX3D(gh, i, y_idx, k);
                rhs[idx] += epsdis[idx] *
                            ((var[idx - 1]           - 2.0 * var[idx] + var[idx + 1]) * inv_dx[0] +
                             (var[idx - 1 * ni * nj] - 2.0 * var[idx] + var[idx + 1 * ni * nj]) * inv_dx[2]);
                }
        break;
    case 3:
#pragma omp parallel for
        for (int k = 2; k < nk - 2; k++)
            for (int i = 2; i < nk - 2; i++) {
                int idx   = CCTK_GFINDEX3D(gh, i, y_idx, k);
                rhs[idx] -= epsdis[idx] / 16.0 *
                            ((var[idx - 2]           - 4.0 * var[idx - 1]            + 6.0 * var[idx] - 4.0 * var[idx + 1]           + var[idx + 2]) * inv_dx[0] +
                             (var[idx - 2 * ni * nj] - 4.0 * var[idx - 1 * ni * nj]  + 6.0 * var[idx] - 4.0 * var[idx + 1 * ni * nj] + var[idx + 2 * ni * nj]) * inv_dx[2]);
                }
        break;
    case 5:
#pragma omp parallel for
        for (int k = 3; k < nk - 3; k++)
            for (int i = 3; i < nk - 3; i++) {
                int idx   = CCTK_GFINDEX3D(gh, i, y_idx, k);
                rhs[idx] += epsdis[idx] / 64 *
                            ((var[idx - 3]           - 6 * var[idx - 2]           + 15 * var[idx - 1]           - 20 * var[idx] + 15 * var[idx + 1]           - 6 * var[idx + 2]           + var[idx + 3]) * inv_dx[0] + 
                             (var[idx - 3 * ni * nj] - 6 * var[idx - 2 * ni * nj] + 15 * var[idx - 1 * ni * nj] - 20 * var[idx] + 15 * var[idx + 1 * ni * nj] - 6 * var[idx + 2 * ni * nj] + var[idx + 3 * ni * nj]) * inv_dx[2]);
//              + (var(i,j-3,k) - 6*var(i,j-2,k) + 15*var(i,j-1,k) - 20*var(i,j,k) + 15*var(i,j+1,k) - 6*var(i,j+2,k) + var(i,j+3,k) ) * idx(2)
                }
        break;
    case 9:
#pragma omp parallel for
        for (int k = 5; k < nk - 5; k++)
            for (int i = 5; i < nk - 5; i++) {
                int idx   = CCTK_GFINDEX3D(gh, i, y_idx, k);
                rhs[idx] += epsdis[idx] / 1024 *
                            ((var[idx - 5] - 10 * var[idx - 4] + 45 * var[idx - 3] - 120 * var[idx - 2] + 210 * var[idx - 1] - 252 * var[idx] +
                              var[idx + 5] - 10 * var[idx + 4] + 45 * var[idx + 3] - 120 * var[idx + 2] + 210 * var[idx + 1]) * inv_dx[0] +
                             (var[idx - 5 * ni * nj] - 10 * var[idx - 4 * ni * nj] + 45 * var[idx - 3 * ni * nj] - 120 * var[idx - 2 * ni * nj] + 210 * var[idx - 1 * ni * nj] - 252 * var[idx] +
                              var[idx + 5 * ni * nj] - 10 * var[idx + 4 * ni * nj] + 45 * var[idx + 3 * ni * nj] - 120 * var[idx + 2 * ni * nj] + 210 * var[idx + 1 * ni * nj]) * inv_dx[2]);
//$              + (var(i,j-5,k) - 10*var(i,j-4,k) + 45*var(i,j-3,k) - 120*var(i,j-2,k) + 210*var(i,j-1,k) - 252*var(i,j,k) + 210*var(i,j+1,k) - 120*var(i,j+2,k) + 45*var(i,j+3,k) - 10*var(i,j+4,k) + var(i,j+5,k)) * idx(2)
            }
        break;
    default:
        CCTK_WARN(0, "internal error");
    }
}