|
|
//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); }
|