//Renata Kopecna

#include "integrals.hh"

#include <Math/Vector4D.h>
#include <Math/Boost.h>
#include <Math/GSLIntegrator.h>

#include <spdlog.h>
#include <funcs.hh>

#include <TMath.h>

double integral_x_to_n(double x, int n){
    return pow(x,n+1)/(n+1.0);
}
double integrate_x_to_n(double a, double b, int n){
    return integral_x_to_n(b,n) - integral_x_to_n(a,n);
}

//this methods are helpful for the phi integration
//calculates int x^n * sin(x) dx
double integral_x_to_n_times_sin_x(double x, int n){
    assert(n >= 0);
    if (n == 0) return -cos(x);
    else        return -pow(x, n)*cos(x) + n*integral_x_to_n_times_cos_x(x, n-1);
}
double integrate_x_to_n_times_sin_x(double a, double b, int n){
    return integral_x_to_n_times_sin_x(b, n) - integral_x_to_n_times_sin_x(a, n);
}

//calculates int x^n * cos(x) dx (recursive using per partes)
double integral_x_to_n_times_cos_x(double x, int n){
    assert(n >= 0);
    if (n == 0) return sin(x);
    else        return pow(x,n)*sin(x) - n*integral_x_to_n_times_sin_x(x, n-1);
}
double integrate_x_to_n_times_cos_x(double a, double b, int n){
    return integral_x_to_n_times_cos_x(b, n)-integral_x_to_n_times_cos_x(a, n);
}

//calculates int x^n * sin(2x) dx (recursive using per partes)
double integral_x_to_n_times_sin_2x(double x, int n){
    assert(n >= 0);
    if (n == 0) return -0.5*cos(2.0*x);
    else        return -pow(x,n)*0.5*cos(2.0*x)
                +0.5*n*integral_x_to_n_times_cos_2x(x,n-1);
}
double integrate_x_to_n_times_sin_2x(double a, double b, int n){
    return integral_x_to_n_times_sin_2x(b, n) - integral_x_to_n_times_sin_2x(a, n);
}

//calculates int x^n * cos(2x) dx (recursive using per partes)
double integral_x_to_n_times_cos_2x(double x, int n){
    assert(n >= 0);
    if (n == 0) return 0.5*sin(2.0*x);
    else        return +0.5*pow(x,n)*sin(2.0*x)
                -0.5*n*integral_x_to_n_times_sin_2x(x,n-1);
}
double integrate_x_to_n_times_cos_2x(double a, double b, int n){
    return integral_x_to_n_times_cos_2x(b, n) - integral_x_to_n_times_cos_2x(a, n);
}

//calculates int x^n * cos(x)^2 dx
double integral_x_to_n_times_cos_x_2(double x, int n){
    assert(n >= 0);
    return +1.0/(1.0 + n)*pow(x,n+1)*cos(x)*cos(x)
            +1.0/(1.0+n)*integral_x_to_n_times_sin_2x(x, n+1);
}

//calculates int x^n * sin(x)^2 dx
double integral_x_to_n_times_sin_x_2(double x, int n){
    assert(n >= 0);
    return +1.0/(1.0 + n)*pow(x,n+1)*sin(x)*sin(x)
            -1.0/(1.0+n)*integral_x_to_n_times_sin_2x(x, n+1);
}

//calculates int x^n * asin(x) dx
double integral_x_to_n_times_asin_x(double x, int n){
    assert(n >= 0);
    if (n == 0){
        return x*asin(x)+sqrt(1-x*x);
    }
    else{
        //    return 1.0/(n+1.0)*pow(x,n)*(x*asin(x)+sqrt(1-x*x))
        //   -n*integral_x_to_n_times_sqrt_1_minus_x2(x, n-1);
        //factor 1/(n+1) missing???
        //fix below!
        return 1.0/(n+1.0)*pow(x,n)*(x*asin(x)+sqrt(1-x*x))
                -n/(n+1.0)*integral_x_to_n_times_sqrt_1_minus_x2(x, n-1);
    }
}

