17 #define NO_ATOMS(ND) (gv.bin[ND]->AvVol*gv.bin[ND]->dustp[0]/ATOMIC_MASS_UNIT/gv.bin[ND]->atomWeight)
93 static const double RELAX = 15.;
121 static const double tlim[5]={0.,50.,150.,500.,DBL_MAX};
122 static const double ppower[4]={2.00,1.30,0.68,0.00};
134 double[],
double[],
double[],
139 double[],
double,
double*,
long,
long,
bool*);
150 double[],
double*,
long*,
178 double *qtemp, *qprob;
181 bool lgNoTdustFailures;
188 const double LIM1 = pow(2.e-6,1./2.);
189 const double LIM2 = pow(6.e-6,1./3.);
196 x = log(0.999*DBL_MAX);
207 qtemp = (
double*)
MALLOC((
size_t)(
NQGRID*
sizeof(double)));
208 qprob = (
double*)
MALLOC((
size_t)(
NQGRID*
sizeof(double)));
210 for( nd=0; nd <
gv.
nBin; nd++ )
229 fprintf(
ioQQQ,
" GrainMakeDiffuse called with unknown storage class: %d\n", scase );
243 qheat(qtemp,qprob,&qnbin,nd);
267 for( j=qnbin-1; j >= 0 && xx > flux*DBL_EPSILON; j-- )
277 bfac = (1. + 0.5*f)*f;
299 bfac = (1. + 0.5*f)*f;
375 lgNoTdustFailures =
true;
376 for( nd=0; nd <
gv.
nBin; nd++ )
380 lgNoTdustFailures =
false;
392 for( nd=0; nd <
gv.
nBin; nd++ )
407 for( nd=0; nd <
gv.
nBin; nd++ )
423 for( nd=0; nd <
gv.
nBin; nd++ )
438 fabs(Comparison1-Comparison2)/Comparison2 <
CONSERV_TOL );
516 dPdlnT = (
double*)
MALLOC((
size_t)(
NQGRID*
sizeof(double)) );
541 PhiDrv[i] = -0.5*(y + oldy);
557 lgNegRate = lgNegRate || ( phiTilde[i] < 0. );
570 sprintf(fnam,
"Lambda_%2.2ld.asc",nd);
571 file = fopen(fnam,
"w");
572 for( i=0; i <
NDEMS; ++i )
573 fprintf(file,
"%e %e %e\n",
579 sprintf(fnam,
"Phi_%2.2ld.asc",nd);
580 file = fopen(fnam,
"w");
582 fprintf(file,
"%e %e\n",
rfield.
anu[i],Phi[i]);
590 Umax =
ufunct(Tmax,nd,&lgBoundErr);
596 deriv = y*c0/(
uderiv(Tmax,nd)*Tmax);
598 fwhm = sqrt(8.*
LN_TWO*c1/deriv);
601 FwhmRatio = fwhm/Umax;
620 ErrorCode = ErrorStart;
628 fprintf(
ioQQQ,
" grain heating: %.4e, integral %.4e, total rate %.4e lgNegRate %c\n",
630 fprintf(
ioQQQ,
" av grain temp %.4e av grain enthalpy (Ryd) %.4e\n",
632 fprintf(
ioQQQ,
" fwhm^2/(4ln2*c2/c0): %.4e fwhm (Ryd) %.4e fwhm/Umax %.4e\n",
633 NumEvents,fwhm,FwhmRatio );
634 fprintf(
ioQQQ,
" HeatingRate1 %.4e HeatingRate2 %.4e lgQHTooWide %c\n",
668 Ulo =
MAX2(Ulo,MinEnth);
673 gv.
bin[nd]->
qtmin = sqrt(minBracket*maxBracket);
679 ErrorCode = ErrorStart;
684 GetProbDistr_LowLimit(nd,rel_tol,Umax,fwhm,Phi,PhiDrv,qtemp,qprob,dPdlnT,qnbin,
685 &new_tmin,&nWideFail,&ErrorCode);
701 GetProbDistr_HighLimit(nd,
STRICT,Umax,fwhm,qtemp,qprob,dPdlnT,&tol,qnbin,&dummy,
714 if( new_tmin < minBracket || new_tmin > maxBracket )
719 if( new_tmin <= minBracket )
720 new_tmin = sqrt(
gv.
bin[nd]->
qtmin*minBracket);
721 if( new_tmin >= maxBracket )
722 new_tmin = sqrt(
gv.
bin[nd]->
qtmin*maxBracket);
736 double help1 =
gv.
bin[nd]->
qtmin*sqrt(DefFac);
737 double help2 = sqrt(
gv.
bin[nd]->
qtmin*maxBracket);
739 new_tmin =
MIN2(help1,help2);
745 double help = sqrt(
gv.
bin[nd]->
qtmin*minBracket);
751 new_tmin = sqrt(minBracket*maxBracket);
756 GetProbDistr_HighLimit(nd,
STRICT,Umax,fwhm,qtemp,qprob,dPdlnT,&tol,qnbin,&new_tmin,
772 fprintf(
ioQQQ,
" GetProbDistr returns code %d\n", ErrorCode );
775 fprintf(
ioQQQ,
" >>>>qheat trying another iteration, qtmin bracket %.4e %.4e",
776 minBracket,maxBracket );
777 fprintf(
ioQQQ,
" nWideFail %ld\n", nWideFail );
805 if( !lgOK && !lgDelta )
807 double Umax2 = Umax*sqrt(tol);
808 double fwhm2 = fwhm*sqrt(tol);
815 GetProbDistr_HighLimit(nd,
RELAX,Umax2,fwhm2,qtemp,qprob,dPdlnT,&tol,qnbin,&dummy,
822 ErrorCode = ErrorCode2;
842 fprintf(
ioQQQ,
" >>>>qheat converged with code: %d\n", ErrorCode );
848 fprintf(
ioQQQ,
" PROBLEM qheat did not converge grain %s in zone %ld, error code = %d\n",
854 fprintf(
gv.
QHPunchFile,
"\nDust Temperature Distribution: grain %s zone %ld\n",
864 for( i=0; i < *qnbin; i++ )
866 fprintf(
gv.
QHPunchFile,
"%.4e %.4e\n", qtemp[i],dPdlnT[i] );
871 fprintf(
gv.
QHPunchFile,
"**** quantum heating was not used\n" );
893 enum {DEBUG_LOC=
false};
920 double NegHeatingRate = 0.;
956 double ehat = gptr->
ehat[i];
957 double cool1,
sign = 1.;
960 if( gptr->
DustZ <= -1 )
991 double integral = 0.;
996 dprintf(
ioQQQ,
" integral test 1: integral %.6e %.6e\n", integral, sum );
1014 double *RateArr,Sum,ESum,DSum,Delta,E_av2,Corr;
1033 double ylo = -exp(-E0/fac);
1038 double yhi = ((E0-
Ehi)/fac-1.)*exp(-Ehi/fac);
1047 Sum = ESum = DSum = 0.;
1056 yhi = ((E0-
Ehi)/fac-1.)*exp(-Ehi/fac);
1058 RateArr[i] = rate*
MAX2(yhi-ylo,0.);
1073 ASSERT( fabs(E_av-E_av2) <= Delta );
1078 phiTilde[i] += RateArr[i]*Corr;
1085 double integral = 0.;
1090 dprintf(
ioQQQ,
" integral test 2: integral %.6e %.6e\n", integral, sum );
1123 const double LIM2 = pow(3.e-6,1./3.);
1124 const double LIM3 = pow(8.e-6,1./4.);
1140 double yhi = -(Ehi/fac+1.)*exp(-Ehi/fac);
1153 yhi = -(x+1.)*exp(-x);
1155 yhi = -(((1./3.)*x - 0.5)*x*x + 1.);
1157 yhi = -(1. - 0.5*x*x);
1160 phiTilde[i] += rate*
MAX2(yhi-ylo,0.);
1168 double integral = 0.;
1173 dprintf(
ioQQQ,
" integral test 3: integral %.6e %.6e\n", integral, sum );
1184 if( NegHeatingRate < 0. )
1186 double scale_fac = (sum+NegHeatingRate)/sum;
1188 phiTilde[i] *= scale_fac;
1190 sum += NegHeatingRate;
1194 double integral = 0.;
1199 dprintf(
ioQQQ,
" integral test 4: integral %.6e %.6e\n", integral, sum );
1269 fprintf(
ioQQQ,
" GetProbDistr_LowLimit called for grain %s\n",
gv.
bin[nd]->
chDstLab );
1276 u1[0] =
ufunct(qtemp[0],nd,&lgBoundErr);
1287 delu[0] = 2.*Lambda[0]/Phi[0];
1289 dPdlnT[0] = p[0]*qtemp[0]*
uderiv(qtemp[0],nd);
1290 RadCooling = 0.5*p[0]*Lambda[0]*delu[0];
1291 NextStep = 0.01*Lambda[0]/Phi[0];
1293 if( NextStep < 0.05*DBL_EPSILON*u1[0] )
1298 else if( NextStep < 5.*DBL_EPSILON*u1[0] )
1309 lgAllNegSlope =
true;
1318 dlnpdlnU = u1[0]*Phi[0]/Lambda[0] - dlnLdlnT*u1[0]/(qtemp[0]*
uderiv(qtemp[0],nd));
1322 (void)
TryDoubleStep(u1,delu,p,qtemp,Lambda,Phi,PhiDrv,NextStep,&dCool,k,nd,&lgBoundErr);
1323 dPdlnT[2] = p[2]*qtemp[2]*
uderiv(qtemp[2],nd);
1325 if( dPdlnT[2] < dPdlnT[0] ) {
1336 double RelError =
TryDoubleStep(u1,delu,p,qtemp,Lambda,Phi,PhiDrv,NextStep,&dCool,k,nd,&lgBoundErr);
1353 NextStep *= sqrt(0.9*rel_tol*QHEAT_TOL/RelError);
1356 if( NextStep < 5.*DBL_EPSILON*u1[k] )
1370 NextStep *=
MIN2(pow(0.9*rel_tol*QHEAT_TOL/
MAX2(RelError,1.e-50),1./3.),4.);
1371 NextStep =
MIN2(NextStep,Lambda[k]/Phi[0]);
1374 dPdlnT[k-1] = p[k-1]*qtemp[k-1]*
uderiv(qtemp[k-1],nd);
1375 dPdlnT[k] = p[k]*qtemp[k]*
uderiv(qtemp[k],nd);
1377 lgAllNegSlope = lgAllNegSlope && ( dPdlnT[k] < dPdlnT[k-2] );
1379 if( dPdlnT[k-1] > maxVal )
1381 maxVal = dPdlnT[k-1];
1384 if( dPdlnT[k] > maxVal )
1390 RadCooling += dCool;
1400 if( p[k] > sqrt(DBL_MAX/100.) )
1408 for( j=0; j <= k; j++ )
1410 p[j] /= DBL_MAX/100.;
1411 dPdlnT[j] /= DBL_MAX/100.;
1413 RadCooling /= DBL_MAX/100.;
1419 ASSERT( p[k] > 0. && dPdlnT[k] > 0. && RadCooling > 0. );
1424 if( k > 0 && k%
NQTEST == 0 )
1426 double wid, avStep, factor;
1439 avStep = log(u1[k]/u1[0])/(double)k;
1443 factor = 1.1 + 3.9*(1.0 - sqrt((
double)k/(
double)
NQGRID));
1444 if( wid/avStep > factor*(
double)
NQGRID )
1481 ScanProbDistr(u1,dPdlnT,nbin,maxVal,nmax,&nstart,&nstart2,&nend,nWideFail,ErrorCode);
1488 for( j=0; j < nbin; j++ )
1497 for( j=nstart; j < nbin; j++ )
1501 *new_tmin = qtemp[j];
1511 if( dPdlnT[nstart2] < 0.999*dPdlnT[nstart2+
NSTARTUP] )
1516 double T1 = qtemp[nstart2];
1517 double T2 = qtemp[nstart2+
NSTARTUP];
1518 double delta_y = log(dPdlnT[nstart2+
NSTARTUP]/dPdlnT[nstart2]);
1519 double c2 = delta_y/(1./
POW3(T1)-1./
POW3(T2));
1520 double help = c2/
POW3(T1) + log(dPdlnT[nstart2]/PROB_CUTOFF_LO);
1521 *new_tmin = pow(c2/help,1./3.);
1525 if( dPdlnT[j] >
SAFETY*PROB_CUTOFF_LO && *new_tmin >= qtmin1 )
1527 double delta_x = log(qtemp[nstart2+
NSTARTUP]/qtemp[nstart2]);
1528 double delta_y = log(dPdlnT[nstart2+
NSTARTUP]/dPdlnT[nstart2]);
1529 delta_x *= log(PROB_CUTOFF_LO/dPdlnT[nstart2])/delta_y;
1530 *new_tmin = qtemp[nstart2]*exp(delta_x);
1531 if( *new_tmin < qtmin1 )
1533 *new_tmin = sqrt( *new_tmin*qtmin1 );
1580 nbin =
RebinQHeatResults(nd,nstart,nend,p,qtemp,qprob,dPdlnT,u1,delu,Lambda,ErrorCode);
1591 for( j=0; j < nbin; j++ )
1604 " zone %ld %s nbin %ld nok %ld nbad %ld total prob %.4e rel_tol %.3e new_tmin %.4e\n",
1651 ASSERT( index >= 0 && index < NQGRID-2 && nd >= 0 && nd < gv.nBin && step > 0. );
1660 for( i=0; i <= index; i++ )
1661 p_max =
MAX2(p_max,p[i]);
1666 for( i=1; i <= 2; i++ )
1673 u1[k] = u1[k-1] + delu[k];
1676 *lgBoundFail = *lgBoundFail || lgErr;
1680 trap1 = trap2 = trap12 = 0.;
1683 for( j=jlo; j < k; j++ )
1685 umin = u1[k] - u1[j];
1695 else if( umin > ulo )
1717 while( ipHi-ipLo > 1 )
1719 long ipMd = (ipLo+ipHi)/2;
1725 bval_jk = Phi[ipLo] + (umin -
rfield.
anu[ipLo])*PhiDrv[ipLo];
1743 trap2 = p[j]*bval_jk;
1745 sum += (trap1 + trap2)*delu[j];
1754 p[k] = (sum + trap1*delu[k])/(2.*Lambda[k] - Phi[0]*delu[k]);
1762 p2k = (sum2 + trap12*step)/(2.*Lambda[k] - Phi[0]*step);
1768 RelErrPk = fabs(p2k-p[k])/p[k];
1772 *cooling =
log_integral(u1[k-2],p[k-2]*Lambda[k-2],u1[k-1],p[k-1]*Lambda[k-1]);
1773 *cooling +=
log_integral(u1[k-1],p[k-1]*Lambda[k-1],u1[k],p[k]*Lambda[k]);
1776 cooling2 =
log_integral(u1[k-2],p[k-2]*Lambda[k-2],u1[k],p2k*Lambda[k]);
1779 RelErrCool = ( index > 0 ) ? fabs(cooling2-(*cooling))/(*cooling) : 0.;
1783 return MAX2(RelErrPk,RelErrCool)/3.;
1800 ASSERT( xx1 > 0. && yy1 > 0. && xx2 > 0. && yy2 > 0. );
1803 eps = log((xx2/xx1)*(yy2/yy1));
1804 if( fabs(eps) < 1.e-4 )
1806 result = xx1*yy1*xx*((eps/6. + 0.5)*eps + 1.);
1810 result = (xx2*yy2 - xx1*yy1)*xx/eps;
1833 const double MIN_FAC_LO = 1.e4;
1834 const double MIN_FAC_HI = 1.e4;
1839 ASSERT( nbin > 0 && nmax >= 0 && nmax < nbin && maxVal > 0. );
1849 for( i=nmax; i >= 0; --i )
1851 if( dPdlnT[i] < minVal )
1859 for( i=nmax; i > *nstart; --i )
1861 double deriv = log(dPdlnT[i]/dPdlnT[i-1])/log(u1[i]/u1[i-1]);
1862 if( deriv > deriv_max )
1871 lgSetLo = ( nmax >= *nend || maxVal/dPdlnT[*nend] < MIN_FAC_HI );
1877 lgSetHi = ( nmax <= *nstart || maxVal/dPdlnT[*nstart] < MIN_FAC_LO );
1879 lgSetHi = ( nmax <= *nstart2 || maxVal/dPdlnT[*nstart2] < MIN_FAC_LO );
1881 if( lgSetLo && lgSetHi )
1901 fprintf(
ioQQQ,
" ScanProbDistr nstart %ld nstart2 %ld nend %ld nmax %ld maxVal %.3e",
1902 *nstart,*nstart2,*nend,nmax,maxVal );
1903 fprintf(
ioQQQ,
" dPdlnT[nstart] %.3e dPdlnT[nstart2] %.3e dPdlnT[nend] %.3e code %d\n",
1904 dPdlnT[*nstart],dPdlnT[*nstart2],dPdlnT[*nend],*ErrorCode );
1956 ASSERT( nstart >= 0 && nstart < nend && nend <
NQGRID );
1969 new_delu = (
double*)
MALLOC((
size_t)(
NQGRID*
sizeof(double)));
1970 new_dPdlnT = (
double*)
MALLOC((
size_t)(
NQGRID*
sizeof(double)));
1971 new_Lambda = (
double*)
MALLOC((
size_t)(
NQGRID*
sizeof(double)));
1972 new_p = (
double*)
MALLOC((
size_t)(
NQGRID*
sizeof(double)));
1973 new_qprob = (
double*)
MALLOC((
size_t)(
NQGRID*
sizeof(double)));
1974 new_qtemp = (
double*)
MALLOC((
size_t)(
NQGRID*
sizeof(double)));
1975 new_u1 = (
double*)
MALLOC((
size_t)(
NQGRID*
sizeof(double)));
1984 help = pow(qtemp[nend]/qtemp[i],1./(1.5*
NQMIN));
2003 double p1,p2,L1,L2,frac,slope;
2004 if( qtemp[i] <= T1 && T1 <= qtemp[i+1] )
2008 frac = log(T1/qtemp[i]);
2009 slope = log(p[i+1]/p[i])/log(qtemp[i+1]/qtemp[i]);
2010 p1 = p[i]*exp(frac*slope);
2011 slope = log(Lambda[i+1]/Lambda[i])/log(qtemp[i+1]/qtemp[i]);
2012 L1 = Lambda[i]*exp(frac*slope);
2021 if( qtemp[i] <= T2 && T2 <= qtemp[i+1] )
2024 help =
ufunct(T2,nd,&lgBoundErr);
2025 uu2 =
MIN2(help,u1[i+1]);
2027 frac = log(T2/qtemp[i]);
2028 slope = log(p[i+1]/p[i])/log(qtemp[i+1]/qtemp[i]);
2029 p2 = p[i]*exp(frac*slope);
2030 slope = log(Lambda[i+1]/Lambda[i])/log(qtemp[i+1]/qtemp[i]);
2031 L2 = Lambda[i]*exp(frac*slope);
2056 }
while( i < nend && ! lgDone );
2067 new_qprob[newnbin] = s0;
2068 new_Lambda[newnbin] = s1/s0;
2072 new_qtemp[newnbin] = exp(y);
2073 new_u1[newnbin] =
ufunct(new_qtemp[newnbin],nd,&lgBoundErr);
2075 new_delu[newnbin] = wid;
2076 new_p[newnbin] = new_qprob[newnbin]/new_delu[newnbin];
2077 new_dPdlnT[newnbin] = new_p[newnbin]*new_qtemp[newnbin]*
uderiv(new_qtemp[newnbin],nd);
2080 RadCooling += new_qprob[newnbin]*new_Lambda[newnbin];
2088 if( newnbin <
NQMIN )
2106 fprintf(
ioQQQ,
" RebinQHeatResults found tol1 %.4e tol2 %.4e\n",
2107 fabs(u1[nend]/Ucheck-1.), fabs(fac-1.) );
2112 ASSERT( fabs(u1[nend]/Ucheck-1.) < 10.*sqrt((
double)(nend-nstart+newnbin))*DBL_EPSILON );
2117 for( i=0; i < newnbin; i++ )
2120 p[i] = new_p[i]/fac;
2121 qtemp[i] = new_qtemp[i];
2122 qprob[i] = new_qprob[i]/fac;
2123 dPdlnT[i] = new_dPdlnT[i]/fac;
2125 delu[i] = new_delu[i];
2126 Lambda[i] = new_Lambda[i];
2129 ASSERT( qtemp[i] > 0. && qprob[i] > 0. );
2194 fprintf(
ioQQQ,
" GetProbDistr_HighLimit called for grain %s\n",
gv.
bin[nd]->
chDstLab );
2202 help1 = Umax*exp(-fac1);
2204 Ulo =
MAX2(help1,help2);
2210 if( fac2 > log(DBL_MAX/10.) )
2215 Uhi = Umax*exp(fac2);
2221 uu1 =
ufunct(T1,nd,&lgErr);
2222 lgBoundErr = lgBoundErr || lgErr;
2223 help1 = log(uu1/Umax);
2224 p1 = c1*exp(-c2*
POW2(help1));
2226 lgBoundErr = lgBoundErr || lgErr;
2230 if( uu1*p1*L1 < 1.e5*DBL_MIN )
2237 help1 = pow(Thi/Tlo,1./(1.2*
NQMIN));
2248 uu2 =
ufunct(T2,nd,&lgErr);
2249 lgBoundErr = lgBoundErr || lgErr;
2250 help1 = log(uu2/Umax);
2251 p2 = c1*exp(-c2*
POW2(help1));
2253 lgBoundErr = lgBoundErr || lgErr;
2261 Lambda[nbin] = s1/s0;
2264 lgBoundErr = lgBoundErr || lgErr;
2265 qtemp[nbin] = exp(y);
2267 p[nbin] = qprob[nbin]/delu[nbin];
2268 dPdlnT[nbin] = p[nbin]*qtemp[nbin]*
uderiv(qtemp[nbin],nd);
2271 RadCooling += qprob[nbin]*Lambda[nbin];
2280 }
while( T2 < Thi && nbin <
NQGRID );
2284 for( i=0; i < nbin; ++i )
2292 *new_tmin = qtemp[0];
2310 fprintf(
ioQQQ,
" GetProbDistr_HighLimit found tol1 %.4e tol2 %.4e\n",
2311 fabs(sum-1.), fabs(sum/fac-1.) );
2312 fprintf(
ioQQQ,
" zone %ld %s nbin %ld total prob %.4e new_tmin %.4e\n",
2329 hok[3] = {1275., 1670., 4359.},
2341 fprintf(
ioQQQ,
" uderiv called with non-positive temperature: %.6e\n" , temp );
2350 numer = (4.15e-22/
EN1RYD)*pow(temp,3.3);
2351 dnumer = (3.3*4.15e-22/
EN1RYD)*pow(temp,2.3);
2352 denom = 1. + 6.51e-03*temp + 1.5e-06*temp*temp + 8.3e-07*pow(temp,2.3);
2353 ddenom = 6.51e-03 + 2.*1.5e-06*temp + 2.3*8.3e-07*pow(temp,1.3);
2356 deriv = (dnumer*denom - numer*ddenom)/
POW2(denom);
2367 for( j = 0; j < 4; j++ )
2369 if( temp >
tlim[j] && temp <=
tlim[j+1] )
2385 x = log10(
MIN2(temp,2000.));
2386 deriv = pow(10.,-21.26+3.1688*x-0.401894*
POW2(x))/
EN1RYD;
2398 else if( N_C <= 100. )
2400 N_H = 2.5*sqrt(N_C);
2407 for( i=0; i < 3; i++ )
2409 double help1 = hok[i]/temp;
2412 double help2 = exp(help1);
2413 double help3 = ( help1 < 1.e-7 ) ? help1*(1.+0.5*help1) : help2-1.;
2420 fprintf(
ioQQQ,
" uderiv called with unknown type for enthalpy function: %d\n", ecase );
2430 fprintf(
ioQQQ,
" uderiv finds non-positive derivative: %.6e, what's up?\n" , deriv );
2449 fprintf(
ioQQQ,
" ufunct called with non-positive temperature: %.6e\n" , temp );
2473 if( enthalpy <= 0. )
2475 fprintf(
ioQQQ,
" inv_ufunct called with non-positive enthalpy: %.6e\n" , enthalpy );
2502 for( nd=0; nd <
gv.
nBin; nd++ )
2505 C_V2 =
uderiv(tdust2,nd);
2511 for( i=1; i <
NDEMS; i++ )
2515 C_V2 =
uderiv(tdust2,nd);
2516 tmid = sqrt(tdust1*tdust2);
2518 for( j=1; j < 4; j++ )
2520 if( tdust1 <
tlim[j] &&
tlim[j] < tdust2 )
2536 for( nd=0; nd <
gv.
nBin; nd++ )
2538 for( i=0; i <
NDEMS; i++ )
2545 for( nd=0; nd <
gv.
nBin; nd++ )
2570 ASSERT( n == 2 || n == 3 );
2577 res = 6.*1.202056903159594*
POW2(x);
2581 res = 24.*1.082323233711138*
POW3(x);
2590 nn = 4*
MAX2(4,2*(
long)(0.05/x));
2591 xx = (
double *)
MALLOC(
sizeof(
double)*(unsigned)nn);
2592 rr = (
double *)
MALLOC(
sizeof(
double)*(unsigned)nn);
2593 aa = (
double *)
MALLOC(
sizeof(
double)*(unsigned)nn);
2594 ww = (
double *)
MALLOC(
sizeof(
double)*(unsigned)nn);
2599 for( i=0; i < nn; i++ )
2601 double help1 = rr[i]/x;
2604 double help2 = exp(help1);
2605 double help3 = ( help1 < 1.e-7 ) ? help1*(1.+0.5*help1) : help2-1.;
2606 res += ww[i]*
powi(rr[i],n+1)*help2/
POW2(help3);
2616 return (
double)n*res;