diff -r -c -b -B cosmomc/camb/cmbmain.f90 cosmomc_wiggles/camb/cmbmain.f90
*** cosmomc/camb/cmbmain.f90	2004-11-29 15:07:14.000000000 +0000
--- cosmomc_wiggles/camb/cmbmain.f90	2004-12-15 13:40:01.000000000 +0000
***************
*** 521,537 ****
  !     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
--- 521,546 ----
  !     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
***************
*** 868,874 ****
        integer lognum
        real(dl)  IntSampleBoost
  
!       IntSampleBoost=AccuracyBoost
  
  !     Fixing the # of k for the integration. 
        
--- 877,885 ----
        integer lognum
        real(dl)  IntSampleBoost
  
! !wig
! !     IntSampleBoost=AccuracyBoost
!       IntSampleBoost=AccuracyBoost/kIntSampleBoost
  
  !     Fixing the # of k for the integration. 
        
diff -r -c -b -B cosmomc/camb/inidriver.F90 cosmomc_wiggles/camb/inidriver.F90
*** cosmomc/camb/inidriver.F90	2004-11-30 13:42:34.000000000 +0000
--- cosmomc_wiggles/camb/inidriver.F90	2004-12-01 12:40:31.000000000 +0000
***************
*** 151,157 ****
--- 151,164 ----
  
             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)
***************
*** 155,165 ****
--- 162,175 ----
             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)//')')
***************
*** 168,175 ****
--- 178,219 ----
  
                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.02d0)
