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_cuba
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
ccc      infile='Plin.z=0.dat'
      infile='nishimichi_Tk.dat'
c
      Omega_m = 0.234d0
      Omega_v = 1.d0 - Omega_m
      sigma8 = 0.76d0
      n_s = 0.961d0
c
ccc      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
ccc         pk(ik) = tf_m
         pk(ik) = ak(ik)**n_s * tf_m**2
      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
c
      integer ndim, ncomp, mineval, maxeval, verbose, last
      double precision epsrel, epsabs
      parameter (ndim = 3)
      parameter (ncomp = 1)
      parameter (epsrel = 1D-4)
      parameter (epsabs = 1D-12)
      parameter (verbose = 0)
      parameter (last = 4)
      parameter (mineval = 0)
      parameter (maxeval = 50000)
      integer key
      parameter (key = 0)
      double precision integral(ncomp), error(ncomp), prob(ncomp)
      integer nregions, neval, fail
c     ---------------------------------------------------
      integer ik, ik_max, ikmax
      parameter(ikmax=5000)
      real*8  ak(ikmax), pk(ikmax), pk13(ikmax), pk22(ikmax)
      real*8  pi, k, kmin, kmax, pklink, k_cut
      external fp13, fp22
      common /pk_data/ ak, pk, ik_max
      common /wave_number/  k, kmin, kmax
      pi = 4.d0 * datan(1.d0)
c     ---------------------------------------------------
c
      write(6,*) ' type k_min ' 
      read(5,*) k_cut
