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

#ifndef TEUKOLSKY_DATA_H
#define TEUKOLSKY_DATA_H

#include <stdarg.h>
#include <stddef.h>
#include <stdint.h>

/**
 * API usage:
 *
 * First, allocate the solver context with td_context_alloc(). All interaction
 * with the solver is done through this context.
 *
 * Fill any fields in the context that are described as settable by the caller.
 * Call td_solve() to solve the equation determined by the option values.
 *
 * The metric/extrinsic curvature can be evaluated with td_eval_metric()/td_eval_curv().
 *
 * Finally, free the solver context with td_context_free().
 */

/**
 * Identifies the specific initial data family to use.
 */
enum TDFamily {
    /**
     * Time-antisymmetric initial data with a cubic radial term multiplied by an
     * exponential and a quadrupole angular term.
     * Similar to that used in Abrahams&Evans PhysRevD v49,n8 (1994).
     * Conformally flat spatial metric.
     *  r       /   x  3      x    \        /    x  2  \
     * K  =  a |  (---)   -  ---    | * exp| - (---)    | sin(2θ)
     *  θ       \   L         L    /        \    L     /
     */
    TD_FAMILY_TIME_ANTISYM_CUBIC = 0,
    /**
     * Time-antisymmetric initial data with a linear radial term multiplied by
     * an exponential and a quadrupole angular term.
     * Conformally flat spatial metric.
     *  r       x       /    x  2  \
     * K  = a  ---  exp| - (---)    | sin(2θ)
     *  θ       L       \    L     /
     */
    TD_FAMILY_TIME_ANTISYM_LINEAR,
};

typedef struct TDContext {
    /**
     *  Solver internals, not to be accessed by the caller
     */
    void *priv;

    /********************************
     *  options, set by the caller  *
     ********************************/

    /**
     * A callback that will be used to print diagnostic messages.
     *
     * Defaults to fprintf(stderr, ...)
     */
    void (*log_callback)(const struct TDContext *td, int level,
                         const char *fmt, va_list);

    /**
     * Arbitrary user data, e.g. to be used by the log callback.
     */
    void *opaque;

    /********************************
     *  initial data parameters     *
     ********************************/

    /**
     * The amplitude parameter for the seed function.
     * Defaults to 1.
     */
    double amplitude;

    /********************************
     *        solver options        *
     ********************************/

    /**
     * The number of basis functions in each direction.
     * [0] - radial, [1] - angular
     */
    unsigned int nb_coeffs[2];

    /**
     * The scaling factor used in the basis functions.
     * [0] - radial, [1] - angular
     */
    double basis_scale_factor[2];

    /**
     * Maximum number of Newton iterations.
     */
    unsigned int max_iter;

    /**
     * Absolute tolerance. The solver is deemed to have converged
     * after maximum difference between iterations is below this bound.
     */
    double atol;

    /**
     * Number of threads to use for parallelization. Set to 0 to automatically
     * detect according to the number of CPU cores.
     */
    unsigned int nb_threads;

    /**
     * Spectral coefficients of the solution are exported here. Mainly useful
     * for debugging, Use td_eval_*() to evaluate the solution.
     *
     * [0]: ψ - 1
     * [1]: K^r_r
     * [2]: K^φ_φ
     */
    double *coeffs[3];

    /**
     * Set to 1 to produce the "upper branch" solution (see thesis for details).
     */
    unsigned int solution_branch;

    /**
     * The family to use, defaults to TD_FAMILY_TIME_ANTISYM_CUBIC.
     */
    enum TDFamily family;
} TDContext;

/**
 * Allocate and initialize the solver.
 */
TDContext *td_context_alloc(void);

/**
 * Free the solver and everything associated with it.
 */
void td_context_free(TDContext **td);

/**
 * Solve the constraint equations and export the expansion coefficients in the
 * context.
 *
 * @return >= 0 on success, a negative error code on failure
 */
int td_solve(TDContext *td, double *coeffs_init[3]);

/**
 * Evaluate the conformal factor ψ at the specified points in spherical
 * coordinates { r, θ, φ } (since the spacetime is axially symmetric, φ is
 * disregarded).
 *
 * @param td the solver context
 * @param nb_coords number of elements in the r and theta arrays
 * @param r values of r coordinates
 * @param theta values of θ coordinates
 * @param diff_order order of the derivatives of the metric to evaluate. The
 *   first element specifies the derivative wrt r, the second
 *   wrt θ. I.e. diff_order[i] = { 0, 0 } evaluates ψ itself,
 *   diff_order = { 0, 1 } evaluates ∂ψ/∂θ etc.
 * @param out a nb_coords-sized array into which the values of the ψ will be
 *   written. out[i] is the value of ψ at the spacetime point { r[i], θ[i], 0 }
 *
 * @return >= 0 on success, a negative error code on failure.
 */
int td_eval_psi(const TDContext *td,
                size_t nb_coords, const double *r, const double *theta,
                const unsigned int diff_order[2],
                double *out);
/**
 * Same as td_eval_psi(), except K_{r}^r ('rr' component of the mixed-index
 * extrinsic curvature tensor) is evaluated.
 */
int td_eval_krr(const TDContext *td,
                size_t nb_coords, const double *r, const double *theta,
                const unsigned int diff_order[2],
                double *out);
/**
 * Same as td_eval_psi(), except K_{φ}^φ ('φφ' component of the mixed-index
 * extrinsic curvature tensor) is evaluated.
 */
int td_eval_kpp(const TDContext *td,
                size_t nb_coords, const double *r, const double *theta,
                const unsigned int diff_order[2],
                double *out);
/**
 * Same as td_eval_psi(), except K_{r}^θ ('rθ' component of the mixed-index
 * extrinsic curvature tensor) is evaluated.
 */
int td_eval_krt(const TDContext *td,
                size_t nb_coords, const double *r, const double *theta,
                const unsigned int diff_order[2],
                double *out);
/**
 * Compute the values of lapse for maximal slicing (i.e. that make the time
 * derivative of K vanish).
 *
 * Parameters are the same as for td_eval_psi().
 */
int td_eval_lapse(const TDContext *td,
                  size_t nb_coords, const double *r, const double *theta,
                  const unsigned int diff_order[2],
                  double *out);

/**
 * Evaluate the residual of the three constraint equations at the specified
 * coordinates.
 */
int td_eval_residual(const TDContext *td,
                     size_t nb_coords_r,     const double *r,
                     size_t nb_coords_theta, const double *theta,
                     double *out[3]);

#endif /* TEUKOLSKY_DATA_H */