+                 P%InitPower%eps2(i)=Ini_Read_Double('eps2('//trim(numstr)//')',0.02d0)
+                 P%InitPower%loginvsigma0(i) = Ini_Read_Double('loginvsigma0('//trim(numstr)//')',0.02d0)
+                 P%InitPower%invsigma0(i) = 10**(P%InitPower%loginvsigma0(i))
+                 P%InitPower%Xsigma0(i) = Ini_Read_Double('Xsigma0('//trim(numstr)//')',0.02d0)
+                 P%InitPower%varphi(i) = Ini_Read_Double('varphi('//trim(numstr)//')',0.02d0)
+                 P%InitPower%k_star= Ini_Read_Double('k_star',0.05d0)
+                 P%InitPower%P_scalar= Ini_Read_Double('P_scalar',1.d-9)
+              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
***************
*** 224,229 ****
--- 268,286 ----
         lAccuracyBoost = Ini_Read_Double('l_accuracy_boost',1.d0)
         lSampleBoost   = Ini_Read_Double('l_sample_boost',1.d0)
  
+ !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'
***************
*** 231,236 ****
--- 288,323 ----
         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 (*,'("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 cosmomc/camb/Makefile cosmomc_wiggles/camb/Makefile
*** cosmomc/camb/Makefile	2004-11-30 13:42:31.000000000 +0000
--- cosmomc_wiggles/camb/Makefile	2004-10-29 12:45:33.000000000 +0100
***************
*** 1,50 ****
! #CAMB Makefile
! #Edit for your compiler
! #Note there are many ifc versions, some of which behave oddly
! 
! #Intel 8 , -openmp toggles mutli-processor:
! #(do "module load dev8" to get v8 compilers on COSMOS)
! F90C     = ifort
! FFLAGS = -openmp -ip -O2 -Vaxlib -W0 -WB -fpp2
! 
! #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 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 and initial power spectrum module
  EQUATIONS     = equations
--- 1,28 ----
! # >>> DESIGNED FOR GMAKE <<<
  
+ #CAMB System unified Makefile
+ 
+ ext=$(shell uname | cut -c1-3)
+ 
+ ifeq ($(ext),IRI)
+ F90C=f90
+ FFLAGS = -O3 -mp -64
+ endif
+ 
+ ifeq ($(ext),Lin)
+ F90C=ifort
+ FFLAGS= -O3 -openmp  -w95 -pc80 -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
  
  #Files containing evolution equations and initial power spectrum module
  EQUATIONS     = equations
***************
*** 52,59 ****
  
  #Driver program
  DRIVER        = inidriver.F90
- #DRIVER        = sigma8.f90
- #DRIVER        = tester.f90
  
  #Settings for building camb_fits
  #Location of FITSIO and name of library 
--- 30,35 ----
***************
*** 64,88 ****
  
  CAMBLIB       = libcamb.a
  
- #Shouldn't need to change anything else...
- 
  F90FLAGS      = $(FFLAGS)
  HEALPIXLD     = -L$(HEALPIXDIR)/lib -lhealpix -L$(FITSDIR) -l$(FITSLIB)
  FC            = $(F90C)
  
  CAMBOBJ       = subroutines.o inifile.o $(POWERSPECTRUM).o recfast.o modules.o bessels.o $(EQUATIONS).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
--- 41,63 ----
  
  CAMBLIB       = libcamb.a
  
  F90FLAGS      = $(FFLAGS)
  HEALPIXLD     = -L$(HEALPIXDIR)/lib -lhealpix -L$(FITSDIR) -l$(FITSLIB)
  FC            = $(F90C)
  
  CAMBOBJ       = subroutines.o inifile.o $(POWERSPECTRUM).o recfast.o modules.o bessels.o $(EQUATIONS).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 cosmomc/camb/modules.f90 cosmomc_wiggles/camb/modules.f90
*** cosmomc/camb/modules.f90	2004-11-29 15:07:14.000000000 +0000
--- cosmomc_wiggles/camb/modules.f90	2004-12-01 13:20:59.000000000 +0000
***************
*** 33,40 ****
          character(LEN=8), parameter :: version = 'Nov_04'
          
          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=8), parameter :: version = 'Nov_04'
          
          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
  
***************
*** 162,175 ****
        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
         
        real(dl), parameter :: OmegaKFlat = 5e-7_dl !Value at which to use flat code
--- 164,183 ----
        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
         
        real(dl), parameter :: OmegaKFlat = 5e-7_dl !Value at which to use flat code
***************
*** 182,188 ****
        integer, parameter:: l0max=4000
  
  !     lmax is max possible number of l's evaluated
!       integer, parameter :: lmax_arr = 100+l0max/7
   
          contains
        
--- 190,203 ----
        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
        
***************
*** 535,540 ****
--- 550,574 ----
          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
***************
*** 636,641 ****
--- 670,678 ----
          end if
          lSet%l0=lind
          lSet%l(1:lind) = ls(1:lind)
+ !wig        
+         endif
+ !end wig
  
        end subroutine initlval
  
***************
*** 655,661 ****
  
        xl = real(lSet%l(1:lSet%l0),dl)
        call spline(xl,iCl(1),max_ind,cllo,clhi,ddCl(1))
!      
              llo=1
              do il=2,lSet%l(max_ind)
                 xi=il
--- 692,699 ----
  
        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=2,lSet%l(max_ind)
                 xi=il
***************
*** 671,676 ****
--- 709,728 ----
                         +(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
  
***************
*** 1698,1704 ****
            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
--- 1750,1759 ----
            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 cosmomc/camb/params.ini cosmomc_wiggles/camb/params.ini
*** cosmomc/camb/params.ini	2005-01-18 18:09:52.000000000 +0000
--- cosmomc_wiggles/camb/params.ini	2005-02-27 15:51:15.000000000 +0000
***************
*** 3,9 ****
  #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}
--- 3,9 ----
  #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}
***************
*** 14,24 ****
  accurate_BB = F
  
  #Files. 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
  
  #Maximum multipole and k*eta. 
--- 14,24 ----
  accurate_BB = F
  
  #Files. 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
  
  #Maximum multipole and k*eta. 
***************
*** 71,76 ****
--- 71,84 ----
  
  #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
***************
*** 79,84 ****
--- 87,117 ----
  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                 = 5e-2
+ 
+ #####################################################################################
+ 
+ #for power_sr power spectra
+ 
+ #slow-roll parameters
+ P_scalar = 2.36978302E-9
+ eps1(1) = 0.104396921E-01
+ eps2(1) = -0.226308387E-01
+ 
+ #wiggles paramters
+ loginvsigma0(1) = 3.47349095
+ Xsigma0(1) = 0.0473405533
+ varphi(1) = 0.
+ 
+ k_star = 1E-2
+ 
+ #####################################################################################
+ 
  #Initial scalar perturbation mode (adiabatic=1, CDM iso=2, Baryon iso=3, 
  # neutrino density iso =4, neutrino velocity iso = 5) 
  initial_condition   = 1
***************
*** 136,145 ****
  
  #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 -
--- 169,178 ----
  
  #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 -
***************
*** 147,152 ****
--- 180,197 ----
  #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 = T
+ 
+ # Remove undersampling of the Cls: if true lSampleBoost has no effect
+ compute_every_l = T
+ 
  
  
  
diff -r -c -b -B cosmomc/camb/power_tilt.f90 cosmomc_wiggles/camb/power_tilt.f90
*** cosmomc/camb/power_tilt.f90	2004-11-29 15:07:14.000000000 +0000
--- cosmomc_wiggles/camb/power_tilt.f90	2004-11-05 14:07:36.000000000 +0000
***************
*** 33,40 ****
  
       private
      
!       character(LEN=*), parameter :: Power_Name = 'power_tilt'
! 
        integer, parameter :: nnmax= 5 
         !Maximum possible number of different power spectra to use
  
--- 33,42 ----
  
       private
      
! !wig      character(LEN=*), parameter :: Power_Name = 'power_tilt'
!       character(80) :: Power_Name = 'power_sr'
!       logical :: Tpl_Phase = .true.
! !end wig
        integer, parameter :: nnmax= 5 
         !Maximum possible number of different power spectra to use
  
***************
*** 51,56 ****
--- 53,103 ----
         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)
+ ! tpl power
+        real(dl) invsigma0(nnmax)
+        real(dl) loginvsigma0(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 ****
--- 117,124 ----
          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 ****
--- 127,152 ----
           AP%k_0_tensor = 0.05
           AP%ScalarPowerAmp = 1
  
+          AP%amposc = 0d0
+          AP%logomega = 2d0
+          AP%omega = 10**(AP%logomega)
+          AP%phase = 0d0
+          AP%logpowosc = 0d0
+          AP%powosc = 10**(AP%logpowosc)
+          AP%k_0_osc = 0.02
+          
+          AP%eps1 = 0.02
+          AP%eps2 = 0.02 
+          AP%P_scalar=1.d-9
+          AP%k_star=0.01d0            
+ 
+          AP%Xsigma0=0d0                
+          AP%loginvsigma0=2. 
+          AP%invsigma0=10**(AP%loginvsigma0)
+          AP%varphi=0.
+         
+ !end wig
+ 
         end subroutine SetDefPowerParams
  
         subroutine InitializePowers(AParamSet,acurv)
***************
*** 86,91 ****
--- 155,164 ----
           !Could read spectra from disk here, do other processing, etc.
  
          real(dl) acurv
+ !wig
+         integer :: i
+         real(dl), parameter :: C_const = -0.7296
+ !end wig
  
          P = AParamSet
          if (P%nn > nnmax) then
***************
*** 96,101 ****
--- 169,262 ----
          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) = 0d0
+               P%as1(i) = P%an(i)-1
+               P%as2(i) = P%n_run(i)/2d0
+               P%amps0(i) = P%amposc(i)
+               P%amps1(i) = 0d0
+               P%amps2(i) = 0d0
+               P%args0(i) = P%phase(i)
+               P%args1(i) = P%omega(i)/P%powosc(i)
+               P%args2(i) = 0d0
+               
+               P%at0(i) = 0d0
+               P%at1(i) = P%ant(i)
+               P%at2(i) = 0d0
+               P%ampt0(i) = P%amps0(i)
+               P%ampt1(i) = 0d0
+               P%ampt2(i) = 0d0
+               P%argt0(i) = P%phase(i)
+               P%argt1(i) = P%omega(i)
+               P%argt2(i) = 0d0
+              
+               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        
+               !First order shape          
+               
+               P%invsigma0(i)=10**(P%loginvsigma0(i))
+               
+               P%as0(i) = 1d0 + P%eps1(i)*(-2d0*(C_const+1d0))-P%eps2(i)*C_const
+               P%as1(i) = -2d0*P%eps1(i) - P%eps2(i)
+               P%as2(i) = 0.d0
+               P%amps0(i) = 2d0*P%Xsigma0(i)
+               P%amps1(i) = 3.141592*P%Xsigma0(i)*(2d0*P%eps1(i) + P%eps2(i))
+               P%amps2(i) = 0d0
+               if (Tpl_Phase) then
+                  P%args0(i) = 2d0*P%invsigma0(i)*(1d0 + P%eps1(i)) + P%varphi(i)
+               else
+                  P%args0(i) = P%varphi(i)
+               endif
+               P%args1(i) = 2d0*P%invsigma0(i)*P%eps1(i)
+               P%args2(i) = 0d0
+               
+               P%at0(i) = 1d0 + P%eps1(i)*(-2d0*(C_const+1d0))
+               P%at1(i) = -2d0*P%eps1(i)
+               P%at2(i) = 0d0
+               P%ampt0(i) = P%amps0(i)
+               P%ampt1(i) = 3.141592*P%Xsigma0(i)*2d0*P%eps1(i)
+               P%ampt2(i) = 0d0
+               P%argt0(i) = P%args0(i)
+               P%argt1(i) = P%args1(i)
+               P%argt2(i) = 0d0
+               
+               !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)*16d0
+ 
+            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
  
--- 281,320 ----
          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)-1d0)
! 
!            ScalarPower = exp(P%as1(in)*lnrat + P%as2(in)*lnrat**2) &
!                 * (1.d0 - 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) * (1d0 - 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 =0d0
!         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)
--- 335,376 ----
          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)-1d0)
