c  Fit a spline to sigma versus M
c
c  Output the fortran data statements defining the spline to
c  sigmacdm_spline.inc.
c  This can then be inserted into the code of sigmacdm_spline.f .
c
c  Compute sigma(M) and the logarithmic slope alpha(m)
c  both from my old trusty fit and directly by integrating P(k).
c
c  This version sets Omega_0=h=Gamma=sigma_8=1. sigma(m) for
c  other values of Gamma and sigma_8 can then be computed from
c  this fit by two simple scaling of the mass and amplitude.
c  See the implementation in subroutine sigmacdm_spline.f .
c
c
c  The numerical integration is more accurate and the trusty fit
c  which is good to a few percent for M>10^10 and is called here
c  as a check.
c 
      program  make_sigma_spline
**************************************variables********************************
      implicit none
      integer N,NN,NT,i,ik,j,iargc,itrans
      real PI
      parameter (N=45,NN=1000,PI=3.1415926,NT=100000)
      real x(N),y(N),y2(N),yp1,ypn,logm,alpha,nspec,
     x err,serr,m,xx,yy,sigmacdm,a(N),a2(N),erra,serra,errm,erram,
     x sigma,lnk,lnkmin,lnkmax,dlnk,sum,pw2k3,rf,Gamma,pwdwk3,alph,suma
      character splinefile*40
*******************************************************************************
      if (iargc().ne.1) stop 'Usage: make_sigma_spline.exe nspec'
      call getreal(1,nspec)  !read the required primodial spectral index
      nspec=0.01*nint(100.0*nspec)  !round to 2 decimal places
      write(0,'(a,f4.2)') 'Adopting primordial spectal index nspec= ',nspec
      write(0,*)
c
      rf=8.0
      Gamma=1.0
      itrans=0 !BBKS transfer function formula
c
      do i=1,N
         logm= 4.0 + 13.0* real(i-1)/real(N-1)
         m=10.0**logm
         x(i)=m
c
c        Integrate  4pi k^3 P(k) W^2(k) and  4pi k^4 P(k) W(k) dW/du(kr)   
         rf=(3.0*m/(4.0*PI*2.7754e+11))**(1.0/3.0)
         lnkmax=5.0-log(rf)
         lnkmin=-9.0-log(rf)
         dlnk=(lnkmax-lnkmin)/real(NT-1)
         sum=0.0
         suma=0.0
         do ik=1,NT
            lnk=lnkmin+dlnk*real(ik-1)
            sum=sum+pw2k3(exp(lnk),rf,Gamma,nspec,itrans)
            suma=suma+pwdwk3(exp(lnk),rf,Gamma,nspec,itrans)
         end do
         lnkmin=lnkmin-dlnk
         lnkmax=lnkmax+dlnk
         sigma=(sum+0.5*pw2k3(exp(lnkmin),rf,Gamma,nspec,itrans)
     x             +0.5*pw2k3(exp(lnkmax),rf,Gamma,nspec) )*dlnk
         alph=(suma+0.5*pwdwk3(exp(lnkmin),rf,Gamma,nspec)
     x             +0.5*pwdwk3(exp(lnkmax),rf,Gamma,nspec) )*dlnk
         alph=alph/(3.0*sigma)
         sigma=68.52*sqrt(sigma)
         y(i)=sigma       !adopt values from
         a(i)=alph        !numerical integration
      end do
      yp1=2.0e+30   !indicates natural spline with zero second
      ypn=2.0e+30   !derivative at the end points.
      call spline(x,y,n,yp1,ypn,y2)
      call spline(x,a,n,yp1,ypn,a2)
      write(0,*)
c
c     Write spline coefficients to a file that later is read by
c     the subroutine in sigmacdm_spline.f
c
      write(splinefile,'(a9,f4.2,a7)') 'sigmacdm_',nspec,'.spline'
      open(10,name=splinefile,status='unknown')
      write(10,*) N  !Nspline
      write(0,*) 'M s=sigma(M) d^2s/dm^2  a=dlns/dlnm  d^2a/dm^2'
      do i=1,N
        write(0,*)  x(i),y(i),y2(i),a(i),a2(i)
        write(10,*) x(i),y(i),y2(i),a(i),a2(i)
      end do
      close(10)
c
      stop
      end
c------------------------------------------------------------------------------
c The CDM-like power spectrum multiplied by 4 PI k^3 W(kr) u.dW/du(kr)