//calculates int x^n * sqrt(1-x^2) dx (recursive using per partes)
double integral_x_to_n_times_sqrt_1_minus_x2(double x, int n){
    assert(n >= 0);
    if (n == 0) return 0.5*asin(x)+0.5*x*sqrt(1-x*x);
    else        return 2.0/(n+2.0)*pow(x, n)*(0.5*asin(x)+0.5*x*sqrt(1-x*x))-n/(n+2.0)*integral_x_to_n_times_asin_x(x, n-1);
}

double integrate_x_to_n_times_sqrt_1_minus_x2(double a, double b, int n){
    return integral_x_to_n_times_sqrt_1_minus_x2(b, n)-integral_x_to_n_times_sqrt_1_minus_x2(a, n);
}

double integral_chebyshev(double x, int n){
    assert(n >= 0);
    if (n == 0)         return x;
    else if (n == 1)    return 0.5*x*x;
    else                return 0.5*(fcnc::chebyshev(x, n+1)/(n+1) - fcnc::chebyshev(x, n-1)/(n-1));
}

double integrate_gauss(double sigma, double mean, double min, double max){
    return 0.5*(TMath::Erf((mean-min)/sqrt(2.0)/sigma) - TMath::Erf((mean-max)/sqrt(2.0)/sigma));
}

double integrate_crystalball(double mean, double sigma, double alpha, double n, double min, double max){

    //implementation taken from RooFit: https://root.cern.ch/doc/master/RooCBShape_8cxx_source.html#l00108

    static const double sqrtPiOver2 = 1.2533141373;
    static const double sqrt2 = 1.4142135624;

    double result = 0.0;
    bool useLog = false;

    if( fabs(n-1.0) < 1.0e-05 )
        useLog = true;

    double sig = fabs((Double_t)sigma);

    //double tmin = (m.min(rangeName)-m0)/sig;
    //double tmax = (m.max(rangeName)-m0)/sig;
    double tmin = (min - mean ) / sigma;
    double tmax = (max - mean ) / sigma;

    if(alpha < 0) {
        double tmp = tmin;
        tmin = -tmax;
        tmax = -tmp;
    }

    double absAlpha = fabs(alpha);

    if( tmin >= -absAlpha ) {
        result += sig*sqrtPiOver2*(   TMath::Erf(tmax/sqrt2)
                                      - TMath::Erf(tmin/sqrt2) );
    }
    else if( tmax <= -absAlpha ) {
        double a = TMath::Power(n/absAlpha,n)*exp(-0.5*absAlpha*absAlpha);
        double b = n/absAlpha - absAlpha;

        if(useLog) {
            result += a*sig*( log(b-tmin) - log(b-tmax) );
        }
        else {
            result += a*sig/(1.0-n)*(   1.0/(TMath::Power(b-tmin,n-1.0))
                                        - 1.0/(TMath::Power(b-tmax,n-1.0)) );
        }
    }
    else {
        double a = TMath::Power(n/absAlpha,n)*exp(-0.5*absAlpha*absAlpha);
        double b = n/absAlpha - absAlpha;

        double term1 = 0.0;
        if(useLog) {
            term1 = a*sig*(  log(b-tmin) - log(n/absAlpha));
        }
        else {
            term1 = a*sig/(1.0-n)*(   1.0/(TMath::Power(b-tmin,n-1.0))
                                      - 1.0/(TMath::Power(n/absAlpha,n-1.0)) );
        }

        double term2 = sig*sqrtPiOver2*(   TMath::Erf(tmax/sqrt2)
                                           - TMath::Erf(-absAlpha/sqrt2) );

        result += term1 + term2;
    }

    return result;

}

