/*
 * Basis sets for pseudospectral methods
 * Copyright (C) 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 "basis.h"
#include "common.h"

/*
 * The basis of even (n = 2 * idx) SB functions (Boyd 2000, Ch 17.9)
 * SB(x, n) = sin((n + 1) arccot(|x| / L))
 * They are symmetric wrt origin and decay as 1/x in infinity.
 */
static double sb_even_eval(double coord, int idx)
{
    double val = atan2(SCALE_FACTOR, coord);

    idx *= 2;   // even only

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

static double sb_even_eval_diff1(double coord, int idx)
{
    double val = atan2(SCALE_FACTOR, coord);

    idx *= 2;   // even only

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

static double sb_even_eval_diff2(double coord, int idx)
{
    double val = atan2(SCALE_FACTOR, coord);

    idx *= 2;   // even only

    return SCALE_FACTOR * (idx + 1) * (2 * 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);
#if QMS_POLAR
    t = (idx + 2) * M_PI / (2 * order + 3);
#else
    t = (idx + 2) * M_PI / (2 * order + 2);
#endif
    return SCALE_FACTOR / tan(t);
}

const BasisSet qms_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 + 2) * M_PI / (2 * order + 4);
    return SCALE_FACTOR / tan(t);
}

const BasisSet qms_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 tl_eval(double coord, int idx)
{
    double t = 2 * atan2(sqrt(SCALE_FACTOR), sqrt(fabs(coord)));

    //idx++;

    return cos(idx * t);// - 1.0;
}

static double tl_eval_diff1(double coord, int idx)
{
    double y = fabs(coord);
    double y12 = sqrt(y);
    double t = 2 * atan2(sqrt(SCALE_FACTOR), y12);

    //idx++;

    if (y < 1e-10)
        return -2.0 * SQR(idx) * pow(-1.0, idx) / SCALE_FACTOR;

    return idx * sqrt(SCALE_FACTOR) * sin(idx * t) / (y12 * (y + SCALE_FACTOR));
}

static double tl_eval_diff2(double coord, int idx)
{
    double y = fabs(coord);
    double y12 = sqrt(y);
    double t = 2 * atan2(sqrt(SCALE_FACTOR), y12);

    //idx++;

    if (y < 1e-10)
        return 4.0 * SQR(idx) * (SQR(idx) + 2) * pow(-1.0, idx) / (3.0 * SQR(SCALE_FACTOR));

    return -(SCALE_FACTOR * SQR(idx) * cos(idx * t) / (y * SQR(y + SCALE_FACTOR)) + sqrt(SCALE_FACTOR) * idx * sin(idx * t) / (y12 * SQR(y + SCALE_FACTOR)) + sqrt(SCALE_FACTOR) * idx * sin(idx * t) / (2 * y * y12 * (y + SCALE_FACTOR)));
}

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

    if (!idx)
        return 0.0;

    idx = order - idx - 1;

    order++;
    idx++;

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

const BasisSet qms_tl_basis = {
    .eval         = tl_eval,
    .eval_diff1   = tl_eval_diff1,
    .eval_diff2   = tl_eval_diff2,
    .colloc_point = tl_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 qms_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)
{
    double x = MIN(1.0, MAX(2.0 * coord / SCALE_FACTOR - 1.0, -1.0));

    //idx += 1;
    return cos(idx * acos(x));
}

static double cheb_eval_diff1(double coord, int idx)
{
    double x = MIN(1.0, MAX(2.0 * coord / SCALE_FACTOR - 1.0, -1.0));
    double t = acos(x);

    //idx += 1;

    if (fabs(x - 1.0) < 1e-8)
        return (2.0 / SCALE_FACTOR) * SQR(idx);
    if (fabs(x + 1.0) < 1e-8)
        return -(2.0 / SCALE_FACTOR) * SQR(idx) * pow(-1.0, idx);
    return (2.0 / SCALE_FACTOR) * idx * sin(idx * t) / sin(t);
}

