/*
 * Copyright 2016 Anton Khirnov <anton@khirnov.net>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

/**
 * @file
 * the template for the loop over the grid points
*/

#define FUNC3(a, b, c) a ## _ ## b ## _ ## c
#define FUNC2(a, b, c) FUNC3(a, b, c)
#define FUNC(x) FUNC2(x, COMP, DIFF)

#define DIFF_RHO (DIFF % 3)
#define DIFF_Z   (DIFF / 3)

static void FUNC(eval_metric)(const BDContext *bd, double *out, ptrdiff_t out_stride,
                              int nb_coords_rho, const double *rho,
                              int nb_coords_z,   const double *z,
                              double *psi[9], double *q[9])
{
    BDPriv *s = bd->priv;

    for (int i = 0; i < nb_coords_z; i++) {
        const double zz = z[i];
        double *dst = out + i * out_stride;

        for (int j = 0; j < nb_coords_rho; j++) {
            const int idx = i * nb_coords_rho + j;
            const double rrho = rho[j];

            const double ppsi = psi[0][idx];
            const double psi2 = SQR(ppsi);
            const double psi3 = psi2 * ppsi;
            const double psi4 = SQR(psi2);

            const double base     = psi4 * exp(2 * q[0][idx]);
            const double dqrho    = DIFF_RHO > 0 ? q[1][idx] + 2 * psi[1][idx] / ppsi : 0;
            const double dqz      = DIFF_Z   > 0 ? q[3][idx] + 2 * psi[3][idx] / ppsi : 0;
            const double dq2rho   = DIFF_RHO > 1 ? q[2][idx] + 2 * psi[2][idx] / ppsi - 2 * SQR(psi[1][idx] / ppsi) : 0;
            const double dq2z     = DIFF_Z > 1 ? q[6][idx] + 2 * psi[6][idx] / ppsi - 2 * SQR(psi[3][idx] / ppsi) : 0;
            const double dq2rho_z = DIFF_RHO && DIFF_Z ? q[4][idx] + 2 * psi[4][idx] / ppsi - 2 * psi[1][idx] * psi[3][idx] / psi2 : 0;

            double val;

            DO_EVAL;

            *dst++ = val;
        }
    }
}

#undef FUNC3
#undef FUNC2
#undef FUNC