double integrate_twotailedcrystalball(double mean, double sigma, double alpha1, double alpha2, double n1, double n2, double min, double max){

    double central = 0.0;
    double left = 0.0;
    double right = 0.0;
    double result = 0.0;

    assert(alpha1 > 0.0);
    assert(alpha2 > 0.0);

    //implementation taken from RooFit: RooDoubleCB.cxx

    static const double rootPiBy2 = sqrt(atan2(0.0,-1.0)/2.0);
    static const double invRoot2 = 1.0/sqrt(2);

    double invwidth = 1.0/sigma;

    double tmin = (min-mean)*invwidth;
    double tmax = (max-mean)*invwidth;
    //spdlog::trace("tmin: {0:f}\t tmax: {0:f} ",tmin, tmax);

    bool isfullrange = (tmin<-1000. && tmax>1000.);

    //compute gaussian (central) contribution
    double central_low  = std::max(min, mean - alpha1*sigma);
    double central_high = std::min(max, mean + alpha2*sigma);

    double tcentral_low  = (central_low-mean)*invwidth;
    double tcentral_high = (central_high-mean)*invwidth;

    //spdlog::trace("central_low: {0:f}\t central_high: {0:f} ",central_low, central_high);
    //spdlog::trace("tcentral_low: {0:f}\t tcentral_high: {0:f} ",tcentral_low, tcentral_high);

    if(central_low < central_high){// is the gaussian part in range?
        central = rootPiBy2*sigma*(TMath::Erf(tcentral_high*invRoot2)-TMath::Erf(tcentral_low*invRoot2));
    }

    //compute left tail;
    if(isfullrange  && (n1-1.0)>1.e-3) {
        left = sigma*TMath::Exp(-0.5*alpha1*alpha1)*n1/(alpha1*(n1-1.));
        spdlog::trace("left: {0:f}", left);
    }
    else{
        double left_low  = min;
        double left_high = std::min(max, mean - alpha1*sigma);
        double thigh     = (left_high-mean)*invwidth;

        if(left_low < left_high){ //is the left tail in range?
            double n1invalpha1 = n1/(fabs(alpha1));
            if(fabs(n1-1.0)>1.e-3){
                double invn1m1 = 1/(n1-1.);
                double leftpow = TMath::Power(n1invalpha1,-n1*invn1m1);
                assert(leftpow > 0.0);
                double left0   = sigma*TMath::Exp(-0.5*alpha1*alpha1)*invn1m1;
                double left1, left2;

                if(max > (mean-alpha1*sigma))
                    left1 = n1invalpha1;
                else
                    left1 = TMath::Power( leftpow*(n1invalpha1 - alpha1 - thigh), 1.-n1);

                if(tmin<-1000.)
                    left2 = 0.;
                else
                    left2 = TMath::Power( leftpow*(n1invalpha1 - alpha1 - tmin ), 1.-n1);

                left = left0*(left1-left2);
            }
            else{
                double A1 = TMath::Power(n1invalpha1,n1)*TMath::Exp(-0.5*alpha1*alpha1);
                double B1 = n1invalpha1-fabs(alpha1);
                left = A1*sigma*(TMath::Log(B1-(left_low-mean)*invwidth) - TMath::Log(B1-(left_high-mean)*invwidth) );
            }
        }
    }

    //compute right tail;
    if(isfullrange && (n2-1.0)>1.e-3) {
        right = sigma*TMath::Exp(-0.5*alpha2*alpha2)*n2/(alpha2*(n2-1.));
    }
    else{
        double right_low  = std::max(min,mean + alpha2*sigma);
        double right_high = max;
        double tlow       = (right_low - mean)*invwidth;

        if(right_low < right_high){ //is the right tail in range?
            double n2invalpha2 = n2/(fabs(alpha2));
            if(fabs(n2-1.0)>1.e-3) {
                double invn2m2 = 1.0/(n2-1.);
                double rightpow = TMath::Power(n2invalpha2,-n2*invn2m2);
                if (rightpow <= 0.0){
                    spdlog::error("rightpow < 0.0! Returning NaN");
                    spdlog::error("TMath::Power(n2invalpha2,-n2*invn2m2) = {0:.30f}",powl(n2/(fabs(alpha2)),-n2*1.0/(n2-1.)));
                    spdlog::error("n2invalpha2: {0:f}\tn2: {1:f}\tinvn2m2: {2:f}",n2invalpha2,n2,invn2m2);
                    spdlog::error("mean: {0:f}\tsigma: {1:f}\talpha1: {2:f}\talpha2: {3:f}\tn1: {4:f}\tn2: {5:f}\tmin: {6:f}\tmax: {7:f}", mean, sigma, alpha1, alpha2, n1, n2, min, max);
                    return nan("");
                }
                double right0 = sigma*TMath::Exp(-0.5*alpha2*alpha2)*invn2m2;
                double right1, right2;

                if(min<(mean+alpha2*sigma)) right1 = n2invalpha2;
                else                        right1 = TMath::Power( rightpow*(n2invalpha2 - alpha2 + tlow), 1.-n2);

                if(tmax>1000.) right2 = 0.;
                else           right2 = TMath::Power( rightpow*(n2invalpha2 - alpha2 + tmax), 1.-n2);

                right = right0*(right1-right2);
            }
            else{
                double A2 = TMath::Power(n2invalpha2,n2)*TMath::Exp(-0.5*alpha2*alpha2);
                double B2 = n2invalpha2-fabs(alpha2);
                right = A2*sigma*(TMath::Log(B2+(right_high-mean)*invwidth) - TMath::Log(B2+(right_low-mean)*invwidth) );
            }
        }
    }

    result = left + central + right;

    if(result <= 0.0){
        spdlog::error("left = {0:e}",left);
        spdlog::error("central = {0:e}",central);
        spdlog::error("right = {0:e}", right);
        spdlog::error("sigma = {0:e}", sigma);
        spdlog::error("n1 = {0:e}", n1);
        spdlog::error("n2 = {0:e}", n2);
        spdlog::error("alpha1 = {0:e}", alpha1);
        spdlog::error("alpha1 = {0:e}", alpha2);
        assert(0);
    }

    if(std::isnan(result) || std::isinf(result)){
        spdlog::critical("Determination of twotailed CB integral failed: {0:f}", result);
        spdlog::critical("mean: {1:f}\tsigma: {2:f}\talpha1: {3:f}\talpha2: {4:f}\tn1: {5:f}\tn2: {6:f}", mean, sigma, alpha1, alpha2, n1, n2);
        spdlog::critical("left: {0:f}\tcentral: {1:f}\tright: {2:f}", left, central, right);

        assert(0);
    }

    return result;
}

