00001
00002
00003 #include "cddefines.h"
00004 #include "physconst.h"
00005 #include "thirdparty.h"
00006 #include "continuum.h"
00007
00008 STATIC double RealF2_1( double alpha, double beta, double gamma, double chi );
00009 STATIC complex<double> Hypergeometric2F1( complex<double> a, complex<double> b, complex<double> c,
00010 double chi, long *NumRenorms, long *NumTerms );
00011 STATIC complex<double> F2_1( complex<double> alpha, complex<double> beta, complex<double> gamma,
00012 double chi, long *NumRenormalizations, long *NumTerms );
00013 STATIC complex<double> HyperGeoInt( double v );
00014 STATIC complex<double> qg32complex( double xl, double xu, complex<double> (*fct)(double) );
00015 STATIC double GauntIntegrand( double y );
00016 STATIC double FreeFreeGaunt( double x );
00017 STATIC double DoBeckert_etal( double etai, double etaf, double chi );
00018 STATIC double DoSutherland( double etai, double etaf, double chi );
00019
00020
00021 static const complex<double> Normalization(1e100, 1e100);
00022 static complex<double> CMinusBMinus1, BMinus1, MinusA;
00023 static double GlobalCHI;
00024 static double Zglobal, HNUglobal, TEglobal;
00025
00026 double cont_gaunt_calc( double temp, double z, double photon )
00027 {
00028 double gaunt, u, gamma2;
00029
00030 Zglobal = z;
00031 HNUglobal = photon;
00032 TEglobal = temp;
00033
00034 u = TE1RYD*photon/temp;
00035 gamma2 = TE1RYD*z*z/temp;
00036
00037 if( log10(u)<-5. )
00038 {
00039
00040 if( log10( gamma2 ) < -0.75187 )
00041 {
00042
00043
00044 gaunt = 0.551329 * ( 0.80888 - log(u) );
00045 }
00046 else
00047 {
00048
00049
00050 gaunt = -0.551329 * (0.5*log(gamma2) + log(u) + 0.056745);
00051 }
00052 }
00053 else
00054 {
00055
00056 gaunt = qg32( 0.01, 1., GauntIntegrand );
00057 gaunt += qg32( 1., 5., GauntIntegrand );
00058 }
00059 ASSERT( gaunt>0. && gaunt<100. );
00060
00061 return gaunt;
00062 }
00063
00064 STATIC double GauntIntegrand( double y )
00065 {
00066 double value;
00067 value = FreeFreeGaunt( y ) * exp(-y);
00068 return value;
00069 }
00070
00071 STATIC double FreeFreeGaunt( double x )
00072 {
00073 double Csum, zeta, etaf, etai, chi, gaunt, z, InitialElectronEnergy, FinalElectronEnergy, photon;
00074 bool lgSutherlandOn = false;
00075 long i;
00076
00077 z = Zglobal;
00078 photon = HNUglobal;
00079 ASSERT( z > 0. );
00080 ASSERT( photon > 0. );
00081
00082
00083 InitialElectronEnergy = sqrt(x) * TEglobal/TE1RYD;
00084 FinalElectronEnergy = photon + InitialElectronEnergy;
00085 ASSERT( InitialElectronEnergy > 0. );
00086
00087
00088 etai = z/sqrt(InitialElectronEnergy);
00089 etaf = z/sqrt(FinalElectronEnergy);
00090 ASSERT( etai > 0. );
00091 ASSERT( etaf > 0. );
00092 chi = -4. * etai * etaf / POW2( etai - etaf );
00093 zeta = etai-etaf;
00094
00095 if( etai>=130.)
00096 {
00097 if( etaf < 1.7 )
00098 {
00099
00100 gaunt = 1.1027 * (1.-exp(-2.*PI*etaf));
00101 }
00102 else if( etaf < 0.1*etai )
00103 {
00104
00105 gaunt = 1. + 0.17282604*pow(etaf,-0.67) - 0.04959570*pow(etaf,-1.33)
00106 - 0.01714286*pow(etaf,-2.) + 0.00204498*pow(etaf,-2.67)
00107 - 0.00243945*pow(etaf,-3.33) - 0.00120387*pow(etaf,-4.)
00108 + 0.00071814*pow(etaf,-4.67) + 0.00026971*pow(etaf,-5.33);
00109 }
00110 else if( zeta > 0.5 )
00111 {
00112
00113 gaunt = 1. + 0.21775*pow(zeta,-0.67) - 0.01312*pow(zeta, -1.33);
00114 }
00115 else
00116 {
00117 double a[10] = {1.47864486, -1.72329012, 0.14420320, 0.05744888, 0.01668957,
00118 0.01580779, 0.00464268, 0.00385156, 0.00116196, 0.00101967};
00119
00120 Csum = 0.;
00121 for( i = 0; i <=9; i++ )
00122 {
00123
00124 Csum += a[i]*RealF2_1( (double)(-i), (double)i, 0.5, 0.5*(1.-zeta) );
00125 }
00126 gaunt = fabs(0.551329*(0.57721 + log(zeta/2.))*exp(PI*zeta)*Csum);
00127 ASSERT( gaunt < 10. );
00128 }
00129 }
00130 else if( lgSutherlandOn )
00131 gaunt = DoSutherland( etai, etaf, chi );
00132 else
00133 gaunt = DoBeckert_etal( etai, etaf, chi );
00134
00135
00136
00137
00138
00139
00142 ASSERT( gaunt > 0. && gaunt<BIGFLOAT );
00143
00144 if( gaunt == 0. )
00145 {
00146 fprintf( ioQQQ, "Uh-Oh! Gaunt is zero! Is this okay?\n");
00147
00148 gaunt = 1e-5;
00149 }
00150
00151 return gaunt;
00152 }
00153
00154
00155 #if defined(__ICC) && defined(__i386) && __INTEL_COMPILER < 910
00156 #pragma optimize("", off)
00157 #endif
00158
00159
00161 STATIC double DoBeckert_etal( double etai, double etaf, double chi )
00162 {
00163 double Delta, BeckertGaunt, MaxFReal, LnBeckertGaunt;
00164 long NumRenorms[2]={0,0}, NumTerms[2]={0,0};
00165 int IndexMinNumRenorms, IndexMaxNumRenorms;
00166 complex<double> a,b,c,F[2];
00167
00168 a = complex<double>( 1., -etai );
00169 b = complex<double>( 0., -etaf );
00170 c = 1.;
00171
00172
00173 F[0] = Hypergeometric2F1( a, b, c, chi, &NumRenorms[0], &NumTerms[0] );
00174
00175 a = complex<double>( 1., -etaf );
00176 b = complex<double>( 0., -etai );
00177
00178
00179 F[1] = Hypergeometric2F1( a, b, c, chi, &NumRenorms[1], &NumTerms[1] );
00180
00181
00182
00183
00184
00185 if( ( MAX2(NumTerms[1],NumTerms[0]) - MIN2(NumTerms[1],NumTerms[0]) >= 2 )
00186 && NumTerms[1]!=-1 && NumTerms[0]!=-1)
00187 {
00188 a = complex<double>( 1., -etai );
00189 b = complex<double>( 0., -etaf );
00190 c = 1.;
00191
00192 NumTerms[0] = MAX2(NumTerms[1],NumTerms[0])+1;
00193 NumTerms[1] = NumTerms[0];
00194 NumRenorms[0] = 0;
00195 NumRenorms[1] = 0;
00196
00197
00198 F[0] = Hypergeometric2F1( a, b, c, chi, &NumRenorms[0], &NumTerms[0] );
00199
00200 a = complex<double>( 1., -etaf );
00201 b = complex<double>( 0., -etai );
00202
00203
00204 F[1] = Hypergeometric2F1( a, b, c, chi, &NumRenorms[1], &NumTerms[1] );
00205
00206 ASSERT( NumTerms[0] == NumTerms[1] );
00207 }
00208
00209
00210 if( log10(abs(F[0])/abs(F[1])) + (NumRenorms[0]-NumRenorms[1])*log10(abs(Normalization)) > 10. )
00211 {
00212 F[1] = 0.;
00213
00214 NumRenorms[1] = NumRenorms[0];
00215 }
00216 else if( log10(abs(F[1])/abs(F[0])) + (NumRenorms[1]-NumRenorms[0])*log10(abs(Normalization)) > 10. )
00217 {
00218 F[0] = 0.;
00219
00220 NumRenorms[0] = NumRenorms[1];
00221 }
00222
00223
00224
00225
00226 MaxFReal = (fabs(F[1].real())>fabs(F[0].real())) ? fabs(F[1].real()):fabs(F[0].real());
00227 while( NumRenorms[0] != NumRenorms[1] )
00228 {
00229
00230 if( MaxFReal > 1e50 )
00231 {
00232 IndexMinNumRenorms = ( NumRenorms[0] > NumRenorms[1] ) ? 1:0;
00233 F[IndexMinNumRenorms] /= Normalization;
00234 ++NumRenorms[IndexMinNumRenorms];
00235 }
00236 else
00237 {
00238 IndexMaxNumRenorms = ( NumRenorms[0] > NumRenorms[1] ) ? 0:1;
00239 F[IndexMaxNumRenorms] = F[IndexMaxNumRenorms]*Normalization;
00240 --NumRenorms[IndexMaxNumRenorms];
00241 }
00242 }
00243
00244 ASSERT( NumRenorms[0] == NumRenorms[1] );
00245
00246
00247
00248
00249 ASSERT( (fabs(F[0].real())<1e+150) && (fabs(F[1].real())<1e+150) &&
00250 (fabs(F[0].imag())<1e+150) && (fabs(F[1].real())<1e+150) );
00251 ASSERT( (fabs(F[0].real())>1e-150) && ((fabs(F[0].imag())>1e-150) || (abs(F[0])==0.)) );
00252 ASSERT( (fabs(F[1].real())>1e-150) && ((fabs(F[1].real())>1e-150) || (abs(F[1])==0.)) );
00253
00254
00255 complex<double> CDelta = F[0]*F[0] - F[1]*F[1];
00256 double renorm = MAX2(fabs(CDelta.real()),fabs(CDelta.imag()));
00257 ASSERT( renorm > 0. );
00258
00259 complex<double> NCDelta( CDelta.real()/renorm, CDelta.imag()/renorm );
00260
00261 Delta = renorm * abs( NCDelta );
00262
00263 ASSERT( Delta > 0. );
00264
00265
00266 if( etaf > 100. )
00267 {
00268
00269 LnBeckertGaunt = 1.6940360 + log(Delta) + log(etaf) + log(etai) - log(fabs(etai-etaf)) - 6.2831853*etaf;
00270 LnBeckertGaunt += 2. * NumRenorms[0] * log(abs(Normalization));
00271 BeckertGaunt = exp( LnBeckertGaunt );
00272 NumRenorms[0] = 0;
00273 }
00274 else
00275 {
00276 BeckertGaunt = Delta*5.4413981*etaf*etai/fabs(etai - etaf)
00277 /(1.-exp(-6.2831853*etai) )/( exp(6.2831853*etaf) - 1.);
00278
00279 while( NumRenorms[0] > 0 )
00280 {
00281 BeckertGaunt *= abs(Normalization);
00282 BeckertGaunt *= abs(Normalization);
00283 ASSERT( BeckertGaunt < BIGDOUBLE );
00284 --NumRenorms[0];
00285 }
00286 }
00287
00288 ASSERT( NumRenorms[0] == 0 );
00289
00290
00291
00292
00293 return BeckertGaunt;
00294 }
00295 #if defined(__ICC) && defined(__i386) && __INTEL_COMPILER < 910
00296 #pragma optimize("", on)
00297 #endif
00298
00299
00300
00301
00303 STATIC double DoSutherland( double etai, double etaf, double chi )
00304 {
00305 double Sgaunt, ICoef, weightI1, weightI0;
00306 long i, NumRenorms[2]={0,0}, NumTerms[2]={0,0};
00307 complex<double> a,b,c,GCoef,kfac,etasum,G[2],I[2],ComplexFactors,GammaProduct;
00308
00309 kfac = complex<double>( fabs((etaf-etai)/(etaf+etai)), 0. );
00310 etasum = complex<double>( 0., etai + etaf );
00311
00312 GCoef = pow(kfac, etasum);
00313
00314
00315 ASSERT( fabs(GCoef.real())<1.0 && fabs(GCoef.imag())<1.0 && ( GCoef.real()!=0. || GCoef.imag()!=0. ) );
00316
00317 for( i = 0; i <= 1; i++ )
00318 {
00319 a = complex<double>( i + 1., -etaf );
00320 b = complex<double>( i + 1., -etai );
00321 c = complex<double>( 2.*i + 2., 0. );
00322
00323
00324 G[i] = Hypergeometric2F1( a, b, c, chi, &NumRenorms[i], &NumTerms[i] );
00325 }
00326
00327
00328
00329
00330
00331 if( MAX2(NumTerms[1],NumTerms[0]) - MIN2(NumTerms[1],NumTerms[0]) > 2
00332 && NumTerms[1]!=-1 && NumTerms[0]!=-1 )
00333 {
00334 NumTerms[0] = MAX2(NumTerms[1],NumTerms[0]);
00335 NumTerms[1] = NumTerms[0];
00336 NumRenorms[0] = 0;
00337 NumRenorms[1] = 0;
00338
00339 for( i = 0; i <= 1; i++ )
00340 {
00341 a = complex<double>( i + 1., -etaf );
00342 b = complex<double>( i + 1., -etai );
00343 c = complex<double>( 2.*i + 2., 0. );
00344
00345 G[i] = Hypergeometric2F1( a, b, c, chi, &NumRenorms[i], &NumTerms[i] );
00346 }
00347
00348 ASSERT( NumTerms[0] == NumTerms[1] );
00349 }
00350
00351 for( i = 0; i <= 1; i++ )
00352 {
00354 ASSERT( fabs(G[i].real())>0. && fabs(G[i].real())<1e100 &&
00355 fabs(G[i].imag())>0. && fabs(G[i].imag())<1e100 );
00356
00357
00358 G[i] *= GCoef;
00359
00360
00361
00362 ICoef = 0.25*pow(-chi, (double)i+1.)*exp( 1.5708*fabs(etai-etaf) )/factorial(2*i);
00363 GammaProduct = cdgamma(complex<double>(i+1.,etai))*cdgamma(complex<double>(i+1.,etaf));
00364 ICoef *= abs(GammaProduct);
00365
00366 ASSERT( ICoef > 0. );
00367
00368 I[i] = ICoef*G[i];
00369
00370 while( NumRenorms[i] > 0 )
00371 {
00372 I[i] *= Normalization;
00373 ASSERT( fabs(I[i].real()) < BIGDOUBLE && fabs(I[i].imag()) < BIGDOUBLE );
00374 --NumRenorms[i];
00375 }
00376
00377 ASSERT( NumRenorms[i] == 0 );
00378 }
00379
00380 weightI0 = POW2(etaf+etai);
00381 weightI1 = 2.*etaf*etai*sqrt(1. + etai*etai)*sqrt(1. + etaf*etaf);
00382
00383 ComplexFactors = I[0] * ( weightI0*I[0] - weightI1*I[1] );
00384
00385
00386 Sgaunt = 1.10266 / etai / etaf * abs( ComplexFactors );
00387
00388 return Sgaunt;
00389 }
00390
00391 #if defined(__ICC) && defined(__i386) && __INTEL_COMPILER < 910
00392 #pragma optimize("", off)
00393 #endif
00394
00395 STATIC complex<double> Hypergeometric2F1( complex<double> a, complex<double> b, complex<double> c,
00396 double chi, long *NumRenorms, long *NumTerms )
00397 {
00398 complex<double> a1, b1, c1, a2, b2, c2, Result, Part[2], F[2];
00399 complex<double> chifac, GammaProduct, Coef, FIntegral;
00401 double Interface1 = 0.4, Interface2 = 10.;
00402 long N_Renorms[2], N_Terms[2], IndexMaxNumRenorms, lgDoIntegral = false;
00403
00404 N_Renorms[0] = *NumRenorms;
00405 N_Renorms[1] = *NumRenorms;
00406 N_Terms[0] = *NumTerms;
00407 N_Terms[1] = *NumTerms;
00408
00409
00410 ASSERT( chi < 0. );
00411
00412
00413
00414
00415
00416 if( fabs(chi) < Interface1 )
00417 {
00418 Result = F2_1( a, b, c, chi, &*NumRenorms, &*NumTerms );
00419 }
00420
00421 else if( fabs(chi) > Interface2 )
00422 {
00423 a1 = a;
00424 b1 = 1.-c+a;
00425 c1 = 1.-b+a;
00426
00427 a2 = b;
00428 b2 = 1.-c+b;
00429 c2 = 1.-a+b;
00430
00431 chifac = -chi;
00432
00433 F[0] = F2_1(a1,b1,c1,1./chi,&N_Renorms[0], &N_Terms[0]);
00434 F[1] = F2_1(a2,b2,c2,1./chi,&N_Renorms[1], &N_Terms[1]);
00435
00436
00437 if( MAX2(N_Terms[1],N_Terms[0]) - MIN2(N_Terms[1],N_Terms[0]) >= 2 )
00438 {
00439 N_Terms[0] = MAX2(N_Terms[1],N_Terms[0]);
00440 N_Terms[1] = N_Terms[0];
00441 N_Renorms[0] = *NumRenorms;
00442 N_Renorms[1] = *NumRenorms;
00443
00444 F[0] = F2_1(a1,b1,c1,1./chi,&N_Renorms[0], &N_Terms[0]);
00445 F[1] = F2_1(a2,b2,c2,1./chi,&N_Renorms[1], &N_Terms[1]);
00446 ASSERT( N_Terms[0] == N_Terms[1] );
00447 }
00448
00449 *NumTerms = MAX2(N_Terms[1],N_Terms[0]);
00450
00451
00452
00453 GammaProduct = (cdgamma(b-a)/cdgamma(b))*(cdgamma(c)/cdgamma(c-a));
00454
00455
00456 Part[0] = F[0]/pow(chifac,a)*GammaProduct;
00457
00458
00459
00460 GammaProduct = (cdgamma(a-b)/cdgamma(a))*(cdgamma(c)/cdgamma(c-b));
00461
00462
00463 Part[1] = F[1]/pow(chifac,b)*GammaProduct;
00464
00465
00466
00467
00468
00469 if( N_Renorms[0] != N_Renorms[1] )
00470 {
00471 IndexMaxNumRenorms = ( N_Renorms[0] > N_Renorms[1] ) ? 0:1;
00472 Part[IndexMaxNumRenorms] *= Normalization;
00473 --N_Renorms[IndexMaxNumRenorms];
00474
00475
00476 ASSERT( N_Renorms[0] == N_Renorms[1] );
00477 }
00478
00479 *NumRenorms = N_Renorms[0];
00480
00481 Result = Part[0] + Part[1];
00482 }
00483
00484 else
00485 {
00486
00487 if( lgDoIntegral )
00488 {
00489
00490
00491 if( abs(b) > abs(a) )
00492 {
00493 complex<double> btemp = b;
00494 b = a;
00495 a = btemp;
00496 }
00497 Coef = cdgamma(c)/(cdgamma(b)*cdgamma(c-b));
00498 CMinusBMinus1 = c-b-1.;
00499 BMinus1 = b-1.;
00500 MinusA = -a;
00501 GlobalCHI = chi;
00502 FIntegral = qg32complex( 0., 0.5, HyperGeoInt );
00503 FIntegral += qg32complex( 0.5, 1., HyperGeoInt );
00504
00505 Result = Coef + FIntegral;
00506 *NumTerms = -1;
00507 *NumRenorms = 0;
00508 }
00509 else
00510 {
00511
00512 a1 = a;
00513 b1 = c-b;
00514 c1 = c;
00515 chifac = 1.-chi;
00516
00517 Result = F2_1(a1,b1,c1,chi/(chi-1.),&*NumRenorms,&*NumTerms)/pow(chifac,a);
00518 }
00519 }
00520
00521
00522 while( fabs(Result.real()) >= 1e50 )
00523 {
00524 Result /= Normalization;
00525 ++*NumRenorms;
00526 }
00527
00528 return Result;
00529 }
00530
00531
00532 STATIC complex<double> F2_1(
00533 complex<double> alpha, complex<double> beta, complex<double> gamma,
00534 double chi, long *NumRenormalizations, long *NumTerms )
00535 {
00536 long i = 3, MinTerms;
00537 bool lgNotConverged = true;
00538 complex<double> LastTerm, Term, Sum;
00539
00540 MinTerms = MAX2( 3, *NumTerms );
00541
00542
00543 Sum = 1./Normalization;
00544 ++*NumRenormalizations;
00545
00546
00547 LastTerm = Sum*alpha*beta/gamma*chi;
00548
00549 Sum += LastTerm;
00550
00551
00552
00553 do{
00554 alpha += 1.;
00555 beta += 1.;
00556 gamma += 1.;
00557
00558
00559 Term = LastTerm*alpha*beta/gamma*chi/(i-1.);
00560
00561 Sum += Term;
00562
00563
00564 if( Sum.real() > 1e100 )
00565 {
00566 Sum /= Normalization;
00567 LastTerm = Term/Normalization;
00568 ++*NumRenormalizations;
00569
00570 fprintf( ioQQQ,"Hypergeometric: Renormalized at term %li. Sum = %.3e %.3e\n",
00571 i, Sum.real(), Sum.imag());
00572 }
00573 else
00574 LastTerm = Term;
00575
00576
00577
00578 if( fabs(LastTerm.real()/Sum.real())<0.001 && fabs(LastTerm.imag()/Sum.imag())<0.001 )
00579 lgNotConverged = false;
00580
00581 if( *NumRenormalizations >= 5 )
00582 {
00583 fprintf( ioQQQ, "We've got too many (%li) renorms!\n",*NumRenormalizations );
00584 }
00585
00586 ++i;
00587
00588 }while ( lgNotConverged || i<MinTerms );
00589
00590 *NumTerms = i;
00591
00592 return Sum;
00593 }
00594 #if defined(__ICC) && defined(__i386) && __INTEL_COMPILER < 910
00595 #pragma optimize("", on)
00596 #endif
00597
00598
00599 STATIC double RealF2_1( double alpha, double beta, double gamma, double chi )
00600 {
00601 long i = 3;
00602 bool lgNotConverged = true;
00603 double LastTerm, Sum;
00604
00605
00606 Sum = 1.;
00607
00608
00609 LastTerm = alpha*beta*chi/gamma;
00610
00611 Sum += LastTerm;
00612
00613
00614
00615 do{
00616 alpha++;
00617 beta++;
00618 gamma++;
00619
00620
00621 LastTerm *= alpha*beta*chi/gamma/(i-1.);
00622
00623 Sum += LastTerm;
00624
00625
00626
00627 if( fabs(LastTerm/Sum)<0.001 )
00628 lgNotConverged = false;
00629
00630 ++i;
00631
00632 }while ( lgNotConverged );
00633
00634 return Sum;
00635 }
00636
00637 STATIC complex<double> HyperGeoInt( double v )
00638 {
00639 return pow(v,BMinus1)*pow(1.-v,CMinusBMinus1)*pow(1.-v*GlobalCHI,MinusA);
00640 }
00641
00642
00643
00644 STATIC complex<double> qg32complex(
00645 double xl,
00646 double xu,
00647
00648 complex<double> (*fct)(double) )
00649 {
00650 double a,
00651 b,
00652 c;
00653 complex<double> y;
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669 a = .5*(xu + xl);
00670 b = xu - xl;
00671 c = .498631930924740780*b;
00672 y = .35093050047350483e-2 * ( (*fct)(a+c) + (*fct)(a-c) );
00673 c = b*.49280575577263417;
00674 y += .8137197365452835e-2 * ( (*fct)(a+c) + (*fct)(a-c) );
00675 c = b*.48238112779375322;
00676 y += .1269603265463103e-1 * ( (*fct)(a+c) + (*fct)(a-c) );
00677 c = b*.46745303796886984;
00678 y += .17136931456510717e-1* ( (*fct)(a+c) + (*fct)(a-c) );
00679 c = b*.44816057788302606;
00680 y += .21417949011113340e-1* ( (*fct)(a+c) + (*fct)(a-c) );
00681 c = b*.42468380686628499;
00682 y += .25499029631188088e-1* ( (*fct)(a+c) + (*fct)(a-c) );
00683 c = b*.3972418979839712;
00684 y += .29342046739267774e-1* ( (*fct)(a+c) + (*fct)(a-c) );
00685 c = b*.36609105937014484;
00686 y += .32911111388180923e-1* ( (*fct)(a+c) + (*fct)(a-c) );
00687 c = b*.3315221334651076;
00688 y += .36172897054424253e-1* ( (*fct)(a+c) + (*fct)(a-c) );
00689 c = b*.29385787862038116;
00690 y += .39096947893535153e-1* ( (*fct)(a+c) + (*fct)(a-c) );
00691 c = b*.2534499544661147;
00692 y += .41655962113473378e-1* ( (*fct)(a+c) + (*fct)(a-c) );
00693 c = b*.21067563806531767;
00694 y += .43826046502201906e-1* ( (*fct)(a+c) + (*fct)(a-c) );
00695 c = b*.16593430114106382;
00696 y += .45586939347881942e-1* ( (*fct)(a+c) + (*fct)(a-c) );
00697 c = b*.11964368112606854;
00698 y += .46922199540402283e-1* ( (*fct)(a+c) + (*fct)(a-c) );
00699 c = b*.7223598079139825e-1;
00700 y += .47819360039637430e-1* ( (*fct)(a+c) + (*fct)(a-c) );
00701 c = b*.24153832843869158e-1;
00702 y += .4827004425736390e-1 * ( (*fct)(a+c) + (*fct)(a-c) );
00703 y *= b;
00704
00705
00706
00707 return( y );
00708 }