<commit_before>#include "gsl/gsl_sf_legendre.h"
#include "gsl/gsl_sf_gegenbauer.h"
#include "gsl/gsl_sf_gamma.h"
#include <math.h>
#include "coeff_helper.h"
#include "bfe_helper.h"
#include <complex.h>

#define SQRT_FOURPI 3.544907701811031

double cc_phi_nlm(double s, double phi, double X, int n, int l, int m) {
    /* Complex conjugate of phi_nlm */
    double tmp = pow(-1.,m) * gsl_sf_fact(l-m) / gsl_sf_fact(l+m);
    return phi_nl(s, n, l) * tmp * gsl_sf_legendre_sphPlm(l, m, X);
}

double STnlm_integrand_help(double phi, double X, double xsi,
                            double density, int n, int l, int m) {
    /*
    Computes the integrand used to compute the expansion
    coefficients, Snlm, Tnlm. The integral is done over:

        * phi: azimuthal angle
        * X: cos(theta), where theta is the colatitude
            (e.g., from spherical coordinates typical to physicists)
        * xsi: (s-1)/(s+1), a radial coordinate mapped to the interval
            [-1,1] rather than [0,inf].
    */
    double _tmp = (1 - xsi);
    double s = (1 + xsi) / _tmp;

    // temporary variables
    double Knl, Anl, krond, _tmp2;

    Knl = 0.5*n*(n + 4*l + 3) + (l + 1)*(2*l + 1);
    if (m == 0) {
        krond = 1.;
    } else {
        krond = 0.;
    }

    _tmp2 = (gsl_sf_gamma(n + 4*l + 3) / (gsl_sf_fact(n) * (n + 2*l + 1.5) * pow(gsl_sf_gamma(2*l + 1.5),2)));
    Anl = 1/(Knl / pow(2., 8*l+6) * _tmp2 * (1 + krond) * M_PI * 2./(2*l+1.) * gsl_sf_fact(l+m) / gsl_sf_fact(l-m));
    // Anl = (-pow(2., 8*l+6) / (Knl * 4*M_PI) *
    //        (gsl_sf_fact(n) * (n + 2*l + 1.5) * pow(gsl_sf_gamma(2*l + 1.5),2)) / gsl_sf_gamma(n+4*l+3));
    return -2. / (_tmp*_tmp) * s*s * cc_phi_nlm(s, phi, X, n, l, m) * Anl * density;
}

extern double c_Snlm_integrand(double phi, double X, double xsi,
                               double density, int n, int l, int m) {
    return STnlm_integrand_help(phi, X, xsi, density, n, l, m) * cos(m*phi);
}

extern double c_Tnlm_integrand(double phi, double X, double xsi,
                               double density, int n, int l, int m) {
    return STnlm_integrand_help(phi, X, xsi, density, n, l, m) * sin(m*phi);
}

extern void c_STnlm_discrete(double *s, double *phi, double *X, double *m_k, int K,
                             int n, int l, int m, double *ST) {
    double _tmp, krond;
    double Knl = 0.5*n*(n + 4*l + 3) + (l + 1)*(2*l + 1);
    double _tmp2 = (gsl_sf_gamma(n + 4*l + 3) / (gsl_sf_fact(n) * (n + 2*l + 1.5) * pow(gsl_sf_gamma(2*l + 1.5),2)));

    if (m == 0) {
        krond = 1.;
    } else {
        krond = 0.;
    }
    double Anl = (2.-krond)/(Knl / pow(2., 8*l+6) * _tmp2 * (1 + krond) * M_PI * 2./(2*l+1.) * gsl_sf_fact(l+m) / gsl_sf_fact(l-m));

    // zero out coeff storage array
    ST[0] = 0.;
    ST[1] = 0.;
    for (int k=0; k<K; k++) {
        _tmp = Anl * m_k[k] * phi_nlm(s[k], phi[k], X[k], n, l, m);
        ST[0] += _tmp * cos(m*phi[k]); // Snlm
        ST[1] += _tmp * sin(m*phi[k]); // Tnlm
    }
}
<commit_msg>fix coefficient comp for discrete points<commit_after>#include "gsl/gsl_sf_legendre.h"
#include "gsl/gsl_sf_gegenbauer.h"
#include "gsl/gsl_sf_gamma.h"
#include <math.h>
#include "coeff_helper.h"
#include "bfe_helper.h"
#include <complex.h>