double integral_x_to_n_times_exp_minus_x(double x, double tau, int n){
    assert(n>=0);
    if (n == 0) return -tau*exp(-x/tau);
    else        return pow(x, n)*(-tau*exp(-x/tau)) + n * tau * integral_x_to_n_times_exp_minus_x(x, tau, n-1);
}



//-----------------------
//        P wave
//-----------------------

//THESE FUNCTIONS NEED TO BE CHECKED
//i correspons to ctl, j to ctk and k to phi


//int[ sin^2(tk) ] for J1s
double integral_f1(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
    //Integrate  sin^2  d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
    return (integral_x_to_n(ctk,n_ctk)-integral_x_to_n(ctk,n_ctk+2))
            *integral_x_to_n(ctl,n_ctl)
            *integral_x_to_n(phi,n_phi);
}
double integrate_f1(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
    return (integrate_x_to_n(ctk_a, ctk_b, n_ctk)-integrate_x_to_n(ctk_a,ctk_b,n_ctk+2))
           *integrate_x_to_n(ctl_a, ctl_b, n_ctl)
           *integrate_x_to_n(phi_a, phi_b, n_phi);
}

//int[ cos^2(tk) ] for J1c
double integral_f2(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
    //Integrate cos^2 thetak d costhetak  d costhetal d phi  times ctk^n1+ctl^n2*phi^n3
    return  integral_x_to_n(ctk,n_ctk+2)
           *integral_x_to_n(ctl,n_ctl)
           *integral_x_to_n(phi,n_phi);
}
double integrate_f2(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
    return  integrate_x_to_n(ctk_a, ctk_b, n_ctk+2)
           *integrate_x_to_n(ctl_a, ctl_b, n_ctl)
           *integrate_x_to_n(phi_a, phi_b, n_phi);
}