+ 
+            TensorPower = exp(P%at1(in)*lnrat) &
+                 * (1.d0 - 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) * (1d0 - 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=0d0
+         endif
+         
+ !end wig
+ 
        end function TensorPower
  
        !Get parameters describing parameterisation (for FITS file)
***************
*** 156,161 ****
--- 381,389 ----
           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 ****
--- 408,460 ----
             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) = '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
diff -r -c -b -B cosmomc/distparams.ini cosmomc_wiggles/distparams.ini
*** cosmomc/distparams.ini	2005-02-23 16:43:30.000000000 +0000
--- cosmomc_wiggles/distparams.ini	2005-02-27 16:01:17.000000000 +0000
***************
*** 31,36 ****
--- 31,37 ----
  
  plot_meanlikes = T
  shade_meanlikes = T
+ td_meanlikes = F
  
  # if non-zero, output _thin file, thinned by thin_factor
  thin_factor = 0
***************
*** 71,82 ****
  #Number of parameters to get covariance matrix for 
  #If you are going to use the output as a proposal density make sure
  #you have map_params = F, and the dimension equal to the number of MCMC parameters
! cov_matrix_dimension = 13
  
! #e.g. colormap('jet')
  matlab_colscheme =
  
! #Parameters to use. If zero use all parameters which have lables.
  plotparams_num = 0
  plotparams = 1 2 3 11
  
--- 72,83 ----
  #Number of parameters to get covariance matrix for 
  #If you are going to use the output as a proposal density make sure
  #you have map_params = F, and the dimension equal to the number of MCMC parameters
! cov_matrix_dimension = 17
  
! #e.g. colormap('jet'), colormap('copper') for td_meanlikes
  matlab_colscheme =
  
! #Parameters to use. If zero use all parameters which have labels.
  plotparams_num = 0
  plotparams = 1 2 3 11
  
***************
*** 92,114 ****
  #w
  lab7 = 
  
  lab8 = n_s
! #n_t
  lab9 =  
  lab10 = n_run
  
! lab11 = log[10^{10} A_s]
  #amp ratio
! lab12 = 
! lab13 =  
  
- lab14 = \Omega_\Lambda
- lab15 = Age/GYr
- lab16 = \Omega_m
- lab17 = \sigma_8
- lab18 = z_{re}
- lab19 = 
- lab20 = H_0
  
  #Need to give limits if prior cuts off distribution where not very small
  limits4 = 0.01 N
--- 93,128 ----
  #w
  lab7 = 
  
+ #n_s or \epsilon_1
  lab8 = n_s
! #n_t or \epsilon_2
  lab9 =  
+ #n_run
  lab10 = n_run
  
! # log(omega) or log(1 / \sigma_0)
! lab11 = log(omega)
! # amp_osc or |x| \sigma_0
! lab12 = amposc
! # phase
! lab13 = phase
! #log(power)
! lab14 = log(power)
! 
! #ln[10^{10} A_s] or ln[10^{10} P_scalar]
! lab15 = ln[10^{10} A_s]
  #amp ratio