#define SQRT_FOURPI 3.544907701811031

double cc_phi_nlm(double s, double phi, double X, int n, int l, int m) {
    /* Complex conjugate of phi_nlm */
    double tmp = pow(-1.,m) * gsl_sf_fact(l-m) / gsl_sf_fact(l+m);
    return phi_nl(s, n, l) * tmp * gsl_sf_legendre_sphPlm(l, m, X);
}

double STnlm_integrand_help(double phi, double X, double xsi,
                            double density, int n, int l, int m) {
    /*
    Computes the integrand used to compute the expansion
    coefficients, Snlm, Tnlm. The integral is done over:

        * phi: azimuthal angle
        * X: cos(theta), where theta is the colatitude
            (e.g., from spherical coordinates typical to physicists)
        * xsi: (s-1)/(s+1), a radial coordinate mapped to the interval
            [-1,1] rather than [0,inf].
    */
    double _tmp = (1 - xsi);
    double s = (1 + xsi) / _tmp;

    // temporary variables
    double Knl, Anl, krond, _tmp2;

    Knl = 0.5*n*(n + 4*l + 3) + (l + 1)*(2*l + 1);
    if (m == 0) {
        krond = 1.;
    } else {
        krond = 0.;
    }

    _tmp2 = (gsl_sf_gamma(n + 4*l + 3) / (gsl_sf_fact(n) * (n + 2*l + 1.5) * pow(gsl_sf_gamma(2*l + 1.5),2)));
    Anl = 1/(Knl / pow(2., 8*l+6) * _tmp2 * (1 + krond) * M_PI * 2./(2*l+1.) * gsl_sf_fact(l+m) / gsl_sf_fact(l-m));
    // Anl = (-pow(2., 8*l+6) / (Knl * 4*M_PI) *
    //        (gsl_sf_fact(n) * (n + 2*l + 1.5) * pow(gsl_sf_gamma(2*l + 1.5),2)) / gsl_sf_gamma(n+4*l+3));
    return -2. / (_tmp*_tmp) * s*s * cc_phi_nlm(s, phi, X, n, l, m) * Anl * density;
}

extern double c_Snlm_integrand(double phi, double X, double xsi,
                               double density, int n, int l, int m) {
    return STnlm_integrand_help(phi, X, xsi, density, n, l, m) * cos(m*phi);
}

extern double c_Tnlm_integrand(double phi, double X, double xsi,
                               double density, int n, int l, int m) {
    return STnlm_integrand_help(phi, X, xsi, density, n, l, m) * sin(m*phi);
}

extern void c_STnlm_discrete(double *s, double *phi, double *X, double *m_k, int K,
                             int n, int l, int m, double *ST) {
    // Taken from Eq. 16 of Lowing et al.
    double Knl = 0.5*n*(n + 4*l + 3) + (l + 1)*(2*l + 1);

    double krond, _tmp;
    if (m == 0) {
        krond = 1.;
    } else {
        krond = 0.;
    }

    double Anl = -pow(2., 8*l+6) / (4*M_PI*Knl) * \
                 (gsl_sf_fact(n) * (n + 2*l + 1.5) * pow(gsl_sf_gamma(2*l + 1.5),2)) / \
                 gsl_sf_gamma(n + 4*l + 3);
    double coeff = Anl * (2.-krond);

    // zero out coeff storage array
    ST[0] = 0.;
    ST[1] = 0.;
    for (int k=0; k<K; k++) {
        _tmp = coeff * m_k[k] * phi_nlm(s[k], phi[k], X[k], n, l, m);
        ST[0] += _tmp * cos(m*phi[k]); // Snlm
        ST[1] += _tmp * sin(m*phi[k]); // Tnlm
    }
}

{'fraction_non_alphanumeric': 0.1742711250205897, 'fraction_numerical': 0.03376708944160764, 'mean_word_length': 2.2557640750670243, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6851214', 'n_tokens_mistral': 2742, 'n_tokens_neox': 2585, 'n_words': 662}