c ******************************************************* c c c c (Standard) Perturbation Theory for c c Non-linver Evolution of Power Spectrum c c using Monte Carlo integration technique c c c c ---------- 2-loop contributions ------------ c c c c Time-stamp: <2008-07-03 00:07:37 ataruya> c c ******************************************************* c c c Note--. c c All the parameter 'ikmax' in this program must be c the same. c program stdPT_2loop c implicit none c integer ik, ikmax, ik_max, itype, ibox, ijob, ispt parameter(ikmax=5000) real*8 ak(ikmax), pk(ikmax) real*8 pk13(ikmax), pk22(ikmax) cc real*8 pk33a(ikmax), pk33b(ikmax), pk24(ikmax), pk15(ikmax) real*8 pi, sigma8, Lbox character infile*50 common /pk_data/ ak, pk, ik_max pi = 4.d0 * atan(1.d0) c ------------------------------------------------- c write(6,*) write(6,*) '**** Perturbation calculation of power spectrum ****' write(6,*) write(6,*) ' --- 2-loop contributions --- ' write(6,*) write(6,*) write(6,*) ' Choose PT correction: 1-loop [1] or 2-loop[2]' read(5,*) ispt if((ispt.ne.1) .and. (ispt.ne.2)) stop write(6,*) 'including finite-volume effect? :type (0 or 1, L_box)' write(6,*) ' (e.g., (0, 300) for boxsize L_box=300Mpc/h ) ' write(6,*) ' ( (1, ???) for ignoring finite-volume effect) ' read(5,*) ibox, Lbox if((ibox.ne.0) .and. (ibox.ne.1)) stop if(ispt.eq.2) then write(6,*) 'restart suspended run[0], start new calcuation[1]' read(5,*) ijob if((ijob.ne.0) .and. (ijob.ne.1)) stop endif c c /////// Load transfer function /////// c call load_transfer_function(sigma8) c write(6,*) ' load transfer function, done ' c c /////// Sigma8 normalization /////// c call normalization_trapez(sigma8) c write(6,*) ' normalization done ' c c //////// Truncation of low-k, high-k in linear P(k) //////// c call truncation_k(ibox, Lbox) c c /////// 1-loop correction of P(k) /////// c if(ispt.eq.1) call calc_one_loop_pk c c /////// 2-loop correction of P(k) /////// c if(ispt.eq.2) call calc_two_loop_pk(ijob) c end c c ******************************************************* c c subroutine load_transfer_function(sigma8) c c ******************************************************* c c c input file is assumed to be the transfer function data c created by CAMB code. c implicit none c integer ik, ikk, ikmax, ik_max parameter(ikmax=5000) real*8 Omega_b, Omega_m, Omega_v real*8 n_s, sigma8 real*8 kk, ak(ikmax), pk(ikmax) real*8 tf_m, pi common /pk_data/ ak, pk, ik_max character infile*50 pi = 4.d0 * datan(1.d0) c -------------------------------------------------------- c infile='Plin.z=0.dat' ccc infile='nishimichi_Tk.dat' c cc Omega_m = 0.234d0 cc Omega_v = 1.d0 - Omega_m cc sigma8 = 0.76d0 cc n_s = 0.961d0 c sigma8 = 0.817d0 c c open(9,file=infile,status='unknown') c do ik=1, ikmax read(9,*,END=10) kk, tf_m ak(ik) = kk pk(ik) = tf_m enddo c 10 continue ik_max = ik - 1 c write(6,*) 'ak(1)=', ak(1) write(6,*) 'ak(ik_max)=', ak(ik_max) write(6,*) 'ik_max=', ik_max pause c close(9) c end c c ******************************************************* c c subroutine normalization_trapez(sigma8) c c ******************************************************* c c implicit none integer ik, ik_max, ikmax parameter(ikmax=5000) real*8 sigma8 real*8 ak(ikmax), pk(ikmax) real*8 r_th, const, x real*8 W_TH, sigma_a, sigma_b, pi common /pk_data/ ak, pk, ik_max pi = 4.d0 * atan(1.d0) r_th = 8.d0 c --------------------------------------------------- c x = ak(1) * r_th if(x.lt.1.d-3) then W_TH = 1.d0 - x*x / 10.d0 + x**4 / 280.d0 else W_TH = 3.d0 * (sin(x) - x * cos(x))/x/x/x endif sigma_a = W_TH * W_TH * pk(1) * ak(1) * ak(1) sigma_a = sigma_a / (2.d0 * pi * pi) c const = 0.d0 do ik=2, ik_max x = ak(ik) * r_th if(x.lt.1.d-3) then W_TH = 1.d0 - x*x / 10.d0 + x**4 / 280.d0 else W_TH = 3.d0 * (sin(x) - x * cos(x))/x/x/x endif sigma_b = W_TH * W_TH * pk(ik) * ak(ik) * ak(ik) sigma_b = sigma_b / (2.d0 * pi * pi) const = const + & (sigma_a + sigma_b) * ( ak(ik) - ak(ik-1) )/ 2.d0 sigma_a = sigma_b enddo c do ik=1, ik_max pk(ik) = sigma8 * sigma8 / const * pk(ik) enddo c end c c ******************************************************* c c subroutine truncation_k(ibox, Lbox) c c ******************************************************* c c implicit none c integer ik, ikk, ikmax, ik_max, ibox parameter(ikmax=5000) real*8 ak(ikmax), pk(ikmax) real*8 akk(ikmax), pkk(ikmax) real*8 kmin, kmax, Lbox, pi common /pk_data/ ak, pk, ik_max pi = 4.d0 * datan(1.d0) c ----------------------------------------------------- c kmin = 5.d-4 ! default value c if(ibox.eq.0) kmin = 2.d0 * pi / Lbox c kmax = 10.d0 ! default value c do ik=1, ik_max akk(ik) = ak(ik) pkk(ik) = pk(ik) enddo c ikk = 1 do ik=1, ik_max if(akk(ik).ge.kmin .and. akk(ik).le.kmax & .and. mod(ik,2).eq.0 ) then ak(ikk) = akk(ik) pk(ikk) = pkk(ik) ikk = ikk + 1 endif enddo c ik_max = ikk -1 c write(6,*) 'ak(1)=', ak(1) write(6,*) 'ak(ik_max)=', ak(ik_max) write(6,*) 'ik_max=', ik_max pause c end c c ******************************************************* c c subroutine find_pk(kk, pklin) c c ******************************************************* c c implicit none integer ik_max, ikmax integer j, jmin, jmax parameter(ikmax=5000) real*8 ak(ikmax), pk(ikmax), kk, s, ds, pklin common /pk_data/ ak, pk, ik_max c ------------------------------------------- c call hunt(ak, ik_max, kk, j) c jmin = j - 2 jmax = j + 2 if(jmin.lt.1) jmin = 1 if(jmax.ge.ik_max) jmax = ik_max c call polint(ak(jmin),pk(jmin),jmax-jmin+1,kk,s,ds) pklin = s c end c c ******************************************************* c c subroutine calc_one_loop_pk c c ******************************************************* c c implicit none integer ik, ik_max, ikmax parameter(ikmax=5000) integer ndim3, ncall13, ncall22, nboost integer idum, nprn, init, itmax real*8 ak(ikmax), pk(ikmax), pk13(ikmax), pk22(ikmax) real*8 pi, k, kmin, kmax, pklink, k_cut cc real*8 avgi, sd, chi2a, region(20) real*8 avgi, region(20) real*8 sd13, sd22, chi2a13, chi2a22 external fp13, fp22 common /pk_data/ ak, pk, ik_max common /wave_number/ k, kmin, kmax COMMON /ranno/ idum pi = 4.d0 * datan(1.d0) c --------------------------------------------------- avgi=0.d0 sd13=0.d0 sd22=0.d0 chi2a13=0.d0 chi2a22=0.d0 ndim3 = 3 itmax = 10 ncall13 = 2000 ncall22 = 2000 nprn = 0 c --------------------------------------------------- c write(6,*) ' input integer idum' read(5,*) idum idum = -abs(idum) write(6,*) ' type k_min ' read(5,*) k_cut c open(10,file='stdPT_2loop.dat',status='unknown', & access='sequential') c kmin = ak(1) kmax = ak(ik_max) c do 10 ik=1, ik_max c k = ak(ik) call find_pk(k, pklink) c if(k.lt.0.1) nboost = 0 if(k.ge.0.1) nboost = 4000 if(k.ge.0.25) nboost = 8000 cc if(k.ge.0.3) nboost = 4000 c pk13(ik)= 0.d0 pk22(ik)= 0.d0 c if(k.le.k_cut) goto 5 c c ////// Monte Carlo integration of p13(k) p22(k) ////// c c region(1) = dlog(kmin) region(1+ndim3) = dlog(kmax) region(2) = 0.d0 region(2+ndim3) = pi region(3) = 0.d0 region(3+ndim3) = 2.d0*pi c init = -1 call vegas(region, ndim3, fp13, init, ncall13+nboost, itmax, & nprn, avgi, sd13, chi2a13) init = 1 call vegas(region, ndim3, fp13, init, ncall13+nboost, itmax, & nprn, avgi, sd13, chi2a13) c pk13(ik) = 6.d0 * pklink * avgi / (2.d0*pi)**3 c c init = -1 call vegas(region, ndim3, fp22, init, ncall22+nboost, itmax, & nprn, avgi, sd22, chi2a22) init = 1 call vegas(region, ndim3, fp22, init, ncall22+nboost, itmax, & nprn, avgi, sd22, chi2a22) c pk22(ik) = 2.d0 * avgi / (2.d0*pi)**3 c 5 write(6,'(A,i4,1p3e16.8)') 'ik, k, pk13, pk22=', & ik, k, pk13(ik), pk22(ik) c c write(10,'(1p8e16.8)') ak(ik), pk(ik), pk13(ik), pk22(ik), & sd13, sd22, chi2a13, chi2a22 c 10 continue c close(10) c end c c ******************************************************* c c subroutine calc_two_loop_pk(ijob) c c ******************************************************* c c implicit none integer ijob, ik, ik_max, ikmax, ik_init integer ix, ix_max, imu, imu_max integer ncall_boost33b, ncall_boost24, ncall_boost15 parameter(ikmax=5000) parameter(ix_max=3500, imu_max=100) integer ndim1, ndim3, ndim4, ndim6 integer ncall33a, ncall33b, ncall24, ncall15 integer idum, nprn, init, itmax real*8 ak(ikmax), pk(ikmax) real*8 pk33a(ikmax), pk33b(ikmax), pk24(ikmax), pk15(ikmax) real*8 k, kmin, kmax, pklink, pi real*8 avgi, sd, chi2a, k_cut, k_cut2 real*8 region(20) external fp33a, fp33b, fp24, fp15, fp33b_4D, fp33b2 common /pk_data/ ak, pk, ik_max common /wave_number/ k, kmin, kmax COMMON /ranno/ idum pi = 4.d0 * datan(1.d0) c --------------------------------------------------- avgi=0.d0 sd=0.d0 chi2a=0.d0 ndim1 = 1 ndim3 = 3 ndim4 = 4 ndim6 = 6 itmax = 10 ncall33a = 2000 ncall33b = 15000 ncall24 = 7500 ncall15 = 15000 nprn = 0 c --------------------------------------------------- c write(6,*) ' input integer idum' read(5,*) idum idum = -abs(idum) write(6,*) ' type k_min ' read(5,*) k_cut write(6,*) ' type k_max ' read(5,*) k_cut2 c if(ijob.eq.0) then call search(ik_init) open(10,file='stdPT_2loop.dat',status='unknown', & access='append') else ik_init = 1 open(10,file='stdPT_2loop.dat',status='unknown', & access='sequential') endif c kmin = ak(1) kmax = ak(ik_max) c do 10 ik=ik_init, ik_max c k = ak(ik) call find_pk(k, pklink) c c ****** boost factor ****** c c if(k.le.0.1) then ncall_boost33b = 0 ncall_boost24 = 0 ncall_boost15 = 0 elseif(k.gt.0.1) then ncall_boost33b = ncall33b*2 ncall_boost24 = ncall24 ncall_boost15 = ncall15*2 elseif(k.gt.0.2) then ncall_boost33b = ncall33b*6 ncall_boost24 = ncall24*4 ncall_boost15 = ncall15*3 elseif(k.gt.0.3) then ncall_boost33b = ncall33b*9 ncall_boost24 = ncall24*8 ncall_boost15 = ncall15*6 endif c pk33a(ik)= 0.0d0 pk33b(ik)= 0.0d0 pk24(ik)= 0.0d0 pk15(ik)= 0.0d0 c if(k.le.k_cut .or. k.ge.k_cut2) goto 5 c c ////// 1D Monte Carlo integration ////// c c region(1) = dlog(kmin) region(1+ndim1) = dlog(kmax) c c */*/*/*/*/* p33(k) A-part */*/*/*/*/*/* c c init = -1 call vegas(region, ndim1, fp33a, init, ncall33a, itmax, & nprn, avgi, sd, chi2a) init = 1 call vegas(region, ndim1, fp33a, init, ncall33a, itmax, & nprn, avgi, sd, chi2a) c pk33a(ik) = avgi / (2.*pi)**2 pk33a(ik) = pklink * pk33a(ik)**2 c c ////// 6D Monte Carlo integration ////// c c region(1) = dlog(kmin) region(1+ndim6) = dlog(kmax) region(2) = 0.d0 region(2+ndim6) = pi region(3) = 0.d0 region(3+ndim6) = 2.d0*pi c region(4) = dlog(kmin) region(4+ndim6) = dlog(kmax) region(5) = 0.d0 region(5+ndim6) = pi region(6) = 0.d0 region(6+ndim6) = 2.d0*pi c c */*/*/*/*/* p33(k) B-part */*/*/*/*/*/* c c init = -1 cc call vegas(region, ndim6, fp33b, init, call vegas(region, ndim6, fp33b2, init, & ncall33b+ncall_boost33b, itmax, nprn, avgi, sd, chi2a) init = 1 cc call vegas(region, ndim6, fp33b, init, call vegas(region, ndim6, fp33b2, init, & ncall33b+ncall_boost33b, itmax, nprn, avgi, sd, chi2a) c pk33b(ik) = 6.d0 * avgi / (2.d0*pi)**6 c c */*/*/*/*/* p24(k) */*/*/*/*/*/* c c init = -1 call vegas(region, ndim6, fp24, init, & ncall24+ncall_boost24, itmax, nprn, avgi, sd, chi2a) init = 1 call vegas(region, ndim6, fp24, init, & ncall24+ncall_boost24, itmax, nprn, avgi, sd, chi2a) c pk24(ik) = 24.d0 * avgi / (2.d0*pi)**6 c c */*/*/*/*/* p15(k) */*/*/*/*/*/* c c init = -1 call vegas(region, ndim6, fp15, init, & ncall15+ncall_boost15, itmax, nprn, avgi, sd, chi2a) init = 1 call vegas(region, ndim6, fp15, init, & ncall15+ncall_boost15, itmax, nprn, avgi, sd, chi2a) c pk15(ik) = 30.d0 * pklink * avgi / (2.d0*pi)**6 c c ////// 6D Monte Carlo integration ////// c c cc region(1) = dlog(kmin/kmax) cc region(1+ndim4) = dlog(kmax/k+1.d0) cc region(2) = -1.d0 cc region(2+ndim4) = 1.d0 cc region(3) = dlog(kmin/kmax) cc region(3+ndim4) = dlog(kmax/kmin) cc region(4) = -1.d0 cc region(4+ndim4) = 1.d0 c c */*/*/*/*/* p33(k) B-part */*/*/*/*/*/* c c cc init = -1 cc call vegas(region, ndim4, fp33b_4D, init, cc & ncall33b, itmax, nprn, avgi, sd, chi2a) cc init = 1 cc call vegas(region, ndim4, fp33b_4D, init, cc & ncall33b, itmax, nprn, avgi, sd, chi2a) c cc pk33b(ik) = 6.d0 * (4.d0 * avgi) * k**6 / (2.d0*pi)**4 c c ////////// Output //////////// c c 5 write(6,'(A,i4,1p5e12.4)') 'ik, k, pk3a, pk33b, pk24, pk15=', & ik, k, pk33a(ik), pk33b(ik), pk24(ik), pk15(ik) c write(10,'(1p6e16.8)') ak(ik), pk(ik), & pk33a(ik), pk33b(ik), pk24(ik), pk15(ik) c 10 continue c close(10) c end c c ******************************************************* c c subroutine search(ik_init) c c ******************************************************* c c c search the number record of last output data c implicit none c integer ik_init, ik, ikmax parameter(ikmax=5000) real*8 ak_dummy, dummy c ------------------------------------------ c open(10,file='stdPT_2loop.dat',status='unknown') c do ik=1, ikmax read(10,*,END=10) ak_dummy, dummy, dummy, dummy, dummy, dummy enddo c 10 ik_init = ik c write(6,*) ik_init, ak_dummy c close(10) c end c c ******************************************************* c c function fp13(pt, wgt) c c ******************************************************* c c c direct 3D integration of fp13 c implicit none integer i real*8 fp13, k, pt(*), wgt, p, pp(3), m_pp(3), kk(3) real*8 kmin, kmax, pklin, F3_sym common /wave_number/ k, kmin, kmax c p = dexp(pt(1)) c pp(1) = p * dsin(pt(2)) * dcos(pt(3)) pp(2) = p * dsin(pt(2)) * dsin(pt(3)) pp(3) = p * dcos(pt(2)) c m_pp(1) = -pp(1) m_pp(2) = -pp(2) m_pp(3) = -pp(3) c kk(1) = 0.d0 kk(2) = 0.d0 kk(3) = k c call find_pk(p, pklin) c fp13 = F3_sym(1, kk, pp, m_pp) * pklin fp13 = fp13 * p**3 * dsin(pt(2)) c end c c ******************************************************* c c function fp22(pt, wgt) c c ******************************************************* c c c direct 3D integration of fp13 c implicit none integer i real*8 fp22, k, pt(*), wgt, p, k_p real*8 pp(3), kp(3), prod_kp real*8 kmin, kmax, pklin_p, pklin_kp, F2_sym common /wave_number/ k, kmin, kmax c p = dexp(pt(1)) c pp(1) = p * dsin(pt(2)) * dcos(pt(3)) pp(2) = p * dsin(pt(2)) * dsin(pt(3)) pp(3) = p * dcos(pt(2)) c kp(1) = -pp(1) kp(2) = -pp(2) kp(3) = k-pp(3) c prod_kp = pp(3) * k c k_p = dsqrt( kp(1)**2 + kp(2)**2 + kp(3)**2 ) c if(k_p.ge.kmin .and. k_p.le.kmax) then c if(prod_kp.le.(0.5d0*k*k)) then call find_pk(p, pklin_p) call find_pk(k_p, pklin_kp) fp22 = F2_sym(1, pp, kp)**2 * pklin_p * pklin_kp fp22 = fp22 * p**3 * dsin(pt(2)) else fp22 = 0.d0 endif c else c fp22 = 0.d0 c endif c fp22= fp22 * 2.d0 c end c c ******************************************************* c c function fp33a(pt, wgt) c c ******************************************************* c c c This is nothing but fp13/2 c implicit none real*8 fp33a, x, k, s, pklin, pt(*) real*8 wgt, kmin, kmax common /wave_number/ k, kmin, kmax c x = dexp(pt(1)) / k fp33a = 0.0 c if(x.lt.0.) then write(6,*) ' warning: fp33 negative !! ' fp33a = 0.0 elseif(x .lt. 5.d-3) then fp33a = -2./3. + 232./315.*x*x - 376./735.*x*x*x*x elseif(x .gt. 5.d2) then s = 1.d0 / x fp33a = -122./315. + 8./105.*s*s -40./1323.*s*s*s*s elseif(x.ge.0.995 .and. x.le.1.005) then fp33a = ( -22. + 2.*(x-1.) -29.*(x-1.)**2 )/63. else fp33a = ( 12./x/x - 158. + 100.*x*x -42.*x*x*x*x + & 3./(x*x*x) * (x*x-1.)**3 * (7.*x*x + 2.) * & dlog(abs((x+1.)/(x-1.))) )/252. endif c call find_pk(k*x, pklin) c fp33a = fp33a * pklin * x * k**3 / 2.d0 c end c c ******************************************************* c c function fp33b(pt, wgt) c c ******************************************************* c c implicit none real*8 fp33b, k, kmin, kmax real*8 wgt, pt(*), F3_sym real*8 pp(3), qq(3), kpq(3), p, q, k_p_q real*8 pklin_p, pklin_q, pklin_kpq common /wave_number/ k, kmin, kmax c p = dexp(pt(1)) q = dexp(pt(4)) c pp(1) = p * dsin(pt(2)) * dcos(pt(3)) pp(2) = p * dsin(pt(2)) * dsin(pt(3)) pp(3) = p * dcos(pt(2)) qq(1) = q * dsin(pt(5)) * dcos(pt(6)) qq(2) = q * dsin(pt(5)) * dsin(pt(6)) qq(3) = q * dcos(pt(5)) kpq(1) = -pp(1)-qq(1) kpq(2) = -pp(2)-qq(2) kpq(3) = k-pp(3)-qq(3) c k_p_q = dsqrt( kpq(1)**2 + kpq(2)**2 + kpq(3)**2 ) c if(k_p_q.ge.kmin .and. k_p_q.le.kmax) then c call find_pk(p, pklin_p) call find_pk(q, pklin_q) call find_pk(k_p_q, pklin_kpq) c fp33b = F3_sym(1, pp, qq, kpq)**2 & * pklin_p * pklin_q * pklin_kpq fp33b = fp33b * p**3 * q**3 * dsin(pt(2)) * dsin(pt(5)) c else c fp33b = 0.d0 c endif c end c c c ******************************************************* c c function fp33b2(pt, wgt) c c ******************************************************* c c implicit none real*8 fp33b2, k, kmin, kmax real*8 wgt, pt(*), F3_sym, kp(3), k_p, prod_q_kp real*8 pp(3), qq(3), kpq(3), p, q, k_p_q real*8 pklin_p, pklin_q, pklin_kpq common /wave_number/ k, kmin, kmax c p = dexp(pt(1)) q = dexp(pt(4)) c pp(1) = p * dsin(pt(2)) * dcos(pt(3)) pp(2) = p * dsin(pt(2)) * dsin(pt(3)) pp(3) = p * dcos(pt(2)) qq(1) = q * dsin(pt(5)) * dcos(pt(6)) qq(2) = q * dsin(pt(5)) * dsin(pt(6)) qq(3) = q * dcos(pt(5)) kpq(1) = -pp(1)-qq(1) kpq(2) = -pp(2)-qq(2) kpq(3) = k-pp(3)-qq(3) c k_p_q = dsqrt( kpq(1)**2 + kpq(2)**2 + kpq(3)**2 ) c kp(1) = -pp(1) kp(2) = -pp(2) kp(3) = k-pp(3) k_p = dsqrt( kp(1)**2 + kp(2)**2 + kp(3)**2 ) prod_q_kp = qq(1)*kp(1) + qq(2)*kp(2) + qq(3)*kp(3) c if(k_p_q.ge.kmin .and. k_p_q.le.kmax) then c if(prod_q_kp .le. 0.5d0*k_p*k_p) then c call find_pk(p, pklin_p) call find_pk(q, pklin_q) call find_pk(k_p_q, pklin_kpq) c fp33b2 = F3_sym(1, pp, qq, kpq)**2 & * pklin_p * pklin_q * pklin_kpq fp33b2 = fp33b2 * p**3 * q**3 * dsin(pt(2)) * dsin(pt(5)) c else fp33b2 = 0.d0 endif c else c fp33b2 = 0.d0 c endif c fp33b2= 2.d0 * fp33b2 c end c c ******************************************************* c c function fp33b_4D(pt, wgt) c c ******************************************************* c c implicit none real*8 fp33b_4D, k, kmin, kmax real*8 wgt, pt(*) real*8 x, mu, y, nu, p, q, pq, kpq, kq real*8 fp33b1, fp33b2, integ_fp33b real*8 pklin_p, pklin_q, pklin_kpq, pklin_pq, pklin_kq common /wave_number/ k, kmin, kmax c x = dexp(pt(1)) mu = pt(2) y = dexp(pt(3)) nu = pt(4) c if( (1.d0 + x*x - 2.d0*x*mu .ge.0.d0).and. & (1.d0 + y*y - 2.d0*y*nu .ge.0.d0) & .and. (x.lt.0.5d0 .or. & (x.ge.0.5d0 .and. (mu.lt.0.5d0/x)) ) ) then c c ///// part I ///// c c q = k * x p = k * y * dsqrt(1.d0 + x*x - 2.d0*x*mu) kpq = k * dsqrt(1.d0 + x*x - 2.d0*x*mu) & * dsqrt(1.d0 + y*y - 2.d0*y*nu) c if( ((q.ge.kmin).and.(q.le.kmax)) .and. & ((p.ge.kmin).and.(p.le.kmax)) .and. & ((kpq.ge.kmin).and.(kpq.le.kmax)) ) then if( y.lt.0.5d0 .or. & (y.ge.0.5d0 .and. (nu.lt.0.5d0/y)) ) then fp33b1 = x**2 * (1.d0+x*x-2.d0*x*mu)**1.5 * y**2 & * integ_fp33b(1, x, mu, y, nu)**2 call find_pk(q, pklin_q) call find_pk(p, pklin_p) call find_pk(kpq, pklin_kpq) c fp33b1 = fp33b1 * pklin_q * pklin_p * pklin_kpq else fp33b1 = 0.d0 endif else fp33b1 = 0.d0 endif c c ///// part II ///// c c p = k * x * y pq = k * x * dsqrt(1.d0 + y*y - 2.d0*y*nu) kq = k * dsqrt(1.d0 + x*x - 2.d0*x*mu) c if( ((p.ge.kmin).and.(p.le.kmax)) .and. & ((pq.ge.kmin).and.(pq.le.kmax)) .and. & ((kq.ge.kmin).and.(kq.le.kmax)) ) then if( y.lt.0.5d0 .or. & (y.ge.0.5d0 .and. (nu.lt.0.5d0/y)) ) then fp33b2 = x**5 * y**2 * integ_fp33b(2, x, mu, y, nu)**2 call find_pk(p, pklin_p) call find_pk(pq, pklin_pq) call find_pk(kq, pklin_kq) c fp33b2 = fp33b2 * pklin_p * pklin_kq * pklin_pq else fp33b2 = 0.d0 endif else fp33b2 = 0.d0 endif else fp33b1 = 0.0d0 fp33b2 = 0.0d0 endif c fp33b_4D = (fp33b1 + fp33b2) * x * y / 2.d0 c end c c ******************************************************* c c function integ_fp33b(a, x, mu, y, nu) c c ******************************************************* c c implicit none integer a real*8 x, mu, y, nu, AA, BB, sigmaab, integ_fp33b real*8 gamma112, gamma121, gamma222, F2_1, F2_2 c F2_1 = 5.d0/7.d0 + 0.5d0 * (nu-y) / dsqrt(1.d0+y*y-2.d0*y*nu) * & ( y/dsqrt(1.d0+y*y-2.d0*y*nu) + dsqrt(1.d0+y*y-2.d0*y*nu)/y ) & + 2.d0/7.d0 * (nu-y)**2 / (1.d0+y*y-2.d0*y*nu) F2_2 = 3.d0/7.d0 + 0.5d0 * (nu-y) / dsqrt(1.d0+y*y-2.d0*y*nu) * & ( y/dsqrt(1.d0+y*y-2.d0*y*nu) + dsqrt(1.d0+y*y-2.d0*y*nu)/y ) & + 4.d0/7.d0 * (nu-y)**2 / (1.d0+y*y-2.d0*y*nu) c AA = 0.5d0 * (1.d0 - x*mu) / (1.d0 + x*x - 2.d0*x*mu) BB = 0.5d0 * mu / x c if(a.eq.1) then gamma112 = AA gamma121 = BB elseif(a.eq.2) then gamma112 = BB gamma121 = AA endif gamma222 = 0.5d0 * (mu - x) / x / (1.d0 + x*x - 2.d0*x*mu) c integ_fp33b = ( sigmaab(3,1,1) * gamma112 & + sigmaab(3,1,2) * gamma222 ) * F2_2 & + sigmaab(3,1,1) * gamma121 * F2_1 c integ_fp33b = integ_fp33b * 2.d0 c end c c ******************************************************* c c function fp24(pt, wgt) c c ******************************************************* c c implicit none real*8 fp24, k, kmin, kmax real*8 wgt, pt(*), F2_sym, F4_sym real*8 pp(3), qq(3), kp(3), m_qq(3), p, q, k_p real*8 pklin_p, pklin_q, pklin_kp common /wave_number/ k, kmin, kmax c p = dexp(pt(1)) q = dexp(pt(4)) c pp(1) = p * dsin(pt(2)) * dcos(pt(3)) pp(2) = p * dsin(pt(2)) * dsin(pt(3)) pp(3) = p * dcos(pt(2)) c qq(1) = q * dsin(pt(5)) * dcos(pt(6)) qq(2) = q * dsin(pt(5)) * dsin(pt(6)) qq(3) = q * dcos(pt(5)) c kp(1) = -pp(1) kp(2) = -pp(2) kp(3) = k-pp(3) c m_qq(1) = -qq(1) m_qq(2) = -qq(2) m_qq(3) = -qq(3) c k_p = dsqrt( kp(1)**2 + kp(2)**2 + kp(3)**2 ) c if(k_p.ge.kmin .and. k_p.le.kmax) then c call find_pk(p, pklin_p) call find_pk(q, pklin_q) call find_pk(k_p, pklin_kp) c fp24 = F2_sym(1, kp, pp) * F4_sym(1, kp, pp, qq, m_qq) & * pklin_p * pklin_q * pklin_kp fp24 = fp24 * p**3 * q**3 * dsin(pt(2)) * dsin(pt(5)) c else c fp24 = 0.d0 c endif c end c c ******************************************************* c c function fp15(pt, wgt) c c ******************************************************* c c implicit none real*8 fp15, k, kmin, kmax real*8 wgt, pt(*), F5_sym real*8 kk(3), pp(3), qq(3), m_pp(3), m_qq(3), p, q real*8 pklin_p, pklin_q common /wave_number/ k, kmin, kmax c p = dexp(pt(1)) q = dexp(pt(4)) c pp(1) = p * dsin(pt(2)) * dcos(pt(3)) pp(2) = p * dsin(pt(2)) * dsin(pt(3)) pp(3) = p * dcos(pt(2)) c qq(1) = q * dsin(pt(5)) * dcos(pt(6)) qq(2) = q * dsin(pt(5)) * dsin(pt(6)) qq(3) = q * dcos(pt(5)) c m_pp(1) = -pp(1) m_pp(2) = -pp(2) m_pp(3) = -pp(3) c m_qq(1) = -qq(1) m_qq(2) = -qq(2) m_qq(3) = -qq(3) c kk(1) = 0.d0 kk(2) = 0.d0 kk(3) = k c call find_pk(p, pklin_p) call find_pk(q, pklin_q) c fp15 = F5_sym(1, pp, qq, kk, m_pp, m_qq) * pklin_p * pklin_q c fp15 = fp15 * p**3 * q**3 * dsin(pt(2)) * dsin(pt(5)) c cc write(6,*) 'in fp15:',F5_sym(1, pp, qq, kk, m_pp, m_qq), cc & F5_sym(1, m_pp, m_qq, qq, pp, kk), cc & F5_sym(1, pp, m_pp, m_qq, kk, qq), cc & F5_sym(1, pp, m_qq, kk, m_pp, qq) c end c c ************************************************************ c SUBROUTINE polint(xa,ya,n,x,y,dy) c c ************************************************************ c implicit none INTEGER n,NMAX REAL*8 dy,x,y,xa(n),ya(n) PARAMETER (NMAX=10) INTEGER i,m,ns REAL*8 den,dif,dift,ho,hp,w,c(NMAX),d(NMAX) ns=1 dif=abs(x-xa(1)) do 11 i=1,n dift=abs(x-xa(i)) if (dift.lt.dif) then ns=i dif=dift endif c(i)=ya(i) d(i)=ya(i) 11 continue y=ya(ns) ns=ns-1 do 13 m=1,n-1 do 12 i=1,n-m ho=xa(i)-x hp=xa(i+m)-x w=c(i+1)-d(i) den=ho-hp if(den.eq.0.)pause 'failure in polint' den=w/den d(i)=hp*den c(i)=ho*den 12 continue if (2*ns.lt.n-m)then dy=c(ns+1) else dy=d(ns) ns=ns-1 endif y=y+dy 13 continue return END C c ************************************************ c c SUBROUTINE hunt(xx,n,x,jlo) c c ************************************************ c c implicit none 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 ******************************************************* c c SUBROUTINE vegas(region,ndim,fxn,init,ncall,itmx,nprn,tgral,sd, *chi2a) c c ******************************************************* c c INTEGER init,itmx,ncall,ndim,nprn,NDMX,MXDIM REAL*8 tgral,chi2a,sd,region(2*ndim),fxn,ALPH,TINY PARAMETER (ALPH=1.5,NDMX=50,MXDIM=10,TINY=1.d-30) EXTERNAL fxn CU USES fxn,ran2,rebin INTEGER i,idum,it,j,k,mds,nd,ndo,ng,npg,ia(MXDIM),kg(MXDIM) REAL*8 calls,dv2g,dxg,f,f2,f2b,fb,rc,ti,tsi,wgt,xjac,xn,xnd,xo, *d(NDMX,MXDIM),di(NDMX,MXDIM),dt(MXDIM),dx(MXDIM),r(NDMX),x(MXDIM), *xi(NDMX,MXDIM),xin(NDMX),ran2 DOUBLE PRECISION schi,si,swgt COMMON /ranno/ idum SAVE if(init.le.0)then mds=1 ndo=1 do 11 j=1,ndim xi(1,j)=1. 11 continue endif if (init.le.1)then si=0. swgt=0. schi=0. endif if (init.le.2)then nd=NDMX ng=1 if(mds.ne.0)then ng=(ncall/2.+0.25)**(1./ndim) mds=1 if((2*ng-NDMX).ge.0)then mds=-1 npg=ng/NDMX+1 nd=ng/npg ng=npg*nd endif endif k=ng**ndim npg=max(ncall/k,2) calls=npg*k dxg=1./ng dv2g=(calls*dxg**ndim)**2/npg/npg/(npg-1.) xnd=nd dxg=dxg*xnd xjac=1./calls do 12 j=1,ndim dx(j)=region(j+ndim)-region(j) xjac=xjac*dx(j) 12 continue if(nd.ne.ndo)then do 13 i=1,nd r(i)=1. 13 continue do 14 j=1,ndim call rebin(ndo/xnd,nd,r,xin,xi(1,j)) 14 continue ndo=nd endif ccc if(nprn.ge.0) write(*,200) ndim,calls,it,itmx,nprn,ALPH,mds,nd, ccc *(j,region(j),j,region(j+ndim),j=1,ndim) endif do 28 it=1,itmx ti=0. tsi=0. do 16 j=1,ndim kg(j)=1 do 15 i=1,nd d(i,j)=0. di(i,j)=0. 15 continue 16 continue 10 continue fb=0. f2b=0. do 19 k=1,npg wgt=xjac do 17 j=1,ndim xn=(kg(j)-ran2(idum))*dxg+1. ia(j)=max(min(int(xn),NDMX),1) if(ia(j).gt.1)then xo=xi(ia(j),j)-xi(ia(j)-1,j) rc=xi(ia(j)-1,j)+(xn-ia(j))*xo else xo=xi(ia(j),j) rc=(xn-ia(j))*xo endif x(j)=region(j)+rc*dx(j) wgt=wgt*xo*xnd 17 continue f=wgt*fxn(x,wgt) f2=f*f fb=fb+f f2b=f2b+f2 do 18 j=1,ndim di(ia(j),j)=di(ia(j),j)+f if(mds.ge.0) d(ia(j),j)=d(ia(j),j)+f2 18 continue 19 continue f2b=sqrt(f2b*npg) f2b=(f2b-fb)*(f2b+fb) if (f2b.le.0.) f2b=TINY ti=ti+fb tsi=tsi+f2b if(mds.lt.0)then do 21 j=1,ndim d(ia(j),j)=d(ia(j),j)+f2b 21 continue endif do 22 k=ndim,1,-1 kg(k)=mod(kg(k),ng)+1 if(kg(k).ne.1) goto 10 22 continue tsi=tsi*dv2g wgt=1./tsi si=si+dble(wgt)*dble(ti) schi=schi+dble(wgt)*dble(ti)**2 swgt=swgt+dble(wgt) tgral=si/swgt chi2a=max((schi-si*tgral)/(it-.99d0),0.d0) sd=sqrt(1./swgt) tsi=sqrt(tsi) if(nprn.ge.0)then cccc write(*,201) it,ti,tsi,tgral,sd,chi2a if(nprn.ne.0)then do 23 j=1,ndim cccc write(*,202) j,(xi(i,j),di(i,j),i=1+nprn/2,nd,nprn) 23 continue endif endif do 25 j=1,ndim xo=d(1,j) xn=d(2,j) d(1,j)=(xo+xn)/2. dt(j)=d(1,j) do 24 i=2,nd-1 rc=xo+xn xo=xn xn=d(i+1,j) d(i,j)=(rc+xn)/3. dt(j)=dt(j)+d(i,j) 24 continue d(nd,j)=(xo+xn)/2. dt(j)=dt(j)+d(nd,j) 25 continue do 27 j=1,ndim rc=0. do 26 i=1,nd if(d(i,j).lt.TINY) d(i,j)=TINY r(i)=((1.-d(i,j)/dt(j))/(log(dt(j))-log(d(i,j))))**ALPH rc=rc+r(i) 26 continue call rebin(rc/xnd,nd,r,xin,xi(1,j)) 27 continue 28 continue return c200 FORMAT(/' input parameters for vegas: ndim=',i3,' ncall=', c *f8.0/28x,' it=',i5,' itmx=',i5/28x,' nprn=',i3,' alph=', c *f5.2/28x,' mds=',i3,' nd=',i4/(30x,'xl(',i2,')= ',g11.4,' xu(', c *i2,')= ',g11.4)) c201 FORMAT(/' iteration no.',I3,': ','integral =',g14.7,'+/- ',g9.2/ c *' all iterations: integral ='g14.7,'+/- ',g9.2,' chi**2/it' c *'n =',g9.2) c202 FORMAT(/' data for axis ',I2/' X delta i ', c *' x delta i ',' x delta i ',/(1x, c *f7.5,1x,g11.4,5x,f7.5,1x,g11.4,5x,f7.5,1x,g11.4)) END c c ******************************************************* c c SUBROUTINE rebin(rc,nd,r,xin,xi) c c ******************************************************* c c INTEGER nd REAL*8 rc,r(*),xi(*),xin(*) INTEGER i,k REAL*8 dr,xn,xo k=0 xn=0. dr=0. do 11 i=1,nd-1 1 if(rc.gt.dr)then k=k+1 dr=dr+r(k) xo=xn xn=xi(k) goto 1 endif dr=dr-rc xin(i)=xn-(xn-xo)*dr/r(k) 11 continue do 12 i=1,nd-1 xi(i)=xin(i) 12 continue xi(nd)=1. return END c c ******************************************************* c c FUNCTION ran2(idum) c c ******************************************************* c c INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV REAL*8 ran2,AM,EPS,RNMX PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1, *IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791, *NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2d-7,RNMX=1.-EPS) INTEGER idum2,j,k,iv(NTAB),iy SAVE iv,iy,idum2 DATA idum2/123456789/, iv/NTAB*0/, iy/0/ if (idum.le.0) then idum=max(-idum,1) idum2=idum do 11 j=NTAB+8,1,-1 k=idum/IQ1 idum=IA1*(idum-k*IQ1)-k*IR1 if (idum.lt.0) idum=idum+IM1 if (j.le.NTAB) iv(j)=idum 11 continue iy=iv(1) endif k=idum/IQ1 idum=IA1*(idum-k*IQ1)-k*IR1 if (idum.lt.0) idum=idum+IM1 k=idum2/IQ2 idum2=IA2*(idum2-k*IQ2)-k*IR2 if (idum2.lt.0) idum2=idum2+IM2 j=1+iy/NDIV iy=iv(j)-idum2 iv(j)=idum if(iy.lt.1)iy=iy+IMM1 ran2=min(AM*iy,RNMX) return END C (C) Copr. 1986-92 Numerical Recipes Software z!0(0.