! lab16 = 
! lab17 =  
! 
! lab18 = \Omega_\Lambda
! lab19 = Age/GYr
! lab20 = \Omega_m
! lab21 = \sigma_8
! lab22 = z_{re}
! lab23 = 
! lab24 = H_0
  
  
  #Need to give limits if prior cuts off distribution where not very small
  limits4 = 0.01 N
diff -r -c -b -B cosmomc/params.ini cosmomc_wiggles/params.ini
*** cosmomc/params.ini	2005-02-23 16:43:37.000000000 +0000
--- cosmomc_wiggles/params.ini	2005-02-27 16:01:17.000000000 +0000
***************
*** 10,16 ****
  samples = 200000
  
  #Feedback level ( 2=lots,1=chatty,0=none)
! feedback = 1
  
  #Temperature at which to Monte-Carlo
  temperature = 1
--- 10,16 ----
  samples = 200000
  
  #Feedback level ( 2=lots,1=chatty,0=none)
! feedback = 3
  
  #Temperature at which to Monte-Carlo
  temperature = 1
***************
*** 45,51 ****
  
  #Can use covariance matrix for proposal density, otherwise use settings below
  #Covariance matrix can be produced using "getdist" program.
! propose_matrix = params_CMB.covmat
  
  #If propose_matrix is blank (first run), can try to use numerical Hessian to 
  #estimate a good propose matrix. As a byproduct you also get an approx best fit point
--- 45,51 ----
  
  #Can use covariance matrix for proposal density, otherwise use settings below
  #Covariance matrix can be produced using "getdist" program.
! #propose_matrix = params_CMB.covmat
  
  #If propose_matrix is blank (first run), can try to use numerical Hessian to 
  #estimate a good propose matrix. As a byproduct you also get an approx best fit point
***************
*** 107,118 ****
  
  #CAMB parameters
  #If we are including tensors
! compute_tensors = F
  #Whether the CMB should be lensed (slows a lot unless also computing matter power)
  CMB_lensing = F
  #increase accuracy_level to run CAMB on higher accuracy
  #(default is about 0.3%, accuracy_level=2 around 0.1% at high l)
! accuracy_level = 1
  
  #If action = 1
  redo_likelihoods = T
--- 107,124 ----
  
  #CAMB parameters
  #If we are including tensors
! compute_tensors = T
  #Whether the CMB should be lensed (slows a lot unless also computing matter power)
  CMB_lensing = F
  #increase accuracy_level to run CAMB on higher accuracy
  #(default is about 0.3%, accuracy_level=2 around 0.1% at high l)
! accuracy_level = 6
! l_accuracy_level = 2
! l_sample_level = 2
! k_src_sample_level = 1.
! k_int_sample_level = 1.
! camb_compute_every_l = T
! 
  
  #If action = 1
  redo_likelihoods = T
***************
*** 143,159 ****
  #w
  param7 = -1 -1 -1 0 0
  
! #n_s
! param8 = 1 0.5 1.5 0.03 0.01
! #n_t
! param9 = 0 0 0 0 0
! #n_run
! param10 = 0 0 0 0 0
! 
! #log[10^10 A_s]
! param11 = 3.1 2.7 4 0.2 0.03
! #amp_ratio
! param12 = 0 0 0 0 0
! #unused fast parameter
! param13 = 0 0 0 0 0
! 
--- 149,210 ----
  #w
  param7 = -1 -1 -1 0 0
  
! ########################################################################################
! #	Decomment ONLY the relevant power spectrum (see power_tilt.f90)                #
! ########################################################################################
! 
! 
! ########################################################################################
! #		        	For Power_Name = power_tilt                            #
! ########################################################################################
! 
! ##n_s
! #param8 = 1 0.5 1.5 0.03 0.02
! ##n_t
! #param9 = 0 0 0 0 0
! ##n_run
! #param10 = 0 0 0 0 0
! 
! ##logfreq
! #param11 = 1.5 1 2.5 0.4 0.1
! ##amp_osc
! #param12 = 0 0 0.8 0.4 0.1
! ##phase
! #param13 = 0 0 3.14 0.4 0.1	
! ##logpower_osc
! #param14 = 0 0 0 0 0
! 
! ##ln[10^10 A_s]
! #param15 = 3.1 2.7 4 0.4 0.1
! ##amp_ratio
! #param16 = 0 0 0 0 0
! ##2dF amp ratio
! ##zero means do it analytically
! #param17 = 0 0 0 0 0
! 
! ########################################################################################
! #		        	For Power_Name = power_sr                              #
! ########################################################################################
! 
! ##eps1
! #param8 = 0.05 0.0001 0.1 0.001 0.004
! ##eps2
! #param9 = -0.04 -0.3 0.3  0.01 0.02
! ##unused
! #param10 = 0 0 0 0 0
! 
! ##loginvsigma0
! #param11 = 3.4 2.5 3.5 0.05 0.1
! ##X times sigma0
! #param12 = 0.05 0 0.45 0.1 0.2
! ##varphi
! #param13 = 0 0 6.28 0.4 0.1
! ##unused
! #param14 = 0 0 0 0 0
! 
! ##ln[P_scalar/cl_norm] (see setting.f90)
! #param15 = 3.1 2.7 4 0.4 0.1
! ##non-zero required for tensors
! #param16 = 1 1 1 0 0
! ##unused
! #param17 = 0 0 0 0 0
diff -r -c -b -B cosmomc/source/CMB_Cls_simple.f90 cosmomc_wiggles/source/CMB_Cls_simple.f90
*** cosmomc/source/CMB_Cls_simple.f90	2004-09-29 21:05:38.000000000 +0100
--- cosmomc_wiggles/source/CMB_Cls_simple.f90	2005-02-27 16:01:17.000000000 +0000
***************
*** 262,267 ****
--- 262,289 ----
     type(CAMBParams)  P 
     integer zix
  