//int[ sin^2(tk) * cos(2tl) ] for J2s
double integral_f3(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
    //Integrate sin^2 thetak cos 2thetal d costhetak d costhetal d phi  times ctk^n1+ctl^n2*phi^n3
    return   (integral_x_to_n(ctk,n_ctk)-integral_x_to_n(ctk,n_ctk+2))
            *(2*integral_x_to_n(ctl,n_ctl+2)-integral_x_to_n(ctl,n_ctl)) //cos(2x) = 2*cos^2(x)-1
            *integral_x_to_n(phi,n_phi);
}
double integrate_f3(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
    return   (integrate_x_to_n(ctk_a, ctk_b, n_ctk)-integrate_x_to_n(ctk_a, ctk_b, n_ctk+2))
            *(2*integrate_x_to_n(ctl_a, ctl_b, n_ctl+2)-integrate_x_to_n(ctl_a, ctl_b, n_ctl)) //cos(2x) = 2*cos^2(x)-1
            *integrate_x_to_n(phi_a,phi_b,n_phi);
}

//int[ cos^2(tk) * cos(2tl)  ] for J2c
double integral_f4(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
    //Integrate cos^2 thetak cos 2thetal d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
    return  integral_x_to_n(ctk,n_ctk+2)
           *(2*integral_x_to_n(ctl,n_ctl+2)-integral_x_to_n(ctl,n_ctl))
           *integral_x_to_n(phi,n_phi);
}
double integrate_f4(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
    //Integrate cos^2 thetak cos 2thetal d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
    return integrate_x_to_n(ctk_a, ctk_b, n_ctk+2)
            *(2*integrate_x_to_n(ctl_a, ctl_b, n_ctl+2)-integrate_x_to_n(ctl_a, ctl_b, n_ctl))
            *integrate_x_to_n(phi_a,phi_b,n_phi);
}

//int[ sin^2(tk) * sin^2(tl) * cos(2phi) ] for J3
double integral_f5(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
    //Integrate sin^2 thetak sin^2 thetal cos 2phi d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
    return  (integral_x_to_n(ctk,n_ctk)-integral_x_to_n(ctk,n_ctk+2))
           *(integral_x_to_n(ctl,n_ctl)-integral_x_to_n(ctl,n_ctl+2))
           *(integral_x_to_n_times_cos_2x(phi,n_phi));
}
double integrate_f5(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
    //Integrate sin^2 thetak sin^2 thetal cos 2phi d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
    return  (integrate_x_to_n(ctk_a, ctk_b, n_ctk)-integrate_x_to_n(ctk_a, ctk_b, n_ctk+2))
           *(integrate_x_to_n(ctl_a, ctl_b, n_ctl)-integrate_x_to_n(ctl_a, ctl_b, n_ctl+2))
           *(integrate_x_to_n_times_cos_2x(phi_a,phi_b,n_phi));
}

//int[ sin(2tk) * sin(2tl) * cos(phi) ] for J4
double integral_f6(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
    //Integrate sin 2thetak sin 2thetal cos phi d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
    return 2*integral_x_to_n_times_sqrt_1_minus_x2(ctk,n_ctk+1)
           *2*integral_x_to_n_times_sqrt_1_minus_x2(ctl,n_ctl+1)
           *integral_x_to_n_times_cos_x(phi,n_phi);
}
double integrate_f6(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
    //Integrate sin 2thetak sin 2thetal cos phi d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
    return 2*integrate_x_to_n_times_sqrt_1_minus_x2(ctk_a, ctk_b, n_ctk+1)
           *2*integrate_x_to_n_times_sqrt_1_minus_x2(ctl_a, ctl_b,n_ctl+1)
           *integrate_x_to_n_times_cos_x(phi_a,phi_b,n_phi);
}

//int [  sin(2tk) * sin(tl) * cos(phi) ] for J5
double integral_f7(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
    //Integrate sin 2thetak sin thetal cos phi d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
    return 2*integral_x_to_n_times_sqrt_1_minus_x2(ctk,n_ctk+1)
           *integral_x_to_n_times_sqrt_1_minus_x2(ctl,n_ctl)
           *integral_x_to_n_times_cos_x(phi,n_phi);
}
double integrate_f7(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
    //Integrate sin 2thetak sin thetal cos phi d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
    return 2*integrate_x_to_n_times_sqrt_1_minus_x2(ctk_a, ctk_b, n_ctk+1)
           *integrate_x_to_n_times_sqrt_1_minus_x2(ctl_a, ctl_b,n_ctl)
           *integrate_x_to_n_times_cos_x(phi_a,phi_b,n_phi);
}

