/*
 * Copyright 2014-2015 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
 * definitions of the q functions in the exponential in the Brill data
 */

#include <math.h>

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

static double q_gundlach(const BDContext *bd, double rho, double z)
{
    return bd->amplitude * SQR(rho) * exp(- (SQR(rho) + SQR(z)));
}

static double dq_rho_gundlach(const BDContext *bd, double rho, double z)
{
    return bd->amplitude * 2 * rho * exp(-SQR(rho) - SQR(z)) * (1 - SQR(rho));
}

static double d2q_rho_gundlach(const BDContext *bd, double rho, double z)
{
    double rho2 = SQR(rho);
    return bd->amplitude * 2 * exp(-rho2 - SQR(z)) * ((1 - rho2) * (1 - 2 * rho2) - 2 * rho2);
}

static double d2q_rho_z_gundlach(const BDContext *bd, double rho, double z)
{
    return -bd->amplitude * 4 * z * rho * exp(-SQR(rho) - SQR(z)) * (1 - SQR(rho));
}

static double dq_z_gundlach(const BDContext *bd, double rho, double z)
{
    return - bd->amplitude * 2 * z * SQR(rho) * exp(-SQR(rho) - SQR(z));
}

static double d2q_z_gundlach(const BDContext *bd, double rho, double z)
{
    return bd->amplitude * 2 * SQR(rho) * exp(-SQR(rho) - SQR(z)) * (2 * SQR(z) - 1);
}

// q function from from PHYSICAL REVIEW D 88, 103009 (2013)
// with σ and ρ_0 hardcoded to 1 for now
const QFunc bdi_q_func_gundlach = {
    .q         = q_gundlach,
    .dq_rho    = dq_rho_gundlach,
    .dq_z      = dq_z_gundlach,
    .d2q_rho   = d2q_rho_gundlach,
    .d2q_z     = d2q_z_gundlach,
    .d2q_rho_z = d2q_rho_z_gundlach,
};

static double q_eppley(const BDContext *bd, double rho, double z)
{
    double r2 = SQR(rho) + SQR(z);
    return bd->amplitude * SQR(rho) / (1.0 + pow(r2, bd->eppley_n / 2.0));
}

static double dq_rho_eppley(const BDContext *bd, double rho, double z)
{
    double A = bd->amplitude;
    double n = bd->eppley_n;
    double rho2 = SQR(rho);
    double z2   = SQR(z);

    double r2 = rho2 + z2;
    double r  = sqrt(r2);
    double rnm2 = pow(r, n - 2);

    return A * rho * (2 * (1 + rnm2 * r2) - n * rho2 * rnm2) / SQR(1 + rnm2 * r2);
}

static double dq_z_eppley(const BDContext *bd, double rho, double z)
{
    double A = bd->amplitude;
    double n = bd->eppley_n;
    double rho2 = SQR(rho);
    double z2   = SQR(z);

    double r2 = rho2 + z2;
    double r  = sqrt(r2);
    double rnm2 = pow(r, n - 2);

    return - A * n * rho2 * z * rnm2 / SQR(1 + rnm2 * r2);
}

static double d2q_rho_eppley(const BDContext *bd, double rho, double z)
{
    double A = bd->amplitude;
    double n = bd->eppley_n;
    double rho2 = SQR(rho);
    double z2   = SQR(z);

    double r2 = rho2 + z2;
    double r  = sqrt(r2);
    double rnm4 = pow(r, n - 4);
    double rn   = rnm4 * SQR(r2);
    double rnp1 = 1 + rn;

    return A * (SQR(n) * SQR(rho2) * rnm4 * (rn - 1) - n * rho2 * rnm4 * (3 * rho2 + 5 * z2) * rnp1 + 2 * SQR(rnp1)) / (rnp1 * SQR(rnp1));
}

static double d2q_z_eppley(const BDContext *bd, double rho, double z)
{
    double A = bd->amplitude;
    double n = bd->eppley_n;
    double rho2 = SQR(rho);
    double z2   = SQR(z);

    double r2 = rho2 + z2;
    double r  = sqrt(r2);
    double rnm4 = pow(r, n - 4);
    double rn   = rnm4 * SQR(r2);
    double rnp1 = 1 + rn;

    return A * n * rho2 * rnm4 * (z2 - rho2 - n * z2 + rn * ((1 + n) * z2 - rho2)) / (rnp1 * SQR(rnp1));
}

static double d2q_rho_z_eppley(const BDContext *bd, double rho, double z)
{
    double A = bd->amplitude;
    double n = bd->eppley_n;
    double rho2 = SQR(rho);
    double z2   = SQR(z);

    double r2 = rho2 + z2;
    double r  = sqrt(r2);
    double rnm4 = pow(r, n - 4);
    double rn   = rnm4 * SQR(r2);
    double rnp1 = 1 + rn;

    return A * n * rho * z * rnm4 * (rn * (n * rho2 - 2 * z2) - 2 * z2 - n * rho2) / (rnp1 * SQR(rnp1));
}

// q function from from PHYSICAL REVIEW D 16, 1609 (1977)
// with λ hardcoded to 1 for now
const QFunc bdi_q_func_eppley = {
    .q         = q_eppley,
    .dq_rho    = dq_rho_eppley,
    .dq_z      = dq_z_eppley,
    .d2q_rho   = d2q_rho_eppley,
    .d2q_z     = d2q_z_eppley,
    .d2q_rho_z = d2q_rho_z_eppley,
};