c....... nr=1 : exact by Suleimanov et al. (2012), Madej et al. (2017) c....... nr=2 : Guilbert (1981), Madej et al. (2017) c....... nr=3 : approximate by Arutyunyan \& Nikogosyan (1980) c....... nr=4 : approximate by Sazonov \& Sunyaev (2000) c....... double precision function profil(nr,eps,eps1,costh,x) implicit none c....... Computation of the probability of scattering a photon, c....... P (\epsilon, s, costh, x), c....... interacting with free electrons in thermal gas. c....... c....... \epsilon = initial energy (in units of mc**2), c....... \epsilon_1 = final energy of a scattered photon, c....... costh = cosine of scattering angle, and c....... c....... s = \epsilon / \epsilon_1 c....... x = m c**2 / kT c....... c....... Thirty-two-point Gauss-Laguerre quadrature is used below. c....... integer i,j,nr,nsize double precision integ,eps,eps1,costh,x,ag,aj,bg,eg,eg2,gammin double precision pi,s,xj,z,bk2,zfg,zfe,zfs,fguilb,fexact,fsazon double precision gam_star,q_up,q_down double precision profil1,profil2,profil3,profil4 parameter ( nsize = 32 ) dimension xj(nsize),aj(nsize) data pi/3.141592654d0/ c....... tablica: XJ - to wezly, AJ - to wagi. data (xj(i),aj(i),i=1,nsize)/0.044489366d0,1.09218342d-1,0.2345261 *1d0,2.10443108d-1,0.576884629d0,2.352132297d-1,1.07244875d0,1.9590 *3336d-1,1.72240878d0,1.29983786d-1,2.528336706d0,7.05786239d-2,3.4 *92213273d0,3.17609125d-2,4.61645677d0,1.191821484d-2,5.9039585d0,3 *.738816295d-3,7.358126733d0,9.8080331d-4,8.98294092d0,2.14864919d- *4,10.78301863d0,3.920342d-5,12.76369799d0,5.9345416d-6,14.93113976 *d0,7.4164046d-7,17.29245434d0,7.6045679d-8,19.85586094d0,6.3506022 *d-9,22.63088901d0,4.28138297d-10,25.62863602d0,2.30589949d-11,28.8 *6210182d0,9.7993793d-13,32.34662915d0,3.2378017d-14,36.10049481d0, *8.17182344d-16,40.1457198d0,1.54213383d-17,44.509208d0,2.11979229d *-19,49.22439499d0,2.05442967d-21,54.3337213d0,1.34698259d-23,59.89 *250916d0,5.6612941d-26,65.97537729d0,1.41856055d-28,72.68762809d0, *1.91337549d-31,80.18744698d0,1.19224876d-34,88.73534042d0,2.671511 *22d-38,98.82954287d0,1.33861694d-42,111.7513981d0,4.5105362d-48/ c....... profil=0.d0 q_down=eps*eps1*(1-costh) q_up=dsqrt((eps-eps1)**2+2*q_down) gam_star=(eps-eps1+q_up*dsqrt(1+2.d0/q_down))/2.d0 integ=0.d0 c....... go to (11,12,13,14), nr return c....... 11 do 16 j=1,nsize zfe=fexact(eps,eps1,costh,xj(j)/x+gam_star) integ=integ+zfe*aj(j) 16 continue integ=integ*dexp((1-gam_star+eps-eps1)*x) profil2=3.d0/32.d0/pi/bk2(x)*integ *eps1/eps profil=profil2 return 12 s=eps/eps1 ag=1.d0-s bg=eps*(1.d0-costh) eg2=1.d0-2.d0*s*costh+s**2 eg=dsqrt(eg2) z=eg2*(1.d0-ag**2/(eg2+bg**2))/(eg2+bg**2) gammin=1.d0/(dsqrt(z)-ag*bg/(eg2+bg**2)) c....... Computations of the integral I (\epsilon, s, \ksi, x) : do 15 j=1,nsize zfg=fguilb(eps,eps1,costh,xj(j)/x+gammin) integ=integ+zfg*aj(j) 15 continue integ=integ*dexp((1.d0-gammin)*x) profil1=3.d0/64.d0/pi**2/eg/eps/bk2(x)*integ profil=profil1 return 13 profil3=dexp((1-gam_star+eps-eps1)*x)/bk2(x)/(8*pi*q_up) profil=profil3 *eps1/eps return 14 zfs=fsazon(eps,eps1,costh,x) profil4=3/(32*pi)*dsqrt(2/pi)*zfs /eps profil=profil4 return end c....... double precision function fsazon (eps,eps1,costh,x) implicit none double precision eps,eps1,costh,x,thetinv,sazon,sum c....... thetinv=1/x sazon=dsqrt(2.d0)*dsqrt(1-costh)*(eps1-eps+eps*eps1*(1-costh)) * /dsqrt(eps**2-2*eps*eps1*costh+eps1**2) sum=(1+costh**2+(1/8.d0-costh-63/8.d0*costh**2+5*costh**3)*thetinv * -costh*(1+costh)/2.d0*sazon**2-3*(1+costh**2)/32.d0/(1-costh)**2 * *sazon**4/thetinv+costh*(1-costh**2)*sazon*eps+(1+costh**2)/8.d0 * /(1-costh)*sazon**3/thetinv*eps+(1-costh)**2*eps*eps1) * *dexp(-sazon**2/4.d0/(1-costh)/thetinv) fsazon=sum*dsqrt(x)*eps1/dsqrt(eps**2-2*eps*eps1*costh+eps1**2) return end c....... double precision function fexact (eps,eps1,costh,gamma) implicit none double precision eps,eps1,costh,gamma,a2_plus,a2_minus, * q_up,q_down,d_plus,d_minus,a_minus,a_plus,a_diff,a2_diff c....... q_down=eps*eps1*(1-costh) q_up=dsqrt((eps-eps1)**2+2*q_down) a2_minus=(gamma-eps)**2+(1+costh)/(1-costh) a2_plus=(gamma+eps1)**2+(1+costh)/(1-costh) d_minus=(a2_plus-a2_minus-q_up**2)/2.d0 d_plus =(a2_plus-a2_minus+q_up**2)/2.d0 c....... a_minus=dsqrt(a2_minus) a_plus=dsqrt(a2_plus) a_diff=(-2*gamma*(eps+eps1)+eps**2-eps1**2)/(a_minus+a_plus) a2_diff=-2*gamma*(eps+eps1)+eps**2-eps1**2 c........ Final formula from Madej et al. (2017): fexact= 2/q_up-(q_down-2)/(q_down*a_plus*a_minus)*a_diff * +((a2_minus+a_minus*a_plus+a2_plus)*q_up**2-a2_diff**2 * -a_minus*a_plus*a_diff**2)/(2*q_down**2 * *a_plus**3*a_minus**3)*a_diff c........ Formula from Suleimanov et al. (2012): c fexact= 2/q_up + (q_down**2-2*q_down-2)/q_down**2 c * *(1/dsqrt(a2_minus)-1/dsqrt(a2_plus)) c * +(d_minus/a2_minus**1.5d0+d_plus/a2_plus**1.5d0)/q_down**2 return end c....... double precision function fguilb (eps,eps1,costh,gamma) implicit none real*8 int1,int2,int3 double precision a,eps,eps1,costh,gamma,b,c,d,pi,s,term,y, * sqab2,sqcd2 data pi/3.141592654d0/ c....... s=eps/eps1 a=1-(1-s+eps*(1-costh)/gamma)*(costh-s)/(1-2*s*costh+s**2) b=-dsqrt((gamma**2-1)/gamma**2-(1-s+eps*(1-costh)/gamma)**2 * /(1-2*s*costh+s**2) )*dsqrt((1-costh**2)/(1-2*s*costh+s**2)) c=1-(1-s+eps*(1-costh)/gamma)*(1-s*costh)/(1-2*s*costh+s**2) d=b*s c....... sqab2=dsqrt(a**2-b**2) sqcd2=dsqrt(c**2-d**2) y=b/sqab2+d/sqcd2 int1=2*pi/sqab2 int2=2*pi/y/(a**2-b**2)**2/(c**2-d**2)*(a*(b*c+a*d)+b*(a*c+b*d) * -a*b*(b*c+a*d)/y/sqab2) int3=pi/y/(a**2-b**2)**3/(c**2-d**2)**2 c....... term=(a*(b*c+a*d)+b*(a*c+b*d)-a*b*(b*c+a*d)/y/dsqrt(a**2-b**2)) * *(2*a*b*c*d/y**2/dsqrt(a**2-b**2)/dsqrt(c**2-d**2)-4*a*c*d/y * /dsqrt(c**2-d**2)-2*a*b*c/y/dsqrt(a**2-b**2)+8*a*c) * +(2*a*d+2*b*c-(b**2*c+2*a*b*d)/y/dsqrt(a**2-b**2) * +a**2*b*(b*c+a*d)/y/(a**2-b**2)**1.5-a**2*b**2 * *(b*c+a*d)/y**2/(a**2-b**2)**2)*c*(a**2-b**2) * *(d/y/dsqrt(c**2-d**2)-2) + (2*a*b-a*b**2/y/dsqrt(a**2-b**2) * -a*b*c*d*(b*c+a*d)/y**2/dsqrt(a**2-b**2)/(c**2-d**2)**1.5) * *(a*b*(c**2-d**2)/y/dsqrt(a**2-b**2)-4*a*(c**2-d**2)) * +2*b*(a**2-b**2)*(c**2-d**2)-b**2/y*dsqrt(a**2-b**2) * *(c**2-d**2)-a**2*b**3*(c**2-d**2)/y**2/(a**2-b**2) * +a**2*b**2*(c**2-d**2)/y/dsqrt(a**2-b**2) * -b*c*d*(b*c+2*a*d)*dsqrt(a**2-b**2)/y**2/dsqrt(c**2-d**2) * -2*a**2*b**2*c*d*(b*c+a*d)/y**3/(a**2-b**2)/dsqrt(c**2-d**2) * +a**2*b*c*d*(b*c+a*d)/y**2/dsqrt(a**2-b**2)/dsqrt(c**2-d**2) int3=int3*term fguilb=2*a*int1 - 2*a/gamma**2*(1-costh)*(1-eps**2*(1 * -costh)/2.d0/s)*int2 + a/gamma**2**2*(1-costh)**2*int3 return end c........ double precision function bk2 (x) implicit none double precision i0,i1,k0,k1,x,t,z c....... Computations of the modified Bessel function of the second c....... order, scaled by exponent c....... c....... K_2 (x) \times exp(x) - for X.GT.0, c....... c....... Method: asymptotic expansions, given by Abramowitz and Stegun c....... are used for computing I_0, I_1, K_0, and K_1. c....... c....... Final step makes use of recurrence relation: c....... K_{n+1} (x) = {2n/x} K_n (x) + K_{n-1} (x) c....... c....... The function has been carefully tested in the range c....... X: between 0.01 and 500. c....... if (2.d0-x.ge.0.d0) go to 42 z=2.d0/x k0 = (((((0.00053208d0*z-0.00251540d0)*z+0.00587872d0)*z * -0.01062446d0)*z+0.02189568d0)*z-0.07832358d0)*z * +1.25331414d0 k1 = (((((-0.00068245d0*z+0.00325614d0)*z-0.00780353d0)*z * +0.01504268d0)*z-0.03655620d0)*z+0.23498619d0)*z * +1.25331414d0 bk2 = (k0+z*k1)/dsqrt(x) return 42 z=(x/2.d0)**2 t=(x/3.75d0)**2 i0 = (((((0.0045813d0*t+0.0360768d0)*t+0.2659732d0)*t * +1.2067492d0)*t+3.0899424d0)*t+3.5156229d0)*t + 1.d0 i1 = (((((.00032411d0*t+.00301532d0)*t+.02658733d0)*t * +.15084934d0)*t+.51498869d0)*t+.87890594d0)*t + 0.5d0 i1=x*i1 k0 = (((((.00000740d0*z+.00010750d0)*z+.00262698d0)*z * +.03488590d0)*z+.23069756d0)*z+.42278420d0)*z * -.57721566d0 - i0*dlog(x/2.d0) k1 = (((((-.00004686d0*z-.00110404d0)*z-.01919402d0)*z * -.18156897d0)*z-.67278579d0)*z+.15443144d0)*z + 1.d0 k1 = k1/x + i1*dlog(x/2.d0) c....... K_{n+1} (x) = {2n/x} K_n (x) + K_{n-1} (x) z = k0 + 2.d0/x*k1 bk2 = z * dexp(x) return end