#include <ctype.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_spline.h>
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"

#define MAX(x,y) (x) > (y) ? (x) : (y)
#define SQR(x) ((x)*(x))

#define MASS 1
/*
 * the constant C from e.g. 0908.1063v4
 */
#define TRUMPET_CONST (3*sqrt(3)*SQR(MASS)/4)
/*
 * small number to avoid r=0 singularities
 */
#define EPS 1E-08
/*
 * isotropic/coordinate radius
 */
#define ISO_R(index) (sqrt(SQR(x[index]) + SQR(y[index]) + SQR(z[index]) + EPS))

static inline double r_areal_to_isotropic(double r, double m)
{
    double par   = sqrt(4*r*r + 4*m*r + 3*m*m);
    double term1 = (2*r + m + par)/4;
    double term2 = (4 + 3*M_SQRT2)*(2*r - 3*m)/(8*r + 6*m + 3*M_SQRT2*par);
    return term1*pow(term2, M_SQRT1_2);
}

/**
 * get maximal isotropic r
 */
static CCTK_REAL get_max_r(CCTK_INT dim)
{
    CCTK_REAL phys_min[dim], phys_max[dim], ext_min[dim], ext_max[dim], int_min[dim], int_max[dim], step[dim];
    CCTK_REAL max = 0;

    GetDomainSpecification(dim, phys_min, phys_max, int_min, int_max, ext_min, ext_max, step);

    for (int i = 0; i < dim; i++) {
        max = MAX(max, abs(phys_min[i]));
        max = MAX(max, abs(phys_max[i]));
    }
    return sqrt(dim)*2*max;
}

static void get_r_tables(double **pr_iso, double **pr_areal, CCTK_INT *size, CCTK_INT dim)
{
#define STEP 0.0001
    double *r_iso, *r_areal, max = get_max_r(dim);
    int count = max*2/STEP, i = 0;

    r_iso   = malloc(count*sizeof(*r_iso));
    r_areal = malloc(count*sizeof(*r_areal));

    for (i = 0;;i++) {
        if (i >= count) {
            count *= 2;
            r_iso = realloc(r_iso, count*sizeof(*r_iso));
            r_areal = realloc(r_areal, count*sizeof(*r_areal));
        }

        r_areal[i] = i*STEP + 1.5;
        r_iso[i]   = r_areal_to_isotropic(r_areal[i], 1);
        if (r_iso[i] > max)
            break;
    }
    *pr_iso   = r_iso;
    *pr_areal = r_areal;
    *size     = i;
}

void trumpet_data(CCTK_ARGUMENTS)
{
    DECLARE_CCTK_ARGUMENTS;
    DECLARE_CCTK_PARAMETERS;

    gsl_interp_accel *acc;
    gsl_spline    *spline;

    double *r_iso, *r_areal;
    CCTK_INT size;

    get_r_tables(&r_iso, &r_areal, &size, cctk_dim);

    spline = gsl_spline_alloc(gsl_interp_cspline, size);
    gsl_spline_init(spline, r_iso, r_areal, size);
    acc = gsl_interp_accel_alloc();

    memset(gxy, 0, sizeof(*gxy)*CCTK_GFINDEX3D(cctkGH, cctk_lsh[0]-1, cctk_lsh[1]-1, cctk_lsh[2]-1));
    memset(gxz, 0, sizeof(*gxy)*CCTK_GFINDEX3D(cctkGH, cctk_lsh[0]-1, cctk_lsh[1]-1, cctk_lsh[2]-1));
    memset(gyz, 0, sizeof(*gxy)*CCTK_GFINDEX3D(cctkGH, cctk_lsh[0]-1, cctk_lsh[1]-1, cctk_lsh[2]-1));

#pragma omp parallel for
    for (int k = 0; k < cctk_lsh[2]; k++)
        for (int j = 0; j < cctk_lsh[1]; j++)
            for (int i = 0; i < cctk_lsh[0]; i++) {
                int   index = CCTK_GFINDEX3D(cctkGH, i, j, k);
                CCTK_REAL r = ISO_R(index);
                CCTK_REAL R = gsl_spline_eval(spline, r, acc);
                CCTK_REAL psi2 = R/r, psi4 = psi2*psi2;
                CCTK_REAL k_fact = TRUMPET_CONST/(r*r*r*r*R);

                gxx[index] = gyy[index] = gzz[index] = psi4;

                kxx[index] = -k_fact*(3*SQR(x[index]) - SQR(r));
                kyy[index] = -k_fact*(3*SQR(y[index]) - SQR(r));
                kzz[index] = -k_fact*(3*SQR(z[index]) - SQR(r));

                kxy[index] = -k_fact*3*x[index]*y[index];
                kxz[index] = -k_fact*3*x[index]*z[index];
                kyz[index] = -k_fact*3*y[index]*z[index];
            }
}

void trumpet_lapse(CCTK_ARGUMENTS)
{
    DECLARE_CCTK_ARGUMENTS;
    DECLARE_CCTK_PARAMETERS;

    gsl_interp_accel *acc;
    gsl_spline    *spline;

    double *r_iso, *r_areal;
    CCTK_INT size;

    get_r_tables(&r_iso, &r_areal, &size, cctk_dim);

    spline = gsl_spline_alloc(gsl_interp_cspline, size);
    gsl_spline_init(spline, r_iso, r_areal, size);
    acc = gsl_interp_accel_alloc();

#pragma omp parallel for
    for (int k = 0; k < cctk_lsh[2]; k++)
        for (int j = 0; j < cctk_lsh[1]; j++)
            for (int i = 0; i < cctk_lsh[0]; i++) {
                int   index = CCTK_GFINDEX3D(cctkGH, i, j, k);
                CCTK_REAL r = sqrt(SQR(x[index]) + SQR(y[index]) + SQR(z[index]) + EPS);
                CCTK_REAL R = gsl_spline_eval(spline, r, acc);

                alp[index] = sqrt(1 - 2*MASS/R + TRUMPET_CONST*TRUMPET_CONST/(R*R*R*R));
            }
}

void trumpet_shift(CCTK_ARGUMENTS)
{
    DECLARE_CCTK_ARGUMENTS;
    DECLARE_CCTK_PARAMETERS;

    gsl_interp_accel *acc;
    gsl_spline    *spline;

    double *r_iso, *r_areal;
    CCTK_INT size;

    get_r_tables(&r_iso, &r_areal, &size, cctk_dim);

    spline = gsl_spline_alloc(gsl_interp_cspline, size);
    gsl_spline_init(spline, r_iso, r_areal, size);
    acc = gsl_interp_accel_alloc();

#pragma omp parallel for
    for (int k = 0; k < cctk_lsh[2]; k++)
        for (int j = 0; j < cctk_lsh[1]; j++)
            for (int i = 0; i < cctk_lsh[0]; i++) {
                int   index = CCTK_GFINDEX3D(cctkGH, i, j, k);
                CCTK_REAL r = sqrt(SQR(x[index]) + SQR(y[index]) + SQR(z[index]) + EPS);
                CCTK_REAL R = gsl_spline_eval(spline, r, acc);

                betax[index] = TRUMPET_CONST*x[index]/(R*R*R);
                betay[index] = TRUMPET_CONST*y[index]/(R*R*R);
                betaz[index] = TRUMPET_CONST*z[index]/(R*R*R);
            }
}

            }
}