//int [ sin^2(tk) * cos(tl) ] for J6s
double integral_f8(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
    //Integrate sin^2 thetak cos thetal d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
    return (integral_x_to_n(ctk,n_ctk)-integral_x_to_n(ctk,n_ctk+2))
            *integral_x_to_n(ctl,n_ctl+1)
            *integral_x_to_n(phi,n_phi);
}
double integrate_f8(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
    //Integrate sin^2 thetak cos thetal d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
    return (integrate_x_to_n(ctk_a, ctk_b, n_ctk)-integrate_x_to_n(ctk_a, ctk_b, n_ctk+2))
            *integrate_x_to_n(ctl_a, ctl_b,n_ctl+1)
            *integrate_x_to_n(phi_a,phi_b,n_phi);
}

//int [ cos^2(tk) * cos(tl) ] for J6c
double integral_f9(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
    //Integrate sin^2 thetak cos thetal d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
    return integral_x_to_n(ctk,n_ctk+2)
           *integral_x_to_n(ctl,n_ctl+1)
           *integral_x_to_n(phi,n_phi);
}
double integrate_f9(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
    //Integrate sin^2 thetak cos thetal d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
    return integrate_x_to_n(ctk_a, ctk_b, n_ctk+2)
           *integrate_x_to_n(ctl_a, ctl_b,n_ctl+1)
           *integrate_x_to_n(phi_a,phi_b,n_phi);
}

//int [ sin(2tk) * sin(tl) * sin(phi) ] for J7
double integral_f10(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
    //Integrate sin 2thetak sin thetal sin phi d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
    return 2*integral_x_to_n_times_sqrt_1_minus_x2(ctk,n_ctk+1)
            *integral_x_to_n_times_sqrt_1_minus_x2(ctl,n_ctl+1)
            *integral_x_to_n_times_sin_x(phi,n_phi);
}
double integrate_f10(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
    //Integrate sin 2thetak sin thetal sin phi d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
    return 2*integrate_x_to_n_times_sqrt_1_minus_x2(ctk_a, ctk_b, n_ctk+1)
            *integrate_x_to_n_times_sqrt_1_minus_x2(ctl_a, ctl_b,n_ctl+1)
            *integrate_x_to_n_times_sin_x(phi_a,phi_b,n_phi);
}

//int [ sin(2tk) * sin(2tl) * sin(phi) ] for J8
double integral_f11(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
    //Integrate sin 2thetak sin 2thetal sin phi d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
    return 2*integral_x_to_n_times_sqrt_1_minus_x2(ctk,n_ctk+1)
           *2*integral_x_to_n_times_sqrt_1_minus_x2(ctl,n_ctl+1)
           *integral_x_to_n_times_sin_x(phi,n_phi);
}
double integrate_f11(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
    //Integrate sin 2thetak sin 2thetal sin phi d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
    return 2*integrate_x_to_n_times_sqrt_1_minus_x2(ctk_a, ctk_b,n_ctk+1)
           *2*integrate_x_to_n_times_sqrt_1_minus_x2(ctl_a, ctl_b,n_ctl+1)
           *integrate_x_to_n_times_sin_x(phi_a,phi_b,n_phi);
}

//int [ sin^2(tk) * sin^2(tl) * sin(2phi) ] for J9
double integral_f12(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
    //Integrate sin^2 thetak sin^2 thetal sin 2phi d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
    return (integral_x_to_n(ctk,n_ctk)-integral_x_to_n(ctk,n_ctk+2))
           *(integral_x_to_n(ctl,n_ctl)-integral_x_to_n(ctl,n_ctl+2))
           *(integral_x_to_n_times_sin_2x(phi,n_phi));
}
double integrate_f12(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
    //Integrate sin^2 thetak sin^2 thetal sin 2phi d costhetak d costhetal d phi times ctk^n1+ctl^n2*phi^n3
    return (integrate_x_to_n(ctk_a, ctk_b,n_ctk)-integrate_x_to_n(ctk_a, ctk_b,n_ctk+2))
           *(integrate_x_to_n(ctl_a, ctl_b,n_ctl)-integrate_x_to_n(ctl_a, ctl_b,n_ctl+2))
           *(integrate_x_to_n_times_sin_2x(phi_a,phi_b,n_phi));
}

