!This module provides the initial power spectra, parameterized as
!
! P \propto (k/k_0_scalar)^(an(in)-1)
!
!for the scalar spectrum, when an(in) is the in'th spectral index. k_0_scalar
!is a pivot scale, fixed here to 0.05/Mpc (change it below as desired).
!
!This module uses the same inputs an(in), ant(in) and rat(in) as CMBFAST, however here
!rat(in) is used to set the ratio of the initial power spectra, so here
!
!** rat(in) is not the Cl quadrupole ratio ***
!
!in general models the quadrupole ratio depends in a complicated way on the ratio of the initial
!power spectra

!The absolute normalization of the Cls is unimportant here, but the relative ratio
!of the tensor and scalar Cls generated with this module will be correct for general models


!The OutputNormalization parameter controls the final output
!Absolute Cls can be obtained by setting OuputNormalization=outNone, otherwise the overall normalization
!of the power spectra doesn't matter

     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? 

     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 initial power spectrum amplitudes
     !Must have UseScalTensRatio=.false. here
     
!Implementation specific variables

      real*8, parameter :: k_0_scalar = 0.05, k_0_tensor = 0.05

    !For the default implementation return power spectra based on spectral indices
      real*8 an(nnmax) !scalar spectral indices
      real*8 ant(nnmax) !tensor spectral indices
      

!Make things visible as neccessary...

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

    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. !We set correct ratio in the initial power spectrum here

       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 ScalarPower,k
        integer in


         ScalarPower=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   
        real*8, parameter :: PiByTwo=3.14159265d0/2.d0
   
        integer in

        TensorPower=rat(in)*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