+ !wig
+    AccuracyBoost = AccuracyLevel
+    lAccuracyBoost = lAccuracyLevel
+    lSampleBoost = lSampleLevel
+    kSrcSampleBoost = kSrcSampleLevel
+    kIntSampleBoost = kIntSampleLevel
+    ComputeEveryl = CambComputeEveryl
+ 
+    write(*,*)
+    write(*,*) '>>>>>>>>>>>>>> accuracy settings <<<<<<<<<<<<<<'
+    write(*,*) 'Power_Name     =    ',Power_Name
+    write(*,*) 'write_all_cls  =   ',write_all_cls
+    write(*,*) 'AccuracyBoost  = ',AccuracyBoost
+    write(*,*) 'lAccuracyBoost = ',lAccuracyBoost          
+    write(*,*) 'DoClsInterpol  =   ',DoClsInterpol
+    write(*,*) 'ComputeEveryl  =   ',ComputeEveryl
+    write(*,*) '>>>>>>>>>>>>>> accuracy settings <<<<<<<<<<<<<<'
+    write(*,*)
+ 
+    if ((Power_Name == 'power_sr').and.(.not.compute_tensors)) stop 'Slow-Roll requires tensors'
+ !end wig
+ 
          Threadnum =num_threads
          w_lam = -1
          call CAMB_SetDefParams(P)
***************
*** 284,292 ****
              P%Transfer%high_precision=.true.
              P%Transfer%k_per_logint=3*nint(AccuracyLevel)
            end if
!           AccuracyBoost = AccuracyLevel
!           lAccuracyBoost = AccuracyLevel
!           lSampleBoost = AccuracyLevel
            P%AccurateReionization = .true.
          end if
                  
--- 306,316 ----
              P%Transfer%high_precision=.true.
              P%Transfer%k_per_logint=3*nint(AccuracyLevel)
            end if
! !wig
! !          AccuracyBoost = AccuracyLevel
! !          lAccuracyBoost = AccuracyLevel
! !          lSampleBoost = AccuracyLevel
! !end wig
            P%AccurateReionization = .true.
          end if
                  
***************
*** 337,343 ****
            write (*,*) 'Computing tensors:', compute_tensors
            write (*,*) 'Doing CMB lensing:',CMB_lensing
            write(*,'(" lmax           = ",1I4)') lmax       
!           if (compute_tensors) write(*,'("lmax_tensor    = ",1I4)') lmax_tensor
            write(*,'(" Number of C_ls = ",1I4)') num_cls
          end if
  
--- 361,367 ----
            write (*,*) 'Computing tensors:', compute_tensors
            write (*,*) 'Doing CMB lensing:',CMB_lensing
            write(*,'(" lmax           = ",1I4)') lmax       
!           if (compute_tensors) write(*,'(" lmax_tensor    = ",1I4)') lmax_tensor
            write(*,'(" Number of C_ls = ",1I4)') num_cls
          end if
  
diff -r -c -b -B cosmomc/source/cmbtypes.f90 cosmomc_wiggles/source/cmbtypes.f90
*** cosmomc/source/cmbtypes.f90	2004-09-29 21:05:38.000000000 +0100
--- cosmomc_wiggles/source/cmbtypes.f90	2005-02-27 16:01:17.000000000 +0000
***************
*** 57,63 ****
    end Type CosmoTheory
  
    logical, parameter ::  Old_format  = .false.
!   logical, parameter :: write_all_Cls = .false. 
     !if false use CAMB's flat interpolation scheme (lossless if models are flat except near lmax when lensed)
  
  contains
--- 57,66 ----
    end Type CosmoTheory
  
    logical, parameter ::  Old_format  = .false.
! !wig
! !  logical, parameter :: write_all_Cls = .false. 
!   logical, parameter :: write_all_Cls = .true.
! !end wig
     !if false use CAMB's flat interpolation scheme (lossless if models are flat except near lmax when lensed)
  
  contains
diff -r -c -b -B cosmomc/source/driver.F90 cosmomc_wiggles/source/driver.F90
*** cosmomc/source/driver.F90	2004-10-14 14:36:03.000000000 +0100
--- cosmomc_wiggles/source/driver.F90	2005-02-27 16:01:17.000000000 +0000
***************
*** 30,35 ****
--- 30,39 ----
          double precision intime
          integer ierror
  
+ !chris
+         integer :: instance_shift
+ !end chris
+ 
   call mpi_init(ierror)
   if (ierror/=MPI_SUCCESS) stop 'MPI fail: rank'
  #endif
***************
*** 52,65 ****
  
          propose_scale = Ini_Read_Real('propose_scale',2.4)
          AccuracyLevel = Ini_Read_Real('accuracy_level',1.)
  
  #ifdef MPI 
          
!         if (instance /= 0) call DoStop('With MPI should not have second parameter')
          call mpi_comm_rank(mpi_comm_world,MPIrank,ierror)
  
          instance = MPIrank +1 !start at 1 for chains
!         write (numstr,*) instance
          rand_inst = instance
          if (ierror/=MPI_SUCCESS) call DoStop('MPI fail')
  
--- 56,84 ----
  
          propose_scale = Ini_Read_Real('propose_scale',2.4)
          AccuracyLevel = Ini_Read_Real('accuracy_level',1.)
+ !wig
+         lAccuracyLevel = Ini_Read_Real('l_accuracy_level',1.)
+         lSampleLevel = Ini_Read_Real('l_sample_level',1.)
+         kSrcSampleLevel = Ini_Read_Real('k_src_sample_level',1.)
+         kIntSampleLevel = Ini_Read_Real('k_int_sample_level',1.)
+         CambComputeEveryl = Ini_Read_Logical('camb_compute_every_l',.false.)
+ !end wig
  
  #ifdef MPI 
          
