#include <math.h>

#include "maximal_slicing_axi.h"

static double sb_even_eval(double coord, int idx)
{
    double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));

    idx *= 2;   // even only

    return sin((idx + 1) * val);
}

static double sb_even_eval_diff1(double coord, int idx)
{
    double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));

    idx *= 2;   // even only

    return - SCALE_FACTOR * (idx + 1) * SGN(coord) * cos((idx + 1) * val) / (SQR(SCALE_FACTOR) + SQR(coord));
}

static double sb_even_eval_diff2(double coord, int idx)
{
    double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));

    idx *= 2;   // even only

    return SCALE_FACTOR * (idx + 1) * SGN(coord) * (2 * fabs(coord) * cos((idx + 1) * val) - SCALE_FACTOR * (idx + 1) * sin((idx + 1) * val)) / SQR(SQR(SCALE_FACTOR) + SQR(coord));
}

static double sb_even_colloc_point(int order, int idx)
{
    double t;

    idx = order - idx - 1;
    //order *= 2;

    //t = (idx + 2) * M_PI / (order + 4);
    t = (idx + 2) * M_PI / (2 * order + 2);
    return SCALE_FACTOR / tan(t);
}

const BasisSet msa_sb_even_basis = {
    .eval         = sb_even_eval,
    .eval_diff1   = sb_even_eval_diff1,
    .eval_diff2   = sb_even_eval_diff2,
    .colloc_point = sb_even_colloc_point,
};

static double tb_even_eval(double coord, int idx)
{
    double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));

    idx++;
    idx *= 2;   // even only

    return cos(idx * val) - 1.0;
}

static double tb_even_eval_diff1(double coord, int idx)
{
    double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));

    idx++;
    idx *= 2;   // even only

    return SCALE_FACTOR * idx * SGN(coord) * sin(idx * val) / (SQR(SCALE_FACTOR) + SQR(coord));
}

static double tb_even_eval_diff2(double coord, int idx)
{
    double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));

    idx++;
    idx *= 2;   // even only

    return -SCALE_FACTOR * idx * SGN(coord) * (2 * fabs(coord) * sin(idx * val) + SCALE_FACTOR * idx * cos(idx * val)) / SQR(SQR(SCALE_FACTOR) + SQR(coord));
}

static double tb_even_colloc_point(int order, int idx)
{
    double t;

    idx = order - idx - 1;
    //order *= 2;

    //t = (idx + 2) * M_PI / (order + 4);
    t = (idx + 3) * M_PI / (2 * (order + 2));
    return SCALE_FACTOR / tan(t);
}

const BasisSet msa_tb_even_basis = {
    .eval         = tb_even_eval,
    .eval_diff1   = tb_even_eval_diff1,
    .eval_diff2   = tb_even_eval_diff2,
    .colloc_point = tb_even_colloc_point,
};

static double full_eval(double coord, int idx)
{
    double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
    int flag = idx & 1;

    idx /= 2;

    if (flag) return sin((2 * idx + 1) * 4 * val);
    else      return cos((2 * idx + 2) * 4 * val) - 1;
}

static double full_eval_diff1(double coord, int idx)
{
    double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
    int flag = idx & 1;

    idx /= 2;

    if (flag) {
        idx = 2 * idx + 1;
        return -4 * idx * SCALE_FACTOR * cos(idx * 4 * val) / (SQR(SCALE_FACTOR) + SQR(coord));
    } else {
        idx = 2 * (idx + 1);
        return 4 * idx * SCALE_FACTOR * sin(idx * 4 * val) / (SQR(SCALE_FACTOR) + SQR(coord));
    }
}

static double full_eval_diff2(double coord, int idx)
{
    double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
    int flag = idx & 1;

    idx /= 2;

    if (flag) {
        idx = 2 * idx + 1;
        return (16 * SQR(idx * SCALE_FACTOR) * cos(idx * 4 * val) - 4 * idx * SCALE_FACTOR * sin(idx * 4 * val) * 2 * coord) / SQR(SQR(SCALE_FACTOR) + SQR(coord));
    } else {
        idx = 2 * (idx + 1);
        return (-16 * SQR(idx * SCALE_FACTOR) * sin(idx * 4 * val) - 4 * idx * SCALE_FACTOR * cos(idx * 4 * val) * 2 * coord) / SQR(SQR(SCALE_FACTOR) + SQR(coord));
    }
}

static double full_colloc_point(int order, int idx)
{
    double t;

    idx = order - idx - 1;

    t = (idx + 0.5) * M_PI / (2 * order);

    return SCALE_FACTOR / tan(t);

}

const BasisSet msa_full_basis = {
    .eval         = full_eval,
    .eval_diff1   = full_eval_diff1,
    .eval_diff2   = full_eval_diff2,
    .colloc_point = full_colloc_point,
};

static double cheb_eval(double coord, int idx)
{
    return cos(2 * idx * acos(coord / SCALE_FACTOR));
}

static double cheb_eval_diff1(double coord, int idx)
{
    return 2 * idx * sin(2 * idx * acos(coord / SCALE_FACTOR)) / sqrt(SQR(SCALE_FACTOR) - SQR(coord));
}

static double cheb_eval_diff2(double coord, int idx)
{
    double t = acos(coord / SCALE_FACTOR);
    return 2 * idx * (cos(2 * idx * t) * 2 * idx / (SQR(SCALE_FACTOR) - SQR(coord)) + sin(2 * idx * t) * coord / pow(SQR(SCALE_FACTOR) - SQR(coord), 1.5));
}

static double cheb_colloc_point(int order, int idx)
{
    return SCALE_FACTOR * cos((idx + 0.01) * M_PI / (4 * order + 4));
}

const BasisSet msa_cheb_basis = {
    .eval       = cheb_eval,
    .eval_diff1 = cheb_eval_diff1,
    .eval_diff2 = cheb_eval_diff2,
    .colloc_point = cheb_colloc_point,
};