c
      open(10,file='stdPT.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
         pk13(ik)= 0.d0
         pk22(ik)= 0.d0
c
         if(k.le.k_cut) goto 5
c
c     ////// Adaptive quadrature of p13(k) p22(k) //////  c
c
         call cuhre(ndim, ncomp, fp13,
     &        epsrel, epsabs, verbose + last, mineval, maxeval,
     &        key,
     &        nregions, neval, fail, integral, error, prob)
c
         pk13(ik) = 6.d0 * pklink * integral(1) / (2.d0*pi)**3 
c
	call cuhre(ndim, ncomp, fp22,
     &    epsrel, epsabs, verbose + last, mineval, maxeval,
     &    key,
     &    nregions, neval, fail, integral, error, prob)
c
         pk22(ik) = 2.d0 * integral(1) / (2.d0*pi)**3 
c
 5       write(6,*) 'ik, k, pk13, pk22=',ik, k, pk13(ik), pk22(ik)
c
         write(10,'(1p4e16.8)') ak(ik), pk(ik), pk13(ik), pk22(ik) 
c
 10   continue
c
      close(10)
c
      end
c
c ******************************************************* c
c
      subroutine calc_two_loop_pk(ijob)
c
c ******************************************************* c
c
      implicit none
c
      integer ndim6, ncomp, mineval, maxeval33b, maxeval24, maxeval15
      integer verbose, last
      double precision epsrel, epsabs
      parameter (ndim6 = 6)
      parameter (ncomp = 1)
      parameter (epsrel = 1D-6)
      parameter (epsabs = 1D-13)
      parameter (verbose = 0)
      parameter (last = 4)
      parameter (mineval = 0)
      parameter (maxeval33b=9000000,maxeval24=6000000,maxeval15=1000000)
      integer nnew33b, nnew24, nnew15
      double precision flatness33b, flatness24, flatness15
      parameter (nnew33b=6000000, nnew24=6000000, nnew15=1000000)
      parameter (flatness33b=5D0, flatness24=15D0, flatness15=15D0)
      integer key
      parameter (key = 0)
      double precision integral(ncomp), error(ncomp), prob(ncomp)
      integer nregions, neval, fail33b, fail24, fail15
c     ---------------------------------------------------
      integer ijob, ik, ik_max, ikmax, ik_init 
      parameter(ikmax=5000)
      real*8  ak(ikmax), pk(ikmax)
      real*8  pk33a(ikmax), pk33b(ikmax), pk24(ikmax), pk15(ikmax)
      real*8  error33b, error24, error15
      real*8  k, kmin, kmax, pklink, pi
      real*8  k_cut, k_cut2 
c
      integer ix, ixmax33a
      parameter(ixmax33a=300)
      real*8 fp33a, w33a(ixmax33a), x33a(ixmax33a)
c
      external fp33b, fp33b2, fp24, fp15
      common /pk_data/ ak, pk, ik_max
      common /wave_number/  k, kmin, kmax
      pi = 4.d0 * datan(1.d0)
c     ---------------------------------------------------
c
      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
         pk33a(ik)= 0.0d0
         pk33b(ik)= 0.0d0
         pk24(ik)= 0.0d0
         pk15(ik)= 0.0d0
         integral(1) = 0.d0
         error(1) = 0.d0
         fail33b = 0 
         fail24 = 0 
         fail15 = 0 
c
         if(k.le.k_cut .or. k.ge.k_cut2) goto 5
c
c     ////// 1D Monte Carlo integration //////  c
c
c     */*/*/*/*/*  p33(k) A-part  */*/*/*/*/*/*   c
c
         call gauleg(log(kmin/k),log(kmax/k),x33a,w33a,ixmax33a)
c
         do ix=1, ixmax33a
            x33a(ix)= dexp(x33a(ix))
            pk33a(ik) = pk33a(ik) + w33a(ix) * fp33a(x33a(ix))
         enddo
c
         pk33a(ik) = pk33a(ik) / (2.*pi)**2 
         pk33a(ik) = pklink * pk33a(ik)**2 
c
c     ////// 6D Monte Carlo integration //////  c
c
c
c     */*/*/*/*/*  p33(k) B-part  */*/*/*/*/*/*   c
c
ccc         call suave(ndim6, ncomp, fp33b,
ccc     &        epsrel, epsabs, verbose + last, mineval, maxeval33b,
ccc     &        nnew33b, flatness33b, 
ccc     &        nregions, neval, fail33b, integral, error, prob)
         call cuhre(ndim6, ncomp, fp33b2,
     &        epsrel, epsabs, verbose + last, mineval, maxeval33b,
     &        key,
     &        nregions, neval, fail33b, integral, error, prob)
c
         pk33b(ik) = 6.d0 * integral(1) / (2.d0*pi)**6
         error33b = 6.d0 * error(1) / (2.d0*pi)**6
c     
c     */*/*/*/*/*  p24(k)  */*/*/*/*/*/*   c
c
ccc         call suave(ndim6, ncomp, fp24,
ccc     &        epsrel, epsabs, verbose + last, mineval, maxeval24,
ccc     &        nnew24, flatness24, 
ccc     &        nregions, neval, fail24, integral, error, prob)
         call cuhre(ndim6, ncomp, fp24,
     &        epsrel, epsabs, verbose + last, mineval, maxeval24,
     &        key,
     &        nregions, neval, fail24, integral, error, prob)
c
         pk24(ik) = 24.d0 * integral(1) / (2.d0*pi)**6
         error24 = 24.d0 * error(1) / (2.d0*pi)**6
c
c     */*/*/*/*/*  p15(k)  */*/*/*/*/*/*   c
c
ccc         call suave(ndim6, ncomp, fp15,
ccc     &        epsrel, epsabs, verbose + last, mineval, maxeval15,
ccc     &        nnew15, flatness15, 
ccc     &        nregions, neval, fail15, integral, error, prob)
         call cuhre(ndim6, ncomp, fp15,
     &        epsrel, epsabs, verbose + last, mineval, maxeval15,
     &        key,
     &        nregions, neval, fail15, integral, error, prob)
c
         pk15(ik) = 30.d0 * pklink * integral(1) / (2.d0*pi)**6
         error15 = 30.d0 * pklink * error(1) / (2.d0*pi)**6
c
c
c     //////////  Output  ////////////   c
c
 5       write(6,'(A,i4,1p5e12.4)') 'ik, k, pk33a, pk33b, pk24, pk15=',
     &        ik, k, pk33a(ik), pk33b(ik), pk24(ik), pk15(ik)
c
         write(10,'(1p9e16.8)') ak(ik), pk(ik), 
     &        pk33a(ik), pk33b(ik), pk24(ik), pk15(ik), 
     &        error33b, error24, error15
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
      subroutine  fp13(ndim, xx, ncomp, ff)
c
c ******************************************************* c
c
c     direct 3D integration of fp13
c
      implicit none
c
      integer ndim, ncomp
      double precision xx(*), ff(*)
      double precision pi
      parameter (pi = 3.14159265358979323846D0)
c --------------------------------------- c      
      real*8  k, p, theta, phi, pp(3), m_pp(3), kk(3)
      real*8  kmin, kmax, pklin, F3_sym, jacobian
      common /wave_number/  k, kmin, kmax
c --------------------------------------- c
c
      p = dexp(dlog(kmin) + (dlog(kmax) - dlog(kmin)) * xx(1))
      theta = xx(2) * pi
      phi = xx(3) * 2.d0 * pi
      jacobian = (dlog(kmax)-dlog(kmin)) * pi * 2.d0 * pi
c
      pp(1) = p * dsin(theta) * dcos(phi)
      pp(2) = p * dsin(theta) * dsin(phi)
      pp(3) = p * dcos(theta) 
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
      ff(1) = F3_sym(1, kk, pp, m_pp) * pklin 
      ff(1) = ff(1) * p**3 * dsin(theta) * jacobian
c
      end
c
c ******************************************************* c
c
      subroutine  fp22(ndim, xx, ncomp, ff)
c
c ******************************************************* c
c
c     direct 3D integration of fp22
c
      implicit none
c
      integer ndim, ncomp
      double precision xx(*), ff(*)
      double precision pi
      parameter (pi = 3.14159265358979323846D0)
c --------------------------------------- c      
      real*8  theta, phi, jacobian, k, 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 --------------------------------------- c      
c
      p = dexp(dlog(kmin) + (dlog(kmax) - dlog(kmin)) * xx(1))
      theta = xx(2) * pi
      phi = xx(3) * 2.d0 * pi
      jacobian = (dlog(kmax)-dlog(kmin)) * pi * 2.d0 * pi
c
      pp(1) = p * dsin(theta) * dcos(phi)
      pp(2) = p * dsin(theta) * dsin(phi)
      pp(3) = p * dcos(theta) 
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)
            ff(1) = F2_sym(1, pp, kp)**2 * pklin_p * pklin_kp 
            ff(1) = ff(1) * p**3 * dsin(theta) * jacobian 
         else
            ff(1) = 0.d0
         endif
c     
      else
c
         ff(1) = 0.d0
c
      endif
c
      ff(1) = ff(1) * 2.d0
c
      end
c
c ******************************************************* c
c
      function  fp33a(x)
c
c ******************************************************* c
c
c     This is nothing but  fp13/2
c
      implicit none
c
      real*8 x, k, s, pklin, kmin, kmax, fp33a
      common /wave_number/  k, kmin, kmax
c --------------------------------------- c      
c
      fp33a = 0.d0
c
      if(x.lt.0.d0) then
         write(6,*) ' warning: fp33a negative !! '
         fp33a = 0.d0
      elseif(x .lt. 5.d-3) then
         fp33a = -2.d0/3.d0 + 232.d0/315.d0*x*x - 376.d0/735.d0*x*x*x*x 
      elseif(x .gt. 5.d2) then
         s = 1.d0 / x
         fp33a = -122.d0/315.d0 + 8.d0/105.d0*s*s -40.d0/1323.d0*s*s*s*s
      elseif(x.ge.0.995 .and. x.le.1.005) then
         fp33a = ( -22.d0 + 2.d0*(x-1.d0) -29.d0*(x-1.d0)**2 )/63.d0
      else 
         fp33a = ( 12.d0/x/x - 158.d0 + 100.d0*x*x -42.d0*x*x*x*x + 
     &        3.d0/(x*x*x) * (x*x-1.d0)**3 * (7.d0*x*x + 2.d0) * 
     &        dlog(abs((x+1.d0)/(x-1.d0))) )/252.d0
      endif
c
      call find_pk(k*x, pklin)
c
      fp33a = fp33a * pklin * x * k**3 / 2.d0
c
      end
c
c ******************************************************* c
c
      subroutine fp33b(ndim, xx, ncomp, ff)
c
c ******************************************************* c
c
      implicit none
c
      integer ndim, ncomp
      double precision xx(*), ff(*)
      double precision pi
      parameter (pi = 3.14159265358979323846D0)
c --------------------------------------- c      
      real*8  theta1, theta2, phi1, phi2
      real*8  k, kmin, kmax, F3_sym, jacobian
      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(dlog(kmin) + (dlog(kmax) - dlog(kmin)) * xx(1))
      theta1 = xx(2) * pi
      phi1 = xx(3) * 2.d0 * pi
      q = dexp(dlog(kmin) + (dlog(kmax) - dlog(kmin)) * xx(4))
      theta2 = xx(5) * pi
      phi2 = xx(6) * 2.d0 * pi
      jacobian = ( (dlog(kmax)-dlog(kmin)) * pi * 2.d0 * pi )**2
c
      pp(1) = p * dsin(theta1) * dcos(phi1)
      pp(2) = p * dsin(theta1) * dsin(phi1)
      pp(3) = p * dcos(theta1) 
      qq(1) = q * dsin(theta2) * dcos(phi2)
      qq(2) = q * dsin(theta2) * dsin(phi2)
      qq(3) = q * dcos(theta2) 
c
      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
         ff(1) = F3_sym(1, pp, qq, kpq)**2 
     &        * pklin_p * pklin_q * pklin_kpq
         ff(1) = ff(1) * p**3 * q**3 * dsin(theta1) * dsin(theta2)
c
      else
c
         ff(1) = 0.d0
c
      endif
c
      ff(1) = ff(1) * jacobian
c
      end
c
c
c ******************************************************* c
c
      subroutine fp33b2(ndim, xx, ncomp, ff)
c
c ******************************************************* c
c
      implicit none
c
      integer ndim, ncomp
      double precision xx(*), ff(*)
      double precision pi
      parameter (pi = 3.14159265358979323846D0)
c --------------------------------------- c      
      real*8  theta1, theta2, phi1, phi2, jacobian
      real*8  k, kmin, kmax 
      real*8  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(dlog(kmin) + (dlog(kmax) - dlog(kmin)) * xx(1))
      theta1 = xx(2) * pi
      phi1 = xx(3) * 2.d0 * pi
      q = dexp(dlog(kmin) + (dlog(kmax) - dlog(kmin)) * xx(4))
      theta2 = xx(5) * pi
      phi2 = xx(6) * 2.d0 * pi
      jacobian = ( (dlog(kmax)-dlog(kmin)) * pi * 2.d0 * pi )**2
c
      pp(1) = p * dsin(theta1) * dcos(phi1)
      pp(2) = p * dsin(theta1) * dsin(phi1)
      pp(3) = p * dcos(theta1) 
      qq(1) = q * dsin(theta2) * dcos(phi2)
      qq(2) = q * dsin(theta2) * dsin(phi2)
      qq(3) = q * dcos(theta2) 
c
      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
            ff(1) = F3_sym(1, pp, qq, kpq)**2 
     &           * pklin_p * pklin_q * pklin_kpq
            ff(1) = ff(1) * p**3 * q**3 * dsin(theta1) * dsin(theta2)
c
         else
            ff(1) = 0.d0
         endif
c
      else
c
         ff(1) = 0.d0
c
      endif
c
      ff(1)= 2.d0 * ff(1) * jacobian
c
      end
c
c ******************************************************* c
c
      subroutine fp24(ndim, xx, ncomp, ff)
c
c ******************************************************* c
c
      implicit none
c
      integer ndim, ncomp
      double precision xx(*), ff(*)
      double precision pi
      parameter (pi = 3.14159265358979323846D0)
c --------------------------------------- c      
      real*8  theta1, theta2, phi1, phi2
      real*8  k, kmin, kmax, jacobian 
      real*8  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(dlog(kmin) + (dlog(kmax) - dlog(kmin)) * xx(1))
      theta1 = xx(2) * pi
      phi1 = xx(3) * 2.d0 * pi
      q = dexp(dlog(kmin) + (dlog(kmax) - dlog(kmin)) * xx(4))
      theta2 = xx(5) * pi
      phi2 = xx(6) * 2.d0 * pi
      jacobian = ( (dlog(kmax)-dlog(kmin)) * pi * 2.d0 * pi )**2
c
      pp(1) = p * dsin(theta1) * dcos(phi1)
      pp(2) = p * dsin(theta1) * dsin(phi1)
      pp(3) = p * dcos(theta1) 
      qq(1) = q * dsin(theta2) * dcos(phi2)
      qq(2) = q * dsin(theta2) * dsin(phi2)
      qq(3) = q * dcos(theta2) 
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
         ff(1) = F2_sym(1, kp, pp) * F4_sym(1, kp, pp, qq, m_qq)
     &        * pklin_p * pklin_q * pklin_kp
         ff(1) = ff(1) * p**3 * q**3 * dsin(theta1) * dsin(theta2)
c
      else
c
         ff(1) = 0.d0
c
      endif
c
      ff(1) = ff(1) * jacobian
c
      end
c
c ******************************************************* c
c
      subroutine fp15(ndim, xx, ncomp, ff)
c
c ******************************************************* c
c
      implicit none
c
      integer ndim, ncomp
      double precision xx(*), ff(*)
      double precision pi
      parameter (pi = 3.14159265358979323846D0)
c --------------------------------------- c      
      real*8  theta1, theta2, phi1, phi2
      real*8  k, kmin, kmax, jacobian 
      real*8  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(dlog(kmin) + (dlog(kmax) - dlog(kmin)) * xx(1))
      theta1 = xx(2) * pi
      phi1 = xx(3) * 2.d0 * pi
      q = dexp(dlog(kmin) + (dlog(kmax) - dlog(kmin)) * xx(4))
      theta2 = xx(5) * pi
      phi2 = xx(6) * 2.d0 * pi
      jacobian = ( (dlog(kmax)-dlog(kmin)) * pi * 2.d0 * pi )**2
c
      pp(1) = p * dsin(theta1) * dcos(phi1)
      pp(2) = p * dsin(theta1) * dsin(phi1)
      pp(3) = p * dcos(theta1) 
      qq(1) = q * dsin(theta2) * dcos(phi2)
      qq(2) = q * dsin(theta2) * dsin(phi2)
      qq(3) = q * dcos(theta2) 
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
      ff(1) = F5_sym(1, pp, qq, kk, m_pp, m_qq) * pklin_p * pklin_q 
c
      ff(1) = ff(1) * p**3 * q**3 * dsin(theta1) * dsin(theta2) 
c
      ff(1) = ff(1) * jacobian
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
      SUBROUTINE gauleg(x1,x2,x,w,n)
c
c ************************************************************
c
      INTEGER n
      REAL*8 x1,x2,x(n),w(n)
      DOUBLE PRECISION EPS
      PARAMETER (EPS=3.d-14)
      INTEGER i,j,m
      DOUBLE PRECISION p1,p2,p3,pp,xl,xm,z,z1
      m=(n+1)/2
      xm=0.5d0*(x2+x1)
      xl=0.5d0*(x2-x1)
      do 12 i=1,m
        z=cos(3.141592654d0*(i-.25d0)/(n+.5d0))
1       continue
          p1=1.d0
          p2=0.d0
          do 11 j=1,n
            p3=p2
            p2=p1
            p1=((2.d0*j-1.d0)*z*p2-(j-1.d0)*p3)/j
11        continue
          pp=n*(z*p1-p2)/(z*z-1.d0)
          z1=z
          z=z1-p1/pp
        if(abs(z-z1).gt.EPS)goto 1
        x(i)=xm-xl*z
        x(n+1-i)=xm+xl*z
        w(i)=2.d0*xl/((1.d0-z*z)*pp*pp)
        w(n+1-i)=w(i)
12    continue
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software z!0(0.