! !chris        
! !        if (instance /= 0) call DoStop('With MPI should not have second parameter')
!         instance_shift=instance
!         write(*,*)'MPI file names shifted by: ',instance_shift
! !end chris
          call mpi_comm_rank(mpi_comm_world,MPIrank,ierror)
  
          instance = MPIrank +1 !start at 1 for chains
! !chris
! !        write (numstr,*) instance
! 	instance_shift = MPIrank + instance_shift
! 	write (numstr,*) instance_shift
! !end chris
          rand_inst = instance
          if (ierror/=MPI_SUCCESS) call DoStop('MPI fail')
  
diff -r -c -b -B cosmomc/source/GetDist.f90 cosmomc_wiggles/source/GetDist.f90
*** cosmomc/source/GetDist.f90	2005-01-18 18:26:36.000000000 +0000
--- cosmomc_wiggles/source/GetDist.f90	2005-02-28 14:31:53.047179482 +0000
***************
*** 65,70 ****
--- 65,73 ----
          logical force_twotail 
          logical smoothing, plot_meanlikes, plot_NDcontours
          logical shade_meanlikes, make_single_samples
+ !3d
+         logical td_meanlikes
+ !end 3d
          integer single_thin,cust2DPlots(max_cols**2)
          integer ix_min(max_cols),ix_max(max_cols)
          real center(max_cols),width(max_cols)   
***************
*** 1141,1153 ****
   !             write (aunit,fmt) (/(I, I=ix_min(j), ix_max(j), 1)/)*width(j)+center(j)
                if (DoShade) then
                if (shade_meanlikes) then
                  write (aunit,'(a)') 'load (fullfile(''plot_data'',''' // trim(plotfile) //'_likes''));'
                  write (aunit,*) 'contourf(x1,x2,'//trim(plotfile)//'_likes,64);'
                else     
                  write (aunit,*) 'contourf(x1,x2,pts,64);'
                end if
! 
                write (aunit,*) 'set(gca,''climmode'',''manual''); shading flat; hold on;'
                end if
                
             
--- 1144,1170 ----
   !             write (aunit,fmt) (/(I, I=ix_min(j), ix_max(j), 1)/)*width(j)+center(j)
                if (DoShade) then
                if (shade_meanlikes) then
+ !3d
                  write (aunit,'(a)') 'load (fullfile(''plot_data'',''' // trim(plotfile) //'_likes''));'
+                  if (td_meanlikes) then
+                     write (aunit,*) 'surfl(x1,x2,'//trim(plotfile)//'_likes);'
+                  else
                  write (aunit,*) 'contourf(x1,x2,'//trim(plotfile)//'_likes,64);'
+                  endif
+ !end 3d
                else     
                  write (aunit,*) 'contourf(x1,x2,pts,64);'
                end if
! !3d
! !              write (aunit,*) 'set(gca,''climmode'',''manual''); shading flat; hold on;'
!               if (td_meanlikes) then
!                  write (aunit,*) 'set(gca,''climmode'',''manual'');shading interp; hold on;'
!                  write (aunit,*) 'grid off;box on;camproj(''orthographic'');'
!                  write (aunit,*) 'set(gca,''climmode'',''manual'');alpha(0.9);axis tight'
!               else
                write (aunit,*) 'set(gca,''climmode'',''manual''); shading flat; hold on;'
+               endif
+ !end 3d
                end if
                
             
***************
*** 1466,1471 ****
--- 1483,1491 ----
            write (*,*) 'WARNING: Mapping params - .covmat file is new params.'
  
      shade_meanlikes = Ini_Read_Logical('shade_meanlikes',.false.)   
+ !3d
+     td_meanlikes = Ini_Read_Logical('td_meanlikes',.false.)  
+ !end 3d   
  
      matlab_col = Ini_Read_String('matlab_colscheme')
      smoothing = Ini_Read_Logical('smoothing',.true.)
***************
*** 1959,1965 ****
--- 1979,1992 ----
        write (51,*)  'set(gcf, ''PaperUnits'',''inches'');'
       
        if (plot_row < plot_col .and. plot_row*plot_col>9) then
+ !3d
+ !            write (51,*) 'set(gcf, ''PaperPosition'',[ 0 0 8 10]);';
+          if (td_meanlikes) then
+             write (51,*) 'set(gcf, ''PaperPosition'',[ 0 0 7 9]);';
+          else
              write (51,*) 'set(gcf, ''PaperPosition'',[ 0 0 8 10]);';
+          endif
+ !end 3d           
            else
              write (51,*) 'set(gcf, ''PaperPosition'',[ 0 0 8 8]);';
        end if
diff -r -c -b -B cosmomc/source/Makefile cosmomc_wiggles/source/Makefile
*** cosmomc/source/Makefile	2004-09-30 20:07:02.000000000 +0100
--- cosmomc_wiggles/source/Makefile	2005-02-27 16:01:17.000000000 +0000
***************
*** 1,90 ****
! #Settings here are for Intel ifc, uncomment others below:
! #You may need to edit the library paths for MKL for Intel
! #Beware of using optmizations which lose accuracy - may give errors when running
! #With Visual Fortan compile open cosmomc.dsw workspace file instead.
! 
! #Intel MPI (assuming mpif77 set to point to ifc)
! #change -lmkl_ia32 to -lmkl_p3 for MKL versions earlier than 6 (6 needed for ifc 8+) 
! F90C     = mpif77
! FFLAGS = -O2 -Vaxlib -ip -W0 -WB -openmp -fpp -DMPI -quiet
! LAPACKL = -L/opt/intel/mkl61/lib/32 -lmkl_lapack -lmkl_ia32 -lguide -lpthread -static
  INCLUDE = -I../camb 
  
