Dear colleague, We have produced two subroutines for the thermal bremsstrahlung Gaunt factors. The first corresponds to our paper "Relativistic Thermal Bremsstrahlung Gaunt Factor for the Intracluster Plasma. II. Analytic Fiting Formulae" (submitted to ApJ; astro-ph/9906342). The second corresponds to our paper "Relativistic Thermal Bremsstrahlung Gaunt Factor for the Intracluster Plasma. III. Analytic Fitting Formula for the Nonrelativistic Exact Gaunt Factor" (submitted to ApJ; astro-ph/9906467). We believe these two subroutines would be sufficient for the analysis of the X-ray data from clusters of galaxies. However, some of you would be interested in still higher-temperature plasmas. In the first subroutine (paper II), we have given fitting formulae for log T (K) < 8.5. For higher temperatures the Gaunt factor rapidly increases, as you see in the figures in our published paper "Relativistic Thermal Bremsstrahlung Gaunt Factor for the Intracluster Plasma" (ApJ 507: 530-557). Therefore it is extremely difficult to produce accurate analytic fitting formulae for log T (K) > 8.5. For log T (K) > 8.5, the relativistic Elwert method described in our paper produces excellent results. The computer program for this calculation is rather simple. Therefore, it would be best to calculate the Gaunt factor directly by using the program for the relativistic Elwert calculation. For this reason we have decided to send a subroutine for the relativistic Elwert calculation to those who would be interested in receiving it. We would be delighted if this subroutine be of any help for your research. When you publish your work by using this subroutine, please quote our ApJ 507: 530-557 paper. We think the present subroutine is a part of this paper. We would be delighted to hear from you whatever comment you might have on this subroutine. Yours very sincerely, Naoki Itoh **************************************************** Naoki Itoh Department of Physics Sophia University 7-1 Kioi-cho, Chiyoda-ku, Tokyo 102 Japan Telephone (81 or 0) 3 3238 3431 FAX (81 or 0) 3 3238 3341 **************************************************** CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Main program: gaunt.f C C Calculates the relativistic Elwert Gaunt factor g_{zj}. C Uses Gauss quadrature integration method. C C Please input the following values. C C zj1 : the minimum of the atomic number C zj2 : the maximum of the atomic number C uu1 : the minimum of log10(u) C uu2 : the maximum of log10(u) C ud : the interval of log10(u) C tt1 : the minimum of log10(T) C tt2 : the maximum of log10(T) C td : the interval of log10(T) C C The Program produces accurate Gaunt factors for C Zj = 1 to 30, log10(u) = -4 to 1, log10(T) = 8 to 11 C CCCCCCCCCCCCCCCCCC 12.July.1999, at Sophia T. Sakamoto CCCCCC implicit real*8(a-z) C C For example C zj1 = 1.d0 zj2 = 30.d0 uu1 = -4.d0 uu2 = 1.d0 ud = 0.1d0 tt1 = 8.d0 tt2 = 11.d0 td = 0.1d0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC dd = 1.D-6 uu2 = uu2 + dd tt2 = tt2 + dd do 45 zj = zj1,zj2 do 5 uu = uu1,uu2,ud do 15 tt = tt1,tt2,td gaunt = gzj(zj,tt,uu) write(6,100) zj,tt,uu,gaunt 100 format('zj = ',f3.0,3x,'log T = ',f6.3,3x, & 'log u = ',f6.3,3x,'gzj = ',f14.7) 15 continue 5 continue 45 continue stop end c The calculation of gzj real*8 function gzj(zj,tt,uu) implicit real*8(a-h,j,o-z) parameter(n=100) real*8 x(n),w(n) n1 = 50 cd = -70.d0 co = 1.579D+5 pai = 3.141592d0 u = 10.d0**uu t = 10.d0**tt rr = zj*zj*co/t q = zj*zj*1.579d0*10.d0**(-4.d0)/(rr*5.93d0) v = cd + 1.d0/q a1 = 1.d0/q + u a2 = 1.d0/q b1 = a1 + 100.d0 b2 = a2 + 100.d0 call legaus(a1,b1,n1,x,w) sum1 = 0.d0 do 10 i = 1,n1 sum1 = sum1 + w(i) * j(x(i),q,v,u,zj) 10 continue call legaus(a2,b2,n1,x,w) sum2 = 0.d0 do 11 i = 1,n1 sum2 = sum2 + w(i) * g0(x(i),q,v,u,z) 11 continue g1 = 3.d0 * dsqrt(6.d0) / (32.d0*dsqrt(pai)) gzj = g1 * q**(7.d0/2.d0) * dexp(u) * sum1 / sum2 return end c----+-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+ c function j(x) c----+-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+ real*8 function j(x,q,v,u,z) implicit real*8(a-z) f = 1.d0/137.036d0 p = 3.141592d0 a = x - u b = x*x - 1.d0/(q*q) c = a*a - 1.d0/(q*q) bf = 2.d0 * dlog((a+dsqrt(c))*q) bi = 2.d0 * dlog((x+dsqrt(b))*q) l = 2.d0 * dlog((a*x+dsqrt(c*b)-1.d0/(q*q))*q/u) j1 = (b/(dexp(x-v)+1.d0))*(a/x)*(1.d0-1.d0/(dexp(a-v)+1.d0))*(1.d0 &-dexp(-2.d0*p*f*z*x/dsqrt(b)))/(1.d0-dexp(-2.d0*p*f*z*a &/dsqrt(c))) j21 = 4.d0/3.d0-2.d0*a*x*(b+c)/(b*c) + (bf*x*c**(-3.d0/2.d0) & + bi*a*b**(-3.d0/2.d0) - bf*bi/dsqrt(b*c))/(q*q) j221 = 8.d0*x*a/(3.d0*dsqrt(b*c)) + u*u*(a*a*x*x+b*c) &*(b*c)**(-3.d0/2.d0) j222 = (a*x+b)*bi*b**(-3.d0/2.d0) + 2.d0*u*a*x/(b*c) & - (a*x+c)*bf*c**(-3.d0/2.d0) j22 = l * (j221 + u/(2.d0*dsqrt(b*c)*q*q)*j222) j = j1 * (j21 + j22) return end c----+-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+ c function g0(x) c----+-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+ real*8 function g0(x,q,v) implicit real*8(a-z) g0 = q*q*q*x*dsqrt(x*x-1.d0/(q*q))/(1.d0+dexp(x-v)) return end c----+-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+ subroutine legaus(a,b,n,x,w) c----+-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+ * gauss-legendre integral for x=[a,b] * x(i) gives the zeros of legendre polynomial of n-th order * w(i) gives associate crystoffel numbers *----+-----+-----+-----+-----+ 4 jun. '85 -----+--- t.cheon -----+ implicit real*8(a-h,o-z) dimension x(n),w(n) data irpmx /1000/ data epsylon /1.d-14/ *----+ +-----+ +-----+ if (n.le.0) then write(6,*)'>> n=',n,' invalid argument (func. legaus)' *----+ +-----+ +-----+ else if (n.eq.1) then x(1)=0.d0 w(1)=2.d0 *----+ +-----+ +-----+ else nc=(n+1)/2 xcan=-1.d0 do 10 i=1,nc do 20 irp=1,irpmx ff=pp(n,xcan,ppd) df=ppd s=0.d0 do 30 j=1,i-1 s=s+1.d0/(xcan-x(j)) 30 continue dif=-ff/(df-ff*s) if(abs(dif).lt.epsylon) go to 99 xcan=xcan+dif 20 continue write(6,*)'### no convergence (subr.legaus)' write(6,*)'### at i=',i write(6,*)'### x-x''=',dif 99 continue x(i)=xcan w(i)=2.d0/(n*pp(n-1,xcan,ddd)*df) rr=0.d0 stepu=df-ff*s stepd1=(2*xcan*df-n*(n+1)*ff)/(1.d0-xcan*xcan) stepd2=-2*df*s+ff*s*s-rr*ff xcan=xcan+stepu/(stepd1+stepd2) 10 continue do 50 i=nc+1,n x(i)=-x(n-i+1) w(i)= w(n-i+1) 50 continue end if *----+ +-----+ +-----+ fa1=0.5d0*(b-a) fa2=0.5d0*(b+a) do 70 i=1,n x(i)=fa1*x(i)+fa2 w(i)=fa1*w(i) 70 continue *----+-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+ end c----+-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+ function pp(n,x,ppd) c----+-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+ * legendre polynomial of order n * ppd returns the 1st derivative of pp(x) *----+-----+-----+-----+-----+ 4 jun. '85 -----+--- t.cheon -----+ implicit real*8(a-h,o-z) p0=1.d0 p1=x pd0=0.d0 pd1=1.d0 if (n.eq.0) then pp = p0 ppd= pd0 else if (n.eq.1) then pp = p1 ppd= pd1 else pi1= p1 pi2= p0 pdi1=pd1 pdi2=pd0 do 10 i=2,n f1=(2.d0*i-1.d0)/(1.d0*i) f2=(1.d0*i-1.d0)/(1.d0*i) pi = f1*x*pi1 - f2*pi2 pdi= f1*(x*pdi1+pi1) - f2*pdi2 pi2 = pi1 pdi2= pdi1 pi1 = pi pdi1= pdi 10 continue pp=pi ppd=pdi end if end