#include <iostream>
#include <complex>
#include <iomanip>

#include <eigen3/Eigen/Dense>
#include <gsl/gsl_sf_legendre.h>

int main(int argc, char* argv[])
{
    std::cout << std::fixed << std::setprecision(16);

    const double pi                                     = 3.141592653589793238;

    std::complex<double> complex_j                      = std::complex<double>(0.0, 1.0);

    int I                                               = 17;
    int J                                               = 9;
    int I_new                                           = 23;
    int J_new                                           = 12;
    int K                                               = 11;

    Eigen::MatrixXcd grid                               = Eigen::MatrixXcd(J, I);
    Eigen::MatrixXcd fourier_coefficients               = Eigen::MatrixXcd(J,       (2*K) + 1);
    Eigen::MatrixXcd spherical_coefficients             = Eigen::MatrixXcd((2*K)+1,    K  + 1);
    Eigen::MatrixXcd filtered_fourier_coefficients      = Eigen::MatrixXcd(J_new,   (2*K) + 1);
    Eigen::MatrixXcd filterd_grid                       = Eigen::MatrixXcd(J_new,    I_new);
    Eigen::MatrixXcd computed_grid                      = Eigen::MatrixXcd(J_new, I_new);

    Eigen::VectorXd gauss_nodes                         = Eigen::VectorXd(J);
    gauss_nodes(0)                                      = -0.9681602395076261;
    gauss_nodes(1)                                      = -0.8360311073266358;
    gauss_nodes(2)                                      = -0.6133714327005904;
    gauss_nodes(3)                                      = -0.3242534234038089;
    gauss_nodes(4)                                      =  0.0000000000000000;
    gauss_nodes(5)                                      =  0.3242534234038089;
    gauss_nodes(6)                                      =  0.6133714327005904;
    gauss_nodes(7)                                      =  0.8360311073266358;
    gauss_nodes(8)                                      =  0.9681602395076261;

    Eigen::VectorXd gauss_weights                       = Eigen::VectorXd(J);
    gauss_weights(0)                                    = 0.0812743883615744;
    gauss_weights(1)                                    = 0.1806481606948574;
    gauss_weights(2)                                    = 0.2606106964029354;
    gauss_weights(3)                                    = 0.3123470770400029;
    gauss_weights(4)                                    = 0.3302393550012598;
    gauss_weights(5)                                    = 0.3123470770400029;
    gauss_weights(6)                                    = 0.2606106964029354;
    gauss_weights(7)                                    = 0.1806481606948574;
    gauss_weights(8)                                    = 0.0812743883615744;

    Eigen::VectorXd gauss_nodes_new                     = Eigen::VectorXd(J_new);
    gauss_nodes_new(0)                                  = -0.9815606342467192;
    gauss_nodes_new(1)                                  = -0.9041172563704749;
    gauss_nodes_new(2)                                  = -0.7699026741943047;
    gauss_nodes_new(3)                                  = -0.5873179542866175;
    gauss_nodes_new(4)                                  = -0.3678314989981802;
    gauss_nodes_new(5)                                  = -0.1252334085114689;
    gauss_nodes_new(6)                                  =  0.1252334085114689;
    gauss_nodes_new(7)                                  =  0.3678314989981802;
    gauss_nodes_new(8)                                  =  0.5873179542866175;
    gauss_nodes_new(9)                                  =  0.7699026741943047;
    gauss_nodes_new(10)                                 =  0.9041172563704749;
    gauss_nodes_new(11)                                 =  0.9815606342467192;

    for (int j=0; j<J_new; ++j)
    {
        double theta                                    = acos(gauss_nodes_new(j));

        for (int i = 0; i<I_new; ++i)
        {
            double phi                                  = (2.0*pi*i) / I_new;
            computed_grid(j, i)                         = cos(3.0 * theta) * sin(2.0 * phi);
        }
    }

    for (int j=0; j<J; ++j)
    {
        double theta                                    = acos(gauss_nodes(j));

        for (int i = 0; i<I; ++i)
        {
            double phi                                  = (2.0*pi*i) / I;
            grid(j, i)                                  = cos(3.0 * theta) * sin(2.0 * phi);
        }
    }

    for (int j=0; j<J; j++)
    {
        for (int m=-K; m<=K; ++m)
        {
            for (int i=0; i<I; i++)
            {
                double phi                              = (2.0*pi*i) / I;
                std::complex<double> exponent           = ( cos(-m*phi) + (complex_j * sin(-m*phi)) );
                fourier_coefficients(j, m+K)            += ( grid(j, i) *  exponent );
            }
            fourier_coefficients(j, m+K)                *= ( sqrt(2.0 * pi) / static_cast<double>(I) );
        }
    }

    for (int m=-K; m<=K; ++m)
    {
        for (int n=abs(m); n<=K; ++n)
        {
            for (int j=0; j<J; j++)
            {
                double theta                            = acos(gauss_nodes(j));
                double w_j                              = gauss_weights(j);
                double na_legendre                      = gsl_sf_legendre_sphPlm(n, abs(m), cos(theta));

                spherical_coefficients(m+K, n)          += (fourier_coefficients(j, m+K) *
                                                            na_legendre                  *
                                                            w_j);
            }
        }
    }

    for (int j=0; j<J_new; j++)
    {
        double theta                                    = acos(gauss_nodes_new(j));

        for (int m=-K; m<=K; ++m)
        {
            for (int n=abs(m); n<=K; ++n)
            {
                double na_legendre                      = gsl_sf_legendre_sphPlm(n, abs(m), cos(theta));

                filtered_fourier_coefficients(j, m+K)   += (spherical_coefficients(m+K, n) *
                                                            na_legendre);
            }
        }
    }

    for (int j=0; j<J_new; j++)
    {
        for (int i=0; i<I_new; i++)
        {
            for (int m=-K; m<=K; ++m)
            {
                double phi                              = (2.0 * pi * i) / I_new;

                std::complex<double> exponent           = ( cos(m*phi) + (complex_j * sin(m*phi)) );

                filterd_grid(j, i)                      += ( sqrt(2.0 * pi) * filtered_fourier_coefficients(j, m+K) * exponent);
            }
        }
    }

    for (int j=0; j<J_new; ++j)
    {
        for (int i=0; i<I_new; ++i)
        {
            std::cout <<filterd_grid(j, i).real()<<"\t"<<filterd_grid(j, i).imag()<<std::endl;
        }
        std::cout << std::endl;
    }

    for (int j=0; j<J_new; ++j)
    {
        for (int i = 0; i<I_new; ++i)
        {
            std::cout <<computed_grid(j, i).real()<<"\t"<<computed_grid(j, i).imag()<<std::endl;
        }
        std::cout << std::endl;
    }
}