/*
 * 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/>.
 */

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

#include "brill_data.h"
#include "internal.h"

#define GRID(arr) (arr[i * nb_coords_rho + j])


/* φφ components */

#define COMP phiphi

#define DIFF 0
#define DO_EVAL val = SQR(rrho) * psi4
#include "eval_metric_template.c"
#undef DIFF
#undef DO_EVAL

#define DIFF 3
#define DO_EVAL val = 4 * SQR(rrho) * psi3 * GRID(psi[3])
#include "eval_metric_template.c"
#undef DIFF
#undef DO_EVAL

#define DIFF 6
#define DO_EVAL val = 4 * SQR(rrho) * (GRID(psi[6]) * psi3 + 3 * psi2 * SQR(GRID(psi[3])))
#include "eval_metric_template.c"
#undef DIFF
#undef DO_EVAL

#define DIFF 1
#define DO_EVAL val = 2 * rrho * psi4 + 4 * SQR(rrho) * psi3 * GRID(psi[1])
#include "eval_metric_template.c"
#undef DIFF
#undef DO_EVAL

#define DIFF 4
#define DO_EVAL val =  8 * rrho * psi3 * GRID(psi[3]) +\
                      12 * SQR(rrho) * psi2 * GRID(psi[3]) * GRID(psi[1]) +\
                       4 * SQR(rrho) * psi3 * GRID(psi[4])
#include "eval_metric_template.c"
#undef DIFF
#undef DO_EVAL

#define DIFF 2
#define DO_EVAL val =  4 * SQR(rrho) * psi3 * GRID(psi[2]) +\
                      12 * SQR(rrho * ppsi * GRID(psi[1])) +\
                      16 * rrho * psi3 * GRID(psi[1]) + 2 * psi4
#include "eval_metric_template.c"
#undef DIFF
#undef DO_EVAL

#undef COMP

/* ρρ components */
#define COMP rhorho

#define DIFF 0
#define DO_EVAL val = base
#include "eval_metric_template.c"
#undef DIFF
#undef DO_EVAL

#define DIFF 3
#define DO_EVAL val = 2 * base * dqz
#include "eval_metric_template.c"
#undef DIFF
#undef DO_EVAL

#define DIFF 6
#define DO_EVAL val = 4 * SQR(dqz) * base + 2 * base * dq2z
#include "eval_metric_template.c"
#undef DIFF
#undef DO_EVAL

#define DIFF 1
#define DO_EVAL val = 2 * base * dqrho
#include "eval_metric_template.c"
#undef DIFF
#undef DO_EVAL

#define DIFF 4
#define DO_EVAL val = 4 * base * dqrho * dqz + 2 * base * dq2rho_z
#include "eval_metric_template.c"
#undef DIFF
#undef DO_EVAL

#define DIFF 2
#define DO_EVAL val = 4 * SQR(dqrho) * base + 2 * base * dq2rho
#include "eval_metric_template.c"
#undef DIFF
#undef DO_EVAL

#undef COMP

typedef void (*EvalMetricFunc)(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]);

static const EvalMetricFunc eval_metric_phiphi[9] = {
    [0] = eval_metric_phiphi_0,
    [1] = eval_metric_phiphi_1,
    [2] = eval_metric_phiphi_2,
    [3] = eval_metric_phiphi_3,
    [4] = eval_metric_phiphi_4,
    [6] = eval_metric_phiphi_6,
};

static const EvalMetricFunc eval_metric_rhorho[9] = {
    [0] = eval_metric_rhorho_0,
    [1] = eval_metric_rhorho_1,
    [2] = eval_metric_rhorho_2,
    [3] = eval_metric_rhorho_3,
    [4] = eval_metric_rhorho_4,
    [6] = eval_metric_rhorho_6,
};

void bdi_eval_metric(const BDContext *bd,
                     const double *rho, int nb_coords_rho,
                     const double *z,   int nb_coords_z,
                     enum BDMetricComponent comp,
                     const unsigned int diff_order[2],
                     double *psi[9], double *q[9],
                     double *out, ptrdiff_t out_stride)
{
    int diff_idx = diff_order[0] * 3 + diff_order[1];
    const EvalMetricFunc *functab;

    if (comp != BD_METRIC_COMPONENT_RHORHO &&
        comp != BD_METRIC_COMPONENT_ZZ     &&
        comp != BD_METRIC_COMPONENT_PHIPHI) {
        memset(out, 0, out_stride * nb_coords_z * sizeof(*out));
        return;
    }

    functab = (comp == BD_METRIC_COMPONENT_PHIPHI) ?
              eval_metric_phiphi : eval_metric_rhorho;

    functab[diff_idx](bd, out, out_stride, nb_coords_rho, rho,
                      nb_coords_z, z, psi, q);
}