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