diff -r -c -b -B CAMB/cmbmain.f90 CAMB_wiggles/cmbmain.f90
*** CAMB/cmbmain.f90	2006-03-16 09:56:42.000000000 +0000
--- CAMB_wiggles/cmbmain.f90	2006-03-24 20:07:34.000000000 +0000
***************
*** 528,544 ****
  !     polarization will be calculated. For low values of k we
  !     use a logarithmic spacing. closed case dealt with by SetClosedkValues
           if (CP%WantScalars .and. CP%Reionization .and. CP%AccuratePolarization) then
!             dlnk0=2._dl/10/AccuracyBoost
              !Need this to get accurate low l polarization
           else
!             dlnk0=5._dl/10/AccuracyBoost
              if (CP%closed) dlnk0=dlnk0/2
           end if
  
           if (CP%AccurateReionization) dlnk0 = dlnk0/2
! 
!          dkn1=0.6_dl/taurst/AccuracyBoost 
!          dkn2=1.1_dl/taurst/AccuracyBoost 
            if (CP%WantTensors .or. CP%WantVectors) then
                dkn1=dkn1  *0.8_dl
                dlnk0=dlnk0/2 !*0.3_dl
--- 528,553 ----
  !     polarization will be calculated. For low values of k we
  !     use a logarithmic spacing. closed case dealt with by SetClosedkValues
           if (CP%WantScalars .and. CP%Reionization .and. CP%AccuratePolarization) then
! !wig
! !            dlnk0=2._dl/10/AccuracyBoost
!             dlnk0=2._dl/10/AccuracyBoost/kSrcSampleBoost
! !end wig
              !Need this to get accurate low l polarization
           else
! !wig
! !            dlnk0=5._dl/10/AccuracyBoost
!             dlnk0=5._dl/10/AccuracyBoost/kSrcSampleBoost
! !end wig
              if (CP%closed) dlnk0=dlnk0/2
           end if
  
           if (CP%AccurateReionization) dlnk0 = dlnk0/2
! !wig
! !         dkn1=0.6_dl/taurst/AccuracyBoost 
! !         dkn2=1.1_dl/taurst/AccuracyBoost 
!          dkn1=0.6_dl/taurst/AccuracyBoost/kSrcSampleBoost
!          dkn2=1.1_dl/taurst/AccuracyBoost/kSrcSampleBoost
! !end wig
            if (CP%WantTensors .or. CP%WantVectors) then
                dkn1=dkn1  *0.8_dl
                dlnk0=dlnk0/2 !*0.3_dl
***************
*** 957,963 ****
        integer lognum
        real(dl)  IntSampleBoost
  
!       IntSampleBoost=AccuracyBoost
  
  !     Fixing the # of k for the integration. 
        
--- 966,974 ----
        integer lognum
        real(dl)  IntSampleBoost
  
! !wig
! !     IntSampleBoost=AccuracyBoost
!       IntSampleBoost=AccuracyBoost/kIntSampleBoost
  
  !     Fixing the # of k for the integration. 
        
diff -r -c -b -B CAMB/inidriver.F90 CAMB_wiggles/inidriver.F90
*** CAMB/inidriver.F90	2005-03-31 00:37:13.000000000 +0100
--- CAMB_wiggles/inidriver.F90	2006-03-27 14:19:45.000000000 +0100
***************
*** 173,179 ****
--- 173,186 ----
  
             P%InitPower%nn = Ini_Read_Int('initial_power_num')
             if (P%InitPower%nn>nnmax) stop 'Too many initial power spectra - increase nnmax in InitialPower'
+ !wig
+ 
+        Power_Name = Ini_Read_String('power_name')
+        
+        if (Power_Name.eq.'power_tilt') then
+ 
             P%InitPower%rat(:) = 1
+           
             do i=1, P%InitPower%nn
                write (numstr,*) i 
                numstr=adjustl(numstr)
***************
*** 177,187 ****
--- 184,197 ----
             do i=1, P%InitPower%nn
                write (numstr,*) i 
                numstr=adjustl(numstr)