static double cheb_eval_diff2(double coord, int idx)
{
    double x = MIN(1.0, MAX(2.0 * coord / SCALE_FACTOR - 1.0, -1.0));
    double t = acos(x);
    double st = sin(t);

    //idx += 1;

    if (fabs(x - 1.0) < 1e-8)
        return SQR(2.0 / SCALE_FACTOR) * SQR(idx) * (SQR(idx) - 1.0) / 3.;
    if (fabs(x + 1.0) < 1e-8)
        return SQR(2.0 / SCALE_FACTOR) * pow(-1.0, idx) * SQR(idx) * (SQR(idx) - 1.0) / 3.;

    return SQR(2.0 / SCALE_FACTOR) * (-SQR(idx) * cos(idx * t) / SQR(st) + idx * cos(t) * sin(idx * t) / (st * SQR(st)));
}

static double cheb_colloc_point(int order, int idx)
{
    double y;

    idx = order - idx - 1;

    //order += 1;
    //idx += 1;
    y = cos(idx * M_PI / (order - 1));
    //y = cos((2 * idx + 1) * M_PI / (2 * order));

    return 0.5 * SCALE_FACTOR * (y + 1.0);
}

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

static double cheb_even_eval(double coord, int idx)
{
    double x = MIN(1.0, MAX(coord / SCALE_FACTOR, -1.0));

    idx *= 2;
    //idx += 1;
    return cos(idx * acos(x));
}

static double cheb_even_eval_diff1(double coord, int idx)
{
    double x = MIN(1.0, MAX(coord / SCALE_FACTOR, -1.0));
    double t = acos(x);

    idx *= 2;
    //idx += 1;

    if (fabs(fabs(x) - 1.0) < 1e-8)
        return (1.0 / SCALE_FACTOR) * SQR(idx);
    return (1.0 / SCALE_FACTOR) * idx * sin(idx * t) / sin(t);
}

static double cheb_even_eval_diff2(double coord, int idx)
{
    double x = MIN(1.0, MAX(coord / SCALE_FACTOR, -1.0));
    double t = acos(x);
    double st = sin(t);

    idx *= 2;
    //idx += 1;

    if (fabs(fabs(x) - 1) < 1e-8)
        return SQR(1.0 / SCALE_FACTOR) * SQR(idx) * (SQR(idx) - 1.0) / 3.;

    return SQR(1.0 / SCALE_FACTOR) * (-SQR(idx) * cos(idx * t) / SQR(st) + idx * cos(t) * sin(idx * t) / (st * SQR(st)));
}

static double cheb_even_colloc_point(int order, int idx)
{
    double y;

    idx = order - idx - 1;

    //order += 1;
    //idx += 1;
    order *= 2;
    y = cos(idx * M_PI / (order - 1));
    y = cos((2 * idx + 2) * M_PI / (2 * order));

    return SCALE_FACTOR * y;
}

const BasisSet qms_cheb_even_basis = {
    .eval         = cheb_even_eval,
    .eval_diff1   = cheb_even_eval_diff1,
    .eval_diff2   = cheb_even_eval_diff2,
    .colloc_point = cheb_even_colloc_point,
};

static double cos_even_eval(double coord, int idx)
{
    return cos(2 * idx * coord);
}

static double cos_even_eval_diff1(double coord, int idx)
{
    return -2 * idx * sin(2 * idx * coord);
}

static double cos_even_eval_diff2(double coord, int idx)
{
    return -4 * SQR(idx) * cos(2 * idx * coord);
}

static double cos_even_colloc_point(int order, int idx)
{
    return M_PI * idx / (2 * order - 0);
}

const BasisSet qms_cos_even_basis = {
    .eval         = cos_even_eval,
    .eval_diff1   = cos_even_eval_diff1,
    .eval_diff2   = cos_even_eval_diff2,
    .colloc_point = cos_even_colloc_point,
};