//-----------------------
//        S wave
//-----------------------

double integral_s_f1(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
    return (integral_x_to_n(ctl,n_ctl)-integral_x_to_n(ctl,n_ctl+2))
            *integral_x_to_n(ctk,n_ctk)
            *integral_x_to_n(phi, n_phi);
}
double integrate_s_f1(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
    return (integrate_x_to_n(ctl_a,ctl_b,n_ctl)-integrate_x_to_n(ctl_a,ctl_b,n_ctl+2))
            *integrate_x_to_n(ctk_a,ctk_b,n_ctk)
            *integrate_x_to_n(phi_a,phi_b,n_phi);
}

double integral_s_f2(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
    return (integral_x_to_n(ctl,n_ctl)-integral_x_to_n(ctl,n_ctl+2))
            *integral_x_to_n(ctk,n_ctk+1)
            *integral_x_to_n(phi, n_phi);
}
double integrate_s_f2(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
    return (integrate_x_to_n(ctl_a,ctl_b,n_ctl)-integrate_x_to_n(ctl_a,ctl_b,n_ctl+2))
            *integrate_x_to_n(ctk_a,ctk_b,n_ctk+1)
            *integrate_x_to_n(phi_a,phi_b,n_phi);
}

double integral_s_f3(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
    return 2*integral_x_to_n_times_sqrt_1_minus_x2(ctl,n_ctl+1)
            *integral_x_to_n_times_sqrt_1_minus_x2(ctk,n_ctk)
            *integral_x_to_n_times_cos_x(phi, n_phi);
}
double integrate_s_f3(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
    return 2*integrate_x_to_n_times_sqrt_1_minus_x2(ctl_a,ctl_b,n_ctl+1)
            *integrate_x_to_n_times_sqrt_1_minus_x2(ctk_a,ctk_b,n_ctk)
            *integrate_x_to_n_times_cos_x(phi_a,phi_b,n_phi);
}

double integral_s_f4(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
    return  integral_x_to_n_times_sqrt_1_minus_x2(ctl,n_ctl)
           *integral_x_to_n_times_sqrt_1_minus_x2(ctk,n_ctk)
           *integral_x_to_n_times_cos_x(phi, n_phi);
}
double integrate_s_f4(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
    return  integrate_x_to_n_times_sqrt_1_minus_x2(ctl_a,ctl_b,n_ctl)
           *integrate_x_to_n_times_sqrt_1_minus_x2(ctk_a,ctk_b,n_ctk)
           *integrate_x_to_n_times_cos_x(phi_a,phi_b,n_phi);
}

double integral_s_f5(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
    return  integral_x_to_n_times_sqrt_1_minus_x2(ctl,n_ctl)
           *integral_x_to_n_times_sqrt_1_minus_x2(ctk,n_ctk)
           *integral_x_to_n_times_sin_x(phi, n_phi);
}
double integrate_s_f5(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
    return  integrate_x_to_n_times_sqrt_1_minus_x2(ctl_a,ctl_b,n_ctl)
           *integrate_x_to_n_times_sqrt_1_minus_x2(ctk_a,ctk_b,n_ctk)
           *integrate_x_to_n_times_sin_x(phi_a,phi_b,n_phi);
}

double integral_s_f6(double ctl, double ctk, double phi, int n_ctl, int n_ctk, int n_phi){
    return 2*integral_x_to_n_times_sqrt_1_minus_x2(ctl,n_ctl+1)
            *integral_x_to_n_times_sqrt_1_minus_x2(ctk,n_ctk)
            *integral_x_to_n_times_sin_x(phi, n_phi);
}
double integrate_s_f6(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, int n_ctl, int n_ctk, int n_phi){
    return 2*integrate_x_to_n_times_sqrt_1_minus_x2(ctl_a,ctl_b,n_ctl+1)
            *integrate_x_to_n_times_sqrt_1_minus_x2(ctk_a,ctk_b,n_ctk)
            *integrate_x_to_n_times_sin_x(phi_a,phi_b,n_phi);
}