/*
 *  libgrpp - a library for the evaluation of integrals over
 *            generalized relativistic pseudopotentials.
 *
 *  Copyright (C) 2021-2023 Alexander Oleynichenko
 */

/*
 * This file contains subroutines for evaluation of angular integrals of the 1st and 2nd type.
 * It also contains functions for construction of matrices of the angular momentum operator
 * in the bases of either real or complex spherical harmonics.
 *
 * For more details on the angular integrals used in RPP integration, see:
 *
 * L. E. McMurchie, E. R. Davidson. Calculation of integrals over ab initio pseudopotentials.
 * J. Comput. Phys. 44(2), 289 (1981). doi: 10.1016/0021-9991(81)90053-x
 */

#include "angular_integrals.h"

#include <math.h>

#include "factorial.h"
#include "libgrpp.h"
#include "spherical_harmonics.h"


double integrate_unitary_sphere_polynomial(int i, int j, int k);


/**
 * Type 1 angular integral.
 * (see MrMurchie & Davidson, formula (28))
 */
double angular_type1_integral(int lambda, int II, int JJ, int KK, double *k)
{
    double sum = 0.0;

    for (int mu = -lambda; mu <= lambda; mu++) {
        double sum2 = 0.0;
        for (int r = 0; r <= lambda; r++) {
            for (int s = 0; s <= lambda; s++) {
                for (int t = 0; t <= lambda; t++) {
                    if (r + s + t == lambda) {
                        double y_lm_rst = spherical_to_cartesian_coef(lambda, mu, r, s, t);
                        double usp_int = integrate_unitary_sphere_polynomial(II + r, JJ + s, KK + t);
                        sum2 += y_lm_rst * usp_int;
                    }
                }
            }
        }
        sum += sum2 * evaluate_real_spherical_harmonic(lambda, mu, k);
    }

    return sum;
}


/**
 * Type 2 angular integral.
 * (see MrMurchie & Davidson, formula (29))
 */
double angular_type2_integral(const int lambda, const int L, const int m, const int a, const int b, const int c,
                              const double *rsh_values)
{
    double sum = 0.0;

    rsh_coef_table_t *rsh_coef_lambda = get_real_spherical_harmonic_table(lambda);
    rsh_coef_table_t *rsh_coef_L = get_real_spherical_harmonic_table(L);

    int ncomb_rst = rsh_coef_lambda->n_cart_comb;
    int ncomb_uvw = rsh_coef_L->n_cart_comb;

    for (int mu = -lambda; mu <= lambda; mu++) {
        double sum2 = 0.0;
        double rsh_value_k = rsh_values[mu + lambda];
        if (fabs(rsh_value_k) < LIBGRPP_ZERO_THRESH) {
            continue;
        }

        for (int icomb_uvw = 0; icomb_uvw < ncomb_uvw; icomb_uvw++) {

            int u = rsh_coef_L->cartesian_comb[3 * icomb_uvw];
            int v = rsh_coef_L->cartesian_comb[3 * icomb_uvw + 1];
            int w = rsh_coef_L->cartesian_comb[3 * icomb_uvw + 2];
            double y_lm_uvw = rsh_coef_L->coeffs[(m + L) * ncomb_uvw + icomb_uvw];
            if (fabs(y_lm_uvw) < LIBGRPP_ZERO_THRESH) {
                continue;
            }

            for (int icomb_rst = 0; icomb_rst < ncomb_rst; icomb_rst++) {

                int r = rsh_coef_lambda->cartesian_comb[3 * icomb_rst];
                int s = rsh_coef_lambda->cartesian_comb[3 * icomb_rst + 1];
                int t = rsh_coef_lambda->cartesian_comb[3 * icomb_rst + 2];
                double y_lam_mu_rst = rsh_coef_lambda->coeffs[(mu + lambda) * ncomb_rst + icomb_rst];
                if (fabs(y_lam_mu_rst) < LIBGRPP_ZERO_THRESH) {
                    continue;
                }

                double usp_int = integrate_unitary_sphere_polynomial(a + r + u, b + s + v, c + t + w);

                sum2 += y_lam_mu_rst * y_lm_uvw * usp_int;
            }
        }

        sum += sum2 * rsh_value_k;
    }

    return sum;
}


/**
 * Integral of the unitary sphere polynomial over full solid angle.
 * (see MrMurchie & Davidson, formula (30))
 */
inline double integrate_unitary_sphere_polynomial(int i, int j, int k)
{
    if ((i % 2 == 0) && (j % 2 == 0) && (k % 2 == 0)) {
        double dfac_i = (double) double_factorial(i - 1);
        double dfac_j = (double) double_factorial(j - 1);
        double dfac_k = (double) double_factorial(k - 1);
        double dfac_ijk = (double) double_factorial(i + j + k + 1);
        return 4 * M_PI * dfac_i * dfac_j * dfac_k / dfac_ijk;
    }
    else {
        return 0.0;
    }
}