! #Intel fortran 8, check you have the latest update from the Intel web pages
! #See Makefile_intel for ifc 7.1 or lower (some versions have problems)
! #F90C     = ifort
! #FFLAGS = -O2 -Vaxlib -ip -W0 -WB -openmp -fpp
! #LAPACKL = -L/opt/intel/mkl61/lib/32 -lmkl_lapack -lmkl_ia32 -lguide -lpthread -static
! #INCLUDE = -I../camb 
! 
! #Intel fortran 8 on COSMOS (do "module load dev8" to get v8 compilers)
! #F90C     = ifort
! #FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp2 -lmpi -lscs  -ldl -DMPI
! #LAPACKL = -lmkl_lapack -lmkl_ipf -limf
! #INCLUDE = -I../camb
! 
! #SGI, -mp toggles multi-processor. Use -O2 if -Ofast gives problems.
! #Not various versions of the compiler are buggy giving erroneous seg faults with -mp.
! #Version 7.3 is OK, currently version 7.4 is bugged, as are some earlier versions.
! #F90C     = f90
! #LAPACKL = -lcomplib.sgimath
! #INCLUDE = -I../camb 
! #FFLAGS  = -Ofast -mp
! 
! #Digital/Compaq fortran, -omp toggles multi-processor
! #F90C    = f90
! #FFLAGS  = -omp -O4 -arch host -math_library fast -tune host -fpe1
! #LAPACKL  = -lcxml
! #INCLUDE = -I../camb 
! 
! #Intel Itanium efc, -openmp toggles multi-processor:
! #F90C     = efc
! #FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp
! #LAPACKL = -L/opt/intel/mkl/lib/64 -lmkl_lapack -lmkl_ipf -lguide -lpthread -static
! #INCLUDE = -I../camb 
! 
! #Intel Itanium efc MPI:
! # mpirun -np 4 ./cosmomc params.ini
! #F90C     = efc
! #FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp -lmpi -lscs  -ldl -DMPI
! #LAPACKL = -L/opt/intel/mkl/lib/64 -lmkl_lapack -lmkl_ipf -limf
! #INCLUDE = -I../camb
! 
! #Absoft ProFortran, single processor, set -cpu:[type] for your local system 
! #F90C     = f95
! #FFLAGS = -O2 -s -cpu:athlon -lU77 -w -YEXT_NAMES="LCS" -YEXT_SFX="_" 
! #LAPACKL =  -llapack -lblas -lg2c 
! #INCLUDE = -I../camb -p../camb
! 
! #NAGF95, single processor:
! #F90C     = f95
! #FFLAGS = -DNAGF95 -O3
! #LAPACKL = -llapack -lblas -lg2c
! #INCLUDE = -I../camb 
! 
! #Sun, single processor:
! #F90C     = f90
! #FFLAGS = -O5
! #LAPACKL = -lsunperf -lfsu
! #INCLUDE = -I../camb -M../camb
! 
! #Sun parallel enterprise:
! #F90C     = f95
! #FFLAGS =  -O4 -xarch=native64 -openmp -ftrap=%none
! #LAPACKL = -lsunperf -lfsu
! #INCLUDE = -I../camb -M../camb
! 
! #IBM XL Fortran, multi-processor (run "module load lapack" then run "gmake")
! #F90C     = xlf90_r $(LAPACK)
! #FFLAGS  = -WF,-DIBMXL -qsmp=omp -qsuffix=f=f90:cpp=F90 -O3 -qstrict -qarch=pwr3 -qtune=pwr3
! #INCLUDE = -I../camb -lessl
! #LAPACKL =
  
  
  PROPOSE = propose.o
  CLSFILE = CMB_Cls_simple.o
  
! #Can use params_H if you prefer more generic parameters
  PARAMETERIZATION = params_CMB.o
  
  CLSLIB  = -L../camb -lcamb $(LAPACKL)
--- 1,49 ----
! # >>> DESIGNED FOR GMAKE <<<
! 
! # Unified Systems makefile for COSMOMC
! # Add FLAGS -DMPI for using MPI
! 
! ext=$(shell uname | cut -c1-3)
! 
! ifeq ($(ext),IRI)
! F90C= f90
! FFLAGS= -O3 -mp -64
! LAPACKL = -lcomplib.sgimath
! INCLUDE = -I../camb
! endif
! 
! ifeq ($(ext),Lin)
! F90C=ifort
! FFLAGS= -O3 -openmp -fpp -w95 -WB -pc80
! LAPACKL = -lmkl_lapack32 -lmkl_p3 
  INCLUDE = -I../camb 
+ endif 
  
! ifeq ($(ext),OSF)
! F90C=f90
! FFLAGS= -omp -O -arch host -math_library fast -tune host -fpe1
! LAPACKL  = -lcxml
! INCLUDE = -I../camb 
! endif
  
+ ifeq ($(ext),Sun)
+ F90C=f90
+ FFLAGS= -O4 -xarch=native64 -openmp -ftrap=%none
+ LAPACKL = -lsunperf -lfsu
+ INCLUDE = -I../camb -M../camb
+ endif
+ 
+ ifeq ($(ext),AIX)
+ F90C     = xlf90_r $(LAPACK)
+ FFLAGS  = -WF,-DIBMXL -qsmp=omp -qsuffix=f=f90:cpp=F90 -O3 -qstrict -qarch=pwr3 -qtune=pwr3
+ INCLUDE = -I../camb -lessl
+ LAPACKL =
+ endif
  
  PROPOSE = propose.o
  CLSFILE = CMB_Cls_simple.o
  
