subroutine calcdfi (eta, beta, $ ff, $ dff, $ ffth, $ dffth) c IN: ETA, BETA c OUT: FF, DFF, FFTH, DFFTH real*8 eta, beta real*8 ff(3), dff(6), ffth(3), dffth(6) c----------------DESCRIPTION----------------------------------- c Evaluation of the Fermi-Dirac integrals c and their derivatives 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 and the derivatives are c dfk_e(eta,beta) = integral (0 to inf.) of c x**k * (1+(1/2)*beta*x)**(1/2) dx / c ( exp((-eta+x)/2) + exp((eta-x)/2 )**2 c and c dfk_b(eta,beta) = integral (0 to inf.) of c (1/4)*x**(k+1)*(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 files c fdints.unf dfdints.unf c if BOTH files exist. Otherwise, The formatted files c fdints.tab dfdints.tab bessel.tab fermi7h.tab c are read and the unformatted files are written. c Initialization will take place if required. Initialization c may be forced by an explicit call to dfrmint_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 DERIVATIVES of THE FERMI-DIRAC INTEGRALS c dF1/2 / deta, dF3/2 / deta, dF5/2 / deta, c dF1/2 / dbeta, dF3/2 / dbeta, dF5/2 / dbeta real df(47,21,6) c ETA VALUES FOR FERMI INTEGRAL TABLES real etta(101) c FERMI INTEGRALS f1/2, f1, ... f3 real fermi(101,6) c DERIVATIVES of THE FERMI INTEGRALS df1/2, df1, ... df3 real dfermi(101,6) c FERMI INTEGRAL f 7/2 --- REQUIRED FOR NR d F 5/2 / d beta real f7half(101) 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, dfa, dfb, dj, dk, $ eb, eb2sr, fa, factor, fb, $ fj, fk, omdj, omdk, phi, $ pi, rtphi, xkki, xkkii, xx, y real*8 ffs(3) real*8 kki, kkii, lgphi, lgbeta real*4 smallbetand, medbetand cfermi real*8 fft, ffthr 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 df, dfermi, f7half save initialized, pi, cnste, cnstf, cnstg, cnsth, $ smallbetand, medbetand c---------------------FIRST EXECUTABLE STATEMENT--------------------- ka = 0 go to 1 entry dfrmint_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) chp open(36,file='dfdints.unf', readonly, open(36,file='dfdints.unf', status='readonly', $ form='unformatted') c read(36) dfh, dfo, dfth, dft, dffh, dfthr read(36) df read(36) dfermi read(36) f7half else c INTIALIZATION FROM FORMATTED FILE open(35,file='fdints.tab', status='readonly') open(37,file='dfdints.tab', status='readonly') open(38,file='bessel.tab', status='readonly') open(39,file='fermi7h.tab', status='readonly') chp open(35,file='fdints.tab', readonly) chp open(37,file='dfdints.tab',readonly) chp open(38,file='bessel.tab', readonly) chp open(39,file='fermi7h.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) read(37,'(2f8.3,6(1pe16.9))') ex, lx, $ (df(kk,j,l),l=1,6) 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,101 read(37,'(f8.3,6(1pe16.9))') ex, $ (dfermi(j,l),l=1,6) read(39,'(f8.3,1pe16.9)') ex, f7half(j) 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 do l = 1,6 df(kk,j,l) = alog(df(kk,j,l)) enddo enddo enddo do j = 1,101 do l = 1,6 fermi(j,l) = alog(fermi(j,l)) dfermi(j,l) = alog(dfermi(j,l)) enddo f7half(j) = alog(f7half(j)) 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) open(36, file='dfdints.unf', status='unknown', $ form='unformatted') write(36) df write(36) dfermi write(36) f7half close(36) 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.8 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) dff(1) = cnste*sqrt(eta)*eb2sr*(1.-cnstg/(eta*eb2sr**2)**2) dffth(1) = -cnstf/(sqrt(eta)*eb2sr)**3 dff(2) = cnste*sqrt(eta**3)*eb2sr*(1. + cnstg*(3. + 6.*eb $ + 2.*eb*eb)/(eta*eb2sr**2)**2) dffth(2) = cnstf*(3. + 6.*eb + 2.*eb*eb)/(sqrt(eta)*eb2sr**3) dff(3) = cnste*sqrt(eta**5)*eb2sr*(1. + cnstg*(15. + 20.*eb $ + 6.*eb*eb)/(eta*eb2sr**2)**2) dffth(3) = cnstf*(15. + 20.*eb + 6.*eb*eb)*sqrt(eta)/eb2sr**3 dff(4) = dff(1)*eta/beta + cnste*(-3./2.)*ay**5*ffs(1) $ + cnstf*(1./2.)*ayy/(sqrt(eta)*beta*eb2sr) dffth(4) = dffth(1)*eta/beta $ + cnstf*(1./2.)*ayy/(sqrt(eta)*beta*eb2sr) dff(5) = dff(2)*eta/beta + cnste*(-5./2.)*ay**7*ffs(2) $ + cnstf*(-1./2.)*sqrt(eta)*(3.+2.*eb)/(beta*eb2sr) dffth(5) = dffth(2)*eta/beta $ + cnstf*(-1./2.)*sqrt(eta)*(3.+2.*eb)/(beta*eb2sr) dff(6) = dff(3)*eta/beta + cnste*(-7./2.)*ay**9*ffs(3) $ + cnstf*(-3./2.)*sqrt(eta**3)*(5.+3.*eb)/(beta*eb2sr) dffth(6) = dffth(3)*eta/beta $ + cnstf*(-3./2.)*sqrt(eta**3)*(5.+3.*eb)/(beta*eb2sr) 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. dff(1) = ff(1) dff(2) = ff(2) dff(3) = ff(3) dff(4) = 1. + beta*(-5./8.+beta*(105./128. $ + beta*(-35./32.))) dff(4) = 3.*factor*dff(4)/8. dff(5) = 1. + beta*(-7./8. + beta*(-1407./2048.)) dff(5) = 15.*factor*dff(5)/16. dff(6) = 1. + beta*(-77/256.) dff(6) = 105.*factor*dff(6)/32. 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) dff(1)=ff(1) dff(2)=ff(2) dff(3)=ff(3) dff(4) = (dee1*rtphi/beta)*(phi*(kkii-kki) - 1.5*kki) dff(5) = (dee1*rtphi/beta**2) $ *(2.*phi*(kki - kkii) + 0.5*(5.*kki + kkii)) dff(6) = (dee1*rtphi/beta**2)*(1.5*kkii $ - 2.*phi*(2.*kki + kkii) $ + 4.*phi**2*(kkii - kki)) 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) dff(4) = (dee*rtphi/beta)*(phi*(kkii-kki) - 1.5*kki) dff(5) = (dee*rtphi/beta**2) $ *(2.*phi*(kki - kkii) + 0.5*(5.*kki + kkii)) dff(6) = (dee*rtphi/beta**2)*(1.5*kkii $ - 2.*phi*(2.*kki + kkii) $ + 4.*phi**2*(kkii - kki)) 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) c dff(4) = (ff(2) - 3.*ff(1)/2.)/beta c dff(5) = (3.*ff(1)-(2.-beta/2.)*ff(2))/beta/beta c dff(6) = ((9.-3./beta)*ff(1)+(9.*beta-4.+2./beta) c $ *ff(2)-(1.+5.*beta/2.)*ff(3))/beta/beta dff(4) = (dee*rtphi/beta)*(phi*(kkii-kki) - 1.5*kki) dff(5) = (dee*rtphi/beta**2) $ *(2.*phi*(kki - kkii) + 0.5*(5.*kki + kkii)) dff(6) = (dee*rtphi/beta**2)*(1.5*kkii $ - 2.*phi*(2.*kki + kkii) $ + 4.*phi**2*(kkii - kki)) c $ + 2.*phi*(2.*phi*(kkii - kki) c $ - 2.*kki - kkii)) endif dff(1) = ff(1) dff(2) = ff(2) dff(3) = ff(3) 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) dff(1) = ff(1) dff(2) = ff(2) dff(3) = ff(3) dff(4) = (1./(2*beta) - 1./(beta*beta))*ff(1) dff(5) = (1./(2*beta) - 1./(beta*beta))*ff(2) dff(6) = (1./(2*beta) - 1./(beta*beta))*ff(3) 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))) dff(l) = exp(ca*(dfermi(jap,2*l-1)) $ + cb*(dfermi(ja, 2*l-1))) if (l.gt.1) dff(l+2) = 0.25*ff(l) enddo c d F 5/2 / d beta = 1/4 f 7/2, dff(6) = 0.25*exp(ca*(f7half(jap)) $ + cb*(f7half(ja))) 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) if (l.gt.1) dff(l+2) = 0.25*ff(l) enddo do l = 1,3 do j = js, je yint(j - js + 1) = dfermi(j,2*l-1) enddo call polynomial_int (etaint, yint, pidimen, $ interp_order, eta, xx, dy) dff(l) = exp(xx) enddo do j = js, je yint(j - js + 1) = f7half(j) enddo call polynomial_int (etaint, yint, pidimen, $ interp_order, eta, xx, dy) dff(6) = 0.25*exp(xx) 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))) dff(l) = btatrt*exp(ca*(dfermi(jap,2*l)) $ + cb*(dfermi(ja, 2*l))) dff(l+3) = ff(l)/(2*beta) 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 do l = 1,3 do j = js, je yint(j - js + 1) = dfermi(j,2*l) enddo call polynomial_int (etaint, yint, pidimen, $ interp_order, eta, xx, dy) dff(l) = btatrt*exp(xx) dff(l+3) = ff(l)/(2*beta) 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 do l = 1,6 dfa = exp(omdj*(df(ja, ka,l)) $ + dj*(df(jap, ka,l))) dfb = exp(omdj*(df(ja,kap,l)) $ + dj*(df(jap,kap,l))) dff(l) = (dk*dfb**(2./3.) $ + omdk*dfa**(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 do l = 1,6 do j = js, je do k = ks, ke ypint(j-js+1,k-ks+1) = df(j,k,l) enddo enddo call TwoDpolyint (etaint, betaint, ypint, pidimen, $ pidimen, interp_order, interp_order, $ eta, betapass, xx, dy) dff(l) = exp(xx) enddo endif endif endif return cfermic ENTRY POINT TO GET VALUES OF FERMI INTEGRALS F2, F3 cfermic INPUT: ETA cfermiC OUTPUT: FFT, FFTHR cfermi cfermi entry fdmsls(eta,fft,ffthr) cfermi cfermi fj = 31. + eta cfermi if (fj.lt.1.) then cfermi fft = 2.*exp(eta) cfermi ffthr = 6.*exp(eta) cfermi elseif (fj.ge.101.) then cfermi fft = (eta**3)/3. cfermi ffthr = (eta**4)/4. cfermi else cfermi ja = int(fj) cfermi if (interp_order.eq.2) then cfermi jap = ja + 1 cfermi dj = fj - dble(ja) cfermi omdj = 1. - dj cfermi fft = exp(omdj*fermi(ja,4) + dj*fermi(jap,4)) cfermi ffthr = exp(omdj*fermi(ja,6) + dj*fermi(jap,6)) cfermi else cfermi je = min (101, ja+interp_order/2) cfermi js = max(1, je - interp_order) cfermi je = js + interp_order cfermi do l = 4,6,2 cfermi do j = js, je cfermi etaint(j - js + 1) = etta(j) cfermi yint(j - js + 1) = fermi(j,l) cfermi enddo cfermi call polynomial_int (etaint, yint, pidimen, cfermi $ interp_order, eta, xx, dy) cfermi if (l.eq.4) then cfermi fft = exp(xx) cfermi elseif (l.eq.4) then cfermi ffthr = exp(xx) cfermi endif cfermi enddo cfermi endif cfermi endif cfermi return end