#include  "OrthPol.h"
// The implementation of the function OrthPol() are from the book "Computation of special functions" by Shanjie Zhang and Jianmming Jin, printing 1996
// chapter 2, page 12


/**
 * Purpose: Compute orthogonal polynomials: Chebyshev, Laguerre, Hermite polynomials, and their derivatives
 * 
 * @param KF
 *      KF = 1 for Chebyshev polynomial of the first kind
 *      KF = 2 for Chebyshev polynomial of the second kind
 *      KF = 3 for Laguerre polynomial
 *      KF = 4 for Hermite polynomial
 * @param n - Integer order of orthogonal polynomial.
 * @param x - Argument of orthogonal polynomial.
 * @return Vector with values of the orthogonal polynomial (the first place) of order n with argument x, and its first derivatives (the second place).
 * @throws None
 */
std::vector<double> OrthPol(const int& KF, const int& n, const double& x)
{
    std::vector<double> result;

    double* PL = new double[n + 1];
    double* DPL = new double[n + 1];
    
    double A = 2.0;
    double B = 0.0;
    double C = 1.0;
    double Y0 = 1.0;
    double Y1 = 2.0 * x;
    double DY0 = 0.0;
    double DY1 = 2.0;
    PL[0] = 1.0;
    PL[1] = 2.0 * x;
    DPL[0] = 0.0;
    DPL[1] = 2.0;

    if (KF == 1)
    {
        Y1 = x;
        DY1 = 1.0;
        PL[1] = x;
        DPL[1] = 1.0;
    }
    else if (KF == 3)
    {
        Y1 = 1.0 - x;
        DY1 = -1.0;
        PL[1] = 1.0 - x;
        DPL[1] = -1.0;
    }
    
    for (int k = 2; k <= n; k++)
    {
        if (KF == 3)
        {   
            A = -1.0 / k;
            B = 2.0 + A;
            C = 1.0 + A;
        }
        else if (KF == 4)
        {
            C = 2.0 * (k - 1.0);
        }
        double YN = (A * x + B) * Y1 - C * Y0;
        double DYN = A * Y1 + (A * x + B) * DY1 - C * DY0;
        PL[k] = YN;
        DPL[k] = DYN;
        Y0 = Y1;
        Y1 = YN;
        DY0 = DY1;
        DY1 = DYN;        
    }
    
    result.push_back(PL[n]);
    result.push_back(DPL[n]);

    return result;
}