******************************************************************************* * * FASTBBN: A Fortran code for computing primordial light element abundances * and chi^2 statistics for comparison with the observed values of * Y(He4), D/H and Li7/H (a default set is encoded below under * "EXPERIMENTAL DATA SET" and can be edited as required). * * References: G. Fiorentini, E. Lisi, S. Sarkar and F.L. Villante, * Phys. Rev. D58 (1998) 063506 [astro-ph/9803177]; * E. Lisi, S. Sarkar and F.L. Villante, * Phys. Rev. D59 (1999) 123520 [hep-ph/9901404]. * * IF YOU USE THIS CODE, PLEASE CITE THE ABOVE PAPERS WHERE IT IS DESCRIBED. * * URL: http://www-thphys.physics.ox.ac.uk/users/SubirSarkar/bbn/fastbbn.f * ******************************************************************************* * * NOTATION: * * effN = effective number of neutrinos, 1 < effN < 5 * eta = nucleon-to-photon ratio (= 2.73 x 10^{-8} \Omega_N h^2) * x = Log(eta)+10, 0 < x < 1 * * Elemental index i: * i=1 : 4He * i=2 : D/H * i=3 : Li/H * i=4 : 3He/H (not used in fits) * * Reaction index k: * k=1 : n decay * k=2 : p(n,gamma)d * k=3 : d(p,gamma)3He * k=4 : d(d,n)3He * k=5 : d(d,p)t * k=6 : t(d,n)4He * k=7 : t(alpha,gamma)7Li * k=8 : 3He(n,p)t * k=9 : 3He(d,p)4He * k=10 : 3He(alpha,gamma)7Be * k=11 : 7Li(p,alpha)4He * k=12 : 7Be(n,p)7Li * ******************************************************************************* * * The most important functions in this code are: * * Y(i,x,effN), which calculates the abundance Y(i) for given * values of x and effN; * * alpha(i,k,x,effN), which calculates the logarithmic derivative of * Y(i) with respect to the reaction rate R(k) for * given values of x and effN. * * CHI2(x,effN,Icorr),which calculates the chi2 including theoretical * and experimental uncertainties (with or without * correlations). * * The subroutine INITIALIZE assigns the relevant input * (including the experimental data set) for calculations and fits * ******************************************************************************* * Typical program structure: * * Set and print relevant input: CALL INITIALIZE(1) * DO loop on x * DO loop on effN * e.g. : * x = 0.7 effN = 3.2 * Calculate chisquare CHI2 c = CHI2(x,effN,1) Print*,' x = ',x,'; N =',effN,'; chisquare = ',c * END DO loop on effN * END DO loop on x STOP END * * ====================================================================== subroutine initialize(Iprint) * ====================================================================== * Sets experimental and ``theoretical'' (reaction rate) input. * Typically, one just needs to update the experimental abundances * Y_ex and their errors dY_ex. They are printed for Iprint \neq 0. * Notation: * ...................................................................... * Experimental input for the abundances * Y_ex (i), i=1,2,3 : experimental values of 4He , D/H, 7Li/H * dY_ex(i), i=1,2,3 : 1 sigma errors of " " " * * (3He is not used in fits) * ...................................................................... * One can also change the following input for the reaction rates: * * R = Ratio of reactions rates 1...12 with respect to the * following default values: * for k=1 the n decay rate is 886.7 sec, from PDG '98 * for k>1 the reaction rates are taken from the paper * ApJ Suppl. 85, 219 (1993) by Smith, Kawano, * and Malaney (SKM) * * [Example: R(1)=1.02 means that neutron lifetime=904.4 sec; * R(2)=0.97 means that p(n,gamma)d rate is 3% smaller * than SKM; etc.] * Changes from the default values R(k)=1 affect the output abundances * and are implemented in the function Y. * * dR = 1 sigma error of R. Present default values: * dR(1) taken from PDB'98; * dR(k), k>1, taken from SKM [see also Table II * in Fiorentini et al., PRD 58, 063506 (1998)]. * * ........................................................................ common/theor/R(12),dR(12) ! shared with CHI2,Y common/exper/Y_ex(3),dY_ex(3) ! shared with CHI2 * ........................................................................ * EXPERIMENTAL DATA SET: * ........................................................................ * 4He D/H 7Li/H * data Y_ex / 0.245 , 3.40e-5 , 1.73e-10 / * data dY_ex/ 0.004 , 0.30e-5 , 0.21e-10 / * * ........................................................................ * INPUT REACTION RATES: * ........................................................................ * k = 1 2 3 4 5 6 7 8 9 10 11 12 * data R/ 1. , 1. ,1. ,1. ,1. ,1. ,1. ,1. ,1. ,1. ,1. ,1. / * data dR/.0021,.07 ,.10 ,.10 ,.10 ,.08 ,.26 ,.10 ,.08 ,.16 ,.08 ,.09 / * ........................................................................ if (Iprint.eq.0) goto 10 print* print*,' -------------------------------------------------------' print*,' 4He mass fraction = ',Y_ex(1),' +/-',dY_ex(1) print*,' D/H = ',Y_ex(2),' +/-',dY_ex(2) print*,' 7Li/H = ',Y_ex(3),' +/-',dY_ex(3) print*,' -------------------------------------------------------' print* 10 return end * ======================================================================== function chi2(x,effN,Icorr) * ======================================================================== * Calculates theoretical uncertainties and total chisquare for given * values of x and of effN, assuming the input assigned in the * subroutine INITIALIZE. * Icorr \neq 0 (default) assumes correlated theoretical errors * Icorr = 0 assumes uncorrelated theoretical errors * ........................................................................ common/theor/R(12),dR(12) ! shared with INITIALIZE, Y common/exper/Y_ex(3),dY_ex(3) ! shared with INITIALIZE * ........................................................................ dimension sigma2(3,3) ! squared (total) error matrix dimension weight(3,3) ! ... and its inverse dimension sigma(3) ! diagonal (total) errors dimension corr(3,3) ! (total) correlation matrix dimension Yth(3),Yex(3),dYex(3) ! abundances multiplied by * ! proper power of 10 (below) * ........................................................................ dimension power(3) ! powers of 10 to make all * ! abundance values of order 1 * ! (avoids ill matrix inversion) data power/ 1.e1, 1.e4, 1.e9/ * ........................................................................ if (x.lt.0.or.x.gt.1.) then print*,'Warning: X value out of range' goto 10 endif if (effN.lt.1.or.effN.gt.5.) then print*,'Warning: Number of neutrinos out of range' goto 10 endif * ........................................................................ do i=1,3 Yth(i) = Y(i,x,effN) *power(i) Yex(i) = Y_ex(i) *power(i) dYex(i) = dY_ex(i) *power(i) enddo * ........................................................................ * Fills sigma2 with theoretical and experimental errors * do i=1,3 do j=1,i sigma2(i,j)=0. do k=1,12 sigma2(i,j) = sigma2(i,j)+ 1 alpha(i,k,x,effN)*alpha(j,k,x,effN)*(dR(k)/R(k))**2 enddo sigma2(i,j) = sigma2(i,j)*Yth(i)*Yth(j) sigma2(j,i) = sigma2(i,j) enddo sigma2(i,i) = sigma2(i,i)+dYex(i)**2 enddo * ........................................................................ * Cancels off-diagonal entries if Icorr=0 if (Icorr.eq.0) then do i=1,3 do j=1,3 if(i.ne.j) sigma2(i,j)=0. enddo enddo endif * ........................................................................ * Auxiliary calculations for possible checks: * * abundance errors sigma = square root of diagonal entries do i=1,3 sigma(i) = sqrt(sigma2(i,i)) enddo * * Correlations between errors (only if Icorr \neq 0) if (icorr.ne.0) then do i=1,3 do j=1,3 corr(i,j)=sigma2(i,j)/sigma(i)/sigma(j) enddo enddo endif * ........................................................................ * Calculation of weight matrix = inverse of the sigma2 matrix * wa=sigma2(1,1) wd=sigma2(2,2) wf=sigma2(3,3) wb=sigma2(1,2) wc=sigma2(1,3) we=sigma2(2,3) det=wa*wd*wf+2.*wb*we*wc-wc**2*wd-we**2*wa-wb**2*wf weight(1,1)=(wd*wf-we**2)/det weight(1,2)=(we*wc-wb*wf)/det weight(1,3)=(wb*we-wc*wd)/det weight(2,2)=(wa*wf-wc**2)/det weight(2,3)=(wb*wc-wa*we)/det weight(3,3)=(wa*wd-wb**2)/det weight(2,1)=weight(1,2) weight(3,1)=weight(1,3) weight(3,2)=weight(2,3) * .................................................................... * Chi2 calculation chi2=0. do i=1,3 do j=1,3 chi2=chi2+(Yth(i)-Yex(i))*weight(i,j)*(Yth(j)-Yex(j)) enddo enddo * ..................................................................... 10 return end * ======================================================================== function Y(i,x,effN) * ======================================================================== * This function calculates the elemental abundance Y for any (reasonable) * effective number of neutrinos effN. The calculation is based on an * approximation which relates the case effN \neq 3 to the case effN = 3, * through an appropriate shift in x, plus a rescaling of the abundances. * * i=1 : 4He mass fraction * i=2 : D/H (by number) * i=3 : 7Li/H (by number) * i=4 : 3He/H (by number) * ....................................................................... common/theor/R(12),dR(12) ! shared with INITIALIZE, CHI2 * Shift vectors: dimension as(4) ! coefficients of (N-3) dimension bs(4) ! coefficients of (N-3)**2 * Rescaling vectors dimension ar(4) ! coefficients of (N-3) dimension br(4) ! coefficients of (N-3)**2 data as/ +1.7012E-002,-3.3683E-002,-3.8980E-002,-3.3890E-002/ data bs/ 0. , 0. ,+2.2750E-003, 0. / data ar/ +5.3574E-002,+2.6214E-002,+9.0458E-002,+3.1001E-004/ data br/ -4.1377E-003,-7.1690E-003,-2.3649E-003,-3.6701E-003/ if (i.le.0.or.i.ge.5) then print*,'Warning: Elemental index out of range' goto 10 endif if (x.lt.0.or.x.gt.1.) then print*,'Warning: X value out of range' goto 10 endif if (effN.lt.1.or.effN.gt.5.) then print*,'Warning: Number of neutrinos out of range' goto 10 endif shift = as(i)*(effN-3.) + bs(i)*(effN-3.)**2 rescl = 1. + ar(i)*(effN-3.) + br(i)*(effN-3.)**2 Y = rescl * Y3nu(i,x+shift) * Includes the effect of possible changes in the reaction rates R(k) * (as assigned in INITIALIZE) Update = 1. do k=1,12 Update = Update + alpha(i,k,x,effN)*(R(k)-1.)/1. enddo Y = Y * Update 10 return end * ======================================================================== function alpha(i,k,x,effN) * ======================================================================== * This function calculates the logarithmic derivatives for the elemental * abundances and for any given (reasonable) effective number of neutrinos * effN. The calculation is based on an approximation which relates the * case effN \neq 3 to the case effN = 3 through an appropriate shift in * x, plus a rescaling of the log derivative. * Shift arrays: dimension as(12,4) ! coefficients of (N-3) dimension bs(12,4) ! coefficients of (N-3)**2 * Rescaling arrays: dimension ar(12,4) ! coefficients of (N-3) dimension br(12,4) ! coefficients of (N-3)**2 data as/ 1 -2.328E-03,-1.600E-02 , 0.,-3.952E-02,-3.534E-02,7*0., ! 4He 2 12*-3.3683E-02, ! D/H 3 12*-3.8980E-02, ! Li 4 12*-3.3890E-02/ ! 3He data bs/ 1 -2.191E-03,+6.409E-04, 0.,+2.440E-03,+1.281E-03,7*0., ! 4He 2 12*0., ! D/H 3 12*2.3649E-03, ! Li 4 12*0./ ! 3He data ar/ 1 -4.77E-02,+6.411E-04, 0.,-7.220E-02,-7.046E-02,7*0., ! 4He 2 12*0., ! D/H 3 12*0., ! Li 4 12*0./ ! 3He data br/ 1 +4.520E-03,+8.939E-05, 0.,+7.582E-03,+7.138E-03,7*0., ! 4He 2 12*0., ! D/H 3 12*0., ! Li 4 12*0./ ! 3He if (i.le.0.or.i.ge.5) then print*,'Warning: Elemental index out of range' goto 10 endif if (k.le.0.or.k.ge.13) then print*,'Warning: Reaction index out of range' goto 10 endif if (x.lt.0.or.x.gt.1.) then print*,'Warning: X value out of range' goto 10 endif if (effN.lt.1.or.effN.gt.5.) then print*,'Warning: Number of neutrinos out of range' goto 10 endif shift = as(k,i)*(effN-3.) + bs(k,i)*(effN-3.)**2 rescl = 1. + ar(k,i)*(effN-3.) + br(k,i)*(effN-3.)**2 alpha = rescl * alpha3nu(i,k,x+shift) 10 return end * ======================================================================== function Y3nu(i,x) * ======================================================================== * This function calculates the elemental abundance Y for N = 3 neutrinos. * The calculation is based on a polynomial approximation, which is * accurate in the extended x range [-0.2,+1.2]. This range is larger * than the default interval [0,1], in order to allow the "shift+rescale" * approximation for number of neutrinos different from 3. * * Warning: this function does not include the effect of changes of * default reaction rates R(k); such changes are implemented in the * function Y(i,x,effN). * * i=1 : 4He mass fraction * i=2 : D/H (by number) * i=3 : 7Li/H (by number) * i=4 : 3He/H (by number) * ........................................................................ * dimension a(8,4) data a/!1,5 2,6 3,7 4,8 1 +0.2233747E+00,+0.5529294E-01,-0.6255666E-01,+0.6843140E-01, * -0.4552996E-01,+0.1339232E-01,+0.1850735E-02,-0.1613542E-02, 2 +0.4809826E-03,-0.1799821E-02,+0.3103564E-02,-0.2568801E-02, * -0.9861666E-04,+0.2185143E-02,-0.1757214E-02,+0.4658893E-03, 3 +0.5391880E-09,-0.2684076E-08,+0.5951011E-08,-0.3852929E-08, * -0.8551272E-08,+0.2127352E-07,-0.1517732E-07,+0.3586446E-08, 4 +0.3444506E-04,-0.6246941E-04,+0.8363856E-04,-0.9718315E-04, * +0.5253057E-04,+0.3667086E-04,-0.6075325E-04,+0.2121180E-04/ * Initial checks ......................................................... if (i.le.0.or.i.ge.5) then print*,'Warning: Elemental index out of range' goto 10 endif if (x.lt.-0.2.or.x.gt.1.2) then print*,'Warning: X value out of range' goto 10 endif * Polynomial expansion ................................................... pol = a(8,i) do n=7,1,-1 pol = pol*x+a(n,i) enddo Y3nu = pol 10 return end * ======================================================================== function alpha3nu(i,k,x) * ======================================================================== * This function calculates the logarithmic derivatives for the elemental * abundances and for N=3 neutrinos. * Initial checks ......................................................... if (i.le.0.or.i.ge.5) then print*,'Warning: Elemental index out of range' goto 10 endif if (k.le.0.or.k.ge.13) then print*,'Warning: Reaction index out of range' goto 10 endif if (x.lt.-0.2.or.x.gt.1.2) then print*,'Warning: X value out of range' goto 10 endif * ........................................................................ if (i.eq.1) then alpha3nu = alpha4He3nu(k,x) elseif (i.eq.2) then alpha3nu = alphaD3nu(k,x) elseif (i.eq.3) then alpha3nu = alphaLi3nu(k,x) elseif (i.eq.4) then alpha3nu = alpha3he3nu(k,x) endif * ........................................................................ 10 return end * ======================================================================== function alpha4He3nu(k,x) * ======================================================================== * This function calculates the logarithmic derivatives for He4 * and for N=3 neutrinos. * The calculation is based on a polynomial approximation. dimension b(8,12) data b/!1,5 2,6 3,7 4,8 1 +0.8129781E+00,-0.1756227E+00,+0.1226484E+00,-0.5153191E-01, * -0.6839997E-01,+0.1199809E+00,-0.7250794E-01,+0.1623198E-01, 2 +0.6321668E-01,-0.2313505E+00,+0.4092393E+00,-0.3741036E+00, * +0.6069960E-01,+0.2198741E+00,-0.2013958E+00,+0.5629734E-01, 3 0. , 0. , 0. , 0. , * 0. , 0. , 0. , 0. , 4 +0.8455370E-02,-0.9458704E-02,+0.1511671E-01,-0.9999229E-02, * -0.8448764E-02,+0.2195439E-01,-0.1624701E-01,+0.4274994E-02, 5 +0.7707554E-02,-0.9410219E-02,+0.1499996E-01,-0.1098785E-01, * -0.4326561E-02,+0.1535783E-01,-0.1139675E-01,+0.2927043E-02, 6 0. , 0. , 0. , 0. , * 0. , 0. , 0. , 0. , 7 0. , 0. , 0. , 0. , * 0. , 0. , 0. , 0. , 8 0. , 0. , 0. , 0. , * 0. , 0. , 0. , 0. , 9 0. , 0. , 0. , 0. , * 0. , 0. , 0. , 0. , O 0. , 0. , 0. , 0. , * 0. , 0. , 0. , 0. , 1 0. , 0. , 0. , 0. , * 0. , 0. , 0. , 0. , 2 0. , 0. , 0. , 0. , * 0. , 0. , 0. , 0. / * Initial checks ......................................................... if (k.le.0.or.k.ge.13) then print*,'Warning: Reaction index out of range' goto 10 elseif (x.lt.-0.2.or.x.gt.1.2) then print*,'Warning: X value out of range' goto 10 endif * Polynomial expansion ................................................... pol = b(8,k) do n=7,1,-1 pol = pol*x+b(n,k) enddo alpha4He3nu = pol 10 return end * ======================================================================== function alphaD3nu(k,x) * ======================================================================== * This function calculates the logarithmic derivatives for D * and for N=3 neutrinos. * The calculation is based on a polynomial approximation. dimension c(8,12) data c/!1,5 2,6 3,7 4,8 1 +0.7093681E+00,-0.7132880E+00,-0.4876382E-01,+0.1201145E+01, * -0.7744275E+00,-0.3913517E+00,+0.7042471E+00,-0.2214912E+00, 2 -0.6871641E+00,+0.1268453E+00,+0.1868096E+01,-0.2187927E+01, * +0.1982162E+01,-0.3224197E+01,+0.2897048E+01,-0.8932150E+00, 3 -0.2394200E-01,-0.8854199E-01,-0.1666217E+00,-0.1872193E+00, * -0.9114170E-01,-0.5496869E-01,+0.3692332E-01,-0.7424205E-01, 4 -0.4271506E+00,-0.1643772E+00,-0.1521546E-01,+0.1439382E+00, * -0.1881865E+00,+0.1371702E+00,-0.3145297E-01,-0.4865135E-02, 5 -0.4181011E+00,-0.1286623E+00,+0.1289934E-01,+0.1850324E+00, * -0.2002014E-01,-0.2932434E+00,+0.3015885E+00,-0.9114459E-01, 6 0. , 0. , 0. , 0. , * 0. , 0. , 0. , 0. , 7 0. , 0. , 0. , 0. , * 0. , 0. , 0. , 0. , 8 -0.7252152E-02,+0.5424960E-02,+0.2910125E-01,+0.8415705E-01, * +0.4563072E-01,-0.4521475E+00,+0.4719815E+00,-0.1449376E+00, 9 -0.2420464E-02,-0.8188533E-02,-0.3052859E-01,-0.3696784E-01, * +0.2034117E+00,-0.1506076E+00,-0.3431342E-01,+0.2941599E-01, O 0. , 0. , 0. , 0. , * 0. , 0. , 0. , 0. , 1 0. , 0. , 0. , 0. , * 0. , 0. , 0. , 0. , 2 0. , 0. , 0. , 0. , * 0. , 0. , 0. , 0. / * Initial checks ......................................................... if (k.le.0.or.k.ge.13) then print*,'Warning: Reaction index out of range' goto 10 elseif (x.lt.-0.2.or.x.gt.1.2) then print*,'Warning: X value out of range' goto 10 endif * Polynomial expansion ................................................... pol = c(8,k) do n=7,1,-1 pol = pol*x+c(n,k) enddo alphaD3nu = pol 10 return end * ======================================================================== function alphaLi3nu(k,x) * ======================================================================== * This function calculates the logarithmic derivatives for 7Li * and for N=3 neutrinos. * The calculation is based on a polynomial approximation. dimension d(8,12) data d/!1,5 2,6 3,7 4,8 1 +0.2018333E+01,-0.5716214E+00,-0.2698362E+01,-0.8717163E+01, * -0.9917103E+01,+0.8180265E+02,-0.9237077E+02,+0.3104329E+02, 2 -0.9568990E+00,-0.1061615E+01,+0.7287105E+01,+0.2500068E+02, * +0.7431018E+01,-0.1754832E+03,+0.2111165E+03,-0.7274162E+02, 3 +0.1854156E-01,+0.1832449E-01,-0.1296782E+01,+0.1878546E+01, * +0.2222105E+02,-0.5231556E+02,+0.4140582E+02,-0.1123666E+02, 4 +0.1177149E+00,-0.7221449E+00,-0.1795135E+01,+0.9849410E+01, * +0.3225521E+02,-0.1120833E+03,+0.1045670E+03,-0.3173481E+02, 5 +0.1853383E+00,-0.1556603E+00,+0.3262168E+00,-0.2790330E+00, * -0.1471552E+01,+0.3228753E+01,-0.2714358E+01,+0.8573074E+00, 6 -0.9993449E+00,-0.3307726E-01,-0.6445227E+00,+0.8243842E+01, * +0.2154546E+02,-0.8552964E+02,+0.8347858E+02,-0.2609343E+02, 7 +0.9693255E+00,-0.3077873E-01,+0.5714311E+00,-0.7386952E+01, * -0.1950403E+02,+0.7694436E+02,-0.7494321E+02,+0.2339704E+02, 8 +0.5350964E-01,+0.2924453E+00,+0.3748815E+00,-0.6103750E+01, * -0.8913786E+01,+0.4736655E+02,-0.4903853E+02,+0.1579026E+02, 9 -0.6525239E-01,-0.2793919E+00,+0.3377858E+00,-0.1345531E+01, * -0.1374399E+02,+0.3570196E+02,-0.2984696E+02,+0.8424692E+01, O +0.3094791E-01,+0.3423862E-01,-0.5821547E+00,+0.7301099E+01, * +0.1973686E+02,-0.7680888E+02,+0.7430356E+02,-0.2307395E+02, 1 -0.1421146E+01,+0.2498440E+00,+0.1106049E+01,+0.3819467E+01, * +0.2665805E+02,-0.8781223E+02,+0.8288752E+02,-0.2551080E+02, 2 -0.8960841E-02,+0.1474491E+00,+0.1186384E+01,-0.9977886E+01, * -0.2606076E+02,+0.1060934E+03,-0.1045214E+03,+0.3277688E+02/ * Initial checks ......................................................... if (k.le.0.or.k.ge.13) then print*,'Warning: Reaction index out of range' goto 10 elseif (x.lt.-0.2.or.x.gt.1.2) then print*,'Warning: X value out of range' goto 10 endif * Polynomial expansion ................................................... pol = d(8,k) do n=7,1,-1 pol = pol*x+d(n,k) enddo alphaLi3nu = pol 10 return end * ======================================================================== function alpha3He3nu(k,x) * ======================================================================== * This function calculates the logarithmic derivatives for He3 * and N=3 neutrinos. * The calculation is based on a polynomial approximation. dimension e(8,12) data e/!1,5 2,6 3,7 4,8 1 +0.1033272E+00,-0.5421601E-01,+0.9241735E+00,-0.5909760E+00, * -0.4146949E+01,+0.8359055E+01,-0.6014316E+01,+0.1540252E+01, 2 +0.1023962E+00,+0.8180470E+00,-0.1327294E+01,-0.4638566E+01, * +0.1605368E+02,-0.1919912E+02,+0.1036359E+02,-0.2117321E+01, 3 +0.6122617E-01,+0.2099102E+00,+0.3215260E+00,+0.1151094E+00, * -0.2309159E+00,+0.8124533E-02,+0.1759614E-01,+0.1301250E-03, 4 +0.3076107E+00,+0.2536181E-01,-0.1072429E+00,-0.4543670E+00, * +0.2133507E+00,+0.6495192E+00,-0.7567540E+00,+0.2437597E+00, 5 -0.5138130E+00,+0.1523988E+00,+0.3166838E+00,-0.6019742E-01, * -0.8584032E-01,+0.1240784E-01,-0.3653218E-01,+0.2873841E-01, 6 -0.3146410E-01,+0.8001962E-01,-0.5958496E-01,-0.1793946E+00, * +0.4860333E+00,-0.4839249E+00,+0.2166142E+00,-0.3631611E-01, 7 0. , 0. , 0. , 0. , * 0. , 0. , 0. , 0. , 8 -0.5556295E+00,-0.7471940E-01,+0.1215732E+01,+0.1078082E+01, * -0.4441796E+01,+0.3535603E+01,-0.7325293E+00,-0.1186327E+00, 9 -0.1182454E+00,-0.4559168E+00,-0.9900842E+00,-0.4210743E+00, * +0.2776241E+01,-0.1801626E+01,-0.5988525E-01,+0.2399508E+00, O 0. , 0. , 0. , 0. , * 0. , 0. , 0. , 0. , 1 0. , 0. , 0. , 0. , * 0. , 0. , 0. , 0. , 2 0. , 0. , 0. , 0. , * 0. , 0. , 0. , 0. / * Initial checks ......................................................... if (k.le.0.or.k.ge.13) then print*,'Warning: Reaction index out of range' goto 10 elseif (x.lt.-0.2.or.x.gt.1.2) then print*,'Warning: X value out of range' goto 10 endif * Polynomial expansion ................................................... pol = e(8,k) do n=7,1,-1 pol = pol*x+e(n,k) enddo alpha3He3nu = pol 10 return end