c     
      real function pwdwk3(k,rf,Gamma,nspec)
*************************************variables*********************************
      implicit none
      real u,q,k,rf,Gamma,PI,trans,win,dwin,nspec
      parameter(PI=3.1415926)
*******************************************************************************
      u=k*rf
      q=k/Gamma
      win = 3.0* (sin(u)/u-cos(u)) / u**2
      dwin = 3.0* (3.0*cos(u)-3.0*sin(u)/u+u*sin(u)) / u**2
      trans =  (log(1.0+2.34*q)/(2.34*q))
     x / ( 1.0+3.89*q + (16.1*q)**2 + (5.46*q)**3 + (6.71*q)**4 ) **0.25
      pwdwk3=  (k**(3.0+nspec)) * trans**2 * win * dwin
      return
      end
c------------------------------------------------------------------------------
c------------------------------------------------------------------------------
c The CDM variance in spheres containing mass m.
c This fits the variance derived from the BBKS CDM power spectrum.
c The units of mass are Msol/h chosen to scale as the mass within a sphere
c of radius r Mpc/h. The scaling of the amplitude with h is such that the
c rms for a sphere of 8 Mpc/h is given by sigma_8 for any h and omega
c In these units if h is varied the characteristic radius goes like
c 1/(omega h) and the characteristic mass like 1/(h^3 omega^2)
c This version also returns the logarithmic slope, alpha, of sigma(m).
c
      real function sigmacdm(m,alpha)
*************************************variables*********************************
      implicit none
      integer ifirst
      real m,ms,a1,a2,a3,p,scla,sclm,m8,sigma,alpha,fac1,fac2,curl,curlp
      save scla,sclm,ifirst
      parameter( a1=8.6594e-12, a2=3.5 ,a3=1.628e+09, p=0.255) 
      data ifirst/1/
*******************************************************************************
      if (ifirst.eq.1) then
        sclm = 1.0
        m8=6.005e+14      !The mass within an 8Mpc/h sphere
        ms=m8
c       The fit to CDM for Gamma=1
        fac1=a2*ms**-0.1
        fac2=a3*ms**-0.63
        curl=(fac1+fac2)
        curlp=curl**p
        sigma = 1.0 / (a1*ms**0.67 * (1.0+curlp )**(1.0/p) )
        scla=1.0/sigma        !scales sigma_8 to required value
c
        ifirst=0   !indicates first call complete and sclm and scla are set
      end if
c        ----------------------------------------------------
      ms=m*sclm
c     The fit to CDM for Gamma=1
      fac1=a2*ms**-0.1
      fac2=a3*ms**-0.63
      curl=(fac1+fac2)
      curlp=curl**p
      sigmacdm=scla  / (a1*ms**0.67 * (1.0+curlp )**(1.0/p) )
      alpha=-0.67 + (0.1*fac1+0.63*fac2)/(curl*(1.0/curlp+1.0))
      return
      end
c------------------------------------------------------------------------------
c CDM-like power spectrum multiplied by k^3 W^2(k)
c
c W(k) is the transform of the top-hat window function
c
c q is the scaled wavenumber k/Gamma and the shape parameter
c Gamma =  Omega_0 h  exp( - Omega_b - Omega_b/Omega_0 )

c     
      real function pw2k3(k,rf,Gamma,nspec,itrans)
*************************************variables*********************************
      implicit none
      integer itrans
      real u,q,k,rf,Gamma,trans,win,nspec,KHORIZON
      parameter(KHORIZON=1.0/3000.0) !defined as c.H_0 in h Mpc^-1
*******************************************************************************
      u=k*rf
      q=k/Gamma
      if (itrans.eq.0) then
c        BBKS transfer function
         trans =  (log(1.0+2.34*q)/(2.34*q**2))
     x / ( (1.0/q)**4+3.89/q**3 + (16.1/q)**2 + 5.46**3/q + 6.71**4 )**0.25
      else
c        Bond & Efstathiou transfer function
         trans=(1.0+(((6.4*q)+((3.0*q)**1.5)+((1.7*q)**2.0))**1.13))**(-1.0/1.13)
      endif
      win=3.0*(sin(u)/u-cos(u))/u**2       ! top hat window function
      pw2k3= k**3 * trans**2 * win**2 * (k/KHORIZON)**nspec
      return
      end