c........ c........ Compton scattering in gas of relativistic thermal electrons. c........ c........ Sample main program for the computation of: c........ 1. angle-averaged Compton redistribution function (fort.3) c........ 2. Compton redistribution function for photons scattered at c........ the right angle (output: fort.4) c........ c........ Algorithms displayed in both output files: 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........ implicit none integer i,idiv,mgi,in,j parameter ( idiv=161, mgi=2000 ) double precision smit(mgi),agt(mgi),xprime(idiv),step,xlow,xhigh double precision x,xdim,z,eps,eps1,temper,slambda,steplg double precision profil,redist1(idiv),redist2(idiv),redist3(idiv) double precision hhh,skkk,ccc,smhy,smel,esu,eee,sigma,pi,eve,evf double precision total1,total2,total3,total4,costh,redist4(idiv) parameter (hhh= 6.62620d-27, skkk= 1.38062d-16, ccc= 2.99793d+10, * smhy= 1.67333d-24, smel= 9.10956d-28, esu= 4.80325d-10, * eee= 4.80325d-10, sigma=5.66961d-05, pi=3.141592654d0) c........ eV - associated energy erg parameter ( eve=1.60219d-12, evf=2.41838d+14 ) c........ step=1/dble(mgi/2.d0) do 35 i=1,mgi agt(i)=step smit(i)= -1.d0-0.5d0/dble(mgi/2.d0)+step*dble(i) 35 continue c........ slambda=1.0d0 ! in Angstr wavelength temper=1.8d7 ! in K temperature steplg=0.001d0 ! auxiliary free parameter xlow=0.1d0 ! auxiliary free parameter xdim=1.d+08*ccc/slambda ! in Hz frequency costh=0.d0 ! sample cosine of scattering angle c read (1,*) slambda,temper,xlow,costh xdim=1.d+08*ccc/slambda ! in Hz frequency c........ x=smel*ccc**2/skkk/temper eps=hhh*xdim/(smel*ccc**2) c....... c....... Computation of the redistribution probability (\nu,\nu^\prime,x) c....... for the initial frequency \nu. Inverted temperature of electrons c....... x = m_e c**2 / kT write (3,'(1hc/1hc,6x,40hCOMPTON FREQUENCY REDISTRIBUTION PROFILE/ . 1hc,4x,6hTEMP =,1pe10.3,13h K and FREQ =,1pe12.5,3h Hz,15x, . 5hmgi =,i7/1hc)') temper,xdim,mgi write (3,'(1hc,15x,10hFREQ PRIME,4x,11hPROBABILITY)') write (3,'(1hc,17x,5hin Hz,10x,6hfor Hz/1hc)') write (3,'(1hc,32x,5hexact,9x,8hGuilbert,7x,6happrox,9x,7hSazonov * )') c....... Integration over final \nu^\prime and outgoing angles: do 81 in=1,idiv xprime(in)=dlog10(xdim)-xlow+steplg*(in-1) xprime(in)=10.d0**xprime(in) eps1=hhh*xprime(in)/(smel*ccc**2) redist1(in)=0.d0 redist2(in)=0.d0 redist3(in)=0.d0 redist4(in)=0.d0 do 82 j=1,mgi redist1(in) = redist1(in) + agt(j)*profil(1,eps,eps1,smit(j),x) redist2(in) = redist2(in) + agt(j)*profil(2,eps,eps1,smit(j),x) redist3(in) = redist3(in) + agt(j)*profil(3,eps,eps1,smit(j),x) redist4(in) = redist4(in) + agt(j)*profil(4,eps,eps1,smit(j),x) 82 continue redist1(in) = 2.d0*pi*redist1(in)*eps/xdim redist2(in) = 2.d0*pi*redist2(in)*eps/xdim redist3(in) = 2.d0*pi*redist3(in)*eps/xdim redist4(in) = 2.d0*pi*redist4(in)*eps/xdim write (3,'(1hc,i10,1p5e15.5)') in,xprime(in),redist1(in), * redist2(in),redist3(in),redist4(in) 81 continue c........ total1=0. total2=0. total3=0. total4=0. do 85 in=1,idiv if (in.eq.1) step=(xprime(2)-xprime(1))/2.d0 if (in.eq.idiv) step=(xprime(idiv)-xprime(idiv-1))/2.d0 if (in.gt.1.and.in.lt.idiv) * step=(xprime(in+1)-xprime(in-1))/2.d0 total1=total1+redist1(in)*step total2=total2+redist2(in)*step total3=total3+redist3(in)*step total4=total4+redist4(in)*step 85 continue write (3,'(1x/25x,0p4f15.5)') total1,total2,total3,total4 c....... c....... Computation of the angle-dependent redistribution probability c....... (\nu,\nu^\prime,costh,x) for the initial frequency \nu and c....... cosine of the scattering angle 'costh'. Inverted temperature c....... of electrons c....... x = m_e c**2 / kT write (4,'(1hc/1hc,6x,40hANGLE-DEPENDENT FREQUENCY REDISTRIBUTION * 8h PROFILE/1hc,4x,6hTEMP =,1pe10.3,10h K, FREQ =,1pe12.5,3h Hz, * 12h and COSTH =,0pf8.4,8x,5hmgi =,i7/1hc)') temper,xdim,costh, * mgi write (4,'(1hc,15x,10hFREQ PRIME,4x,11hPROBABILITY)') write (4,'(1hc,17x,5hin Hz,10x,6hfor Hz/1hc)') write (4,'(1hc,32x,5hexact,9x,8hGuilbert,7x,6happrox,9x,7hSazonov * )') c........ Scattering profile for a single scattering cosine: do 71 in=1,idiv xprime(in)=dlog10(xdim)-xlow+steplg*(in-1) xprime(in)=10.d0**xprime(in) eps1=hhh*xprime(in)/(smel*ccc**2) redist1(in)=0.5d0*profil(1,eps,eps1,costh,x)*eps/xdim redist2(in)=0.5d0*profil(2,eps,eps1,costh,x)*eps/xdim redist3(in)=0.5d0*profil(3,eps,eps1,costh,x)*eps/xdim redist4(in)=0.5d0*profil(4,eps,eps1,costh,x)*eps/xdim write (4,'(1hc,i10,1p5e15.5)') in,xprime(in),redist1(in), * redist2(in),redist3(in),redist4(in) 71 continue c........ stop end