+              
+               if (P%WantScalars .or. P%WantTransfer) then
                P%InitPower%an(i) = &
                     Ini_Read_Double('scalar_spectral_index('//trim(numstr)//')')
  
                P%InitPower%n_run(i) = &
                     Ini_Read_Double('scalar_nrun('//trim(numstr)//')',0._dl)
+               end if
      
                if (P%WantTensors) then
                   P%InitPower%ant(i) = Ini_Read_Double('tensor_spectral_index('//trim(numstr)//')')
***************
*** 190,197 ****
--- 200,245 ----
  
                P%InitPower%ScalarPowerAmp(i) = Ini_Read_Double('scalar_amp('//trim(numstr)//')',1.d0) 
                !Always need this as may want to set tensor amplitude even if scalars not computed
+              
+              P%InitPower%logomega(i) = Ini_Read_Double('osc_logomega('//trim(numstr)//')',2.d0)
+              P%InitPower%amposc(i) = Ini_Read_Double('osc_amp('//trim(numstr)//')',0.d0)             
+              P%InitPower%omega(i) = 10**(P%InitPower%logomega(i))
+              P%InitPower%phase(i) = Ini_Read_Double('osc_phase('//trim(numstr)//')',0.d0)
+              P%InitPower%logpowosc(i) = Ini_Read_Double('osc_logpower('//trim(numstr)//')',0.d0)
+              P%InitPower%powosc(i) = 10**(P%InitPower%logpowosc(i))
             end do
        
+           P%InitPower%k_0_osc = Ini_Read_Double('osc_pivot',1d-2)          
+ 
+        elseif (Power_Name.eq.'power_sr') then
+           
+           do i=1, P%InitPower%nn
+              write (numstr,*) i 
+              numstr=adjustl(numstr)
+              
+              if (P%WantScalars .or. P%WantTransfer) then
+                 P%InitPower%eps1(i)= Ini_Read_Double('eps1('//trim(numstr)//')',0.01d0)
+                 P%InitPower%eps2(i)=Ini_Read_Double('eps2('//trim(numstr)//')',0.01d0)
+ 		P%InitPower%eps3(i)=Ini_Read_Double('eps3('//trim(numstr)//')',0.00d0)
+                 P%InitPower%logeps1invsigma0(i) &
+ 		= Ini_Read_Double('logeps1invsigma0('//trim(numstr)//')',1d0)
+                 P%InitPower%loginvsigma0(i) = P%InitPower%logeps1invsigma0(i) &
+                      - log10(P%InitPower%eps1(i))
+                 P%InitPower%invsigma0(i) = 10**(P%InitPower%loginvsigma0(i))
+                 P%InitPower%Xsigma0(i) = Ini_Read_Double('Xsigma0('//trim(numstr)//')',0d0)
+                 P%InitPower%varphi(i) = Ini_Read_Double('varphi('//trim(numstr)//')',0.d0)
+                 P%InitPower%k_star= Ini_Read_Double('k_star',0.05d0)
+                 P%InitPower%P_scalar= 1e-10*exp(Ini_Read_Double('ln[10^{10} P_scalar]',1d0))
+              endif
+              
+           enddo
+           
+        else
+           stop 'No such initial power spectrum'
+        endif
+ 
+ !end wig
+ 
              if (P%WantScalars .or. P%WantTransfer) then
              P%Scalar_initial_condition = Ini_Read_Int('initial_condition',initial_adiabatic)
              if (P%Scalar_initial_condition == initial_vector) then
***************
*** 252,257 ****
--- 300,318 ----
           call Ini_SaveReadValues(trim(outroot) //'params.ini',1)
         end if
  
+ !wig
+        kSrcSampleBoost= Ini_Read_Double('k_src_sample_boost',1.d0)
+        kIntSampleBoost= Ini_Read_Double('k_int_sample_boost',1.d0)
+        DoClsInterpol  = Ini_Read_Logical('do_cls_interpol',.true.)
+        ComputeEveryl  = Ini_Read_Logical('compute_every_l',.true.)
+ 
+        if (FeedbackLevel > 0) then
+           write (*,*)'do_cls_interpol is',DoClsInterpol
+           write (*,*)'compute_every_l is',ComputeEveryl
+           write (*,*)'primordial power spectra is: ',Power_Name
+        endif
+ !end wig
+ 
         call Ini_Close
  
         if (.not. CAMB_ValidateParams(P)) stop 'Stopped due to parameter error'
***************
*** 259,264 ****
--- 320,356 ----
         if (FeedbackLevel > 0) then
           Age = CAMB_GetAge(P) 
           write (*,'("Age of universe/GYr  = ",f7.3)') Age  
+ !wig
+          if (Power_Name.eq.'power_tilt') then
+ 
+              do i=1,P%InitPower%nn
+                 write (*,'("Scalar tilt          = ",f12.8)') P%InitPower%an(i)
+                 write (*,'("Scalar running       = ",f12.8)') P%InitPower%n_run(i)
+                 if (P%WantTensors) then
+                    write (*,'("primordial T/S       = ",f12.8)') P%InitPower%rat(i)
+                    write (*,'("tensor tilt          = ",f12.8)') P%InitPower%ant(i)
+                 endif                
+                 write (*,'("Wiggles amplitude    = ",f12.8)') P%InitPower%amposc(i)
+                 write (*,'("Wiggles frequency    = ",f12.8)') P%InitPower%logomega(i)
+                 write (*,'("Wiggles phase        = ",f12.8)') P%InitPower%phase(i)
+                 write (*,'("Wiggles power        = ",f12.8)') P%InitPower%powosc(i)
+              enddo
+              write (*,'("Wiggles pivot        = ",f12.8)') P%InitPower%k_0_osc
+ 
+           elseif (Power_Name.eq.'power_sr') then
+ 
+              do i=1,P%InitPower%nn
+                 write (*,'("Slow-roll epsilon1   = ",f12.8)') P%InitPower%eps1(i)
+                 write (*,'("Slow-roll epsilon2   = ",f12.8)') P%InitPower%eps2(i)
+ 		write (*,'("Slow-roll epsilon3   = ",f12.8)') P%InitPower%eps3(i)
+                 write (*,'("TPL amplitude        = ",f12.8)') P%InitPower%Xsigma0(i)
+                 write (*,'("TPL scale            = ",f12.8)') P%InitPower%loginvsigma0(i)
+                 write (*,'("TPL phase            = ",f12.8)') P%InitPower%varphi(i)
+              enddo
+           else
+              stop 'Power_Name unknown'
+           endif
+ !end wig
         end if 
  
         call CAMB_GetResults(P)
diff -r -c -b -B CAMB/Makefile CAMB_wiggles/Makefile
*** CAMB/Makefile	2006-03-14 11:07:24.000000000 +0000
--- CAMB_wiggles/Makefile	2006-03-26 19:13:30.000000000 +0100
***************
*** 1,57 ****
! #CAMB Makefile
! #Edit for your compiler
! #Note there are many ifc versions, some of which behave oddly
! 
! #Intel 8 , -openmp toggles mutli-processor:
! F90C     = ifort
! FFLAGS = -openmp -ip -O2 -Vaxlib -W0 -WB -fpp2
! 
! # Intel 9 on IA-64 (eg. COSMOS)
! # (do "module load icomp90" before compiling)
! #F90C = ifort
! #FFLAGS = -openmp -fpp2 -w -O3 -ip -mP2OPT_hlo_prefetch=F
! 
! #Intel ifc, add -openmp for multi-processor (some have bugs):
! #F90C     = ifc
! #FFLAGS = -O2 -Vaxlib -ip -W0 -WB -quiet -fpp2
! #some systems can can also add e.g. -tpp7 -xW
! 
! #SGI, -mp toggles multi-processor. Use -O2 if -Ofast gives problems.
! #F90C     = f90
! #FFLAGS  = -Ofast -mp
! 
! #Digital/Compaq fortran, -omp toggles multi-processor
! #F90C    = f90
! #FFLAGS  = -omp -O4 -arch host -math_library fast -tune host -fpe1
! 
! #Absoft ProFortran, single processor:
! #F90C     = f95
! #FFLAGS = -O2 -cpu:athlon -s -lU77 -w -YEXT_NAMES="LCS" -YEXT_SFX="_"
! 
! #NAGF95, single processor:
! #F90C     = f95
! #FFLAGS = -DNAGF95 -O3
! 
! #G95 compiler
! #F90C   = g95
! #FFLAGS = -O2
! 
! #Sun, single processor:
! #F90C     = f90
! #FFLAGS = -O5
! 
! #Sun V880
! #F90C = mpf90
! #FFLAGS =  -O4 -openmp -ftrap=%none -dalign -DMPI
! 
! #Sun parallel enterprise:
! #F90C     = f95
! #FFLAGS =  -O2 -xarch=native64 -openmp -ftrap=%none
! #try removing -openmp if get bus errors. -03, -04 etc are dodgy.
! 
! #IBM XL Fortran, multi-processor (run gmake)
! #F90C     = xlf90_r
! #FFLAGS  = -qsmp=omp -qsuffix=f=f90:cpp=F90 -O3 -qstrict -qarch=pwr3 -qtune=pwr3
  
  
  #Files containing evolution equations initial power spectrum module
--- 1,33 ----
! # >>> DESIGNED FOR GMAKE <<<
! 
! #CAMB System unified Makefile
! 
! ext=$(shell uname | cut -c1-3)
! 
! ifeq ($(ext),IRI)
! F90C=f90
! FFLAGS = -Ofast -mp -n32 -LANG:recursive=ON
! endif
! 
! ifeq ($(ext),Lin)
! F90C=ifort
! FFLAGS= -O  -openmp -fpp
! endif 
! 
! ifeq ($(ext),OSF)
! F90C=f90
! FFLAGS= -omp -O -arch host -math_library fast -tune host -fpe1
! endif
! 
! ifeq ($(ext),Sun)
! F90C=f90
! FFLAGS= -O4 -xarch=native64 -openmp -ftrap=%none
! endif
! 
! ifeq ($(ext),AIX)
! F90C=xlf90_r
! FFLAGS= -O4 -qsmp=omp -qmaxmem=-1 -qstrict -qfree=f90 -qsuffix=f=f90:cpp=F90
! endif
  
  
  #Files containing evolution equations initial power spectrum module
***************
*** 83,99 ****
  CAMBOBJ       = subroutines.o inifile.o $(POWERSPECTRUM).o recfast.o modules.o \
  	bessels.o $(EQUATIONS).o $(NONLINEAR).o lensing.o cmbmain.o camb.o
  
! default: camb
  
! all: camb $(CAMBLIB)
  
! camb: $(CAMBOBJ) $(DRIVER)
  	$(F90C) $(F90FLAGS) $(CAMBOBJ) $(DRIVER) -o $@
  
  $(CAMBLIB): $(CAMBOBJ)
  	ar -r $@ $?
  
! camb_fits: writefits.f90 $(CAMBOBJ) $(DRIVER)
  	$(F90C) $(F90FLAGS) -I$(HEALPIXDIR)/include $(CAMBOBJ) writefits.f90 $(DRIVER) $(HEALPIXLD) -DWRITE_FITS -o $@
  
  %.o: %.f90
--- 59,75 ----
  CAMBOBJ       = subroutines.o inifile.o $(POWERSPECTRUM).o recfast.o modules.o \
  	bessels.o $(EQUATIONS).o $(NONLINEAR).o lensing.o cmbmain.o camb.o
  
! default: camb.$(ext)
  
! all: camb.$(ext) $(CAMBLIB)
  
! camb.$(ext): $(CAMBOBJ) $(DRIVER)
  	$(F90C) $(F90FLAGS) $(CAMBOBJ) $(DRIVER) -o $@
  
  $(CAMBLIB): $(CAMBOBJ)
  	ar -r $@ $?
  
! camb_fits.$(ext): writefits.f90 $(CAMBOBJ) $(DRIVER)
  	$(F90C) $(F90FLAGS) -I$(HEALPIXDIR)/include $(CAMBOBJ) writefits.f90 $(DRIVER) $(HEALPIXLD) -DWRITE_FITS -o $@
  
  %.o: %.f90
diff -r -c -b -B CAMB/modules.f90 CAMB_wiggles/modules.f90
*** CAMB/modules.f90	2006-03-16 11:17:06.000000000 +0000
--- CAMB_wiggles/modules.f90	2006-03-27 15:50:09.000000000 +0100
***************
*** 33,40 ****
          character(LEN=*), parameter :: version = 'Mar_06'
          
          integer :: FeedbackLevel = 0 !if >0 print out useful information about the model
! 
!         logical, parameter :: DebugMsgs=.false. !Set to true to view progress and timing
  
          logical, parameter :: DebugEvolution = .false. !Set to true to do all the evolution for all k
  
--- 33,42 ----
          character(LEN=*), parameter :: version = 'Mar_06'
          
          integer :: FeedbackLevel = 0 !if >0 print out useful information about the model
! !wig
! !        logical, parameter :: DebugMsgs=.false. !Set to true to view progress and timing
!         logical, parameter :: DebugMsgs=.true.
! !end wig
  
          logical, parameter :: DebugEvolution = .false. !Set to true to do all the evolution for all k
  
***************
*** 163,176 ****
        real(dl) :: lSampleBoost=1._dl
            !Increase lSampleBoost to increase sampling in lSamp%l for Cl interpolation
            
!       real(dl) :: AccuracyBoost =1._dl  
  
            !Decrease step sizes, etc. by this parameter. Useful for checking accuracy.
            !Can also be used to improve speed significantly if less accuracy is required.              
            !or improving accuracy for extreme models. 
            !Note this does not increase lSamp%l sampling or massive neutrino q-sampling
  
!       real(sp) :: lAccuracyBoost=1. 
            !Boost number of multipoles integrated in Boltzman heirarchy
  
        integer, parameter :: lmin = 2
--- 165,184 ----
        real(dl) :: lSampleBoost=1._dl
            !Increase lSampleBoost to increase sampling in lSamp%l for Cl interpolation
            
! !wig
! !      real(dl) :: AccuracyBoost =1._dl  
!       real(dl) :: AccuracyBoost =8._dl
! !end wig
  
            !Decrease step sizes, etc. by this parameter. Useful for checking accuracy.
            !Can also be used to improve speed significantly if less accuracy is required.              
            !or improving accuracy for extreme models. 
            !Note this does not increase lSamp%l sampling or massive neutrino q-sampling
  
! !wig
! !      real(sp) :: lAccuracyBoost=1. 
!       real(sp) :: lAccuracyBoost=2.
! !end wig
            !Boost number of multipoles integrated in Boltzman heirarchy
  
        integer, parameter :: lmin = 2
***************
*** 186,192 ****
        integer, parameter:: l0max=4000
  
  !     lmax is max possible number of l's evaluated
!       integer, parameter :: lmax_arr = 100+l0max/7
   
          contains
        
--- 194,207 ----
        integer, parameter:: l0max=4000
  
  !     lmax is max possible number of l's evaluated
! !wig
! !      integer, parameter :: lmax_arr = 100+l0max/7
!       integer, parameter :: lmax_arr = 100+l0max
!       real(dl) :: kSrcSampleBoost=1._dl
!       real(dl) :: kIntSampleBoost=1._dl
!       logical :: DoClsInterpol=.true.
!       logical :: ComputeEveryl=.true.
! !end wig
   
          contains
        
***************
*** 539,544 ****
--- 554,578 ----
          integer lind, lvar, step,top,bot,ls(lmax_arr)
          real(dl) AScale
      
+ !wig        
+         if (ComputeEveryl) then
+            write(*,*) 'Warning: computing every l'           
+            if (max_l.gt.lmax_arr) then
+               write(*,*)'lmax_arr maxl',lmax_arr,max_l
+               stop
+            endif
+            lind=0
+            do lvar=2, max_l
+               lind=lind+1
+               ls(lind)=lvar 
+            end do
+ 
+            lSet%l0=lind
+            lSet%l = ls(1:lind)
+ 
+         else
+ !end wig
+ 
          Ascale=scale/lSampleBoost       
        
          lind=0
***************
*** 640,645 ****
--- 674,682 ----
          end if
          lSet%l0=lind
          lSet%l(1:lind) = ls(1:lind)
+ !wig        
+         endif
+ !end wig
          
        end subroutine initlval
  
***************
*** 659,665 ****
  
        xl = real(lSet%l(1:lSet%l0),dl)
        call spline(xl,iCl(1),max_ind,cllo,clhi,ddCl(1))
!      
              llo=1
              do il=lmin,lSet%l(max_ind)
                 xi=il
--- 696,703 ----
  
        xl = real(lSet%l(1:lSet%l0),dl)
        call spline(xl,iCl(1),max_ind,cllo,clhi,ddCl(1))
! !wig
!       if (DoClsInterpol) then
              llo=1
              do il=lmin,lSet%l(max_ind)
                 xi=il
***************
*** 675,680 ****
--- 713,732 ----
                         +(b0**3-b0)*ddCl(lhi))*ho**2/6
                
              end do
+       else
+         llo=1
+          do il=2,lSet%l(max_ind)
+             xi=il
+             if ((xi > lSet%l(llo+1)).and.(llo < max_ind)) then
+                llo=llo+1
+             end if
+ 
+             all_Cl(il) = iCl(llo+1)
+ 
+          end do
+          all_Cl(2)=iCl(1)
+       endif
+ !end wig 
  
        end subroutine InterpolateClArr
  
***************
*** 1872,1878 ****
            Type(MatterTransferData), intent(in) :: MTrans
            character(LEN=Ini_max_string_len), intent(IN) :: FileNames(*)
            integer itf,in,i
!           integer, parameter :: points = 500
            real outpower(points,CP%InitPower%nn)
            character(LEN=80) fmt
            real minkh,dlnkh
--- 1924,1933 ----
            Type(MatterTransferData), intent(in) :: MTrans
            character(LEN=Ini_max_string_len), intent(IN) :: FileNames(*)
            integer itf,in,i
! !wig
! !          integer, parameter :: points = 500
!           integer, parameter :: points = 10000
! !end wig
            real outpower(points,CP%InitPower%nn)
            character(LEN=80) fmt
            real minkh,dlnkh
diff -r -c -b -B CAMB/params.ini CAMB_wiggles/params.ini
*** CAMB/params.ini	2005-03-30 16:05:28.000000000 +0100
--- CAMB_wiggles/params.ini	2006-03-27 14:17:20.000000000 +0100
***************
*** 6,12 ****
  #What to do
  get_scalar_cls = T
  get_vector_cls = F
! get_tensor_cls = F
  get_transfer   = F
  
  #if do_lensing then scalar_output_file contains additional columns of l^4 C_l^{pp} and l^3 C_l^{pT}
--- 6,12 ----
  #What to do
  get_scalar_cls = T
  get_vector_cls = F
! get_tensor_cls = T
  get_transfer   = F
  
  #if do_lensing then scalar_output_file contains additional columns of l^4 C_l^{pp} and l^3 C_l^{pT}
***************
*** 64,69 ****
--- 64,77 ----
  
  #Initial power spectrum, amplitude, spectral index and running
  initial_power_num         = 1
+ 
+ #Type of initial power spectrum: power_tilt or power_sr
+ power_name = power_sr
+ 
+ ######################################################################################
+ 
+ #for power_tilt power spectra
+ 
  scalar_amp(1)             = 2.3e-9
  scalar_spectral_index(1)  = 1
  scalar_nrun(1)            = 0
***************
*** 72,77 ****
--- 80,112 ----
  initial_ratio(1)          = 1
  #note vector modes use the scalar settings above
  
+ #primordial power_tilted wiggle
+ osc_logomega(1)           = 1.
+ osc_amp(1)                = 0.1
+ osc_phase(1)              = 0.
+ osc_logpower(1)           = 0.
+ osc_pivot                 = 0.05
+ 
+ #####################################################################################
+ 
+ #for power_sr power spectra
+ 
+ #slow-roll parameters
+ ln[10**10 P_scalar] = 1.
+ eps1(1) = 0.104396921E-01
+ eps2(1) = -0.226308387E-01
+ eps3(1) = 0.
+ 
+ #wiggles paramters
+ logeps1invsigma0(1) = 1.4921786
+ Xsigma0(1) = 0.0473405533
+ varphi(1) = 0.
+ 
+ #tpl pivot
+ k_star = 0.05
+ 
+ #####################################################################################
+ 
  #Initial scalar perturbation mode (adiabatic=1, CDM iso=2, Baryon iso=3, 
  # neutrino density iso =4, neutrino velocity iso = 5) 
  initial_condition   = 1
***************
*** 102,112 ****
  
  
  #Output files not produced if blank. make camb_fits to use use the FITS setting.
! scalar_output_file = scalCls.dat
! vector_output_file = vecCls.dat
! tensor_output_file = tensCls.dat
! total_output_file  = totCls.dat
! lensed_output_file = lensedCls.dat
  FITS_filename      = scalCls.fits
  
  ##Optional parameters to control the computation speed,accuracy and feedback
--- 137,147 ----
  
  
  #Output files not produced if blank. make camb_fits to use use the FITS setting.
! scalar_output_file = scalar_cls.dat
! vector_output_file = vector_cls.dat
! tensor_output_file = tensor_cls.dat
! total_output_file  = total_cls.dat
! lensed_output_file = lensed_cls.dat
  FITS_filename      = scalCls.fits
  
  ##Optional parameters to control the computation speed,accuracy and feedback
***************
*** 147,159 ****
  
  #Increase accuracy_boost to decrease time steps, use more k values,  etc.
  #Decrease to speed up at cost of worse accuracy. Suggest 0.8 to 3.
! accuracy_boost          = 1
  
  #Larger to keep more terms in the hierarchy evolution. 
! l_accuracy_boost        = 1
  
  #Increase to use more C_l values for interpolation.
  #Increasing a bit will improve the polarization accuracy at l up to 200 -
  #interpolation errors may be up to 3%
  #Decrease to speed up non-flat models a bit
  l_sample_boost          = 1
--- 182,206 ----
  
  #Increase accuracy_boost to decrease time steps, use more k values,  etc.
  #Decrease to speed up at cost of worse accuracy. Suggest 0.8 to 3.
! accuracy_boost          = 6
  
  #Larger to keep more terms in the hierarchy evolution. 
! l_accuracy_boost        = 2
  
  #Increase to use more C_l values for interpolation.
  #Increasing a bit will improve the polarization accuracy at l up to 200 -
  #interpolation errors may be up to 3%
  #Decrease to speed up non-flat models a bit
  l_sample_boost          = 1
+ 
+ # Boost the sampling on k sources only, in unit of AccuracyBoost
+ k_src_sample_boost      = 1.
+ 
+ # Boost the sampling on k line of sight integral only, in unit of AccuracyBoost
+ k_int_sample_boost      = 1.
+ 
+ #Use spline interpolation for the Cls
+ do_cls_interpol         = F
+  
+ # Remove undersampling of the Cls: if true lSampleBoost has no effect
+ compute_every_l         = T
diff -r -c -b -B CAMB/power_tilt.f90 CAMB_wiggles/power_tilt.f90
*** CAMB/power_tilt.f90	2005-03-31 00:37:13.000000000 +0100
--- CAMB_wiggles/power_tilt.f90	2006-03-27 14:41:10.000000000 +0100
***************
*** 33,40 ****
  
       private
      
!       character(LEN=*), parameter :: Power_Name = 'power_tilt'
! 
        integer, parameter :: nnmax= 5 
         !Maximum possible number of different power spectra to use
  
--- 33,43 ----
  
       private
      
! !wig      character(LEN=*), parameter :: Power_Name = 'power_tilt'
!       character(80) :: Power_Name = 'power_sr'
!       logical, parameter :: useTplPhase = .false.      
!       logical :: useSecondOrder = .false.
! !end wig
        integer, parameter :: nnmax= 5 
         !Maximum possible number of different power spectra to use
  
***************
*** 51,56 ****
--- 54,106 ----
         real(dl) k_0_scalar, k_0_tensor
         real(dl) ScalarPowerAmp(nnmax)
   
+ !wig
+ ! slow-roll power
+        real(dl) P_scalar
+        real(dl) k_star
+        real(dl) eps1(nnmax)
+        real(dl) eps2(nnmax)
+        real(dl) eps3(nnmax)
+ ! tpl power
+        real(dl) invsigma0(nnmax)
+        real(dl) loginvsigma0(nnmax)
+        real(dl) logeps1invsigma0(nnmax)
+        real(dl) Xsigma0(nnmax)
+        real(dl) varphi(nnmax) 
+ ! osc power
+        real(dl) amposc(nnmax)
+        real(dl) omega(nnmax)
+        real(dl) logomega(nnmax)
+        real(dl) phase(nnmax)
+        real(dl) powosc(nnmax)
+        real(dl) logpowosc(nnmax)
+        real(dl) k_0_osc
+ !Internal parameters              
+        real(dl) ks_pivot          !scalar pivot wave number
+        real(dl) kt_pivot          !tensor pivot wave number
+        real(dl) ko_pivot          !wiggles pivot wave number
+        real(dl) P_0_scalar(nnmax) !scalar normalisation during inflation
+        real(dl) P_0_tensor(nnmax) !tensor normalisation during inflation
+        real(dl) as0(nnmax)        !scalar power spectrum coefficients: order 0
+        real(dl) as1(nnmax)        !order 1
+        real(dl) as2(nnmax)        !order 2
+        real(dl) amps0(nnmax)      !amplitudes of tpl effect in scalar power: order 0
+        real(dl) amps1(nnmax)      !order 1
+        real(dl) amps2(nnmax)      !order 2
+        real(dl) args0(nnmax)      !argument of the modulations: order 0
+        real(dl) args1(nnmax)      !order 1
+        real(dl) args2(nnmax)      !order 2
+        real(dl) at0(nnmax)        !Tensor power spectrum coefficients: same as scalars
+        real(dl) at1(nnmax)
+        real(dl) at2(nnmax)
+        real(dl) ampt0(nnmax) 
+        real(dl) ampt1(nnmax)
+        real(dl) ampt2(nnmax)
+        real(dl) argt0(nnmax)
+        real(dl) argt1(nnmax) 
+        real(dl) argt2(nnmax)
+ !end wig
+         
        end Type InitialPowerParams
  
        real(dl) curv  !Curvature contant, set in InitializePowers     
***************
*** 70,75 ****
--- 120,127 ----
          Type (InitialPowerParams) :: AP
  
           AP%nn     = 1 !number of initial power spectra
+ !wig
+            
           AP%an     = 1 !scalar spectral index
           AP%n_run   = 0 !running of scalar spectral index
           AP%ant    = 0 !tensor spectra index
***************
*** 78,83 ****
--- 130,157 ----
           AP%k_0_tensor = 0.05
           AP%ScalarPowerAmp = 1
  
+          AP%amposc = 0.
+          AP%logomega = 2.
+          AP%omega = 10**(AP%logomega)
+          AP%phase = 0.
+          AP%logpowosc = 0.
+          AP%powosc = 10**(AP%logpowosc)
+          AP%k_0_osc = 0.05
+          
+          AP%eps1 = 0.01
+          AP%eps2 = 0.01 
+          AP%eps3 = 0.
+          AP%P_scalar=1.d-9
+          AP%k_star = 0.05            
+ 
+          AP%Xsigma0 = 0.                
+          AP%logeps1invsigma0 = 0. 
+          AP%loginvsigma0 = AP%logeps1invsigma0 - log10(AP%eps1) 
+          AP%invsigma0=10**(AP%loginvsigma0)
+          AP%varphi=0.
+         
+ !end wig
+ 
         end subroutine SetDefPowerParams
  
         subroutine InitializePowers(AParamSet,acurv)
***************
*** 86,91 ****
--- 160,170 ----
           !Could read spectra from disk here, do other processing, etc.
  
          real(dl) acurv
+ !wig
+         integer :: i
+         real(dl), parameter :: C_const = -0.7296
+         real(dl), parameter :: pi = 3.14159265
+ !end wig
  
          P = AParamSet
          if (P%nn > nnmax) then
***************
*** 96,101 ****
--- 175,319 ----
          curv=acurv         
  
  !Write implementation specific code here...
+ !wig
+         if (Power_Name.eq.'power_tilt') then
+            
+            P%ks_pivot = P%k_0_scalar
+            P%kt_pivot = P%k_0_tensor
+            P%ko_pivot = P%k_0_osc
+ 
+            do i=1, P%nn        
+                                                 
+               P%omega(i) = 10**(P%logomega(i))
+               P%powosc(i) = 10**(P%logpowosc(i))
+ 
+               P%as0(i) = 0.
+               P%as1(i) = P%an(i)-1
+               P%as2(i) = P%n_run(i)/2.
+               P%amps0(i) = P%amposc(i)
+               P%amps1(i) = 0.
+               P%amps2(i) = 0.
+               P%args0(i) = P%phase(i)
+               P%args1(i) = P%omega(i)/P%powosc(i)
+               P%args2(i) = 0.
+               
+               P%at0(i) = 0.
+               P%at1(i) = P%ant(i)
+               P%at2(i) = 0.
+               P%ampt0(i) = P%amps0(i)
+               P%ampt1(i) = 0.
+               P%ampt2(i) = 0.
+               P%argt0(i) = P%phase(i)
+               P%argt1(i) = P%omega(i)
+               P%argt2(i) = 0.
+              
+               P%P_0_scalar(i) = P%ScalarPowerAmp(i)
+               P%P_0_tensor(i) = P%rat(i)*P%ScalarPowerAmp(i)
+ 
+            enddo
+ 
+         elseif (Power_Name.eq.'power_sr') then
+            if(abs(curv) .gt. 1.d-3) then
+               write(*,*)'Non flat cosmology, curv = ',curv
+               write(*,*)'Cannot use SR or TPL spectrum.'
+               stop
+            endif
+ 
+            P%ks_pivot = P%k_star
+            P%kt_pivot = P%k_star
+            P%ko_pivot = P%k_star
+ 
+            do i=1, P%nn        
+ 
+ !TODO: derive second order tpl corrections?
+               if (useSecondOrder) then
+                  if (P%Xsigma0(i).ne.0.) then
+                     useSecondOrder=.false.
+                     write(*,*)'PowerIni: SR second order disabled!'
+                  endif
+               endif
+ 
+ !first order tpl corrections
+               
+               P%loginvsigma0(i) = P%logeps1invsigma0(i) - log10(P%eps1(i))
+               P%invsigma0(i)=10**(P%loginvsigma0(i))
+               
+               P%amps0(i) = 2.*P%Xsigma0(i)
+               P%amps1(i) = pi*P%Xsigma0(i)*(2.*P%eps1(i) + P%eps2(i))
+               P%amps2(i) = 0.
+               if (useTplPhase) then
+                  P%args0(i) = 2.*P%invsigma0(i)*(1. + P%eps1(i)) + P%varphi(i)
+               else
+                  P%args0(i) = P%varphi(i)
+               endif
+               P%args1(i) = 2.*P%invsigma0(i)*P%eps1(i)
+               P%args2(i) = 0.
+ 
+               P%ampt0(i) = P%amps0(i)
+               P%ampt1(i) = pi*P%Xsigma0(i)*2.*P%eps1(i)
+               P%ampt2(i) = 0.
+               P%argt0(i) = P%args0(i)
+               P%argt1(i) = P%args1(i)
+               P%argt2(i) = 0.
+ 
+ 
+ !first and second order slowroll
+ 
+               P%as0(i) = 1. + P%eps1(i)*(-2.*(C_const+1.))-P%eps2(i)*C_const
+               P%as1(i) = -2.*P%eps1(i) - P%eps2(i)
+               P%as2(i) = 0.
+                             
+               P%at0(i) = 1. + P%eps1(i)*(-2.*(C_const+1.))
+               P%at1(i) = -2.*P%eps1(i)
+               P%at2(i) = 0.
+ 
+               if (useSecondOrder) then
+                                   
+                  P%as0(i) = P%as0(i) &
+                       + P%eps1(i)*P%eps1(i) * (2.*C_const**2 + 2.*C_const + pi**2/2. - 5.) &
+                       + P%eps1(i)*P%eps2(i) * (C_const**2 - C_const + 7.*pi**2/12. - 7.) &
+                       + P%eps2(i)*P%eps2(i) * (0.5*C_const**2 + pi**2/8. - 1.) &
+                       + P%eps2(i)*P%eps3(i) * (-0.5*C_const**2 + pi**2/24.)
+ 
+                  P%as1(i) = P%as1(i) &
+                       + P%eps1(i)*P%eps1(i) * 2.*(2.*C_const + 1.) &
+                       + P%eps1(i)*P%eps2(i) * (2.*C_const - 1.) &
+                       + P%eps2(i)*P%eps2(i) * C_const &
+                       + P%eps2(i)*P%eps3(i) * (-C_const)
+ 
+                  P%as2(i) = P%as2(i) &
+                       + P%eps1(i)*P%eps1(i) * 4. &
+                       + P%eps1(i)*P%eps2(i) * 2. &
+                       + P%eps2(i)*P%eps2(i) * 1. &
+                       + P%eps2(i)*P%eps3(i) * (-1.)
+ 
+                  
+                  P%at0(i) = P%at0(i) &
+                       +  P%eps1(i)*P%eps1(i) * (2.*C_const**2 + 2.*C_const + pi**2/2. - 5.) &
+                       +  P%eps1(i)*P%eps2(i) * (-C_const**2 - 2.*C_const + pi**2/12. - 2.)
+                  
+                  P%at1(i) = P%at1(i) &
+                       + P%eps1(i)*P%eps1(i) * 2.*(2.*C_const + 1.) &
+                       + P%eps1(i)*P%eps2(i) * 2.*(-C_const - 1.)
+ 
+                  P%at2(i) = P%at2(i)  &
+                       + P%eps1(i)*P%eps1(i) * 4. &
+                       + P%eps1(i)*P%eps2(i) * (-2.)
+                       
+               endif
+ 
+               
+               !Relative normalisation of scalar and tensor spectra
+                  
+               P%P_0_scalar(i)=P%P_scalar/(P%as0(i))
+               P%P_0_tensor(i)=P%P_scalar/(P%as0(i))*P%eps1(i)*16.
+ 
+            enddo
+ 
+         else
+            stop 'No such initial power spectra'
+         endif
+ !end wig
  
         end subroutine InitializePowers
         
***************
*** 120,127 ****
          real(dl) ScalarPower,k, lnrat
          integer in
  
!           lnrat = log(k/P%k_0_scalar)
!           ScalarPower=P%ScalarPowerAmp(in)*exp((P%an(in)-1)*lnrat + P%n_run(in)/2*lnrat**2)   
   
        end function ScalarPower
  
--- 338,377 ----
          real(dl) ScalarPower,k, lnrat
          integer in
  
! !wig
!         real :: angle
!         
! !          lnrat = log(k/P%k_0_scalar)
! !          ScalarPower=P%ScalarPowerAmp(in)*exp((P%an(in)-1)*lnrat + P%n_run(in)/2*lnrat**2)   
!                  
!         lnrat = log(k/P%ks_pivot)
! 
!         if (Power_Name.eq.'power_tilt') then
! 
!            angle = P%args0(in) + P%args1(in) * ((k/P%ko_pivot)**P%powosc(in)-1.)
! 
!            ScalarPower = exp(P%as1(in)*lnrat + P%as2(in)*lnrat**2) &
!                 * (1. - P%amps0(in)*cos(angle))
! 
!         elseif (Power_Name.eq.'power_sr') then
! 
!            angle = P%args0(in) + P%args1(in) * log(k/P%ko_pivot) 
! 
!            ScalarPower = (P%as0(in) + P%as1(in)*lnrat + 0.5*P%as2(in)*lnrat**2) &
!                 * (1. - P%amps0(in)*cos(angle)) - P%amps1(in)*sin(angle)
! 
!         else
!            stop 'No such initial power spectra'
!         endif
! 
!         ScalarPower = ScalarPower*P%P_0_scalar(in)
! 
!         if(ScalarPower .lt. 0.0) then
!            ScalarPower =0.
!         end if
! 
! 
! !end wig
   
        end function ScalarPower
  
***************
*** 142,151 ****
          real(dl), parameter :: PiByTwo=3.14159265d0/2._dl
     
          integer in
  
-         TensorPower=P%rat(in)*P%ScalarPowerAmp(in)*exp(P%ant(in)*log(k/P%k_0_tensor))
          if (curv < 0) TensorPower=TensorPower*tanh(PiByTwo*sqrt(-k**2/curv-3)) 
  
        end function TensorPower
  
        !Get parameters describing parameterisation (for FITS file)
--- 392,433 ----
          real(dl), parameter :: PiByTwo=3.14159265d0/2._dl
     
          integer in
+ !wig
+         real :: lnrat,angle
+ 
+ !        TensorPower=P%rat(in)*P%ScalarPowerAmp(in)*exp(P%ant(in)*log(k/P%k_0_tensor))
+ !        if (curv < 0) TensorPower=TensorPower*tanh(PiByTwo*sqrt(-k**2/curv-3)) 
+         
+         lnrat = log(k/P%kt_pivot)
+ 
+         if (Power_Name.eq.'power_tilt') then
+ 
+            angle = P%argt0(in) + P%argt1(in) * ((k/P%ko_pivot)**P%powosc(in)-1.)
+ 
+            TensorPower = exp(P%at1(in)*lnrat) &
+                 * (1. - P%ampt0(in)*cos(angle))
  
          if (curv < 0) TensorPower=TensorPower*tanh(PiByTwo*sqrt(-k**2/curv-3)) 
  
+         elseif (Power_Name.eq.'power_sr') then
+ 
+            angle = P%argt0(in) + P%argt1(in) * log(k/P%ko_pivot)
+ 
+            TensorPower = (P%at0(in) + P%at1(in)*lnrat + 0.5*P%at2(in)*lnrat**2) &
+                 * (1. - P%ampt0(in)*cos(angle)) - P%ampt1(in)*sin(angle)
+ 
+         else
+            stop 'No such initial power spectra'
+         endif
+ 
+         TensorPower = TensorPower*P%P_0_tensor(in)
+         
+         if(TensorPower .lt. 0.0) then
+            TensorPower=0.
+         endif
+         
+ !end wig
+ 
        end function TensorPower
  
        !Get parameters describing parameterisation (for FITS file)
***************
*** 156,161 ****
--- 438,446 ----
           logical, intent(IN) :: Scal, Tens
           integer num, Power_Descript
           num=0
+ !wig
+          if (Power_Name.eq.'power_tilt') then
+ 
           if (Scal) then
           num=num+1
           Keys(num) = 'n_s'
***************
*** 180,185 ****
--- 465,520 ----
             Vals(num) = P%rat(in)
           end if
           end if
+             num=num+1
+             Keys(num) = 'amposc'
+             Vals(num) = P%amposc(in)
+             num=num+1
+             Keys(num) = 'logomega'
+             Vals(num) = P%logomega(in)
+             num=num+1
+             Keys(num) = 'phase'
+             Vals(num) = P%phase(in)
+             num=num+1
+             Keys(num) = 'o_logpow'
+             Vals(num) = P%logpowosc(in)
+             num=num+1
+             Keys(num) = 'o_pivot'
+             Vals(num) = P%k_0_osc
+ 
+          elseif (Power_Name.eq.'power_sr') then
+ 
+             num=num+1
+             Keys(num) = 'k_star'
+             Vals(num) = P%k_star
+             num=num+1 
+             Keys(num) = 'P_scalar'
+             Vals(num) = P%P_scalar     
+             num=num+1 
+             Keys(num) = 'eps1'
+             Vals(num) = P%eps1(in)
+             num=num+1 
+             Keys(num) = 'eps2'
+             Vals(num) = P%eps2(in)
+             num=num+1 
+             Keys(num) = 'eps3'
+             Vals(num) = P%eps2(in)
+             num=num+1 
+             Keys(num) = 'isigma0'
+             Vals(num) = P%invsigma0(in)
+             num=num+1
+             Keys(num) = 'lisigma0'
+             Vals(num) = P%loginvsigma0(in)
+             num=num+1 
+             Keys(num) = 'Xsigma0'
+             Vals(num) = P%Xsigma0(in)
+             num=num+1 
+             Keys(num) = 'varphi'
+             Vals(num) = P%varphi(in)
+ 
+          else
+             stop 'No such initial power spectra'
+          endif
+ !end wig
           Power_Descript = num
  
         end  function Power_Descript
