subroutine calcfi (eta, beta, ff) c IN: ETA, BETA c OUT: FF real*8 eta, beta real*8 ff(3) c----------------DESCRIPTION----------------------------------- c Evaluation of the Fermi-Dirac integrals c by interpolation from table for c eta = mu/kt between -30 and 70, c beta = kT/mcc between 10**-6 and 10**4, c and by a combination of assymptotic formulae and c interpolation for eta and beta outside these limits. c The integrals are c c fk(eta,beta) = integral (0 to inf.) of c x**k * ( 1+(1/2)*beta*x )**(1/2) dx / c ( exp(x-eta) + 1 ) c---------------------HP-UX NOTE------------------------------- c On HP-UX, the readonly specification on OPEN statement c uses nonstandard syntax. See lines beginning with comment c chp c------------------------AUTHORS:------------------------------ c Ken Van Riper, Los Alamos National Laboratory c kvr@lanl.gov c Juan Antonia Miralles, University of Valencia c miralles@kerr.ific.uv.es c (addresses as of September, 1995) c------------REFERENCES:-------------------------------------- c "Accurate Evaluation of Fermi-Dirac Integrals and c Their Derivatives for Arbitrary Degeneracy and Relativity" c by J. A. Miralles and K. A. Van Riper, 1995, c submitted to Astrophysical Journal c "Equation of State of an Ideal Fermi Gas", 1977, c by S. A. Bludman and K. A. Van Riper, c Astrophysical Journal, 212}, p859. c------------------------INTIALIZATION------------------------ c Tables are read from the (FORTRAN) unformatted file c fdints.unf c if the file exists. Otherwise, The formatted files c fdints.tab bessel.tab c are read and the unformatted file is written. c Initialization will take place if required. Initialization c may be forced by an explicit call to frmint_itl c------------------VARIABLE DECLARATIONS-------------------------- c ETA AND LOG10 BETA GRIDS FOR CENTRAL TABLE : real etat(47), lgbtat(21), ex, lx c FERMI-DIRAC INTEGRALS F1/2, F3/2, F5/2 CENTRAL TABLE real f(47,21,3) c ETA VALUES FOR FERMI INTEGRAL TABLES real etta(101) c FERMI INTEGRALS f1/2, f1, ... f3 real fermi(101,6) c GRID FOR BESSEL FUNCTIONS real lgzt(63) c BESSEL FUNCTIONS real bessel(63,2) real*8 algy, ay, ayy, btatrt, ca, cb, $ cnste, cnstf, cnstg, cnsth, da, db, $ dee, dee1, dj, dk, $ eb, eb2sr, fa, factor, fb, $ fj, fk, omdj, omdk, phi, $ pi, rtphi, xkki, xkkii, xx, y real*8 ffs(3), ffth(3) real*8 kki, kkii, lgphi, lgbeta real*4 smallbetand, medbetand integer j, ja, jap, je, js, k, ka, kap, $ ke, kk, ks, l c VARIABLES FOR POLYNOMIAL INTERPOLATION integer pidimen parameter (pidimen=7) real*8 ypint (pidimen,pidimen), etaint(pidimen), $ betaint(pidimen), yint(pidimen), dy, betapass logical initialized, lv, lvd c EVAL value on return is which table or interp`n block was used c INTERP_ORDER Order of polynomial interpolation c If not explicity set, --> 6 on initialization integer eval, interp_order common /evalt/ eval, interp_order data interp_order /-1/ data initialized /.false./ save etat, lgbtat, etta, lgzt, f, fermi, bessel save initialized, pi, cnste, cnstf, cnstg, cnsth, $ smallbetand, medbetand c---------------------FIRST EXECUTABLE STATEMENT--------------------- ka = 0 go to 1 entry frmint_itl ka = 1 1 continue if (.not.initialized) then c-------------------------INITIALIZATION BLOCK------------------------- inquire (file = 'fdints.unf', exist=lv) inquire (file = 'dfdints.unf', exist=lvd) if (lv.and.lvd) then c INTIALIZATION FROM UNFORMATTED FILE chp open(34,file='fdints.unf', readonly, open(34,file='fdints.unf', status='readonly', $ form='unformatted') read(34) etat,lgbtat read(34) f read(34) etta read(34) fermi read(34) lgzt, bessel close(34) else c INTIALIZATION FROM FORMATTED FILE open(35,file='fdints.tab', status='readonly') open(38,file='bessel.tab', status='readonly') chp open(35,file='fdints.tab', readonly) chp open(38,file='bessel.tab', readonly) do kk = 1,47 do j = 1,21 read(35,'(2f8.3,3(1pe16.9))') etat(kk),lgbtat(j), $ (f(kk,j,l),l=1,3) enddo enddo do j = 1,101 read(35,'(f8.3,6(1pe16.9))') $ etta(j), (fermi(j,l),l=1,6) enddo do j = 1,59,1 read(38,'(f6.2,2(1pe17.10))') lgzt(j), $ bessel(j,1), bessel(j,2) enddo close(35) c STORE AS LOGARITHMS do kk = 1,47 do j = 1,21 do l = 1,3 f(kk,j,l) = alog(f(kk,j,l)) enddo enddo enddo do j = 1,101 do l = 1,6 fermi(j,l) = alog(fermi(j,l)) enddo enddo do j = 1,59 do l = 1,2 bessel(j,l) = alog(bessel(j,l)) enddo enddo c CREATE BINARY FILE WITH LOG OF VALUES open(34, file='fdints.unf', status='unknown', $ form='unformatted') write(34) etat,lgbtat write(34) f write(34) etta write(34) fermi write(34) lgzt, bessel close(34) endif cnste = 1./sqrt(2.) pi = 2.*abs(asin(1.)) cnstg = pi**2/6. cnstf = cnstg*cnste cnsth = sqrt(pi)/2. if (interp_order.le.0) interp_order = 6 c UNCOMMENT FOLLOWING AS DESIRED: cw write (6,'(a,i3)') ' Interpolation Order = ', interp_order c Bondaries for ND treatments smallbetand = -2.0 c use -2.0 for smallbetand when only c calculating functions w/o derivatives medbetand = 0.3 climits write (6,'(a)') '2 limits for Chebyshev use' climits read (5,*) smallbetand, medbetand cw write (6,'(a,f6.2)') cw $ ' Small Beta expansion used for Log10 beta < ', smallbetand cw write (6,'(a,f6.2,a,f6.2)') cw $ ' Chebyshev expansion used between ', smallbetand, cw $ ' < Log10 beta <= ', medbetand cw write (6,'(a,f6.2,a,f6.2)') cw $ ' Bessel Function Table used between ', medbetand, cw $ ' < Log10 beta <= ', 4.0 initialized = .true. c--------------------------------END INITIALIZATION BLOCK-------------- endif if (ka.eq.1) return lgbeta = log10(beta) phi = 1.0/beta if (eta .gt. 70.0) then c LARGE ETA eval = 1 ayy = 1. + eta*beta y = sqrt(eta*beta*(2.+eta*beta)) ay = 1.0/(sqrt(beta)) algy = log(y + ayy) eb = eta*beta eb2sr = sqrt(2.+eb) if (y.gt.0.05) then ffs(1) = 0.5*(y*ayy-algy) ffs(2) = (1./3.)*y**3 - 0.5*(y*ayy-algy) else ffs(1) = y**3*(1./3. + y*y*(-0.1 + y*y*(3./56. $ + y*y*(-5./144. + 35./1408.*y*y)))) ffs(2) = y**5*(0.1 + y*y*(-3./56. $ + y*y*(5./144. - 35./1408.*y*y))) endif if (y.gt.0.1) then ffs(3) = 0.625*y*ayy*(1.0+0.4*y**2) $ - (2./3.)*y**3 - 0.625*algy else ffs(3) = y**7*(1./28. + y*y*(-1./36. + 15./704.*y*y)) endif ffth(1) = cnstf*ayy/(sqrt(eta)*eb2sr) ff(1) = cnste*ay**3*ffs(1)+ffth(1) ffth(2) = cnstf*sqrt(eta)*(3.+2.*eb)/eb2sr ff(2) = cnste*ay**5*ffs(2)+ffth(2) ffth(3) = cnstf*(sqrt(eta))**3*(5.+3.*eb)/eb2sr ff(3) = cnste*ay**7*ffs(3)+ffth(3) elseif (eta .lt. -30.0) then if (lgbeta .lt. smallbetand) then eval = 2 c SMALL ETA, SMALL BETA factor = cnsth*exp(eta) ff(1) = 1. + beta*(3./8.+beta*(-15./128.+beta*(105./1024. $ + beta*(-105./1024.)))) ff(1) = factor*ff(1) ff(2) = 1. + beta*(5./8.+beta*(-35./128. $ + beta*(-2345./16384.))) ff(2) = 3.*factor*ff(2)/2. ff(3) = 1. + beta*(7./8.+beta*(-539./4090.)) ff(3) = 15.*factor*ff(3)/4. elseif(lgbeta. le. medbetand) then c Clenshaw and Tooper's Chebyshev expansion c Note that the factor exp(phi) is not c contained in the functions when evaluated c with the expansion. eval = 7 call besscheb (lgbeta, kki, kkii) dee1 = cnste*exp(eta) rtphi = sqrt(phi) ff(1) = dee1 * rtphi*kki ff(2) = dee1*rtphi**3*(kkii-kki) ff(3) = dee1*rtphi**5*(2.0*kki + (3.0*beta-2.0)*kkii) elseif (lgbeta .le. 4.0) then eval = 3 c SMALL ETA CENTRAL BETA lgphi = log10(phi) ja = int(10.0*lgphi + 41.0) if (ja.eq.59) then jap = 60 da = 0. db = 1. else jap = ja + 1 da = 10.0*(lgphi-lgzt(ja)) db = 10.*(lgzt(jap) - lgphi) endif if (interp_order.eq.2) then kki = exp(da *bessel(jap,1) + db *bessel(ja,1)) kkii = exp(da*bessel(jap,2) + db*bessel(ja,2)) dee = cnste*exp(eta+phi) rtphi = sqrt(phi) ff(1) = dee * rtphi*kki ff(2) = dee*rtphi**3*(kkii-kki) ff(3) = dee*rtphi**5*(2.0*kki + (3.0*beta-2.0)*kkii) else je = min (59, ja+interp_order/2) js = max(1, je - interp_order) je = js + interp_order do j = js, je betaint(j - js + 1) = lgzt(j) yint(j - js + 1) = bessel(j,1) enddo call polynomial_int (betaint, yint, pidimen, $ interp_order, lgphi, xkki, dy) kki = exp(xkki) do j = js, je yint(j - js + 1) = bessel(j,2) enddo call polynomial_int (betaint, yint, pidimen, $ interp_order, lgphi, xkkii, dy) kkii = exp(xkkii) dee = cnste*exp(eta+phi) rtphi = sqrt(phi) ff(1) = dee * rtphi*kki ff(2) = dee*rtphi**3*(kkii-kki) ff(3) = dee*rtphi**5*(2.0*kki + (3.0*beta-2.0)*kkii) endif else eval = 4 c SMALL ETA, LARGE BETA ff(1) = exp(eta+phi)*sqrt(beta*0.5) ff(2) = 2.0*ff(1) ff(3) = 6.0*ff(1) endif else if (lgbeta .lt. -6.0) then eval = 5 c SMALL BETA, CENTRAL ETA ja = int( 31.0 + eta) if (interp_order.eq.2) then if (ja.eq.101) then jap = 101 ca = 0. cb = 1. else jap = ja + 1 ca = eta - etta(ja) cb = etta(jap) - eta endif do l = 1,3 ff(l) = exp(ca*( fermi(jap,2*l-1)) $ + cb*(fermi(ja, 2*l-1))) enddo else je = min (101, ja+interp_order/2) js = max(1, je - interp_order) je = js + interp_order do j = js, je etaint(j - js + 1) = etta(j) enddo do l = 1,3 do j = js, je etaint(j - js + 1) = etta(j) yint(j - js + 1) = fermi(j,2*l-1) enddo call polynomial_int (etaint, yint, pidimen, $ interp_order, eta, xx, dy) ff(l) = exp(xx) enddo endif elseif (lgbeta .gt. 4.0) then eval = 6 c LARGE BETA, CENTRAL ETA ja = int( 31.0 + eta) btatrt = sqrt (beta*0.5) if (interp_order.eq.2) then jap = ja + 1 if (ja.eq.101) then jap = 101 ca = 0. cb = 1. else ca = eta - etta(ja) cb = etta(jap) - eta endif do l = 1,3 ff(l) = btatrt*exp(ca*( fermi(jap,2*l)) $ + cb*(fermi(ja, 2*l))) enddo else je = min (101, ja+interp_order/2) js = max(1, je - interp_order) je = js + interp_order do j = js, je etaint(j - js + 1) = etta(j) enddo do l = 1,3 do j = js, je etaint(j - js + 1) = etta(j) yint(j - js + 1) = fermi(j,2*l) enddo call polynomial_int (etaint, yint, pidimen, $ interp_order, eta, xx, dy) ff(l) = btatrt*exp(xx) enddo endif else eval = 0 c CENTRAL BLOCK c write(6,*) ' CENTRAL BLOCK (5)' if (eta.le.-5.0) then fj = 0.2*eta + 7.0 elseif (eta.le. -1.0) then fj = eta + 11.0 elseif (eta.le.1.0) then fj = 10.0*eta + 20.0 elseif (eta.le.5.0) then fj = eta + 29.0 else fj = 0.2*eta + 33.0 if (fj.eq.47.0) then ja = 47 jap = 47 dj = 0.0 go to 70 endif endif ja = int(fj) jap = ja + 1 dj = fj - dble(ja) 70 continue fk = 2.0*lgbeta + 13.0 ka = int (fk) if (ka.eq.21) then kap = 21 dk = 0.0 else kap = ka + 1 dk = fk - dble(ka) endif omdk = 1.0 - dk omdj = 1.0 - dj if (interp_order.eq.2) then c K in beta direction, J in eta direction do l = 1,3 fa = exp(omdj*(f(ja, ka,l)) $ + dj*(f(jap, ka,l))) fb = exp(omdj*(f(ja,kap,l)) $ + dj*(f(jap,kap,l))) ff(l) = (dk*fb**(2./3.) $ + omdk*fa**(2./3.))**1.5 enddo else c POLYNOMIAL INTERPOLATION je = min (47, ja+interp_order/2) js = max(1, je - interp_order) je = js + interp_order ke = min (21, ka+interp_order/2) ks = max(1, ke - interp_order) ke = ks + interp_order do j = js, je etaint(j - js + 1) = etat(j) enddo do k = ks, ke betapass = log10(beta) betaint(k - ks + 1) = lgbtat(k) enddo do l = 1,3 do j = js, je do k = ks, ke ypint(j-js+1,k-ks+1) = f(j,k,l) enddo enddo call TwoDpolyint (etaint, betaint, ypint, pidimen, $ pidimen, interp_order, interp_order, $ eta, betapass, xx, dy) ff(l) = exp(xx) enddo endif endif endif return end