PROGRAM Generate_Table ! ! Integrate a dimensionless generalized NFW profile to relate the integral, f(x), ! to a dimensionless angular radius, x=theta/theta500. This is needed to find ! the integral value given theta [i.e., inverse function of f(x)]. ! ! Output file "x_fx.txt" contains ! 1st: x ! 2nd: f(x), where x=theta/theta500 ! ! May 17, 2016: E.Komatsu IMPLICIT none double precision :: xin=0d0,xout=6d0 ! pressure profile is cut at r=xout*r500 double precision :: y,fx,tol=1d-6 double precision :: pgnfw,rombint integer :: i external pgnfw,rombint ! calculate profiles open(1,file='x_fx.txt') do i=1,100000 y=xin + 1d-5*dble(i)*(xout-xin) ! y=theta/theta500 fx=rombint(pgnfw,-dsqrt(xout**2d0-y**2d0),dsqrt(xout**2d0-y**2d0),tol,y) write(1,'(2E18.5)')y,fx enddo close(1) END PROGRAM Generate_Table !----------------------------------------------------------------------------- DOUBLE PRECISION FUNCTION pgnfw(l,y) ! x=r/r500 IMPLICIT none double precision :: l,x,y double precision :: c=1.81d0,g=0.31d0,a=1.33d0,b=4.13d0,P0=6.41d0 ! Planck 2013 x=dsqrt(l**2d0+y**2d0) pgnfw=P0/(c*x)**g/(1d0+(c*x)**a)**((b-g)/a) return END FUNCTION pgnfw