CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C This program calculates the emissivity Wee(tau) C for the electron-electron thermal bremsstrahlung defined as C follows. C d Eee C Wee(tau) = ------- = mc^2 integral_{0}^{infinity} dk k Pee(k,tau) C dt dV C is the total energy emitted per unit time, per unit volume C by electron gas of uniform number density ne at temperature T, C where tau =kBT/mc^2 is the electron thermal energy in units of C the electron rest mass. C This program returns the analytic fitting function Wee^fit(tau) C defined by eq. (4.3) in Ref.1. C C Ref.1: S. Nozawa, K. Takahashi, Y. Kohyama and N. Itoh, C Analytic fitting formulae for relativistic electron-electron C thermal bremsstrahlung, accepted to A&A (2009) C C INPUTS C This program requires a dimensionless variable tau. C tau = kBT/mc^2 ----- real*8, 50eV < kBT C C OUTPUTS C This program returns the emissivity Wee. C Wee ----- real*8 C Wee = 1.192x10^{-22} [ne(cm^-3)]^2 tau^3/2*Gfit(tau) C erg cm^-3 sec^-1 C Please refer to eq. (4.3) in Ref.1 for detail. C C 2009.01 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PROGRAM MAIN implicit real*8 (a-h, o-z) data pi, emass/ 3.141592654d0, 510.998902d0/ write(*,*)'-----------------------------------------------------------' write(*,*)'Input: electron thermal energy kBT (in keV)' write(*,*)' kBT > 50 eV' write(*,*)'-----------------------------------------------------------' 100 write(*,*)'kBT (in keV) = ? kBT=0 will quit this program' read(*,*) T if (T .ne. 0.d0) then tau = T/emass CALL FUNCFx(TAU,Weefit) write(*,*)'Weefit =',Weefit write(*,*)'------------------------------------------------------' GOTO 100 else end if END SUBROUTINE FUNCFx(TAU,Weefit) IMPLICIT REAL*8 (A-Z) data pi, emass/ 3.141592654d0, 510.998902d0/ IF (TAU.ge.(0.05D0/emass).and.TAU.lt.(1D0/emass)) THEN G=GI(tau) ELSE IF (TAU.ge.(1D0/emass).and.TAU.lt.(300D0/emass)) THEN G=GII(tau) ELSE IF (TAU.ge.(300D0/emass).and.TAU.le.(7000D0/emass)) THEN G=GIII(tau) ELSE IF (TAU.gt.(7000D0/emass)) THEN G=GIV(tau) END IF Weefit=1.192D-22*DSQRT(TAU)**3*G RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Function GI(tau) returns the fitting function of Wee(x,tau). C C GI(x, tau) = DSQRT( 8D0 / (3D0*PI) ) * Jfit C C inputs: 50eV < tau < 1keV (tau = KbT/mc^2) C CThis analytic fitting formula was obtained by Itoh, Kawana & Nozawa (2002a). C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC REAL FUNCTION GI*8(TAU) IMPLICIT REAL*8 (A-Z) real*8 b(1:11) data b/ $ 2.21564D0 , 1.83879D-1,-1.33575D-1, 5.89871D-2, $ -1.45904D-2,-7.10244D-4, 2.80940D-3,-1.70485D-3, $ 5.26075D-4, 9.94159D-5,-1.06851D-4/ data pi, emass/ 3.141592654d0, 510.998902d0/ logth=LOG10(TAU) theta=(logth+2.65d0)/1.35d0 jfit= b( 1) * +b( 2)*theta * +b( 3)*theta** 2d0 * +b( 4)*theta** 3d0 * +b( 5)*theta** 4d0 * +b( 6)*theta** 5d0 * +b( 7)*theta** 6d0 * +b( 8)*theta** 7d0 * +b( 9)*theta** 8d0 * +b(10)*theta** 9d0 * +b(11)*theta**10d0 GI=DSQRT(8D0/(3D0*PI))*jfit return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Function GII(tau) returns the fitting function of Wee(x,tau). C C GII(tau) = G0+G1+G2 C C G0=2D0*AA(2)+AA(1)+AA(0)+0.5D0*BB(1)+BB(0) C G1=G1+GAM *AA(i)*CC(j) sum for i= 0 to 8, j= 2 to 6 C G2=G2+GAM/Z*BB(i)*CC(j) sum for i= 0 to 8, j= 2 to 6 C GAM is FUNCTION of GAMMA : GAMMA(Z), Z=i+j/8D0+1D0 C C AA2(tau) = A2(0) + A2(1)*tau^(1/8) + A2(2)*tau^(2/8) + A2(3)*tau^(3/8) C + A2(4)*tau^(4/8) + A2(5)*tau^(5/8) + A2(6)*tau^(6/8) C + A2(7)*tau^(7/8) + A2(8)*tau^(8/8) C C Same function forms for AA1(tau), AA0(tau), BB1(tau) and BB0(tau). C C inputs: 1keV < tau < 300keV (tau = KbT/mc^2) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC real function GII*8(tau) implicit real*8 (a-h, o-z) real*8 A2(0:8), A1(0:8), A0(0:8), B1(0:8), B0(0:8) real*8 C6(0:8), C5(0:8), C4(0:8), C3(0:8), C2(0:8) real*8 AA(0:2), BB(0:2), CC(2:6) data A2/ 0.9217D0, -13.4988D0, 76.4539D0, -217.8301D0, & 320.9753D0, -188.0667D0, -82.4161D0, 163.7191D0, & -60.0248D0/ data A1/ -9.3647D0, 95.9186D0, -397.0172D0, 842.9376D0, & -907.3076D0, 306.8802D0, 291.2983D0, -299.0253D0, & 76.3461D0/ data A0/ -37.3698D0, 380.3659D0, -1489.8014D0, 2861.4150D0, & -2326.3704D0, -691.6118D0, 2853.7893D0, -2040.7952D0, & 492.5981D0/ data B1/ -8.6991D0, 63.3830D0, -128.8939D0, -135.0312D0, & 977.5838D0, -1649.9529D0, 1258.6812D0, -404.7461D0, & 27.3354D0/ data B0/ -11.6281D0, 125.6066D0, -532.7489D0, 1142.3873D0, & -1156.8545D0, 75.0102D0, 996.8114D0, -888.1895D0, & 250.1386D0/ data C2/ -5.7752D0, 46.2097D0, -160.7280D0, 305.0070D0, & -329.5420D0, 191.0770D0, -46.2718D0, 0.0000D0, & 0.0000D0/ data C3/ 30.5586D0, -248.2177D0, 874.1964D0, -1676.9028D0, & 1828.8677D0, -1068.9366D0, 260.5656D0, 0.0000D0, & 0.0000D0/ data C4/ -54.3272D0, 450.9676D0, -1616.5987D0, 3148.1061D0, & -3478.3930D0, 2055.6693D0, -505.6789D0, 0.0000D0, & 0.0000D0/ data C5/ 36.2625D0, -310.0972D0, 1138.0531D0, -2260.8347D0, & 2541.9361D0, -1525.2058D0, 380.0852D0, 0.0000D0, & 0.0000D0/ data C6/ -8.4082D0, 74.7925D0, -282.9540D0, 576.3930D0, & -661.9390D0, 404.2930D0, -102.2330D0, 0.0000D0, & 0.0000D0/ AA=0D0 BB=0D0 CC=0D0 do 100 i = 0, 8 AA(2) = AA(2) + A2(i)*tau**(i/8.d0) AA(1) = AA(1) + A1(i)*tau**(i/8.d0) AA(0) = AA(0) + A0(i)*tau**(i/8.d0) BB(1) = BB(1) + B1(i)*tau**(i/8.d0) BB(0) = BB(0) + B0(i)*tau**(i/8.d0) CC(6) = CC(6) + C6(i)*tau**(i/6.d0) CC(5) = CC(5) + C5(i)*tau**(i/6.d0) CC(4) = CC(4) + C4(i)*tau**(i/6.d0) CC(3) = CC(3) + C3(i)*tau**(i/6.d0) CC(2) = CC(2) + C2(i)*tau**(i/6.d0) 100 continue G0=0D0 G1=0D0 G2=0D0 do 200 Xi = 0D0, 2D0 i=idint(Xi) do 210 Xj = 2D0, 6D0 j=idint(Xj) Z=Xi+Xj/8D0+1D0 CALL GAMMAFUNC(Z,GAM) G1=G1+GAM *AA(i)*CC(j) G2=G2+GAM/Z*BB(i)*CC(j) 210 continue 200 continue G0=2D0*AA(2)+AA(1)+AA(0)+0.5D0*BB(1)+BB(0) GII=G0+G1+G2 return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Function GIII(tau) returns the fitting function of Wee(x,tau). C C GIII(tau) = 2D0*AA(2)+AA(1)+AA(0)+0.5D0*BB(1)+BB(0) C C AA(2) = A2(0) + A2(1)*tau^(1/8) + A2(2)*tau^(2/8) + A2(3)*tau^(3/8) C + A2(4)*tau^(4/8) + A2(5)*tau^(5/8) + A2(6)*tau^(6/8) C + A2(7)*tau^(7/8) + A2(8)*tau^(8/8) C C Same function forms for AA(1), AA(0), BB(1) and BB(0). C C inputs: 300keV < tau < 7MeV (tau = kBT/mc^2) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC real function GIII*8(tau) implicit real*8 (a-h, o-z) real*8 A2(0:8), A1(0:8), A0(0:8), B1(0:8), B0(0:8) real*8 AA(0:2), BB(0:1) data A2/ 64.7512D0, -213.8956D0, 174.1432D0, 136.5088D0, & -271.4899D0, 89.3210D0, 58.2584D0, -46.0807D0, & 8.7301D0/ data A1/ 49.7139D0, -189.7746D0, 271.0298D0, -269.7807D0, & 420.4812D0, -576.6247D0, 432.7790D0, -160.5365D0, & 23.3925D0/ data A0/ 52.1633D0, -257.0313D0, 446.8161D0, -293.0585D0, & 0.0000D0, 77.0474D0, -23.8718D0, 0.0000D0, & 0.1997D0/ data B1/ 376.4322D0, -1223.3635D0, 628.6787D0, 2237.3946D0, & -3828.8387D0, 2121.7933D0, -55.1667D0, -349.4321D0, & 92.2059D0/ data B0/ -8.5862D0, 34.1348D0, -116.3287D0, 296.5451D0, & -393.4207D0, 237.5497D0, -30.6000D0, -27.6170D0, & 8.8453D0/ AA(2)=0D0 AA(1)=0D0 AA(0)=0D0 BB(1)=0D0 BB(0)=0D0 do 100 i = 0, 8 AA(2) = AA(2) + A2(i)*tau**(i/8.d0) AA(1) = AA(1) + A1(i)*tau**(i/8.d0) AA(0) = AA(0) + A0(i)*tau**(i/8.d0) BB(1) = BB(1) + B1(i)*tau**(i/8.d0) BB(0) = BB(0) + B0(i)*tau**(i/8.d0) 100 continue G0=2D0*AA(2)+AA(1)+AA(0)+0.5D0*BB(1)+BB(0) GIII=G0 return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Function GIV(tau) returns C the extreme-relativistic approximation of Wee(x,tau). C C inputs: 7MeV < tau (tau = kBT/mc^2) C CThis extreme-relativistic approximation was obtained by Alexanian (1968). C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC REAL FUNCTION GIV*8(TAU) IMPLICIT REAL*8 (A-Z) data pi, emass/ 3.141592654d0, 510.998902d0/ GUM=DEXP(0.5772156649D0) CON1=9D0/PI/DSQRT(TAU) Gt=DLOG(2D0*TAU/GUM)+5D0/4D0 GIV=CON1*Gt return end c----------------------------------------------------------------------c c c c gamma function (input) c c c c----------------------------------------------------------------------c SUBROUTINE GAMMAFUNC(Z,gamma) REAL*8 Z,gamma INTEGER i CALL re_gamma_init CALL re_gamma_func(gamma , z ) RETURN END subroutine re_gamma_init c----------------------------------------------------------------------c c c c gamma function (initialization) c c c c----------------------------------------------------------------------c c implicit none c* c* real*8 b2n , pi integer*4 ic c* real*8 cc , ln2pi2, pi2 integer*4 n_08 parameter (n_08 = 8) common / re_gam / cc(n_08), ln2pi2, pi2 c* c* dimension b2n(n_08) c* c* c* ----- setting the Bernoulli´s number ----- c* b2n(1) = 1.d0/ 6.d0 b2n(2) = -1.d0/ 30.d0 b2n(3) = 1.d0/ 42.d0 b2n(4) = -1.d0/ 30.d0 b2n(5) = 5.d0/ 66.d0 b2n(6) = -691.d0/2730.d0 b2n(7) = 7.d0/ 6.d0 b2n(8) = -3617.d0/ 510.d0 c* c* do ic = 1, n_08 cc(ic) = b2n(ic)/dble(2*ic*(2*ic-1)) end do c* pi = acos(-1.d0) pi2 = 2.d0*pi ln2pi2 = log(pi2)/2.d0 c* c* return end subroutine re_gamma_func(gamma , z ) c----------------------------------------------------------------------c c c c gamma function c c c c----------------------------------------------------------------------c c implicit none c* c* real*8 gamma , z real*8 loggam, zz , norm , one c* c* one = 1.d0 norm = one zz = z c* c* ----- setting Re(zz) > 10 ----- call re_gamma_norm (zz , norm , one ) c* c* ----- setting log(gamma) ----- call re_gamma_log (loggam, zz ) c* c* ----- setting gamma(z) ----- call re_gamma_value(gamma , loggam, norm ) c* c* return end subroutine re_gamma_norm(zz , norm , one ) c----------------------------------------------------------------------c c c c normalization factor c c c c----------------------------------------------------------------------c c implicit none c* c* real*8 zz , norm , one real*8 rr , cut_vl parameter (cut_vl = 10.d0) integer*4 icount, mxloop parameter (mxloop = 1250) c* c* icount = 0 10 continue icount = icount + 1 rr = zz if (rr .lt. cut_vl) 1 then norm = norm/zz zz = zz + one if (icount .lt. mxloop) go to 10 write(6,9000) icount, zz stop 99 end if c* c* return 9000 format('program error at gamma_norm(icount, zz) :',/, 1 i6,1p6e12.4) end subroutine re_gamma_log(ln , z ) c----------------------------------------------------------------------c c c c log(gamma) c c c c----------------------------------------------------------------------c c implicit none c* c* real*8 ln , z real*8 z1 , z2 , v1 , v2 c* real*8 cc , ln2pi2, pi2 integer*4 n_08 parameter (n_08 = 8) common / re_gam / cc(n_08), ln2pi2, pi2 c* c* z1 = 1.d0/z z2 = z1*z1 v2 = (((((((cc(8)*z2 + cc(7))*z2 + cc(6))*z2 + cc(5))*z2 1 + cc(4))*z2 + cc(3))*z2 + cc(2))*z2 + cc(1))*z1 v1 = (z - 0.5d0)*log(z) - z + ln2pi2 ln = v1 + v2 c* c* return end subroutine re_gamma_value(gm , ln , norm ) c----------------------------------------------------------------------c c c c log(gamma) c c c c----------------------------------------------------------------------c c implicit none c* c* real*8 gm , ln , norm real*8 z1 , z2 c* real*8 cc , ln2pi2, pi2 integer*4 n_08 parameter (n_08 = 8) common / re_gam / cc(n_08), ln2pi2, pi2 c* c* z1 = ln z2 = exp(z1) gm = norm*z2 c* c* return end subroutine re_log_gamma_func(gamma , z ) c----------------------------------------------------------------------c c c c gamma function c c c c----------------------------------------------------------------------c c implicit none c* c* real*8 gamma , z real*8 loggam, zz , norm , one c* c* one = 1.d0 norm = one zz = z c* c* ----- setting Re(zz) > 10 ----- call re_gamma_norm (zz , norm , one ) c* c* ----- setting log(gamma) ----- call re_gamma_log (loggam, zz ) c* c* ----- setting gamma(z) ----- call re_log_gamma_value(gamma , loggam, norm ) c* c* return end subroutine re_log_gamma_value(gm , ln , norm ) c----------------------------------------------------------------------c c c c log(gamma) c c c c----------------------------------------------------------------------c c implicit none c* c* real*8 gm , ln , norm c* c* gm = ln + log(norm) c* c* return end