CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C This program calculates the photon production rate Pee(k,tau) C for the electron-electron thermal bremsstrahlung defined as C follows. C d Nee C Pee(k,tau) = -------- C dt dV dk C is the number of photons emitted per unit time, per unit volume, C per unit unit energy interval by electron gas of uniform number C density ne at temperature T, where k is the photon energy in C units of the electron rest mass and tau =kBT/mc^2 is the electron C thermal energy in units of the electron rest mass. C This program returns the analytic fitting function Pee^fit(x,tau) C defined by eq. (4.1) 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 two dimensionless variables x and tau. C x = k/tau ----- real*8, 10^-4 < x < 10 C tau = kBT/mc^2 ----- real*8, kBT > 50eV C C OUTPUTS C This program returns the photon production rate Pee. C Pee ----- real*8 C Pee = 1.455x10^{-16} [ne(cm^-3)]^2 exp(-x)/x*1/sqrt(tau) C *Gfit(x,tau) cm^-3 sec^-1 C Please refer to eq. (4.1) in Ref.1 for detail. C C 2009.01 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PROGRAM MAIN implicit real*8 (a-h, o-z) data emass/ 510.998902d0/ write(*,*)'----------------------------------------------------------' write(*,*)'Inputs:' write(*,*)'kBT (kBT > 50eV): the electron thermal energy' write(*,*)' x = k/tau (10^-4 < x < 10) : k is the photon energy' write(*,*)' in units of the electron rest mass' write(*,*)'----------------------------------------------------------' 100 write(*,*)'kBT = ? kBT=0 will quit this program.' read(*,*) T if (T .ne. 0.d0) then tau = T/emass write(*,*)'x = ?' read(*,*) x CALL FUNCFx(X,TAU,Peefit) write(*,*)'x, Peefit(x,tau) =', x, Peefit write(*,*)'------------------------------------------------------' GOTO 100 else end if END SUBROUTINE FUNCFx(X,TAU,Peefit) 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(x,tau) ELSE IF (TAU.ge.(1D0/emass).and.TAU.lt.(300D0/emass)) THEN G=GII(x, tau)*FCC(x, tau) ELSE IF (TAU.ge.(300D0/emass).and.TAU.le.(7000D0/emass)) THEN G=GIII(x, tau) ELSE IF (TAU.gt.(7000D0/emass)) THEN G=GIV(x, tau) END IF Peefit=1.455D-16*DEXP(-X)/X/DSQRT(TAU)*G CCC Peefit=DEXP(-X)/X/DSQRT(TAU)*G RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Function GI(x, tau) returns the fitting function of Pee(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(X,TAU) implicit real*8(a-z) real a(1:121) data a/ $ 3.15847D0 , 2.46819D-2, -2.11118D-2, 1.24009D-2, $ -5.41633D-3, 1.70070D-3, -3.05111D-4, -1.21721D-4, $ 1.77611D-4, -2.05480D-5, -3.58754D-5, $ -2.52430D0 , 1.03924D-1, -8.53821D-2, 4.73623D-2, $ -1.91406D-2, 5.39773D-3, -7.26681D-4, -7.47266D-4, $ 8.73517D-4, -6.92284D-5, -1.80305D-4, $ 4.04877D-1, 1.98935D-1, -1.52444D-1, 7.51656D-2, $ -2.58034D-2, 4.13361D-3, 4.67015D-3, -2.20675D-3, $ -2.67582D-3, 2.95254D-5, 1.40751D-3, $ 6.13466D-1, 2.18843D-1, -1.45660D-1, 5.07201D-2, $ -2.23048D-3, -1.14273D-2, 1.24789D-2, -2.74351D-3, $ -4.57871D-3, -1.70374D-4, 2.06757D-3, $ 6.28867D-1, 1.20482D-1, -4.63705D-2, -2.25247D-2, $ 5.07325D-2, -3.23280D-2, -1.16976D-2, -1.00402D-3, $ 2.96622D-2, -5.43191D-4, -1.23098D-2, $ 3.29441D-1, -4.82390D-2, 8.16592D-2, -8.17151D-2, $ 5.94414D-2, -2.19399D-2, -1.13488D-2, -2.38863D-3, $ 1.89850D-2, 2.50978D-3, -8.81767D-3, $ -1.71486D-1, -1.20811D-1, 9.87296D-2, -4.59297D-2, $ -2.11247D-2, 1.76310D-2, 6.31446D-2, -2.28987D-3, $ -8.84093D-2, 4.45570D-3, 3.46210D-2, $ -3.68685D-1, -4.46133D-4, -3.24743D-2, 5.05096D-2, $ -5.05387D-2, 2.23352D-2, 1.33830D-2, 7.79323D-3, $ -2.93629D-2, -2.80083D-3, 1.23727D-2, $ -7.59200D-2, 8.88749D-2, -8.82637D-2, 5.58818D-2, $ 9.20453D-3, -4.59817D-3, -8.54735D-2, 7.98332D-3, $ 1.02966D-1, -5.68093D-3, -4.04801D-2, $ 1.60187D-1, 2.50320D-2, -7.52221D-3, -9.11885D-3, $ 1.67321D-2, -8.24286D-3, -6.47349D-3, -3.80435D-3, $ 1.38957D-2, 1.10618D-3, -5.68689D-3, $ 8.37729D-2, -1.28900D-2, 1.99419D-2, -1.71348D-2, $ -3.47663D-3, -3.90032D-4, 3.72266D-2, -4.25035D-3, $ -4.22093D-2, 2.33625D-3, 1.66733D-2/ data pi, emass/ 3.141592654d0, 510.998902d0/ logth=LOG10(TAU) log u=LOG10(X) the = (logth+2.65)/1.35 u = (logu+1.50)/2.50 jfit =(A(1)+A(2)*the+A(3)*the**2+A(4)*the**3+A(5)*the**4 * +A(6)*the**5+A(7)*the**6+A(8)*the**7 * +A(9)*the**8+A(10)*the**9+A(11)*the**10) * +(A(12)+A(13)*the+A(14)*the**2+A(15)*the**3+A(16)*the**4 * +A(17)*the**5+A(18)*the**6+A(19)*the**7 * +A(20)*the**8+A(21)*the**9+A(22)*the**10)*u * +(A(23)+A(24)*the+A(25)*the**2+A(26)*the**3+A(27)*the**4 * +A(28)*the**5+A(29)*the**6+A(30)*the**7 * +A(31)*the**8+A(32)*the**9+A(33)*the**10)*u**(+2) * +(A(34)+A(35)*the+A(36)*the**2+A(37)*the**3+A(38)*the**4 * +A(39)*the**5+A(40)*the**6+A(41)*the**7 * +A(42)*the**8+A(43)*the**9+A(44)*the**10)*u**(+3) * +(A(45)+A(46)*the+A(47)*the**2+A(48)*the**3+A(49)*the**4 * +A(50)*the**5+A(51)*the**6+A(52)*the**7 * +A(53)*the**8+A(54)*the**9+A(55)*the**10)*u**(+4) * +(A(56)+A(57)*the+A(58)*the**2+A(59)*the**3+A(60)*the**4 * +A(61)*the**5+A(62)*the**6+A(63)*the**7 * +A(64)*the**8+A(65)*the**9+A(66)*the**10)*u**(+5) * +(A(67)+A(68)*the+A(69)*the**2+A(70)*the**3+A(71)*the**4 * +A(72)*the**5+A(73)*the**6+A(74)*the**7 * +A(75)*the**8+A(76)*the**9+A(77)*the**10)*u**(+6) * +(A(78)+A(79)*the+A(80)*the**2+A(81)*the**3+A(82)*the**4 * +A(83)*the**5+A(84)*the**6+A(85)*the**7 * +A(86)*the**8+A(87)*the**9+A(88)*the**10)*u**(+7) * +(A(89)+A(90)*the+A(91)*the**2+A(92)*the**3+A(93)*the**4 * +A(94)*the**5+A(95)*the**6+A(96)*the**7 * +A(97)*the**8+A(98)*the**9+A(99)*the**10)*u**(+8) * +(A(100)+A(101)*the+A(102)*the**2+A(103)*the**3+A(104)*the**4 * +A(105)*the**5+A(106)*the**6+A(107)*the**7 * +A(108)*the**8+A(109)*the**9+A(110)*the**10)*u**(+9) * +(A(111)+A(112)*the+A(113)*the**2+A(114)*the**3+A(115)*the**4 * +A(116)*the**5+A(117)*the**6+A(118)*the**7 * +A(119)*the**8+A(120)*the**9+A(121)*the**10)*u**(+10) GI=DSQRT(8D0/(3D0*PI))*jfit return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Function GII(x, tau) returns the fitting function of Pee(x,tau). C C GII(x, tau) = AA2*x^2 + AA1*x + AA0*x C - EXP(x) * Ei(-x) * ( BB1*x + BB0*x ) 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(x, 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) data pi, emass/ 3.141592654d0, 510.998902d0/ 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/ AA2=0D0 AA1=0D0 AA0=0D0 BB1=0D0 BB0=0D0 CALL FUNCE(X,Ei) do 10 i = 0, 8 AA2 = AA2 + A2(i)*tau**(i/8.d0) AA1 = AA1 + A1(i)*tau**(i/8.d0) AA0 = AA0 + A0(i)*tau**(i/8.d0) BB1 = BB1 + B1(i)*tau**(i/8.d0) BB0 = BB0 + B0(i)*tau**(i/8.d0) 10 continue GII=AA2*X**2+AA1*X+AA0-DEXP(X)*Ei*(BB1*X+BB0) return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Function GIII(x, tau) returns the fitting function of Pee(x,tau). C C GIII(x, tau) = AA2*x^2 + AA1*x + AA0*x C - EXP(x) * Ei(-x) * ( BB1*x + BB0*x ) 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: 300keV < tau < 7MeV (tau = kBT/mc^2) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC real function GIII*8(x, 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) data pi, emass/ 3.141592654d0, 510.998902d0/ 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/ AA2=0D0 AA1=0D0 AA0=0D0 BB1=0D0 BB0=0D0 CALL FUNCE(X,Ei) do 10 i = 0, 8 AA2 = AA2 + A2(i)*tau**(i/8.d0) AA1 = AA1 + A1(i)*tau**(i/8.d0) AA0 = AA0 + A0(i)*tau**(i/8.d0) BB1 = BB1 + B1(i)*tau**(i/8.d0) BB0 = BB0 + B0(i)*tau**(i/8.d0) 10 continue GIII=AA2*X**2+AA1*X+AA0-DEXP(X)*Ei*(BB1*X+BB0) return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Function GIV(x, tau) returns C the extreme-relativistic approximation of Pee(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(X,TAU) IMPLICIT REAL*8 (A-Z) data pi, emass/ 3.141592654d0, 510.998902d0/ GAM=DEXP(0.5772156649D0) CON1=3D0/(4D0*PI*SQRT(TAU)) A1=28D0/3D0+2D0*X +X**2/2D0 A2= 8D0/3D0+4D0/3D0*X+X**2 A3= 8D0/3D0-4D0/3D0*X+X**2 CALL FUNCE(X,Ei) Fx=A1+2D0*A2*DLOG(2D0*TAU/GAM)-DEXP(X)*Ei*A3 GIV=CON1*Fx RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Function FCC(x, tau) returns the fitting function of C the Colomb correction with the use of the Elwert factor. C C FCC(x, tau) = 1 + p(tau)*x^(2/8) + q(tau)*x^(3/8) C + r(tau)*x^(4/8) + s(tau)*x^(5/8) C + t(tau)*x^(6/8) C C p(tau) = p0/sqrt(tau) + p1*tau^(1/6) + p2*tau^(2/6) + p3*tau^(3/6) C + p4*tau^(4/6) + p5*tau^(5/6) C C Same function forms for q(tau), r(tau), s(tau) and t(tau). C C inputs: 1keV < tau < 300keV (tau = kBT/mc^2) C 10^-4 < x < 10 (x = kc/tau) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC real function FCC*8(x, tau) implicit real*8 (a-h, o-z) real*8 p(0:6), q(0:6), r(0:6), s(0:6), t(0:6) data pi, emass/ 3.141592654d0, 510.998902d0/ data p/ -5.7752D0, 46.2097D0, -160.7280D0, * 305.0070D0, -329.5420D0, 191.0770D0, -46.2718D0/ data q/ 30.5586D0, -248.2177D0, 874.1964D0, * -1676.9028D0, 1828.8677D0, -1068.9366D0, 260.5656D0/ data r/ -54.3272D0, 450.9676D0, -1616.5987D0, * 3148.1061D0, -3478.3930D0, 2055.6693D0, -505.6789D0/ data s/ 36.2625D0, -310.0972D0, 1138.0531D0, * -2260.8347D0, 2541.9361D0, -1525.2058D0, 380.0852D0/ data t/ -8.4082D0, 74.7925D0, -282.9540D0, * 576.3930D0, -661.9390D0, 404.2930D0, -102.2330D0/ ptau =0D0 qtau =0D0 rtau =0D0 stau =0D0 ttau =0D0 do 10 i = 0, 6 ptau = ptau + p(i)*tau**(i/6.d0) qtau = qtau + q(i)*tau**(i/6.d0) rtau = rtau + r(i)*tau**(i/6.d0) stau = stau + s(i)*tau**(i/6.d0) ttau = ttau + t(i)*tau**(i/6.d0) 10 continue FCC = 1.d0 + ptau*x**(2.d0/8.d0) + qtau*x**(3.d0/8.d0) * + rtau*x**(4.d0/8.d0) + stau*x**(5.d0/8.d0) * + ttau*x**(6.d0/8.d0) return end SUBROUTINE FUNCF(Ifin,N) IMPLICIT NONE REAL*8 N,I,Ifin N=1D0 DO 10 I=1D0,Ifin N=N*I 10 CONTINUE RETURN END SUBROUTINE FUNCE(X,Ei) IMPLICIT NONE REAL*8 ANS,Ei,GAM REAL*8 T1,T2 REAL*8 A1,A2,A3 REAL*8 X,T,F REAL*8 M,N REAL*8 XT,WT INTEGER iNT,IXT DIMENSION XT(20),WT(20) DATA GAM / 0.57721566490153286060D0 / DATA XT /0.005396838D0,0.028303937D0,0.069040301D0,0.127034037D0, & 0.201908993D0,0.293890939D0,0.404249262D0,0.535818423D0, & 0.693715292D0,0.886478749D0,1.128058627D0,1.44151356D0 , & 1.866303878D0,2.473721274D0,3.402622768D0,4.952726389D0, & 7.871906009D0,14.48429371D0,35.33077428D0,185.2936858D0/ DATA WT /0.013834412D0,0.031913834D0,0.049457133D0,0.066460888D0, & 0.083319379D0,0.100847649D0,0.120329957D0,0.143643116D0, & 0.173542463D0,0.214251752D0,0.272638839D0,0.360614529D0, & 0.500321925D0,0.736334738D0,1.167598113D0,2.043782594D0, & 4.118375485D0,10.37584757D0,39.8368777D0 ,474.9872281D0/ ANS=0D0 A1=0D0 A3=0D0 IF (X.le.1D0) THEN A1=GAM A2=DLOG(X) DO 110 M=1D0,18D0 CALL FUNCF(M,N) A3=A3+(-X)**M/(N*M) 110 CONTINUE ANS=A1+A2+A3 Ei=ANS ELSE iNT=20 T1=0D0 T2=1D30 DO 100 IXT=1,iNT T=XT(IXT) F=DEXP(-T)/(T+X) ANS=ANS+F*WT(IXT) 100 CONTINUE Ei=-DEXP(-X)*ANS END IF RETURN END