ypdfdks/ 000755 000767 000024 00000000000 12726614445 015101 5 ustar 00eiichirokomatsu staff 000000 000000 ypdfdks/chebyshev.f 000644 000767 000024 00000003022 12577165440 017225 0 ustar 00eiichirokomatsu staff 000000 000000 SUBROUTINE chebft(a,b,c,n,func) INTEGER n,NMAX REAL a,b,c(n),func,PI EXTERNAL func PARAMETER (NMAX=50, PI=3.141592653589793d0) INTEGER j,k REAL bma,bpa,fac,y,f(NMAX) DOUBLE PRECISION sum bma=0.5*(b-a) bpa=0.5*(b+a) do 11 k=1,n y=cos(PI*(k-0.5)/n) f(k)=func(y*bma+bpa) 11 continue fac=2./n do 13 j=1,n sum=0.d0 do 12 k=1,n sum=sum+f(k)*cos((PI*(j-1))*((k-0.5d0)/n)) 12 continue c(j)=fac*sum 13 continue return END C (C) Copr. 1986-92 Numerical Recipes Software D041&0(9p#3. SUBROUTINE chder(a,b,c,cder,n) INTEGER n REAL a,b,c(n),cder(n) INTEGER j REAL con cder(n)=0. cder(n-1)=2*(n-1)*c(n) if(n.ge.3)then do 11 j=n-2,1,-1 cder(j)=cder(j+2)+2*j*c(j+1) 11 continue endif con=2./(b-a) do 12 j=1,n cder(j)=cder(j)*con 12 continue return END C (C) Copr. 1986-92 Numerical Recipes Software D041&0(9p#3. FUNCTION chebev(a,b,c,m,x) INTEGER m REAL chebev,a,b,x,c(m) INTEGER j REAL d,dd,sv,y,y2 if ((x-a)*(x-b).gt.0.) then write(*,*) a,b,x pause 'x not in range in chebev' endif d=0. dd=0. y=(2.*x-a-b)/(b-a) y2=2.*y do 11 j=m,2,-1 sv=d d=y2*d-dd+c(j) dd=sv 11 continue chebev=y*d-dd+0.5*c(1) return END C (C) Copr. 1986-92 Numerical Recipes Software D041&0(9p#3. ypdfdks/cosmo.f90 000644 000767 000024 00000000143 12577165440 016537 0 ustar 00eiichirokomatsu staff 000000 000000 MODULE cosmo ! cosmological parameters double precision :: om0,ode0,w,obh2,h0 END MODULE cosmo ypdfdks/da.f90 000644 000767 000024 00000005313 12577165440 016007 0 ustar 00eiichirokomatsu staff 000000 000000 !----------------------------------------------------------------------------- ! ! Proper angular diameter distance, returned in units of h^-1 Mpc. ! Note h^-1! ! ! *** Radiation density is ignored in the calculation, so ! this routine cannot be used to calculate the distance to ! the last scattering surface! *** ! ! << sample code >> ! ! USE cosmo ! USE angular_distance ! double precision :: da,redshift ! external da ! ! Specify three cosmological parameters in double precision. ! ! The data type has been defined in MODULE cosmo. ! ode0=0.723d0 ! om0=0.277d0 ! w = -1d0 ! CALL setup_da ! redshift=10d0 ! print*,'omega matter=',om0 ! print*,'omega de=',ode0 ! print*,'w=',w ! print*,'Da(z) at z=',redshift,'is',da(redshift),' h^-1 Mpc' ! end ! ! August 23, 2008 ! E. Komatsu ! !----------------------------------------------------------------------------- MODULE angular_distance INTEGER, parameter :: nw=501 DOUBLE PRECISION, dimension(nw) :: xa,za,z2a contains SUBROUTINE Setup_Da USE cosmo ! tabulate da(x), where x=ln(a) IMPLICIT none INTEGER :: i DOUBLE PRECISION :: x,eta,tol=1d-7,rombint,ok0 EXTERNAL one_over_h ok0=1d0-om0-ode0 do i=1,nw x = (i-nw)*0.01 ! x=ln(a)=-5.00,-4.99,...,0 xa(i) = x eta=rombint(one_over_h,0d0,dexp(-x)-1d0,tol) if(abs(ok0)<1.d-5)then za(i)=2998d0*exp(x)*eta ! flat model else if(ok0>0)then za(i)=2998d0*exp(x)*sinh(sqrt(ok0)*eta)/sqrt(ok0) ! open model else za(i)=2998d0*exp(x)*sin(sqrt(-ok0)*eta)/sqrt(-ok0) ! closed model endif enddo CALL spline(xa,za,nw,1.d30,1.d30,z2a) ! evaluate z2a(i)=d^2[da(i)]/dy^2 return END SUBROUTINE Setup_Da END MODULE angular_distance !----------------------------------------------------------------------------- DOUBLE PRECISION FUNCTION da(z) ! proper distance in units of h^-1 Mpc. USE angular_distance ! Returns da(z) by interpolating the table generated by Setup_Da. IMPLICIT none DOUBLE PRECISION, intent(IN) :: z DOUBLE PRECISION :: a,b,hh,x INTEGER :: jlo x=-dlog(1.+z) ! x=ln(a)=-ln(1+z) CALL hunt(xa,nw,x,jlo) hh=xa(jlo+1)-xa(jlo) a=(xa(jlo+1)-x)/hh b=(x-xa(jlo))/hh da=(a*za(jlo)+b*za(jlo+1)+((a**3-a)*z2a(jlo)+(b**3-b)*z2a(jlo+1))*(hh**2)/6.) return END FUNCTION Da !----------------------------------------------------------------------------- DOUBLE PRECISION FUNCTION one_over_h(z) USE cosmo ! one_over_h(z) = H_0/H(z) (dimensionless) ! x=a IMPLICIT none DOUBLE PRECISION, intent(IN) :: z double precision :: x,ok0,or0=0d0 ! radiation density is ignored! ok0=1d0-om0-ode0 x=1d0/(1d0+z) one_over_h = 1d0/sqrt(om0/x**3d0+ok0/x**2d0+or0/x**4d0+ode0/x**(3d0+3d0*w)) return END FUNCTION one_over_h ypdfdks/dverk.f 000644 000767 000024 00000121335 12577165440 016370 0 ustar 00eiichirokomatsu staff 000000 000000 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine dverk (n, fcn, x, y, xend, tol, ind, c, nw, w) integer n, ind, nw, k double precision x, y(n), xend, tol, c(*), w(nw,9), temp c c*********************************************************************** c * c note added 11/14/85. * c * c if you discover any errors in this subroutine, please contact * c * c kenneth r. jackson * c department of computer science * c university of toronto * c toronto, ontario, * c canada m5s 1a4 * c * c phone: 416-978-7075 * c * c electronic mail: * c uucp: {cornell,decvax,ihnp4,linus,uw-beaver}!utcsri!krj * c csnet: krj@toronto * c arpa: krj.toronto@csnet-relay * c bitnet: krj%toronto@csnet-relay.arpa * c * c dverk is written in fortran 66. * c * c the constants dwarf and rreb -- c(10) and c(11), respectively -- are * c set for a vax in double precision. they should be reset, as * c described below, if this program is run on another machine. * c * c the c array is declared in this subroutine to have one element only, * c although more elements are referenced in this subroutine. this * c causes some compilers to issue warning messages. there is, though, * c no error provided c is declared sufficiently large in the calling * c program, as described below. * c * c the following external statement for fcn was added to avoid a * c warning message from the unix f77 compiler. the original dverk * c comments and code follow it. * c * c*********************************************************************** c external fcn c c*********************************************************************** c * c purpose - this is a runge-kutta subroutine based on verner's * c fifth and sixth order pair of formulas for finding approximations to * c the solution of a system of first order ordinary differential * c equations with initial conditions. it attempts to keep the global * c error proportional to a tolerance specified by the user. (the * c proportionality depends on the kind of error control that is used, * c as well as the differential equation and the range of integration.) * c * c various options are available to the user, including different * c kinds of error control, restrictions on step sizes, and interrupts * c which permit the user to examine the state of the calculation (and * c perhaps make modifications) during intermediate stages. * c * c the program is efficient for non-stiff systems. however, a good * c variable-order-adams method will probably be more efficient if the * c function evaluations are very costly. such a method would also be * c more suitable if one wanted to obtain a large number of intermediate * c solution values by interpolation, as might be the case for example * c with graphical output. * c * c hull-enright-jackson 1/10/76 * c * c*********************************************************************** c * c use - the user must specify each of the following * c * c n number of equations * c * c fcn name of subroutine for evaluating functions - the subroutine * c itself must also be provided by the user - it should be of * c the following form * c subroutine fcn(n, x, y, yprime) * c integer n * c double precision x, y(n), yprime(n) * c *** etc *** * c and it should evaluate yprime, given n, x and y * c * c x independent variable - initial value supplied by user * c * c y dependent variable - initial values of components y(1), y(2), * c ..., y(n) supplied by user * c * c xend value of x to which integration is to be carried out - it may * c be less than the initial value of x * c * c tol tolerance - the subroutine attempts to control a norm of the * c local error in such a way that the global error is * c proportional to tol. in some problems there will be enough * c damping of errors, as well as some cancellation, so that * c the global error will be less than tol. alternatively, the * c control can be viewed as attempting to provide a * c calculated value of y at xend which is the exact solution * c to the problem y' = f(x,y) + e(x) where the norm of e(x) * c is proportional to tol. (the norm is a max norm with * c weights that depend on the error control strategy chosen * c by the user. the default weight for the k-th component is * c 1/max(1,abs(y(k))), which therefore provides a mixture of * c absolute and relative error control.) * c * c ind indicator - on initial entry ind must be set equal to either * c 1 or 2. if the user does not wish to use any options, he * c should set ind to 1 - all that remains for the user to do * c then is to declare c and w, and to specify nw. the user * c may also select various options on initial entry by * c setting ind = 2 and initializing the first 9 components of * c c as described in the next section. he may also re-enter * c the subroutine with ind = 3 as mentioned again below. in * c any event, the subroutine returns with ind equal to * c 3 after a normal return * c 4, 5, or 6 after an interrupt (see options c(8), c(9)) * c -1, -2, or -3 after an error condition (see below) * c * c c communications vector - the dimension must be greater than or * c equal to 24, unless option c(1) = 4 or 5 is used, in which * c case the dimension must be greater than or equal to n+30 * c * c nw first dimension of workspace w - must be greater than or * c equal to n * c * c w workspace matrix - first dimension must be nw and second must * c be greater than or equal to 9 * c * c the subroutine will normally return with ind = 3, having * c replaced the initial values of x and y with, respectively, the value * c of xend and an approximation to y at xend. the subroutine can be * c called repeatedly with new values of xend without having to change * c any other argument. however, changes in tol, or any of the options * c described below, may also be made on such a re-entry if desired. * c * c three error returns are also possible, in which case x and y * c will be the most recently accepted values - * c with ind = -3 the subroutine was unable to satisfy the error * c requirement with a particular step-size that is less than or * c equal to hmin, which may mean that tol is too small * c with ind = -2 the value of hmin is greater than hmax, which * c probably means that the requested tol (which is used in the * c calculation of hmin) is too small * c with ind = -1 the allowed maximum number of fcn evaluations has * c been exceeded, but this can only occur if option c(7), as * c described in the next section, has been used * c * c there are several circumstances that will cause the calculations * c to be terminated, along with output of information that will help * c the user determine the cause of the trouble. these circumstances * c involve entry with illegal or inconsistent values of the arguments, * c such as attempting a normal re-entry without first changing the * c value of xend, or attempting to re-enter with ind less than zero. * c * c*********************************************************************** c * c options - if the subroutine is entered with ind = 1, the first 9 * c components of the communications vector are initialized to zero, and * c the subroutine uses only default values for each option. if the * c subroutine is entered with ind = 2, the user must specify each of * c these 9 components - normally he would first set them all to zero, * c and then make non-zero those that correspond to the particular * c options he wishes to select. in any event, options may be changed on * c re-entry to the subroutine - but if the user changes any of the * c options, or tol, in the course of a calculation he should be careful * c about how such changes affect the subroutine - it may be better to * c restart with ind = 1 or 2. (components 10 to 24 of c are used by the * c program - the information is available to the user, but should not * c normally be changed by him.) * c * c c(1) error control indicator - the norm of the local error is the * c max norm of the weighted error estimate vector, the * c weights being determined according to the value of c(1) - * c if c(1)=1 the weights are 1 (absolute error control) * c if c(1)=2 the weights are 1/abs(y(k)) (relative error * c control) * c if c(1)=3 the weights are 1/max(abs(c(2)),abs(y(k))) * c (relative error control, unless abs(y(k)) is less * c than the floor value, abs(c(2)) ) * c if c(1)=4 the weights are 1/max(abs(c(k+30)),abs(y(k))) * c (here individual floor values are used) * c if c(1)=5 the weights are 1/abs(c(k+30)) * c for all other values of c(1), including c(1) = 0, the * c default values of the weights are taken to be * c 1/max(1,abs(y(k))), as mentioned earlier * c (in the two cases c(1) = 4 or 5 the user must declare the * c dimension of c to be at least n+30 and must initialize the * c components c(31), c(32), ..., c(n+30).) * c * c c(2) floor value - used when the indicator c(1) has the value 3 * c * c c(3) hmin specification - if not zero, the subroutine chooses hmin * c to be abs(c(3)) - otherwise it uses the default value * c 10*max(dwarf,rreb*max(weighted norm y/tol,abs(x))), * c where dwarf is a very small positive machine number and * c rreb is the relative roundoff error bound * c * c c(4) hstart specification - if not zero, the subroutine will use * c an initial hmag equal to abs(c(4)), except of course for * c the restrictions imposed by hmin and hmax - otherwise it * c uses the default value of hmax*(tol)**(1/6) * c * c c(5) scale specification - this is intended to be a measure of the * c scale of the problem - larger values of scale tend to make * c the method more reliable, first by possibly restricting * c hmax (as described below) and second, by tightening the * c acceptance requirement - if c(5) is zero, a default value * c of 1 is used. for linear homogeneous problems with * c constant coefficients, an appropriate value for scale is a * c norm of the associated matrix. for other problems, an * c approximation to an average value of a norm of the * c jacobian along the trajectory may be appropriate * c * c c(6) hmax specification - four cases are possible * c if c(6).ne.0 and c(5).ne.0, hmax is taken to be * c min(abs(c(6)),2/abs(c(5))) * c if c(6).ne.0 and c(5).eq.0, hmax is taken to be abs(c(6)) * c if c(6).eq.0 and c(5).ne.0, hmax is taken to be * c 2/abs(c(5)) * c if c(6).eq.0 and c(5).eq.0, hmax is given a default value * c of 2 * c * c c(7) maximum number of function evaluations - if not zero, an * c error return with ind = -1 will be caused when the number * c of function evaluations exceeds abs(c(7)) * c * c c(8) interrupt number 1 - if not zero, the subroutine will * c interrupt the calculations after it has chosen its * c preliminary value of hmag, and just before choosing htrial * c and xtrial in preparation for taking a step (htrial may * c differ from hmag in sign, and may require adjustment if * c xend is near) - the subroutine returns with ind = 4, and * c will resume calculation at the point of interruption if * c re-entered with ind = 4 * c * c c(9) interrupt number 2 - if not zero, the subroutine will * c interrupt the calculations immediately after it has * c decided whether or not to accept the result of the most * c recent trial step, with ind = 5 if it plans to accept, or * c ind = 6 if it plans to reject - y(*) is the previously * c accepted result, while w(*,9) is the newly computed trial * c value, and w(*,2) is the unweighted error estimate vector. * c the subroutine will resume calculations at the point of * c interruption on re-entry with ind = 5 or 6. (the user may * c change ind in this case if he wishes, for example to force * c acceptance of a step that would otherwise be rejected, or * c vice versa. he can also restart with ind = 1 or 2.) * c * c*********************************************************************** c * c summary of the components of the communications vector * c * c prescribed at the option determined by the program * c of the user * c * c c(10) rreb(rel roundoff err bnd) * c c(1) error control indicator c(11) dwarf (very small mach no) * c c(2) floor value c(12) weighted norm y * c c(3) hmin specification c(13) hmin * c c(4) hstart specification c(14) hmag * c c(5) scale specification c(15) scale * c c(6) hmax specification c(16) hmax * c c(7) max no of fcn evals c(17) xtrial * c c(8) interrupt no 1 c(18) htrial * c c(9) interrupt no 2 c(19) est * c c(20) previous xend * c c(21) flag for xend * c c(22) no of successful steps * c c(23) no of successive failures * c c(24) no of fcn evals * c * c if c(1) = 4 or 5, c(31), c(32), ... c(n+30) are floor values * c * c*********************************************************************** c * c an overview of the program * c * c begin initialization, parameter checking, interrupt re-entries * c ......abort if ind out of range 1 to 6 * c . cases - initial entry, normal re-entry, interrupt re-entries * c . case 1 - initial entry (ind .eq. 1 or 2) * c v........abort if n.gt.nw or tol.le.0 * c . if initial entry without options (ind .eq. 1) * c . set c(1) to c(9) equal to zero * c . else initial entry with options (ind .eq. 2) * c . make c(1) to c(9) non-negative * c . make floor values non-negative if they are to be used * c . end if * c . initialize rreb, dwarf, prev xend, flag, counts * c . case 2 - normal re-entry (ind .eq. 3) * c .........abort if xend reached, and either x changed or xend not * c . re-initialize flag * c . case 3 - re-entry following an interrupt (ind .eq. 4 to 6) * c v transfer control to the appropriate re-entry point....... * c . end cases . * c . end initialization, etc. . * c . v * c . loop through the following 4 stages, once for each trial step . * c . stage 1 - prepare . * c***********error return (with ind=-1) if no of fcn evals too great . * c . calc slope (adding 1 to no of fcn evals) if ind .ne. 6 . * c . calc hmin, scale, hmax . * c***********error return (with ind=-2) if hmin .gt. hmax . * c . calc preliminary hmag . * c***********interrupt no 1 (with ind=4) if requested.......re-entry.v * c . calc hmag, xtrial and htrial . * c . end stage 1 . * c v stage 2 - calc ytrial (adding 7 to no of fcn evals) . * c . stage 3 - calc the error estimate . * c . stage 4 - make decisions . * c . set ind=5 if step acceptable, else set ind=6 . * c***********interrupt no 2 if requested....................re-entry.v * c . if step accepted (ind .eq. 5) * c . update x, y from xtrial, ytrial * c . add 1 to no of successful steps * c . set no of successive failures to zero * c**************return(with ind=3, xend saved, flag set) if x .eq. xend * c . else step not accepted (ind .eq. 6) * c . add 1 to no of successive failures * c**************error return (with ind=-3) if hmag .le. hmin * c . end if * c . end stage 4 * c . end loop * c . * c begin abort action * c output appropriate message about stopping the calculations, * c along with values of ind, n, nw, tol, hmin, hmax, x, xend, * c previous xend, no of successful steps, no of successive * c failures, no of fcn evals, and the components of y * c stop * c end abort action * c * c*********************************************************************** c c ****************************************************************** c * begin initialization, parameter checking, interrupt re-entries * c ****************************************************************** c c ......abort if ind out of range 1 to 6 if (ind.lt.1 .or. ind.gt.6) go to 500 c c cases - initial entry, normal re-entry, interrupt re-entries go to (5, 5, 45, 1111, 2222, 2222), ind c case 1 - initial entry (ind .eq. 1 or 2) c .........abort if n.gt.nw or tol.le.0 5 if (n.gt.nw .or. tol.le.0.d0) go to 500 if (ind.eq. 2) go to 15 c initial entry without options (ind .eq. 1) c set c(1) to c(9) equal to 0 do 10 k = 1, 9 c(k) = 0.d0 10 continue go to 35 15 continue c initial entry with options (ind .eq. 2) c make c(1) to c(9) non-negative do 20 k = 1, 9 c(k) = dabs(c(k)) 20 continue c make floor values non-negative if they are to be used if (c(1).ne.4.d0 .and. c(1).ne.5.d0) go to 30 do 25 k = 1, n c(k+30) = dabs(c(k+30)) 25 continue 30 continue 35 continue c initialize rreb, dwarf, prev xend, flag, counts c(10) = 2.d0**(-56) c(11) = 1.d-35 c set previous xend initially to initial value of x c(20) = x do 40 k = 21, 24 c(k) = 0.d0 40 continue go to 50 c case 2 - normal re-entry (ind .eq. 3) c .........abort if xend reached, and either x changed or xend not 45 if (c(21).ne.0.d0 .and. + (x.ne.c(20) .or. xend.eq.c(20))) go to 500 c re-initialize flag c(21) = 0.d0 go to 50 c case 3 - re-entry following an interrupt (ind .eq. 4 to 6) c transfer control to the appropriate re-entry point.......... c this has already been handled by the computed go to . c end cases v 50 continue c c end initialization, etc. c c ****************************************************************** c * loop through the following 4 stages, once for each trial step * c * until the occurrence of one of the following * c * (a) the normal return (with ind .eq. 3) on reaching xend in * c * stage 4 * c * (b) an error return (with ind .lt. 0) in stage 1 or stage 4 * c * (c) an interrupt return (with ind .eq. 4, 5 or 6), if * c * requested, in stage 1 or stage 4 * c ****************************************************************** c 99999 continue c c *************************************************************** c * stage 1 - prepare - do calculations of hmin, hmax, etc., * c * and some parameter checking, and end up with suitable * c * values of hmag, xtrial and htrial in preparation for taking * c * an integration step. * c *************************************************************** c c***********error return (with ind=-1) if no of fcn evals too great if (c(7).eq.0.d0 .or. c(24).lt.c(7)) go to 100 ind = -1 return 100 continue c c calculate slope (adding 1 to no of fcn evals) if ind .ne. 6 if (ind .eq. 6) go to 105 call fcn(n, x, y, w(1,1)) c(24) = c(24) + 1.d0 105 continue c c calculate hmin - use default unless value prescribed c(13) = c(3) if (c(3) .ne. 0.d0) go to 165 c calculate default value of hmin c first calculate weighted norm y - c(12) - as specified c by the error control indicator c(1) temp = 0.d0 if (c(1) .ne. 1.d0) go to 115 c absolute error control - weights are 1 do 110 k = 1, n temp = dmax1(temp, dabs(y(k))) 110 continue c(12) = temp go to 160 115 if (c(1) .ne. 2.d0) go to 120 c relative error control - weights are 1/dabs(y(k)) so c weighted norm y is 1 c(12) = 1.d0 go to 160 120 if (c(1) .ne. 3.d0) go to 130 c weights are 1/max(c(2),abs(y(k))) do 125 k = 1, n temp = dmax1(temp, dabs(y(k))/c(2)) 125 continue c(12) = dmin1(temp, 1.d0) go to 160 130 if (c(1) .ne. 4.d0) go to 140 c weights are 1/max(c(k+30),abs(y(k))) do 135 k = 1, n temp = dmax1(temp, dabs(y(k))/c(k+30)) 135 continue c(12) = dmin1(temp, 1.d0) go to 160 140 if (c(1) .ne. 5.d0) go to 150 c weights are 1/c(k+30) do 145 k = 1, n temp = dmax1(temp, dabs(y(k))/c(k+30)) 145 continue c(12) = temp go to 160 150 continue c default case - weights are 1/max(1,abs(y(k))) do 155 k = 1, n temp = dmax1(temp, dabs(y(k))) 155 continue c(12) = dmin1(temp, 1.d0) 160 continue c(13) = 10.d0*dmax1(c(11),c(10)*dmax1(c(12)/tol,dabs(x))) 165 continue c c calculate scale - use default unless value prescribed c(15) = c(5) if (c(5) .eq. 0.d0) c(15) = 1.d0 c c calculate hmax - consider 4 cases c case 1 both hmax and scale prescribed if (c(6).ne.0.d0 .and. c(5).ne.0.d0) + c(16) = dmin1(c(6), 2.d0/c(5)) c case 2 - hmax prescribed, but scale not if (c(6).ne.0.d0 .and. c(5).eq.0.d0) c(16) = c(6) c case 3 - hmax not prescribed, but scale is if (c(6).eq.0.d0 .and. c(5).ne.0.d0) c(16) = 2.d0/c(5) c case 4 - neither hmax nor scale is provided if (c(6).eq.0.d0 .and. c(5).eq.0.d0) c(16) = 2.d0 c c***********error return (with ind=-2) if hmin .gt. hmax if (c(13) .le. c(16)) go to 170 ind = -2 return 170 continue c c calculate preliminary hmag - consider 3 cases if (ind .gt. 2) go to 175 c case 1 - initial entry - use prescribed value of hstart, if c any, else default c(14) = c(4) if (c(4) .eq. 0.d0) c(14) = c(16)*tol**(1.0d0/6.0d0) go to 185 175 if (c(23) .gt. 1.d0) go to 180 c case 2 - after a successful step, or at most one failure, c use min(2, .9*(tol/est)**(1/6))*hmag, but avoid possible c overflow. then avoid reduction by more than half. temp = 2.d0*c(14) if (tol .lt. (2.d0/.9d0)**6*c(19)) + temp = .9d0*(tol/c(19))**(1.0d0/6.0d0)*c(14) c(14) = dmax1(temp, .5d0*c(14)) go to 185 180 continue c case 3 - after two or more successive failures c(14) = .5d0*c(14) 185 continue c c check against hmax c(14) = dmin1(c(14), c(16)) c c check against hmin c(14) = dmax1(c(14), c(13)) c c***********interrupt no 1 (with ind=4) if requested if (c(8) .eq. 0.d0) go to 1111 ind = 4 return c resume here on re-entry with ind .eq. 4 ........re-entry.. 1111 continue c c calculate hmag, xtrial - depending on preliminary hmag, xend if (c(14) .ge. dabs(xend - x)) go to 190 c do not step more than half way to xend c(14) = dmin1(c(14), .5d0*dabs(xend - x)) c(17) = x + dsign(c(14), xend - x) go to 195 190 continue c hit xend exactly c(14) = dabs(xend - x) c(17) = xend 195 continue c c calculate htrial c(18) = c(17) - x c c end stage 1 c c *************************************************************** c * stage 2 - calculate ytrial (adding 7 to no of fcn evals). * c * w(*,2), ... w(*,8) hold intermediate results needed in * c * stage 3. w(*,9) is temporary storage until finally it holds * c * ytrial. * c *************************************************************** c temp = c(18)/1398169080000.d0 c do 200 k = 1, n w(k,9) = y(k) + temp*w(k,1)*233028180000.d0 200 continue call fcn(n, x + c(18)/6.d0, w(1,9), w(1,2)) c do 205 k = 1, n w(k,9) = y(k) + temp*( w(k,1)*74569017600.d0 + + w(k,2)*298276070400.d0 ) 205 continue call fcn(n, x + c(18)*(4.d0/15.d0), w(1,9), w(1,3)) c do 210 k = 1, n w(k,9) = y(k) + temp*( w(k,1)*1165140900000.d0 + - w(k,2)*3728450880000.d0 + + w(k,3)*3495422700000.d0 ) 210 continue call fcn(n, x + c(18)*(2.d0/3.d0), w(1,9), w(1,4)) c do 215 k = 1, n w(k,9) = y(k) + temp*( - w(k,1)*3604654659375.d0 + + w(k,2)*12816549900000.d0 + - w(k,3)*9284716546875.d0 + + w(k,4)*1237962206250.d0 ) 215 continue call fcn(n, x + c(18)*(5.d0/6.d0), w(1,9), w(1,5)) c do 220 k = 1, n w(k,9) = y(k) + temp*( w(k,1)*3355605792000.d0 + - w(k,2)*11185352640000.d0 + + w(k,3)*9172628850000.d0 + - w(k,4)*427218330000.d0 + + w(k,5)*482505408000.d0 ) 220 continue call fcn(n, x + c(18), w(1,9), w(1,6)) c do 225 k = 1, n w(k,9) = y(k) + temp*( - w(k,1)*770204740536.d0 + + w(k,2)*2311639545600.d0 + - w(k,3)*1322092233000.d0 + - w(k,4)*453006781920.d0 + + w(k,5)*326875481856.d0 ) 225 continue call fcn(n, x + c(18)/15.d0, w(1,9), w(1,7)) c do 230 k = 1, n w(k,9) = y(k) + temp*( w(k,1)*2845924389000.d0 + - w(k,2)*9754668000000.d0 + + w(k,3)*7897110375000.d0 + - w(k,4)*192082660000.d0 + + w(k,5)*400298976000.d0 + + w(k,7)*201586000000.d0 ) 230 continue call fcn(n, x + c(18), w(1,9), w(1,8)) c c calculate ytrial, the extrapolated approximation and store c in w(*,9) do 235 k = 1, n w(k,9) = y(k) + temp*( w(k,1)*104862681000.d0 + + w(k,3)*545186250000.d0 + + w(k,4)*446637345000.d0 + + w(k,5)*188806464000.d0 + + w(k,7)*15076875000.d0 + + w(k,8)*97599465000.d0 ) 235 continue c c add 7 to the no of fcn evals c(24) = c(24) + 7.d0 c c end stage 2 c c *************************************************************** c * stage 3 - calculate the error estimate est. first calculate * c * the unweighted absolute error estimate vector (per unit * c * step) for the unextrapolated approximation and store it in * c * w(*,2). then calculate the weighted max norm of w(*,2) as * c * specified by the error control indicator c(1). finally, * c * modify this result to produce est, the error estimate (per * c * unit step) for the extrapolated approximation ytrial. * c *************************************************************** c c calculate the unweighted absolute error estimate vector do 300 k = 1, n w(k,2) = ( w(k,1)*8738556750.d0 + + w(k,3)*9735468750.d0 + - w(k,4)*9709507500.d0 + + w(k,5)*8582112000.d0 + + w(k,6)*95329710000.d0 + - w(k,7)*15076875000.d0 + - w(k,8)*97599465000.d0)/1398169080000.d0 300 continue c c calculate the weighted max norm of w(*,2) as specified by c the error control indicator c(1) temp = 0.d0 if (c(1) .ne. 1.d0) go to 310 c absolute error control do 305 k = 1, n temp = dmax1(temp,dabs(w(k,2))) 305 continue go to 360 310 if (c(1) .ne. 2.d0) go to 320 c relative error control do 315 k = 1, n temp = dmax1(temp, dabs(w(k,2)/y(k))) 315 continue go to 360 320 if (c(1) .ne. 3.d0) go to 330 c weights are 1/max(c(2),abs(y(k))) do 325 k = 1, n temp = dmax1(temp, dabs(w(k,2)) + / dmax1(c(2), dabs(y(k))) ) 325 continue go to 360 330 if (c(1) .ne. 4.d0) go to 340 c weights are 1/max(c(k+30),abs(y(k))) do 335 k = 1, n temp = dmax1(temp, dabs(w(k,2)) + / dmax1(c(k+30), dabs(y(k))) ) 335 continue go to 360 340 if (c(1) .ne. 5.d0) go to 350 c weights are 1/c(k+30) do 345 k = 1, n temp = dmax1(temp, dabs(w(k,2)/c(k+30))) 345 continue go to 360 350 continue c default case - weights are 1/max(1,abs(y(k))) do 355 k = 1, n temp = dmax1(temp, dabs(w(k,2)) + / dmax1(1.d0, dabs(y(k))) ) 355 continue 360 continue c c calculate est - (the weighted max norm of w(*,2))*hmag*scale c - est is intended to be a measure of the error per unit c step in ytrial c(19) = temp*c(14)*c(15) c c end stage 3 c c *************************************************************** c * stage 4 - make decisions. * c *************************************************************** c c set ind=5 if step acceptable, else set ind=6 ind = 5 if (c(19) .gt. tol) ind = 6 c c***********interrupt no 2 if requested if (c(9) .eq. 0.d0) go to 2222 return c resume here on re-entry with ind .eq. 5 or 6 ...re-entry.. 2222 continue c if (ind .eq. 6) go to 410 c step accepted (ind .eq. 5), so update x, y from xtrial, c ytrial, add 1 to the no of successful steps, and set c the no of successive failures to zero x = c(17) do 400 k = 1, n y(k) = w(k,9) 400 continue c(22) = c(22) + 1.d0 c(23) = 0.d0 c**************return(with ind=3, xend saved, flag set) if x .eq. xend if (x .ne. xend) go to 405 ind = 3 c(20) = xend c(21) = 1.d0 return 405 continue go to 420 410 continue c step not accepted (ind .eq. 6), so add 1 to the no of c successive failures c(23) = c(23) + 1.d0 c**************error return (with ind=-3) if hmag .le. hmin if (c(14) .gt. c(13)) go to 415 ind = -3 return 415 continue 420 continue c c end stage 4 c go to 99999 c end loop c c begin abort action 500 continue c write(6,505) ind, tol, x, n, c(13), xend, nw, c(16), c(20), + c(22), c(23), c(24), (y(k), k = 1, n) 505 format( /// 1h0, 58hcomputation stopped in dverk with the followin +g values - + / 1h0, 5hind =, i4, 5x, 6htol =, 1pd13.6, 5x, 11hx =, + 1pd22.15 + / 1h , 5hn =, i4, 5x, 6hhmin =, 1pd13.6, 5x, 11hxend =, + 1pd22.15 + / 1h , 5hnw =, i4, 5x, 6hhmax =, 1pd13.6, 5x, 11hprev xend =, + 1pd22.15 + / 1h0, 14x, 27hno of successful steps =, 0pf8.0 + / 1h , 14x, 27hno of successive failures =, 0pf8.0 + / 1h , 14x, 27hno of function evals =, 0pf8.0 + / 1h0, 23hthe components of y are + // (1h , 1p5d24.15) ) c stop c c end abort action c end ypdfdks/fx.f90 000644 000767 000024 00000002702 12720033162 016020 0 ustar 00eiichirokomatsu staff 000000 000000 !----------------------------------------------------------------------------- ! ! Read in, tabulate, and interpolate a function, and find an inverse; i.e., ! given a function of x, f=f(x), we find a value of x for a given value of f. ! ! [1st column] x ! [2nd column] f(x) ! ! May 18, 2016 ! E. Komatsu ! !----------------------------------------------------------------------------- MODULE Fx IMPLICIT none INTEGER :: jlo,ndata DOUBLE PRECISION, dimension(:), allocatable :: xa,ya,y2a contains SUBROUTINE Open_File(filename,n) IMPLICIT none DOUBLE PRECISION, dimension(:), allocatable :: x,f CHARACTER(len=128) :: filename INTEGER, intent(IN) :: n INTEGER :: i ndata=n ALLOCATE(xa(ndata),ya(ndata),y2a(ndata)) ALLOCATE(x(ndata),f(ndata)) open(1,file=filename,status='old',form='formatted') print*,'read in '//trim(filename) do i=1,ndata read(1,*)x(i),f(i) enddo close(1) xa=f ya=x DEALLOCATE(x,f) CALL spline(xa,ya,ndata,1.d30,1.d30,y2a) return END SUBROUTINE Open_File SUBROUTINE Close_File DEALLOCATE(xa,ya,y2a) return END SUBROUTINE Close_File END MODULE Fx DOUBLE PRECISION FUNCTION Find_x(f) Use Fx IMPLICIT none DOUBLE PRECISION :: a,b,h,x,y,f x = f CALL hunt(xa,ndata,x,jlo) h=xa(jlo+1)-xa(jlo) a=(xa(jlo+1)-x)/h b=(x-xa(jlo))/h y=a*ya(jlo)+b*ya(jlo+1)+((a**3-a)*y2a(jlo)+(b**3-b)*y2a(jlo+1))*(h**2)/6. Find_x = y return END FUNCTION Find_x ypdfdks/generate_table.f90 000644 000767 000024 00000002253 12726541544 020363 0 ustar 00eiichirokomatsu staff 000000 000000 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 ypdfdks/growth.f90 000644 000767 000024 00000007121 12577165440 016734 0 ustar 00eiichirokomatsu staff 000000 000000 !----------------------------------------------------------------------------- ! ! Linear Growth Factor, g(z)=D(z)(1+z), which is constant and equal to ! unity during the matter dominated era. Also computes its derivative ! with respect to lna, dg/dlna. ! ! *** Radiation density is ignored in the calculation *** ! ! << sample code >> ! ! USE cosmo ! USE growth ! double precision :: g,dgdlna,redshift ! external g,dgdlna ! ! Specify three cosmological parameters in double precision. ! ! The data type has been defined in MODULE cosmo. ! ode0=0.723d0 ! om0=0.277d0 ! w = -1d0 ! CALL setup_growth ! redshift=10d0 ! print*,'omega matter=',om0 ! print*,'omega de=',ode0 ! print*,'w=',w ! print*,'g(z) at z=',redshift,'is',g(redshift) ! print*,'dg/dlna(z) at z=',redshift,'is',dgdlna(redshift) ! end ! ! August 23, 2008 ! E. Komatsu ! !----------------------------------------------------------------------------- MODULE Growth INTEGER, parameter :: nw=501 DOUBLE PRECISION, dimension(nw) :: xa,ga,g2a,dga,dg2a contains SUBROUTINE Setup_Growth ! tabulate g(x), where x=ln(a) IMPLICIT none INTEGER :: i,ind DOUBLE PRECISION :: x,xini,y(2),c(24),work(2,9),tol=1.d-8 EXTERNAL fderiv ind=1 do i=1,nw x = (i-nw)*0.01 ! x=ln(a)=-5.00,-4.99,...,0 xa(i) = x if(i.eq.1)then xini = xa(i) y(1)=1d0 ! initial condition, g=1, at z=-5.00 (z=147.413) y(2)=0d0 ! initial condition, g'=1, at z=-5.00 (z=147.413) else call dverk(2,fderiv,xini,y,x,tol,ind,c,2,work) xini = x endif ga(i) = y(1) dga(i)= y(2) enddo CALL spline(xa,ga,nw,1.d30,1.d30,g2a) ! evaluate g2a(i)=d^2[g(i)]/dy^2 CALL spline(xa,dga,nw,1.d30,1.d30,dg2a)! evaluate g2a(i)=d^2[dg/dlna(i)]/dy^2 return END SUBROUTINE Setup_Growth END MODULE Growth !----------------------------------------------------------------------------- SUBROUTINE fderiv(n,x,y,yprime) USE cosmo ! x = ln(a) = -ln(1+z) ! y(1) = g ! y(2) = g' ! See Eq.(76) of Komatsu et al. (2008) [WMAP 5-year] IMPLICIT none INTEGER, intent(IN) :: n DOUBLE PRECISION :: x,y(n),yprime(n) DOUBLE PRECISION :: ok,ode,ok0,h2,z,a DOUBLE PRECISION, parameter :: or0 = 0d0 ! ignore radiation density z=dexp(-x)-1d0 a=dexp(x) ok0=1d0-om0-ode0 h2=om0/a**3d0+or0/a**4d0+ok0/a**2d0+ode0/a**(3d0*(1.+w)) ok=ok0/h2*dexp(-2d0*x) ode=ode0/h2*dexp(-3d0*x*(1d0+w)) yprime(1)= y(2) yprime(2)= -(2.5d0+0.5d0*(ok-3d0*w*ode))*y(2) & -(2d0*ok+1.5d0*(1d0-w)*ode)*y(1) return END SUBROUTINE fderiv !----------------------------------------------------------------------------- DOUBLE PRECISION FUNCTION g(z) USE growth ! Returns g(z) = D(z)(1+z) by interpolating the table generated by ! Setup_Growth. IMPLICIT none DOUBLE PRECISION, intent(IN) :: z DOUBLE PRECISION :: a,b,hh,x INTEGER :: jlo x=-dlog(1d0+z) ! x=ln(a)=-ln(1+z) CALL hunt(xa,nw,x,jlo) hh=xa(jlo+1)-xa(jlo) a=(xa(jlo+1)-x)/hh b=(x-xa(jlo))/hh g=a*ga(jlo)+b*ga(jlo+1)+((a**3-a)*g2a(jlo)+(b**3-b)*g2a(jlo+1))*(hh**2)/6. return END FUNCTION g !----------------------------------------------------------------------------- DOUBLE PRECISION FUNCTION dgdlna(z) USE growth ! Returns dg/dlna(z) by interpolating the table generated by ! Setup_Growth. IMPLICIT none DOUBLE PRECISION, intent(IN) :: z DOUBLE PRECISION :: a,b,hh,x INTEGER :: jlo x=-dlog(1d0+z) ! x=ln(a)=-ln(1+z) CALL hunt(xa,nw,x,jlo) hh=xa(jlo+1)-xa(jlo) a=(xa(jlo+1)-x)/hh b=(x-xa(jlo))/hh dgdlna=a*dga(jlo)+b*dga(jlo+1)+((a**3-a)*dg2a(jlo)+(b**3-b)*dg2a(jlo+1))*(hh**2)/6. return END FUNCTION dgdlna ypdfdks/hunt.f 000644 000767 000024 00000001610 12577165440 016224 0 ustar 00eiichirokomatsu staff 000000 000000 SUBROUTINE hunt(xx,n,x,jlo) INTEGER jlo,n REAL*8 x,xx(n) INTEGER inc,jhi,jm LOGICAL ascnd ascnd=xx(n).gt.xx(1) if(jlo.le.0.or.jlo.gt.n)then jlo=0 jhi=n+1 goto 3 endif inc=1 if(x.ge.xx(jlo).eqv.ascnd)then 1 jhi=jlo+inc if(jhi.gt.n)then jhi=n+1 else if(x.ge.xx(jhi).eqv.ascnd)then jlo=jhi inc=inc+inc goto 1 endif else jhi=jlo 2 jlo=jhi-inc if(jlo.lt.1)then jlo=0 else if(x.lt.xx(jlo).eqv.ascnd)then jhi=jlo inc=inc+inc goto 2 endif endif 3 if(jhi-jlo.eq.1)return jm=(jhi+jlo)/2 if(x.gt.xx(jm).eqv.ascnd)then jlo=jm else jhi=jm endif goto 3 END C (C) Copr. 1986-92 Numerical Recipes Software D041&0(9p#3. ypdfdks/integrand.f90 000644 000767 000024 00000010356 12726613276 017402 0 ustar 00eiichirokomatsu staff 000000 000000 double precision function integrand(lnM,z,y) USE cosmo USE sigma USE growth USE angular_distance USE nfw USE fx ! A bin-integrated PDF of Compton Y is given by !
(Y) = int dz dV/dz int dlnM dn/dlnM*pi*[theta^2(Y)-theta^2(Y+dY)]
! This routine returns the first term of integrand, dV/dz*dn/dlnM*pi*theta^2(Y)
! The pressure profile is taken from the Planck 2013 generalized NFW (gNFW)
! pressure profile fit, with the hydrostatic mass bias given by "HSEbias" (default 1.2)
! May 21, 2016: E.Komatsu
IMPLICIT none
double precision, intent(IN) :: lnM,z,y
double precision :: HSEbias=1.2d0
double precision :: xin=1d-5,xout=6d0 ! pressure profile is cut at r=xout*r500
double precision :: omega,tol=1d-6,pi=3.1415926535d0
double precision :: zin=1d-5 ! redshift of the input matter power spectrum
double precision :: mvir,rvir,rs,delc,m500,r500,theta500
double precision :: rhoc,Ez,da,dvdz,norm
double precision :: deltac=1.6865d0,mf,dndlnRh,dndlnMh,lnnu,dlnnudlnRh,Rh
double precision :: g,scalefactor,m180d,m200d
double precision :: find_x,theta
real :: chebev,lnsigma2,lnRh,dlnsigma2dlnRh
integer :: i
external find_x
mvir=dexp(lnM)
! compute omega's, rhoc, and delc
omega=om0*(1d0+z)**3d0/Ez(z)**2d0 ! Omega(z); E(z)=H(z)/H0
rhoc=2.775d11*Ez(z)**2d0 ! critical density in units of h^2 M_sun/Mpc^3
! Eq.(6) of Bryan & Norman (1998), valid only for a flat universe!
delc=18d0*pi*pi+82d0*(omega-1d0)-39d0*(omega-1d0)**2d0
! One may also use Eq.(C19) of Nakamura & Suto (1997):
!!$ delc=omega*18d0*pi*pi*(1d0+0.4093d0*(1d0/omega-1d0)**0.9052d0)
! compute virial and NFW parameters
rvir=(3d0*mvir/4d0/pi/delc/rhoc)**(1d0/3d0) ! virial radius, h^-1 Mpc
! Concentration parameter from Duffy et al. (2008)
cvir=7.85d0*(mvir/2d12)**(-0.081)/(1d0+z)**0.71
! One may also use the concentration parameter from Seljak (2000)
!!$ cvir=10d0*(mvir/3.42d12)**(-0.2)/(1d0+z)
rs=rvir/cvir ! NFW scale radius, h^-1 Mpc
! convert to m500
CALL mvir2mdel(mvir,rs,cvir,500d0*rhoc,m500)
m500=m500/HSEbias
r500=(3d0*m500/4d0/pi/500d0/rhoc)**(1d0/3d0) ! h^-1 Mpc
theta500=r500/da(z)
! to find an angular radius at which the Compton Y value is equal to a given valule
norm=1.65d0*(h0/0.7d0)**2d0*Ez(z)**(8d0/3d0) &
*(m500/3d14/0.7d0)**(2d0/3d0+0.12d0)/0.5176d0*(r500/h0) &
*283d0/50d0/2.725d6*(0.7d0/h0)**1.5d0
if(Y/norm<6.d0)then ! if Y/norm>6, this halo does not contribute because it is not bright enough
theta=find_x(Y/norm)*theta500
else
theta=0d0
endif
! calculate mass function
! Sheth&Tormen, Tinker et al., and Bocquet et al.'s mass functions are given for the
! overdensity mass M200d (with respect to the mean mass density rather than
! the critical density).
CALL mvir2mdel(mvir,rs,cvir,200d0*omega*rhoc,m200d)
Rh=(3d0*m200d/4d0/pi/om0/2.775d11)**(1d0/3d0) ! h^-1 Mpc
! Alternatively, one may wish to use Jenkins et al.'s mass function,
! which is given for the overdensity mass M180d (with respect to the mean
! mass density rather than the critical density).
!!$ CALL mvir2mdel(mvir,rs,cvir,180d0*omega*rhoc,m180d)
!!$ Rh=(3d0*m180d/4d0/pi/om0/2.775d11)**(1d0/3d0) ! h^-1 Mpc
lnRh=real(dlog(Rh))
lnsigma2=CHEBEV(lnR1,lnR2,c,ndim,lnRh) ! ln(sigma^2)
dlnsigma2dlnRh=CHEBEV(lnR1,lnR2,cder,ndim,lnRh) ! dln(sigma^2)/dlnRh
scalefactor=g(z)/g(zin)*(1d0+zin)/(1d0+z)
lnnu=2d0*dlog(deltac)-dble(lnsigma2)-2d0*dlog(scalefactor) ! ln(nu)
dlnnudlnRh=-dble(dlnsigma2dlnRh) ! dln(nu)/dlnRh
dndlnRh=(3d0/4d0/pi)*dlnnudlnRh*mf(lnnu,z)/Rh**3d0
!!$ dndlnRh=(3d0/4d0/pi)*dlnnudlnRh*mf(lnnu)/Rh**3d0 ! if using ST or Jenkins
dndlnMh=dndlnRh/3d0 ! in units of h^3 Mpc^-3
! volume element
dvdz=(1d0+z)**2d0*da(z)**2d0*2998d0/Ez(z) ! h^-3 Mpc^3
! This is the integrand: dV/dz*dn/dlnM*pi*theta^2, to be integrated over lnM and z
integrand=dvdz*dndlnMh*pi*theta**2d0
return
end function integrand
!-----------------------------------------------------------------------------
DOUBLE PRECISION FUNCTION Ez(z)
USE cosmo
! Ez = H(z)/H0 (dimensionless)
! x=a
IMPLICIT none
DOUBLE PRECISION, intent(IN) :: z
double precision :: x,ok0,or0=0d0 ! radiation density is ignored!
ok0=1d0-om0-ode0
x=1d0/(1d0+z)
Ez=dsqrt(om0/x**3d0+ok0/x**2d0+or0/x**4d0+ode0/x**(3d0+3d0*w))
return
END FUNCTION Ez
ypdfdks/linearpk.f90 000644 000767 000024 00000004533 12577165440 017233 0 ustar 00eiichirokomatsu staff 000000 000000 !-----------------------------------------------------------------------------
!
! Read in, tabulate, and interpolate the linear matter power spectrum.
! The input file is in the following format (same as CAMB format):
!
! [1st column] k/h (i.e., wavenumber in units of h Mpc^-1)
! [2nd column] h^3 P(k) (i.e., P(k) in units of h^-3 Mpc^3)
!
! << sample code >>
!
! USE linearpk
! double precision :: linear_pk, k_ov_h
! character(len=128) :: filename
! integer :: n
! external linear_pk
! filename='wmap5baosn_max_likelihood_matterpower.dat'
! n=896 ! # of lines in the file
! CALL open_linearpk(filename,n)
! k_ov_h=1d0
! print*,'P(k) at k=',k_ov_h,' h Mpc^-1 is',linear_pk(k_ov_h),' h^3 Mpc^-3'
! end
!
! August 23, 2008
! E. Komatsu
!
!-----------------------------------------------------------------------------
MODULE LinearPk
IMPLICIT none
INTEGER :: jlo,ndata
DOUBLE PRECISION, dimension(:), allocatable :: xa,ya,y2a
contains
SUBROUTINE Open_LinearPk(filename,n)
IMPLICIT none
DOUBLE PRECISION, dimension(:), allocatable :: aklin,pk_lin
CHARACTER(len=128) :: filename
INTEGER, intent(IN) :: n
INTEGER :: i
ndata=n
ALLOCATE(xa(ndata),ya(ndata),y2a(ndata))
ALLOCATE(aklin(ndata),pk_lin(ndata))
open(1,file=filename,status='old',form='formatted')
print*,'read in '//trim(filename)
do i=1,ndata
read(1,*)aklin(i),pk_lin(i) ! k/h & h^3 P(k) in CAMB format!
enddo
close(1)
xa=dlog(aklin)
ya=dlog(pk_lin)
DEALLOCATE(aklin,pk_lin)
CALL spline(xa,ya,ndata,1.d30,1.d30,y2a)
return
END SUBROUTINE Open_LinearPk
SUBROUTINE Close_LinearPk
DEALLOCATE(xa,ya,y2a)
return
END SUBROUTINE Close_LinearPk
END MODULE LinearPk
DOUBLE PRECISION FUNCTION Linear_Pk(ak)
Use LinearPk
IMPLICIT none
DOUBLE PRECISION :: a,b,h,x,y,ak
x = dlog(ak)
CALL hunt(xa,ndata,x,jlo)
h=xa(jlo+1)-xa(jlo)
a=(xa(jlo+1)-x)/h
b=(x-xa(jlo))/h
y=a*ya(jlo)+b*ya(jlo+1)+((a**3-a)*y2a(jlo)+(b**3-b)*y2a(jlo+1))*(h**2)/6.
Linear_Pk = dexp(y)
return
END FUNCTION Linear_Pk
DOUBLE PRECISION FUNCTION dlnPkdlnk(ak)
Use LinearPk
IMPLICIT none
DOUBLE PRECISION :: a,b,h,x,y,ak
x = dlog(ak)
CALL hunt(xa,ndata,x,jlo)
h=xa(jlo+1)-xa(jlo)
a=(xa(jlo+1)-x)/h
b=(x-xa(jlo))/h
y=(ya(jlo+1)-ya(jlo))/h+(-(3.*a**2-1.)*y2a(jlo)+(3.*b**2-1.)*y2a(jlo+1))*h/6.
dlnPkdlnk = y
return
END FUNCTION DlnPkdlnk
ypdfdks/Makefile 000644 000767 000024 00000002331 12726614266 016541 0 ustar 00eiichirokomatsu staff 000000 000000 FC = g95
FLAGS = -O3
### Use Bocquet et al. (2015) "Hydro" mass function
OBJS = fx.o cosmo.o nfw.o da.o linearpk.o sigma.o spline.o hunt.o chebyshev.o growth.o dverk.o integrand.o qromb.o qgaus2.o rombint.o mvir2mdel.o zbrent.o mf_magneticum.o
### Use Bocquet et al. (2015) "DM Only" mass function
#OBJS = fx.o cosmo.o nfw.o da.o linearpk.o sigma.o spline.o hunt.o chebyshev.o growth.o dverk.o integrand.o qromb.o qgaus2.o rombint.o mvir2mdel.o zbrent.o mf_magneticum_dm.o
### Use Tinker et al. (2008) mass function
#OBJS = fx.o cosmo.o nfw.o da.o linearpk.o sigma.o spline.o hunt.o chebyshev.o growth.o dverk.o integrand.o qromb.o qgaus2.o rombint.o mvir2mdel.o zbrent.o mf_tinker_redshift.o
### Use Tinker et al. (2010) mass function
#OBJS = fx.o cosmo.o nfw.o da.o linearpk.o sigma.o spline.o hunt.o chebyshev.o growth.o dverk.o integrand.o qromb.o qgaus2.o rombint.o mvir2mdel.o zbrent.o mf_tinker10.o
.SUFFIXES: .f90
.f90.o:
$(FC) $(FFLAGS) -c $<
default: ypdf generate_table
ypdf: $(OBJS) ypdf.o
$(FC) $(FFLAGS) -o $@ $(OBJS) $@.o $(LDFLAGS)
generate_table: generate_table.o rombint1.o
$(FC) $(FFLAGS) -o $@ rombint1.o $@.o $(LDFLAGS)
clean:
-rm -f *.o *.mod *~
tidy: clean
-rm -f ypdf generate_table
-rm -f y_pdf.txt x_fx.txt
ypdfdks/mf_jenkins.f90 000644 000767 000024 00000007260 12577260252 017546 0 ustar 00eiichirokomatsu staff 000000 000000 !-----------------------------------------------------------------------------
!
! Halo multiplicity function in double precision, mf(lnnu),
! where lnnu=ln(nu) and nu is called the "threshold,"
! given by nu = [delta_c/sigma(R,z)]^2. Here, delta_c=1.6865, sigma(R,z)
! is the r.m.s. mass fluctuation within a top-hat smoothing of scale R
! at a redshift z.
!
! Ref. Jenkins et al., MNRAS, 321, 372 (2001)
! Note that "mf(lnnu)" here corresponds to f(sigma)/2 in this reference,
! where sigma=sqrt(1.6865/nu).
!
! IMPORTANT: Jenkins et al.'s mass function is NOT normalized, i.e.,
!
! int_{-infinity}^{+infinity} dln(nu) mf(lnnu) = 1
!
! is NOT satisfied! Also, we use the fitting function Eq.(B3) for Delta=180.
! (That is, the mass is defined within a spherical region whose mean
! density is equal to 200 times the MEAN MASS density of the universe.)
!
! The halo mass function can be computed in the following way.
!
! (1) First, the comoving number density of halos per ln(nu) is related
! to the mass function, dn/dM, as
!
! dM M dn/dM = rho_m dln(nu) mf(lnnu)
!
! where rho_m = (average mass density of the universe today)
!
! Dividing both sides by dM, one finds
!
! dn/dlnM = rho_m dln(nu)/dlnM mf(lnnu)/M
!
! (2) nu is more directly related to R, as nu=[delta_c/sigma(R,z)]^2.
! So, it is more convenient to write dn/dM as
!
! dn/dlnM = dlnR/dlnM dn/dlnR
!
! and compute dn/dlnR first, via
!
! dn/dlnR = rho_m dln(nu)/dlnR mf(lnnu)/M(R)
!
! Now, we can use M(R)=(4pi/3)rho_m R^3 to obtain
!
! dn/dlnR = (3/4pi) dln(nu)/dlnR mf(lnnu)/R^3
!
! with ln(nu) and dln(nu)/dlnR related to lnR via
! - ln(nu) = 2*ln(1.6865)-ln[sigma^2(R,z)]
! - dln(nu)/dlnR = -dln[sigma^2(R)]/dlnR (which is indep. of z)
!
! (3) Once dn/dlnR is obtained as a function of R, one can
! compute dn/dlnM using dlnR/dlnM=1/3, and obtain
!
! dn/dlnM = (1/3) dn/dlnR
!
! with M(R)=(4pi/3)rho_m R^3,
! and rho_m=2.775d11 (omega_matter h^2) M_solar Mpc^-3.
!
!
! << sample code >>
!
! USE linearpk
! USE sigma
! double precision :: deltac=1.6865d0,mf,dndlnRh,dndlnMh
! double precision :: lnnu,dlnnudlnRh,Mh
! REAL :: chebev
! REAL :: lnsigma2,lnRh,Rh,dlnsigma2dlnRh
! character(len=128) :: filename
! integer :: n
! filename='wmap5baosn_max_likelihood_matterpower.dat'
! n=896 ! # of lines in the file
! CALL get_linearpk(filename,n)
! CALL compute_sigma2
! Rh=8. ! h^-1 Mpc
! lnRh = alog(Rh)
! if(lnRh>lnR2.or.lnRh