!This module provides the initial power spectra given by sets of slow-roll parameters. AL 29/3/01.
!Based on astro-ph/9911225

!Fill arrays slow_roll_epsilon(nnmax), slow_roll_delta(nnmax) with nn sets of parameters
!The output is not COBE normalized and the scalar and tensor results are in the ratio predicted
!by slow-roll inflation

     module InitialPower   
   
     implicit none   

     private
    
      real*8 curv  !Curvature contant, set in InitializePowers     
     
!How we want to compute the Cls
!Essential variables

   logical UseScalTensRatio  
      !Normalize Cl to fixed quadrupole ratio? Definitely not- here we set the initial power spectra amplitudes.

     integer, parameter :: outCOBE=0, outNone=1
     integer OutputNormalization 

!Change the OutputNormalization parameter to determine the normalization of the output Cls

  !If outCOBE then the output is l(l+1)Cl/2pi normalized to COBE like CMBFAST

  !If outNone the initial power spectra (ScalarPower and TensorPower) should return absolute values
  !and no normalization takes place. Only makes sense if  UseScalTensRatio = false
  !The output in this case is l(l+1)Cl/2pi.  
  !If you are only interested in the relative normalization of the input scalar and tensor power spectra you can use
  !UseScalTensRatio=false with outCOBE to fix the total overal normalization

  !For any other value the scalar temperature output is normalized to 1 at l=OutputNormalization
    
     real*8, parameter ::  OutputDenominator =2*3.1415927d0 
      !When using outNone the output is l(l+1)Cl/OutputDenominator

     integer, parameter :: nnmax= 10 
       !Maximum possible number of different power spectra to use

     integer nn 
      !The actual number of power spectra to use

     real*8 rat(nnmax) !ratio of scalar to tensor quadrupole, ignore
     real*8 an(nnmax),ant(nnmax) !After computation contains spectral indices
     !ignore
         
!Implementation specific variables

! Slow roll parameters
! See asto-ph/9911225

      real*8 slow_roll_epsilon(nnmax), slow_roll_delta(nnmax)

! Pivot scales for the scalar and tensor initial power spectra
! Best the same for this setup
      real*8, parameter ::  k_0_scalar = 0.05d0
      real*8, parameter ::  k_0_tensor = 0.05d0

      real*8, parameter ::  C_const = -0.7296
    
!Make things visible as neccessary...

      public  InitializePowers, ScalarPower, TensorPower, OutputNormalization, UseScalTensRatio 
      public outNone,outCOBE, OutputDenominator
      public nn,nnmax,rat,an,ant,slow_roll_epsilon, slow_roll_delta

    contains
       
    
      subroutine InitializePowers(acurv)
         !Called before computing final Cls in cmbmain.f90
         !Could read spectra from disk here, do other processing, etc.

        real*8 acurv

        curv=acurv         

!Write implementation specific code here...

        OutputNormalization = outNone
        UseScalTensRatio = .false.

      end subroutine InitializePowers
       

      function ScalarPower(k,in)

       !"in" gives the index of the power to return for this k
       !ScalarPower = const for scale invariant spectrum
       !The normalization is defined so that for adiabatic perturbations the gradient of the 3-Ricci 
       !scalar on co-moving hypersurfaces receives power
       ! < |D_a R^{(3)}|^2 > = int dk/k 16 k^6/S^6 (1-3K/k^2)^2 ScalarPower(k) 
       !In other words ScalarPower is the power spectrum of chi, the conserved quantity equal to
       !Phi + 2/3 \Omega^{-1} \frac{H^{-1}\Phi' + \Phi}{1+w} during inflation (w=p/\rho)
       !Near the end of inflation < |\chi(x)|^2 > = \int dk/k ScalarPower(k)
       !and chi is also equal to 3/2 psi.
       !Here nu^2 = (k^2 + curv)/|curv| 

        real*8 ep,delta,ScalarPower,k,ScalarPowerAmp
        integer in

        ep=slow_roll_epsilon(in)
        delta=slow_roll_delta(in)
        
        an(in) = 1.d0 - 4*ep + 2*delta
        ScalarPowerAmp = (1-2*ep -2*C_const*(2*ep-delta))/ep
        
        ScalarPower=ScalarPowerAmp*exp((an(in)-1)*log(k/k_0_scalar))    

      end function ScalarPower

      
      function TensorPower(k,in)
      
       !TensorPower= const for scale invariant spectrum
       !The normalization is defined so that
       ! < h_{ij}(x) h^{ij}(x) > = \sum_nu nu /(nu^2-1) (nu^2-4)/nu^2 TensorPower(k)
       !for a closed model
       ! < h_{ij}(x) h^{ij}(x) > = int d nu /(nu^2+1) (nu^2+4)/nu^2 TensorPower(k)
       !for an open model
       !"in" gives the index of the power spectrum to return 
       !Here nu^2 = (k^2 + 3*curv)/|curv| 


        real*8 TensorPower,k, TensorPowerAmp,ep,delta
        real*8, parameter :: PiByTwo=3.14159265d0/2.d0
   
        integer in

        ep=slow_roll_epsilon(in)
        delta=slow_roll_delta(in)
        TensorPowerAmp = 16.d0*(1-2*(C_const+1.d0)*ep)
        ant(in) = -2*ep
       
        TensorPower=TensorPowerAmp*exp(ant(in)*log(k/k_0_tensor))
        if (curv < 0) TensorPower=TensorPower*tanh(PiByTwo*sqrt(-k**2/curv-3)) 

      end function TensorPower



     end module InitialPower