! #PARAMETRIZATION = params_H.o
  PARAMETERIZATION = params_CMB.o
  
  CLSLIB  = -L../camb -lcamb $(LAPACKL)
***************
*** 99,107 ****
  	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
  
  
! default: cosmomc
  
! all : cosmomc getdist
  
  .f.o:
  	f77 $(F90FLAGS) -c $<
--- 58,66 ----
  	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
  
  
! default: cosmomc.$(ext)
  
! all : cosmomc.$(ext) getdist.$(ext)
  
  .f.o:
  	f77 $(F90FLAGS) -c $<
***************
*** 112,124 ****
  %.o: %.F90
  	$(F90C) $(F90FLAGS) -c $*.F90
  
! 
! cosmomc: $(OBJFILES) 	
! 	$(F90C) -o ../cosmomc $(OBJFILES) $(CLSLIB) $(F90FLAGS)
! 
  
  clean:
  	rm -f *.o *.mod *.d *.pc *.obj ../core
  
! getdist: $(DISTFILES)
! 	$(F90C) -o ../getdist $(DISTFILES) $(CLSLIB) $(F90FLAGS) 
--- 71,81 ----
  %.o: %.F90
  	$(F90C) $(F90FLAGS) -c $*.F90
  
! cosmomc.$(ext): $(OBJFILES) 	
! 	$(F90C) -o ../$@ $(OBJFILES) $(CLSLIB) $(F90FLAGS)
  
  clean:
  	rm -f *.o *.mod *.d *.pc *.obj ../core
  
! getdist.$(ext): $(DISTFILES)
! 	$(F90C) -o ../$@ $(DISTFILES) $(CLSLIB) $(F90FLAGS) 
diff -r -c -b -B cosmomc/source/params_CMB.f90 cosmomc_wiggles/source/params_CMB.f90
*** cosmomc/source/params_CMB.f90	2004-09-29 21:05:38.000000000 +0100
--- cosmomc_wiggles/source/params_CMB.f90	2005-02-27 16:01:17.000000000 +0000
***************
*** 80,86 ****
  
         integer, intent(in) :: in
  
! 
         if (Power_Name == 'power_tilt') then
  
         P%InitPower%ScalarPowerAmp(in) = cl_norm*CMB%norm(norm_As)
--- 80,86 ----
  
         integer, intent(in) :: in
  
! !wig
         if (Power_Name == 'power_tilt') then
  
         P%InitPower%ScalarPowerAmp(in) = cl_norm*CMB%norm(norm_As)
***************
*** 90,99 ****
--- 90,119 ----
         P%InitPower%ant(in) = CMB%InitPower(2)
         P%InitPower%n_run(in) = CMB%InitPower(3)
  
+           P%InitPower%logomega(in) = CMB%InitPower(4)
+           P%InitPower%omega(in) = 10**(P%InitPower%logomega(in))
+           P%InitPower%amposc(in) = CMB%InitPower(5)          
+           P%InitPower%phase(in) = CMB%InitPower(6)
+           P%InitPower%logpowosc(in) = CMB%InitPower(7)
+           P%InitPower%powosc(in) = 10**(P%InitPower%logpowosc(in))
+           
+        elseif (Power_Name == 'power_sr') then
+          
+           P%InitPower%P_scalar =  cl_norm*CMB%norm(norm_As)
+           P%InitPower%eps1(in) =  CMB%InitPower(1)
+           P%InitPower%eps2(in) =  CMB%InitPower(2)
+           !CMB%InitPower(3) unused. Slot for eps3
+ 
+           P%InitPower%loginvsigma0(in) = CMB%InitPower(4)
+           P%InitPower%invsigma0(in) = 10**(P%InitPower%loginvsigma0(in))
+           P%InitPower%Xsigma0(in) = CMB%InitPower(5)
+           P%InitPower%varphi(in)  = CMB%InitPower(6)
+           !CMB%InitPower(7) unused.
  
         else
           stop 'params_CMB:Wrong initial power spectrum'
         end if
+ !end wig
  
      end subroutine SetCAMBInitPower
   
diff -r -c -b -B cosmomc/source/settings.f90 cosmomc_wiggles/source/settings.f90
*** cosmomc/source/settings.f90	2004-09-29 21:05:38.000000000 +0100
--- cosmomc_wiggles/source/settings.f90	2005-02-27 16:01:17.000000000 +0000
***************
*** 5,10 ****
--- 5,17 ----
    implicit none
  
    real :: AccuracyLevel = 1.
+ !wig
+   real :: lAccuracyLevel = 1.
+   real :: lSampleLevel = 1.
+   real :: kSrcSampleLevel = 1.
+   real :: kIntSampleLevel = 1.
+   logical :: CambComputeEveryl = .false.
+ !end wig
    !Set to >1 to use CAMB etc on higher accuracy settings. 
    !Does not affect MCMC (except making it all slower)
  
***************
*** 13,19 ****
  !num_hard is the number of 'hard' parameters
  
    integer, parameter :: num_hard =7 
!   integer, parameter :: num_initpower = 3 
    integer, parameter :: num_norm = 3 
       !Should be 2 or larger (scale for the scalar Cl, and the T/S ratio
  
--- 20,29 ----
  !num_hard is the number of 'hard' parameters
  
    integer, parameter :: num_hard =7 
! !wig
! !  integer, parameter :: num_initpower = 3 
!   integer, parameter :: num_initpower = 7
! !end wig 
    integer, parameter :: num_norm = 3 
       !Should be 2 or larger (scale for the scalar Cl, and the T/S ratio
  
