/* With the functions in this file 1-dimensional definite integrals can be calculated using either - a Newton-Cotes algorithm : * trapezoidal rule * Simpson's rule - a Gauss-Legendre algorithm. */ #include #include #include #include #define MAXIT 20 /* this gives the max ( <= 32 !) number of iterations that is allowed */ #define FRACACC 1.0e-5 #define EPS 3.0e-11 /* relative precision */ #define FUNC(x) ((*func)(x)) /* func is a pointer to the function to be integrated */ #define FUNC2(x) ((*func2)(x)) /* func2 is a pointer to the function to be integrated */ #define FUNC3(x) ((*func3)(x)) /* func3 is a pointer to the function to be integrated */ #define FUNC4(x) ((*func4)(x)) /* func4 is a pointer to the function to be integrated */ /* The function int_trapzd() - which makes use of the function trapzd() - calculates the integral using the extended trapezoidal rule. See : Press et al., Numerical Recipes in C ,Cambridge 1992, 4.1 and 4.2 */ double int_trapzd(double (*func)(double x), double a, double b) { double previous = 1.0E-30; double integral = 1.0E-25; double fracacc; int j=1; fracacc = fabs( (previous - integral) / previous ); while( (fracacc > FRACACC) && (j <= MAXIT) ) { previous = integral; integral = trapzd(func, a, b, j); fracacc = fabs( (previous - integral) / previous ); j++; } return integral; } /* The function int_simp() - which makes use of the function trapzd() - calculates the integral using the extended simpson rule. See : Press et al., Numerical Recipes in C ,Cambridge 1992, 4.1 and 4.2 */ double int_simp(double (*func)(double x), double a, double b) { double previous = 1.0e-15; double integral = 1.0e-10; double integralterm1 , integralterm2 = 0.0; double fracacc; int j=1; fracacc = fabs( (previous - integral) / previous ); while( (fracacc > FRACACC) && (j <= MAXIT) ) { previous = integral; integralterm1 = trapzd(func, a, b, j++); integral = ((1.0 / 3.0) * (4.0 * integralterm1 - integralterm2)); integralterm2 = integralterm1; fracacc = fabs( (previous - integral) / previous ); } return integral; } double trapzd(double (*func)(double x), double a, double b, int n) { double coord, tnm, sum, del; static double s; long j,it; if (n == 1) { s = .5 * ( b - a) * (FUNC(a) + FUNC(b)); return s; } else { for (it = 1, j = 1; j < (n - 1); j++) it = it << 1; tnm = it; del = (b - a) / tnm; coord = a + .5 * del; for (sum=0.0,j=1; j<=it; j++,coord += del) sum += FUNC(coord); s = .5 * (s + (b - a) *sum / tnm); return s; } } /* The function gauss_leg() calculates the roots and weights for the GAUSS-LEGENDRE n-point quadrature formula. See : Press et al., Numerical Recipes in C, Cambridge UP 1992, section 4.5.*/ void gauss_leg(double x1, double x2, double *x, double *w, int n) { int m, j, i; double z1, z, xm, xl, pp, p3, p2, p1; m = (n + 1) / 2; xm = 0.5 * (x2 + x1); xl = 0.5 * (x2 - x1); for (i = 1; i <= m; i++) /* loop over desired roots */ { z = cos(3.141592654 * (i - 0.25)/(n + 0.5) ); /* Starting approximation to ith root */ do { p1 = 1.0; p2 = 0.0; for (j = 1; j <= n; j++) { p3 = p2; p2 = p1; p1 = ( (2.0 * j - 1.0) * z * p2 - (j - 1.0) * p3) / j ; } pp = n * (z * p1 - p2) / (z * z - 1.0); z1 = z; z = z1 - p1 / pp; /* Newton-Raphson interpolation formula */ } while (fabs(z - z1) > EPS); x[i] = xm - xl * z; /* Roots scaled to desired interval */ x[n + 1 - i] = xm + xl * z; /* symmetric counterparts of roots */ w[i] = 2.0 * xl / ( (1.0 - z * z) * pp * pp ); /* weights */ w[n + 1 - i] = w[i]; } } /* The function int_gauss() calculates a one-dimensional definite integral using a gaussian quadrature formula (in the present case we use the GAUSS-LEGENDRE quadrature formula. - a and b are the lower and upper integration limits respectively, - n stands for the n-point quadrature formula (all polynomials up to degree (2 * n - 1) are integrated exactly), - x[1...n] and w[1...n] contain the abscissa and weights of the n-point Gauss-Legendre quadrature formula and are calculated by the function gauss_leg(). */ double int_gauss(double (*func)(double), double a, double b, int n) { int j; double sum = 0.0; static double *x, *w; x = dvector(1, n); w = dvector(1, n); gauss_leg(a, b, x, w, n); for (j = 1; j <= n; j++) { sum += w[j] * FUNC(x[j]); } free_dvector(x, 1, n); free_dvector(w, 1, n); return sum; } double int_gauss2(double (*func2)(double), double a, double b, int n) { int j; double sum = 0.0; static double *x2, *w2; x2 = dvector(1, n); w2 = dvector(1, n); gauss_leg(a, b, x2, w2, n); for (j = 1; j <= n; j++) { sum += w2[j] * FUNC2(x2[j]); } free_dvector(x2, 1, n); free_dvector(w2, 1, n); return sum; } double int_gauss3(double (*func3)(double), double a, double b, int n) { int j; double sum = 0.0; static double *x3, *w3; x3 = dvector(1, n); w3 = dvector(1, n); gauss_leg(a, b, x3, w3, n); for (j = 1; j <= n; j++) { sum += w3[j] * FUNC3(x3[j]); } free_dvector(x3, 1, n); free_dvector(w3, 1, n); return sum; } double int_gauss4(double (*func4)(double), double a, double b, int n) { int j; double sum = 0.0; static double *x4, *w4; x4 = dvector(1, n); w4 = dvector(1, n); gauss_leg(a, b, x4, w4, n); for (j = 1; j <= n; j++) { sum += w4[j] * FUNC4(x4[j]); } free_dvector(x4, 1, n); free_dvector(w4, 1, n); return sum; } /* The function int_gauss_infinit() calculates a one-dimensional definite integral from zero to infinity using the GAUSS-LEGENDRE quadrature formula. - c is an arbitrary transformation factor that can be used to check the routine (the result should be independent of c), - n stands for the n-point quadrature formula (all polynomials up to degree (2 * n - 1) are integrated exactly), - x[1...n] and w[1...n] contain the abscissa and weights of the n-point Gauss-Legendre quadrature formula and are calculated by the function gauss_leg(). */ double int_gauss_infinit(double (*func)(double), double c, int n) { int j; double sum = 0.0, *q, *W; static double *x, *w; x = dvector(1, n); w = dvector(1, n); q = dvector(1, n); W = dvector(1, n); gauss_leg(- 0.99999, 0.99999, x, w, n); for (j = 1; j <= n; j++) { q[j] = c * tan(3.141592654 / 4.0 * (1.0 + x[j]) ); W[j] = c * 3.141592654 / 4.0 * w[j] * 1.0 / pow( cos(3.141592654 / 4.0 * (1.0 + x[j]) ), 2.0); sum += W[j] * FUNC(q[j]); } free_dvector(x, 1, n); free_dvector(w, 1, n); free_dvector(q, 1, n); free_dvector(W, 1, n); return sum; } /* The function gauss_leg_infinit() calculates the roots and weights for the GAUSS-LEGENDRE n-point quadrature formula for integration limits zero and infinity. */ void gauss_leg_infinit(double *q, double *W, double c, int n) { int j; static double *x, *w; x = dvector(1, n); w = dvector(1, n); gauss_leg(- 0.99999, 0.99999, x, w, n); for (j = 1; j <= n; j++) { q[j] = c * tan(3.141592654 / 4.0 * (1.0 + x[j]) ); W[j] = c * 3.141592654 / 4.0 * w[j] * 1.0 / pow( cos(3.141592654 / 4.0 * (1.0 + x[j]) ), 2.0); } free_dvector(x, 1, n); free_dvector(w, 1, n); } /* FONCTIONS DE POIDS A UTILISER DANS MULTIINT.C*/ /* a=(Ebeam[MeV]/melectron)^2 et b=1 */ /* - W1(x) = 1/(ax+b) */ /* - W2(x) = 1/(ax+b)^2 */ /* - W3(x) = 1/(ax+b)^3 */ double W1(double x, double a, double b ) { double W; W= 1./( a*x + b ); return W; } double W2(double x, double a, double b ) { double W; W= 1./pow(( a*x + b ),2.); return W; } double W3(double x, double a, double b ) { double W; W= 1./pow(( a*x + b ),3.); return W; } void MOMENTS1(int nmax, double a, double b,double x1, double x2, double *ct) { int i,j,idif; double t0, t1, sum, phas, den, s1, y; double *yt, *tt; yt=dvector(1,nmax); tt=dvector(1,nmax); t0=log( fabs(a*x2+b)/fabs(a*x1+b) ); ct[0]=t0/a; yt[1]=b/a; tt[1]=x2-x1; for ( i=2; i<=nmax; i++) { yt[i]=yt[i-1]*b/a; tt[i]=pow(x2,i) - pow(x1,i); } y=b/a; t1=x2-x1; ct[1]=t1-y*t0; ct[1]=ct[1]/a; for (i=2; i<=nmax; i++) { sum=0.; phas=-1.; for (j=1; j<=i; j++) { idif=i-j; den=(double)idif; if (idif > 0) s1=tt[i-j]*yt[j]/den; else s1=yt[i]*t0; sum=sum+s1*phas; phas=-phas; } sum=sum+tt[i]/((double)i); ct[i]=sum/a; } free_dvector(yt, 1, nmax); free_dvector(tt, 1, nmax); } void MOMENTS2(int nmax, double a, double b,double x1, double x2, double *ct1, double *ct2) { int i; double di; for (i=1; i<=nmax; i++) { di=(double)i; ct2[i]= pow(x1,di)/(a*x1+b) -pow(x2,di)/(a*x2+b) +di*ct1[i-1]; ct2[i]=ct2[i]/a; } ct2[0]= 1./(a*x1+b) -1./(a*x2+b); ct2[0]=ct2[0]/a; } void MOMENTS3(int nmax, double a, double b, double c, double x1, double x2, double *ct3) /* fonction de poids f(x)=1/(a*x^2+b*x+c) */ { double t0, delta,deltap,rnum2,rden2,rnum1,rden1,t2,t1,y2,y1; int m; t0= log( fabs((a*x2*x2+b*x2+c) / (a*x1*x1+b*x1+c)) ); delta=b*b -4.*a*c; deltap=-delta; if (delta > 0.) { rnum2= sqrt(delta) - ( b + 2.*a*x2); rden2= sqrt(delta) + ( b + 2.*a*x2); rnum1= sqrt(delta) - ( b + 2.*a*x1); rden1= sqrt(delta) + ( b + 2.*a*x1); t2=rnum2/rden2; t1=rnum1/rden1; } if (delta == 0.) { y2=b+2.*a*x2; y1=b+2.*a*x1; ct3[0]=-2.*(1./y2 - 1./y1); printf(" ct30 = %e\n", ct3[0]); } if (delta < 0.) { y2= atan( (b+2.*a*x2)/sqrt(deltap) ); y1= atan( (b+2.*a*x1)/sqrt(deltap) ); ct3[0]=2./sqrt(deltap) * (y2-y1); } ct3[1]=1./(2.*a)*t0 -b/(2.*a) *ct3[0]; ct3[2]= 1./a * (x2-x1) - b/(2.*a*a) * t0 +(b*b -2.*a*c)/(2.*a*a) * ct3[0]; for (m=3; m<=nmax; m++) { ct3[m]= -1./(a*(1.-(double)m)) * ( pow(x2,(double)m-1.) - pow(x1,(double)m-1.) ) - b/a * ct3[m-1] - c/a * ct3[m-2]; } } void MOMENTS3bis(int nmax, double a, double b, double c, double x1, double x2, double *ct3) /* fonction de poids f(x)=1/(a*x^2+b*x+c) */ { double t0, delta,deltap,rnum2,rden2,rnum1,rden1,t2,t1,y2,y1; int m; t0= log( fabs((a*x2*x2+b*x2+c) / (a*x1*x1+b*x1+c)) ); delta=b*b -4.*a*c; deltap=-delta; if (delta > 0.) { rnum2= sqrt(delta) - ( b + 2.*a*x2); rden2= sqrt(delta) + ( b + 2.*a*x2); rnum1= sqrt(delta) - ( b + 2.*a*x1); rden1= sqrt(delta) + ( b + 2.*a*x1); t2=rnum2/rden2; t1=rnum1/rden1; ct3[0]=log( fabs(t2/t1) ); ct3[0]= ct3[0]/sqrt(delta); } if (delta == 0.) { y2=b+2.*a*x2; y1=b+2.*a*x1; ct3[0]=-2.*(1./y2 - 1./y1); } if (delta < 0.) { y2= atan( (b+2.*a*x2)/sqrt(deltap) ); y1= atan( (b+2.*a*x1)/sqrt(deltap) ); ct3[0]=2./sqrt(deltap) * (y2-y1); } ct3[1]=1./(2.*a)*t0 -b/(2.*a) *ct3[0]; ct3[2]= 1./a * (x2-x1) - b/(2.*a*a) * t0 +(b*b -2.*a*c)/(2.*a*a) * ct3[0]; for (m=3; m<=nmax; m++) { ct3[m]= -1./(a*(1.-(double)m)) * ( pow(x2,(double)m-1.) - pow(x1,(double)m-1.) ) - b/a * ct3[m-1] - c/a * ct3[m-2]; } } void MOMENTS4(int nmax, double a, double b, double c, double x1, double x2, double *ct3, double *ct4) /* fonction de poids f(x)=1/(a*x^2+b*x+c)^2 */ { double delta,deltap,t2,t1,u2,u1,t0,y; int m,k; delta= b*b -4.*a*c; deltap=-delta; t2= ( b + 2.*a*x2 ) / ( a*x2*x2 + b*x2 + c ); t1= ( b + 2.*a*x1 ) / ( a*x1*x1 + b*x1 + c ); ct4[0]= t2-t1 + 2.*a*ct3[0]; ct4[0]= ct4[0]/deltap; t2=( 2.*c + b*x2 )/ ( a*x2*x2 + b*x2 +c ); t1=( 2.*c + b*x1 )/ ( a*x1*x1 + b*x1 +c ); ct4[1]= -(t2-t1) - b*ct3[0]; ct4[1]= ct4[1]/deltap; t2= 1. / ( a*x2*x2 + b*x2 +c ); t1= 1. / ( a*x1*x1 + b*x1 +c ); u2= x2 / ( a*x2*x2 + b*x2 +c ); u1= x1 / ( a*x1*x1 + b*x1 +c ); ct4[2]= c*b/a * (t2-t1) + ( b*b - 2.*a*c )/a * (u2-u1) + 2.*c*ct3[0]; ct4[2]= ct4[2]/deltap; t0= log( fabs( (a*x2*x2 + b*x2 + c) / (a*x1*x1 + b*x1 + c) ) ); y= c*( 2.*a*c - b*b ) / (a*a) * (t2-t1) + b*( 3.*a*c - b*b ) / (a*a) * (u2-u1) - b*( 6.*a*c -b*b) / (2.*a*a) * ct3[0]; y=y/deltap; ct4[3]= 1./ (2.*a*a) * t0 +y; for (m=4; m<=nmax; m++) { ct4[m]= -1./ ((double)(3-m)*a) * ( t2*pow(x2,(double)(m-1)) -t1*pow(x1,(double)(m-1)) ) -(double)(2-m)/(double)(3-m) * b/a * ct4[m-1] +(double)(m-1)/(double)(3-m) * c/a * ct4[m-2]; } } double fact(int n) { int i; double f; f=1.; for (i=n; i>=1; --i) f=f*i; return f; } double producttermes(int m, int n) { int i,p,init; init=1; for (i=0; i<=n-1; i++) init=init*(m+i); p=init; return p; } void gl_f3(int m_max, double a, double b, double c, int nb_tranche, double *xmin, double *xmax, double *anu3_se2) { int i, ntr, npt, m; int nb_point=10; double x1, w1; double *xgauss3, *wgauss3; wgauss3 = dvector(1,nb_point); xgauss3 = dvector(1,nb_point); for (i=0; i<=m_max; i++) anu3_se2[i]=0.0; gauss_leg(-1.,1.,xgauss3,wgauss3,nb_point); for (ntr=1; ntr<=nb_tranche; ntr++) { for (npt=1; npt<=nb_point; npt++) { x1 = ( (xmax[ntr]-xmin[ntr])/2. )*xgauss3[npt] + (xmax[ntr]+xmin[ntr])/2. ; w1 = wgauss3[npt]*(xmax[ntr]-xmin[ntr])/2.; for (m=0; m<=m_max; m++) anu3_se2[m] = anu3_se2[m] + w1 * pow(x1,m)/(a*x1*x1 + b*x1 + c); } } free_dvector(wgauss3,1,nb_point); free_dvector(xgauss3,1,nb_point); } void gl_f4(int m_max, double a, double b, double c, int nb_tranche, double *xmin, double *xmax, double *anu4_se2) { int i, ntr, npt, m; int nb_point=10; double x1, w1; double *xgauss4, *wgauss4; wgauss4 = dvector(1,nb_point); xgauss4 = dvector(1,nb_point); for (i=0; i<=m_max; i++) anu4_se2[i]=0.0; gauss_leg(-1.,1.,xgauss4,wgauss4,nb_point); for (ntr=1; ntr<=nb_tranche; ntr++) { for (npt=1; npt<=nb_point; npt++) { x1 = ( (xmax[ntr]-xmin[ntr])/2. )*xgauss4[npt] + (xmax[ntr]+xmin[ntr])/2. ; w1 = wgauss4[npt]*(xmax[ntr]-xmin[ntr])/2.; for (m=0; m<=m_max; m++) anu4_se2[m] = anu4_se2[m] + w1 * pow(x1,m)/pow((a*x1*x1 + b*x1 + c),2.); } } free_dvector(wgauss4,1,nb_point); free_dvector(xgauss4,1,nb_point); } /* The function int_cmn() calculates the integration over x of the ratio x**m/(ax+b)**n. According to the value of a/b, the integral is evaluated with threee different methods: - a/b<0.02 Expansion in powers of a/b. - 0.021 Analytic calculation (recurrence relation) . */ void int_cmn(int nmax, double a, double b, double x1, double x2, int *mmax, double **cmn) { int m, n, s; double eps1=1.0, eps2=0.02, ratio; double num1, num2, denom, sum; double pwr, func; double *anu1; ratio = a/b; if (fabs(ratio)>=eps1) { anu1 = dvector(0,mmax[1]); MOMENTS1(mmax[1], a, b, x1, x2, anu1); for (m=0; m<=mmax[0]; m++) cmn[m][0] = (double)1./(double)(m+1) * (pow(x2,m+1) - pow(x1,m+1)); for (m=0; m<=mmax[1]; m++) cmn[m][1] = anu1[m]; for (n=2; n<=nmax; n++) { cmn[0][n] = (double)1./pow(a*x1+b,n-1) - (double)1./pow(a*x2+b,n-1); cmn[0][n] = cmn[0][n]/a/(double)(n-1); for (m=1; m<=mmax[n]; m++) { cmn[m][n] = pow(x1,(double)m) / pow(a*x1+b,(double)(n-1)) - pow(x2,(double)m) / pow(a*x2+b,(double)(n-1)) + (double)m * cmn[m-1][n-1]; cmn[m][n] = cmn[m][n]/a/(double)(n-1); } } free_dvector(anu1, 0, mmax[1]); } if (fabs(ratio)eps2) { int ngl=40; int igl; double y; double *xgauss, *wgauss; xgauss = dvector(1, ngl); wgauss = dvector(1, ngl); gauss_leg(x1, x2, xgauss, wgauss, ngl); for (n=0; n<=nmax; n++) { for (m=0; m<=mmax[n]; m++) { sum = 0.0; for (igl=1; igl<=ngl; igl++) { y = xgauss[igl]; func = pow(y, m) / pow( a*y+b, n); sum = sum + wgauss[igl]*func; } cmn[m][n] = sum; } } free_dvector(xgauss, 1, ngl); free_dvector(wgauss, 1, ngl); } if (fabs(ratio)<=eps2) { int jflag, kflag, lflag; double rho; double r2, r3, r4, r5, r6, r7, r8, r9; double c1, c2, c3, c4, c5, c6, c7, c8, c9; double t0, t1, t2, t3, t4, t5, t6, t7, t8, t9; double f0, f1, f2, f3, f4, f5, f6, f7, f8, f9; double term0, term1, term2, term3, term4, term5, term6, term7, term8, term9; for (m=0; m<=mmax[0]; m++) cmn[m][0] = (double)1./(double)(m+1) * (pow(x2,m+1) - pow(x1,m+1)); for (n=1; n<=nmax; n++) { for (m=0; m<=mmax[n]; m++) { t0 = (double)1./(double)(m+1); f0 = pow(x2,m+1) - pow(x1,m+1); term0 = t0 * f0; c1 = (double)n; t1 = - c1/(double)(m+2) * ratio; f1 = pow(x2,m+2) - pow(x1,m+2); term1 = t1 * f1; r2 = ratio * ratio; c2 = c1 * (double)(n+1); t2 = c2 * r2 / fact(2) / (double)(m+3); f2 = pow(x2,m+3) - pow(x1,m+3); term2 = t2 * f2; r3 = r2 * ratio; c3 = c2 * (double)(n+2); t3 = -c3 * r3 / fact(3) / (double)(m+4); f3 = pow(x2,m+4) - pow(x1,m+4); term3 = t3 * f3; jflag = 0; kflag = 0; lflag = 0; rho = fabs(term3/term0); if (rho<1.E-20) { cmn[m][n] = pow(b, -n) * (term0+term1+term2+term3); jflag = 1; kflag = 1; lflag = 1; } if (jflag==0) { r4 = r3 * ratio; c4 = c3 * (double)(n+3); t4 = c4 * r4 / fact(4) / (double)(m+5); f4 = pow(x2,m+5) - pow(x1,m+5); term4 = t4 * f4; r5 = r4 * ratio; c5 = c4 * (double)(n+4); t5 = -c5 * r5 / fact(5) / (double)(m+6); f5 = pow(x2,m+6) - pow(x1,m+6); term5 = t5 * f5; rho = fabs(term5/term0); if (rho<1.E-20) { cmn[m][n] = pow(b, -n) * (term0+term1+term2+term3+term4+term5); kflag = 1; lflag = 1; } } if (kflag==0) { r6 = r5 * ratio; c6 = c5 * (double)(n+5); t6 = c6 * r6 / fact(6) / (double)(m+7); f6 = pow(x2,m+7) - pow(x1,m+7); term6 = t6 * f6; r7 = r6 * ratio; c7 = c6 * (double)(n+6); t7 = -c7 * r7 / fact(7) / (double)(m+8); f7 = pow(x2,m+8) - pow(x1,m+8); term7 = t7 * f7; rho = fabs(term7/term0); if (rho<1.E-20) { cmn[m][n] = pow(b, -n) * (term0+term1+term2+term3+term4+term5+term6+term7); lflag = 1; } } if (lflag==0) { r8 = r7 * ratio; c8 = c7 * (double)(n+7); t8 = c8 * r8 / fact(8) / (double)(m+9); f8 = pow(x2,m+9) - pow(x1,m+9); term8 = t8 * f8; r9 = r8 * ratio; c9 = c8 * (double)(n+8); t9 = -c9 * r9 / fact(9) / (double)(m+10); f9 = pow(x2,m+10) - pow(x1,m+10); term9 = t9 * f9; cmn[m][n] = pow(b, -n) * (term0+term1+term2+term3+term4+term5+term6+term7+term8+term9); } } } } } /* The function int_dmn() calculates the integration over x of the ratio x**m/(ax**2+bx+c)**n (n=1 or 2). According to the value of a/b, the integral is evaluated with threee different methods: - a/c<0.02 Expansion in powers of a/b. - 0.021 Analytic calculation (recurrence relation) . */ void int_dmn(double a, double b, double c, double x1, double x2, int *mmax, double **dmn) { int i, m, n, nmax=2; double eps1= 1.0, eps2= 0.02, ratio, temp; double delta, t0, t1, t2, u1, u2; double num1, num2, denom1, denom2; double sum, func; ratio = a/c; if (fabs(ratio)>=eps1) { delta = b*b - (double)4.*a*c; t0 = log( fabs( (a*x2*x2+b*x2+c)/(a*x1*x1+b*x1+c) ) ); /* n=1 */ if (delta>(double)0.) { num1 = sqrt(delta) - (b + (double)2.*a*x1); denom1 = sqrt(delta) + (b + (double)2.*a*x1); num2 = sqrt(delta) - (b + (double)2.*a*x2); denom2 = sqrt(delta) + (b + (double)2.*a*x2); t1 = num1/denom1; t2 = num2/denom2; dmn[0][1] = log(fabs(t2/t1)) / sqrt(delta); } if (delta==(double)0.) { t1 = b + (double)2.*a*x1; t2 = b + (double)2.*a*x2; dmn[0][1] = -(double)2. * ((double)1./t2 - (double)1./t1); } if (delta<(double)0.) { t1 = atan( (b+(double)2.*a*x1)/sqrt(-delta) ); t2 = atan( (b+(double)2.*a*x2)/sqrt(-delta) ); dmn[0][1] = (double)2./sqrt(-delta) * (t2-t1); } if (mmax[1]>=1) { dmn[1][1] = ( (double)1.*t0 - b*dmn[0][1]) / ((double)2.*a); if (mmax[1]>=2) { dmn[2][1] = (double)1./a * (x2-x1) + ( - b*t0 + (b*b - (double)2.*a*c) * dmn[0][1] ) / ((double)2.*a*a); for (m=3; m<=mmax[1]; m++) { dmn[m][1] = (double)1./(a*(double)(m-1)) * (pow(x2,m-1) - pow(x1,m-1)) - b/a * dmn[m-1][1] - c/a * dmn[m-2][1]; } } } /* n=2 */ temp = dmn[0][1]; t1 = (b+(double)2.*a*x1) / (a*x1*x1+b*x1+c); t2 = (b+(double)2.*a*x2) / (a*x2*x2+b*x2+c); dmn[0][2] = ( t2-t1 + (double)2.*a*temp ) / (-delta); if (mmax[2]>=1) { t1 = (b*x1+(double)2.*c) / (a*x1*x1+b*x1+c); t2 = (b*x2+(double)2.*c) / (a*x2*x2+b*x2+c); dmn[1][2] = ( t2-t1 + b*temp ) / delta; if (mmax[2]>=2) { t1 = (double)1. / (a*x1*x1+b*x1+c); t2 = (double)1. / (a*x2*x2+b*x2+c); u1 = x1 / (a*x1*x1+b*x1+c); u2 = x2 / (a*x2*x2+b*x2+c); dmn[2][2] = ( c*b/a*(t2-t1) + (b*b-(double)2.*a*c)/a*(u2-u1) + (double)2.*c*temp ) / (-delta); if (mmax[2]>=3) { dmn[3][2] = c*((double)2.*a*c-b*b)/(a*a) * (t2-t1) + b*((double)3.*a*c-b*b)/(a*a) * (u2-u1) - b*((double)6.*a*c-b*b)/((double)2.*a*a) * temp; dmn[3][2] = dmn[3][2]/(-delta) + (double)1./((double)2.*a*a) * t0; for (m=4; m<=mmax[2]; m++) { dmn[m][2] = - (double)1./((double)(3-m)*a) * ( t2*pow(x2,m-1) - t1*pow(x1,m-1) ) - (double)(2-m)/(double)(3-m) * b/a * dmn[m-1][2] + (double)(m-1)/(double)(3-m) * c/a * dmn[m-2][2]; } } } } } if (fabs(ratio)eps2) { int ngl=40; int igl; double y; double *xgauss, *wgauss; xgauss = dvector(1, ngl); wgauss = dvector(1, ngl); gauss_leg(x1, x2, xgauss, wgauss, ngl); for (n=1; n<=nmax; n++) { for (m=0; m<=mmax[n]; m++) { sum = 0.0; for (igl=1; igl<=ngl; igl++) { y = xgauss[igl]; func = pow(y, m) / pow( a*y*y+b*y+c, n); sum = sum + wgauss[igl]*func; } dmn[m][n] = sum; } } free_dvector(xgauss, 1, ngl); free_dvector(wgauss, 1, ngl); } if (fabs(ratio)<=eps2) { int m_cmn, n_cmn, is, l; int mmax_dmn, order=6; int *mmax_cmn; double s1, s2, t1, t2; double **cmn; mmax_dmn = mmax[1]; for (l=2; l<=nmax; l++) if (mmax[l]>mmax_dmn) mmax_dmn=mmax[l]; n_cmn = order; m_cmn = 2*(order-1) + mmax_dmn; mmax_cmn = ivector(0,n_cmn); mmax_cmn[n_cmn] = 2*(n_cmn-1) + mmax_dmn; for (l=n_cmn-1; l>=0; l--) mmax_cmn[l] = mmax_cmn[l+1] -1; cmn = dmatrix(0, m_cmn, 0, n_cmn); int_cmn(n_cmn, b/c, 1., x1, x2, mmax_cmn, cmn); for (m=0; m<=mmax[1]; m++) { s1 = cmn[m][1]/c; for (is=2; is<=order; is++) { t1 = pow(-a/c, is-1); s1 = s1 + (double)1./c * t1 * cmn[2*(is-1)+m][is]; } dmn[m][1] = s1; } for (m=0; m<=mmax[2]; m++) { s2 = cmn[m][2]/pow(c,2); for (is=3; is<=order; is++) { t2 = pow(-a/c, is-2) * (double)(is-1); s2 = s2 + pow((double)1./c, 2) * t2 * cmn[2*(is-2)+m][is]; } dmn[m][2] = s2; } free_ivector(mmax_cmn, 0, n_cmn); free_dmatrix(cmn, 0, m_cmn, 0, n_cmn); } } /* The function int_cmn_anal() gives the analytic result of the integration over x of the ratio x**m/(ax+b)**n (n=1 or 2) when the denominateur has NOT any zeros in the domain of integration ([0,1]). */ void int_cmn_anal(double a, double b, int *mmax, double **cmn) { double a2, a3, a4, b2, b3, b4; a2 = pow(a,2); a3 = pow(a,3); a4 = pow(a,4); b2 = pow(b,2); b3 = pow(b,3); b4 = pow(b,4); cmn[0][1] = log(fabs(1+a/b)) / a; cmn[1][1] = (1 - b*cmn[0][1]) / a; if ( mmax[1]>=2 ) cmn[2][1] = (0.5 - b/a + b2/a*cmn[0][1]) / a; if ( mmax[1]>=3 ) cmn[3][1] = (1.0/3.0 - b/(2*a) + b2/a2 - b3/a2*cmn[0][1]) / a; if ( mmax[1]>=4 ) cmn[4][1] = (1.0/4.0 - b/(3*a) + b2/(2*a2) - b3/a3 + b4/a3*cmn[0][1]) / a; cmn[0][2] = 1 / (b*(a+b)); cmn[1][2] = ( cmn[0][1] - 1/(a+b) ) / a; if ( mmax[2]>=2 ) cmn[2][2] = ( (a+2*b)/(a+b) - 2*b*cmn[0][1] ) / a2; if ( mmax[2]>=3 ) cmn[3][2] = ( (a2-3*a*b-6*b2)/(2*(a+b)) + 3*b2*cmn[0][1] ) / a3; if ( mmax[2]>=4 ) cmn[4][2] = ( (a3-2*a2*b+6*a*b2+10*b3)/(3*(a+b)) - 4*b3*cmn[0][1] ) / a4; } /* The function int_dmn_anal() gives the analytic result of the integration over x of the ratio x**m/(ax**2+bx+c)**n (n=1 or 2) when the denominator has NOT any zeros in the domain of integration ([0,1]). */ void int_dmn_anal(double a, double b, double c, int *mmax, double **dmn) { double delta, r1, r2; delta = b*b - (double)4.*a*c; if (delta>0) { r1 = (-b+sqrt(delta)) / (2*a); r2 = (-b-sqrt(delta)) / (2*a); dmn[0][1] = log( fabs( (1-1/r1)/(1-1/r2) ) ) / ( a*(r1-r2) ); } else { dmn[0][1] = (double)2./sqrt(-delta) * ( atan((b+(double)2.*a)/sqrt(-delta)) - atan(b/sqrt(-delta)) ); } dmn[1][1] = ( r1*log(fabs(1-1/r1)) - r2*log(fabs(1-1/r2)) ) / ( a*(r1-r2) ); if ( mmax[1]>=2 ) dmn[2][1] = ( log(fabs((a+b+c)/c)) - b*dmn[0][1] ) / (2*a) ; dmn[0][2] = ( a*b + b*b - 2*a*c ) / ( c*delta*(a+b+c) ) - 2*a/delta * dmn[0][1]; dmn[1][2] = -(2*a + b) / ( delta*(a+b+c) ) + b/delta * dmn[0][1]; if ( mmax[2]>=2 ) dmn[2][2] = (2*c + b) / ( delta*(a+b+c) ) - 2*c/delta * dmn[0][1]; if ( mmax[2]>=3 ) dmn[3][2] = log( fabs((a+b+c)/c) ) / (2*a*a) + ( 2*a*c - b*c - b*b ) / ( a*delta*(a+b+c) ) + b*(6*a*c - b*b) / (2*a*a*delta) * dmn[0][1]; if ( mmax[2]>=4 ) dmn[4][2] = (1/(a+b+c) - 2*b*dmn[3][2] - 3*c*dmn[2][2])/a; } /* The function int_cmn_cp() calculates in the complex plane the integration over x [0,1] of the ratio x**m/(ax+b)**n (n=1 or 2) when the possible zeros of the denominateur are out of the domain of integration. (Residus of the poles are not treated.) Domain of integration is an half-circle of radius 1/2 and centered at (1/2,0i). */ void int_cmn_cp(double a, double b, dcomplex **res, int m_max, int n_max, double epsilon) { int i, ngl_theta=150; int alpha, beta, degre; double *theta_gauss, *wgauss, theta_max=PI; dcomplex x, sqr_x, arg_x; dcomplex jacob, Bc, d1, d2; for (alpha=0; alpha<=m_max; alpha++) { for (beta=0; beta<=n_max; beta++) { res[alpha][beta] = Complex(0,0); } } theta_gauss = dvector(1,ngl_theta); wgauss = dvector(1,ngl_theta); gauss_leg(0, PI, theta_gauss, wgauss, ngl_theta); for (i=1; i<=ngl_theta; i++) { arg_x = Complex( cos(theta_gauss[i]), - epsilon*sin(theta_gauss[i]) ); jacob = Cmul( arg_x, Complex(0, epsilon*0.5) ); x = CRmul( 0.5, Cadd( Complex(1,0), arg_x ) ); sqr_x = Cmul( x, x ); Bc = Cadd( CRmul(a,x), Complex(b,0) ); d1 = Cmul( jacob, Cinv(Bc) ); d2 = Cmul( jacob, Cinv( Cmul(Bc,Bc) ) ); res[0][1] = Cadd( res[0][1], CRmul(wgauss[i],d1) ); res[0][2] = Cadd( res[0][2], CRmul(wgauss[i],d2) ); for (degre=1; degre<=m_max; degre++) { res[degre][1] = Cadd( res[degre][1], CRmul( wgauss[i],Cmul(Cpow(x,degre),d1) ) ); if (n_max>=2) { res[degre][2] = Cadd( res[degre][2], CRmul( wgauss[i],Cmul(Cpow(x,degre),d2) ) ); } } } /* End loop Gauss */ } /* The function int_dmn_cp() calculates in the complex plane the integration over x [0,1] of the ratio x**m/(ax**2+bx+c)**n (n=1 or 2) when the possible zeros of the denominateur are out of the domain of integration. (Residus of the poles are not treated.) Domain of integration is an half-circle of radius 1/2 and centered at (1/2,0i). */ void int_dmn_cp(double a, double b, double c, dcomplex **res, int m_max, double epsilon) { int i, ngl_theta=150; int alpha, beta, degre; double *theta_gauss, *wgauss, theta_max=PI; dcomplex x, sqr_x, arg_x; dcomplex jacob, Ac, d1, d2; for (alpha=0; alpha<=m_max; alpha++) { for (beta=0; beta<=2; beta++) { res[alpha][beta] = Complex(0,0); } } theta_gauss = dvector(1,ngl_theta); wgauss = dvector(1,ngl_theta); gauss_leg(0, PI, theta_gauss, wgauss, ngl_theta); for (i=1; i<=ngl_theta; i++) { arg_x = Complex( cos(theta_gauss[i]), - epsilon*sin(theta_gauss[i]) ); jacob = Cmul( arg_x, Complex(0, epsilon*0.5) ); x = CRmul( 0.5, Cadd( Complex(1,0), arg_x ) ); sqr_x = Cmul( x, x ); Ac = Cadd( CRmul(a,sqr_x), Cadd( CRmul(b,x), Complex(c,0) ) ); d1 = Cmul( jacob, Cinv(Ac) ); d2 = Cmul( jacob, Cinv( Cmul(Ac,Ac) ) ); res[0][1] = Cadd( res[0][1], CRmul(wgauss[i],d1) ); res[1][1] = Cadd( res[1][1], CRmul( wgauss[i],Cmul(x,d1) ) ); res[0][2] = Cadd( res[0][2], CRmul(wgauss[i],d2) ); res[1][2] = Cadd( res[1][2], CRmul( wgauss[i],Cmul(x,d2) ) ); for (degre=2; degre<=m_max; degre++) { res[degre][1] = Cadd( res[degre][1], CRmul( wgauss[i],Cmul(Cpow(x,degre),d1) ) ); res[degre][2] = Cadd( res[degre][2], CRmul( wgauss[i],Cmul(Cpow(x,degre),d2) ) ); } } /* printf("Methode Plan Complexe : "); Cecho(res[1][1]); printf("\n"); */ } /* The function int_cmn_cp_anal() calculates analytically the real and imaginary parts of the integral over x [0,1]: x**m/(ax+b)**n (n=1 or 2) NB: EPSILON = +1 (-1) IF THE POLE IS ABOVE (BELOW) THE REAL AXIS !!*/ void int_cmn_cp_anal(double xmin, double xmax, double a, double b, dcomplex **res, int m_max, int n_max, double epsilon) { int alpha, beta, i, m, n; double *fact, *J; double pole, sign_a, pascal; pole = -b/a; sign_a = 1.; if(a<0) sign_a = -1.; for (alpha=0; alpha<=m_max; alpha++) { for (beta=0; beta<=n_max; beta++) { res[alpha][beta] = Complex(0,0); } } /*----------------------------------------------------*/ fact = dvector(0, m_max); fact[0] = 1; for(n=1; n<=m_max; n++) { fact[n] = fact[n-1] * (double)n; } /*_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/*/ /*_/ n=1 _/*/ /*_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/*/ J = dvector(1, m_max+1); J[1] = log( (xmax-pole)/(pole-xmin) ); for(n=2; n<=m_max+1; n++) { J[n] = (pow(xmax-pole, n-1) - pow(xmin-pole, n-1)) / (n-1); } res[0][1].re = J[1]; for(m=1; m<=m_max; m++) { for(i=m; i>=0; i--) { pascal = fact[m]/(fact[i]*fact[m-i]); res[m][1].re = res[m][1].re + J[i+1] * pascal * pow(pole, m-i); } } for (m=0; m<=m_max; m++) { res[m][1].re = res[m][1].re / a; res[m][1].im = epsilon * PI * pow(pole,m) / a; } free_dvector(J, 1, m_max+1); /*_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/*/ /*_/ n=2 _/*/ /*_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/*/ if(n_max>=2) { J = dvector(0, m_max); J[0] = 1/(xmin-pole) - 1/(xmax-pole); if(m_max>=1) J[1] = log( (xmax-pole)/(pole-xmin) ); if(m_max>=2) J[2] = xmax - xmin; for(n=3; n<=m_max; n++) { J[n] = (pow(xmax-pole, n-1) - pow(xmin-pole, n-1)) / (n-1); } res[0][2].re = J[0]/pow(a,2); for(m=1; m<=m_max; m++) { for(i=m; i>=0; i--) { pascal = fact[m]/(fact[i]*fact[m-i]); res[m][2].re = res[m][2].re + J[i] * pascal * pow(pole, m-i); } } for (m=1; m<=m_max; m++) { res[m][2].re = res[m][2].re / pow(a,2); res[m][2].im = epsilon * PI * m * pow(pole,m-1) / pow(a,2.) ; } free_dvector(J, 0, m_max); } /* End n_max>=2 */ free_dvector(fact, 0, m_max); } /* The function int_dmn_cp_anal() calculates analytically the real and imaginary parts of the integral over x [0,1]: x**m/(ax**2+bx+c)**n (n=1 or 2) NB: THE CASE OF TWO POLES IN THE DOMAIN OF INTEGRATION WITH CONJUGATE IMAGINARY PARTS IS TREATED!!*/ void int_dmn_cp_anal(double xmin, double xmax, double a, double b, double c, dcomplex **res, int m_max, int n_max, double epsilon) { int i, m, n, m_cmn, n_cmn; int *mmax_cmn; double discri, delta, r1, r2, signe_a=1.; double **cmn; dcomplex **cmn_cp; /*----------------------------------------------------*/ if(a<0) signe_a = -1.; discri= pow(b,2) - 4*a*c; delta = discri / pow(a,2); r1 = (-b + sqrt(discri)) / (2*a); r2 = (-b - sqrt(discri)) / (2*a); m_cmn = m_max; n_cmn = 2; /*----------------------------------------------------*/ mmax_cmn = ivector(0, n_cmn); mmax_cmn[n_cmn] = m_cmn; for (i=n_cmn-1; i>=0; i--) mmax_cmn[i] = m_cmn; cmn = dmatrix(0, m_cmn, 0, n_cmn); cmn_cp = dcomp_matrix(0, m_cmn, 0, n_cmn); if(r1>xmin && r1=2) res[m][2] = Csub( CRmul(1/delta,cmn_cp[m][2]), CRmul(2*pow(delta,-1.5),cmn_cp[m][1]) ); } if(r2>xmin && r2=2) res[m][2] = Cadd( res[m][2], Cadd( CRmul(1/delta,cmn_cp[m][2]), CRmul(2*pow(delta,-1.5),cmn_cp[m][1]) ) ); } for (m=0; m<=m_max; m++) { res[m][1].re = res[m][1].re / a; res[m][1].im = epsilon * res[m][1].im / a; if(n_max>=2) { res[m][2].re = res[m][2].re / pow(a,2.); res[m][2].im = epsilon * res[m][2].im / pow(a,2.); } } free_dcomp_matrix(cmn_cp, 0, m_cmn, 0, n_cmn); free_dmatrix(cmn, 0, m_cmn, 0, n_cmn); }