diff -r -c -b -B cosmomc/camb/cmbmain.f90 cosmomc_wiggles/camb/cmbmain.f90
*** cosmomc/camb/cmbmain.f90	2010-01-27 11:54:17.000000000 +0100
--- cosmomc_wiggles/camb/cmbmain.f90	2011-03-23 17:23:54.856627807 +0100
***************
*** 635,651 ****
  !     use a logarithmic spacing. closed case dealt with by SetClosedkValues
  
           if (CP%WantScalars .and. CP%Reion%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
--- 635,661 ----
  !     use a logarithmic spacing. closed case dealt with by SetClosedkValues
  
           if (CP%WantScalars .and. CP%Reion%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
***************
*** 1029,1034 ****
--- 1039,1048 ----
         qmax_int = min(qmax,max_bessels_etak/CP%tau0)
  
         IntSampleBoost=AccuracyBoost 
+ !wig
+ !     IntSampleBoost=AccuracyBoost
+       IntSampleBoost=AccuracyBoost*kIntSampleBoost
+ 
         if (do_bispectrum) IntSampleBoost = IntSampleBoost * 2  
         if (hard_bispectrum) IntSampleBoost = IntSampleBoost * 2  
  
diff -r -c -b -B cosmomc/camb/inidriver.F90 cosmomc_wiggles/camb/inidriver.F90
*** cosmomc/camb/inidriver.F90	2009-12-15 16:34:40.000000000 +0100
--- cosmomc_wiggles/camb/inidriver.F90	2010-05-27 15:13:25.893017186 +0200
***************
*** 256,261 ****
--- 256,275 ----
           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'
***************
*** 267,272 ****
--- 281,319 ----
         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').or.(Power_Name.eq.'power_dbi')) 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 (*,'("DBI       delta1     = ",f12.8)') P%InitPower%del1(i)
+                 write (*,'("DBI       gamma      = ",f12.8)') P%InitPower%gam(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	2009-10-22 15:59:31.000000000 +0200
--- cosmomc_wiggles/camb/Makefile	2010-05-27 15:13:25.893017186 +0200
***************
*** 1,69 ****
! #CAMB Makefile
  
! #Edit for your compiler
! #Note there are many ifc versions, some of which behave oddly
  
- #Intel , -openmp toggles mutli-processor:
- #note version 10.0 gives wrong result for lensed when compiled with -openmp [fixed in 10.1]
- F90C     = ifort
- FFLAGS = -openmp -O2 -ip -W0 -WB -fpp2 -vec_report0
- 
- # 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
- 
- #G95 compiler
- #F90C   = g95
- #FFLAGS = -O2
- 
- #Gfortran compiler: if pre v4.3 add -D__GFORTRAN__
- #F90C     = gfc
- #FFLAGS =  -O2 
- 
- #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
- 
- #PGF90
- #F90C = pgf90
- #FFLAGS = -O2 -DESCAPEBACKSLASH
- 
- #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  = -DESCAPEBACKSLASH -DIBMXL -qsmp=omp -qsuffix=f=f90:cpp=F90 -O3 -qstrict -qarch=pwr3 -qtune=pwr3
- 
- #Settings for building camb_fits
- #Location of FITSIO and name of library
- FITSDIR       = /home/cpac/cpac-tools/lib
- FITSLIB       = cfitsio
- #Location of HEALPIX for building camb_fits
- HEALPIXDIR    = /home/cpac/cpac-tools/healpix
  
- include ./Makefile_main
\ No newline at end of file
--- 1,91 ----
!  # >>> 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=gfortran
! FFLAGS= -O -fopenmp #-D__GFORTRAN__
! 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
! EQUATIONS     = equations
! POWERSPECTRUM = power_tilt
! REIONIZATION = reionization
! RECOMBINATION = recfast
! 
! #Module doing non-linear scaling
! NONLINEAR     = halofit
! 
! #Driver program
! DRIVER        = inidriver.F90
! #DRIVER        = sigma8.f90
! #DRIVER        = tester.f90
! 
! 
! CAMBLIB       = libcamb.a
! 
! #Shouldn't need to change anything else...
! 
! F90FLAGS      = $(FFLAGS)
! HEALPIXLD     = -L$(HEALPIXDIR)/lib -lhealpix -L$(FITSDIR) -l$(FITSLIB)
! 
! CAMBOBJ       = constants.o utils.o subroutines.o inifile.o $(POWERSPECTRUM).o $(RECOMBINATION).o $(REIONIZATION).o modules.o \
!         bessels.o $(EQUATIONS).o $(NONLINEAR).o lensing.o cmbmain.o camb.o
! 
! default: camb.$(ext)
! 
! all: camb.$(ext) $(CAMBLIB)
! 
! subroutines.o: constants.o utils.o
! $(POWERSPECTRUM): subroutines.o  inifile.o
! $(RECOMBINATION).o: subroutines.o inifile.o
! $(REIONIZATION).o: constants.o inifile.o
! modules.o: $(REIONIZATION).o
! bessels.o: modules.o
! $(EQUATIONS): bessels.o
! $(NONLINEAR).o:  modules.o
! lensing.o: bessels.o
! cmbmain.o: lensing.o
! camb.o: cmbmain.o
! 
! 
! camb.$(ext): $(CAMBOBJ) $(DRIVER)
! 	$(F90C) $(F90FLAGS) $(CAMBOBJ) $(DRIVER) -o $@
! 
! $(CAMBLIB): $(CAMBOBJ)
! 	ar -r $@ $?
! 
! camb_fits: writefits.f90 $(CAMBOBJ) $(DRIVER)
! 	$(F90C) $(F90FLAGS) -I$(HEALPIXDIR)/include tils.o $(CAMBOBJ) writefits.f90 $(DRIVER) $(HEALPIXLD) -DWRITE_FITS -o $@
! 
! %.o: %.f90
! 	$(F90C) $(F90FLAGS) -c $*.f90
! 
! utils.o:
! 	$(F90C) $(F90FLAGS) -c utils.F90	
! 
! clean:
! 	-rm -f *.o *.a *.d core *.mod
  
  
diff -r -c -b -B cosmomc/camb/modules.f90 cosmomc_wiggles/camb/modules.f90
*** cosmomc/camb/modules.f90	2010-05-19 15:55:22.000000000 +0200
--- cosmomc_wiggles/camb/modules.f90	2010-05-27 15:13:25.893017186 +0200
***************
*** 37,44 ****
          character(LEN=*), parameter :: version = 'May_10'
          
          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
            
--- 37,46 ----
          character(LEN=*), parameter :: version = 'May_10'
          
          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=.false.
! !end wig
  
          logical, parameter :: DebugEvolution = .false. !Set to true to do all the evolution for all k
            
***************
*** 185,200 ****
  
        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  !!!!!
--- 187,204 ----
  
        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  !!!!!
***************
*** 211,216 ****
--- 215,226 ----
  
  !     lmax is max possible number of l's evaluated
        integer, parameter :: lmax_arr = l0max !!!100+l0max/7
+ !wig
+       real(dl) :: kSrcSampleBoost=1._dl
+       real(dl) :: kIntSampleBoost=1._dl
+       logical :: DoClsInterpol=.true.
+       logical :: ComputeEveryl=.true.
+ !end wig
   
          contains
        
***************
*** 562,567 ****
--- 572,596 ----
          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       
         
          if (lSampleBoost >=50) then
***************
*** 697,702 ****
--- 726,734 ----
          end if
          lSet%l0=lind
          lSet%l(1:lind) = ls(1:lind)
+ !wig        
+         endif
+ !end wig
          
        end subroutine initlval
  
***************
*** 716,722 ****
  
        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
--- 748,755 ----
  
        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
***************
*** 732,737 ****
--- 765,784 ----
                         +(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
  
***************
*** 2073,2078 ****
--- 2120,2129 ----
               dlnkh = 0.02
               points = log(MTrans%TransferData(Transfer_kh,MTrans%num_q_trans,itf)/minkh)/dlnkh+1
  !             dlnkh = log(MTrans%TransferData(Transfer_kh,MTrans%num_q_trans,itf)/minkh)/(points-0.999)
+ !wig
+              write(*,*)'Transfert_SaveToFiles: points= ',points
+              write(*,*)'you may need to increase this to sample HF oscillations!'
+ !endwig
               allocate(outpower(points,CP%InitPower%nn))
               do in = 1, CP%InitPower%nn
                call Transfer_GetMatterPower(MTrans,outpower(1,in), itf, in, minkh,dlnkh, points)
diff -r -c -b -B cosmomc/camb/params.ini cosmomc_wiggles/camb/params.ini
*** cosmomc/camb/params.ini	2010-01-26 19:42:44.000000000 +0100
--- cosmomc_wiggles/camb/params.ini	2010-05-27 15:13:25.893017186 +0200
***************
*** 1,12 ****
  #Parameters for CAMB
  
  #output_root is prefixed to output file names
! output_root = test
  
  #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}
--- 1,12 ----
  #Parameters for CAMB
  
  #output_root is prefixed to output file names
! output_root = tpl
  
  #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}
***************
*** 63,68 ****
--- 63,76 ----
  #Fraction of total omega_nu h^2 accounted for by each eigenstate, eg. 0.5 0.5
  nu_mass_fractions = 1
  
+ #Type of initial power spectrum: power_tilt or power_sr
+ power_name = power_sr
+  
+ ######################################################################################
+ 
+ #for power_tilt power spectra
+ 
+ 
  #Initial power spectrum, amplitude, spectral index and running. Pivot k in Mpc^{-1}.
  initial_power_num         = 1
  pivot_scalar              = 0.05
***************
*** 96,101 ****
--- 104,140 ----
  RECFAST_Heswitch = 6
  RECFAST_Hswitch  = T
  
+ #primordial power_tilted wiggle
+ osc_logomega(1)           = 1.
+ osc_amp(1)                = 0.
+ 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] = 0.303161E+01
+ eps1(1) = 0.210391E-02
+ eps2(1) = 0.421584E-01
+ eps3(1) = 0
+ 
+ #dbi parameter
+ gam(1) = 1.
+ del1(1) = 0.
+ 
+ #wiggles parameters
+ logeps1invsigma0(1) = 1.
+ Xsigma0(1) = 0.
+ 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
***************
*** 187,189 ****
--- 226,240 ----
  #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         = T
+  
+ # Remove undersampling of the Cls: if true lSampleBoost has no effect
+ compute_every_l         = F
diff -r -c -b -B cosmomc/camb/power_tilt.f90 cosmomc_wiggles/camb/power_tilt.f90
*** cosmomc/camb/power_tilt.f90	2008-03-13 15:40:20.000000000 +0100
--- cosmomc_wiggles/camb/power_tilt.f90	2010-05-27 15:13:25.893017186 +0200
***************
*** 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, save :: useSecondOrder = .false.
! !end wig
        integer, parameter :: nnmax= 5 
         !Maximum possible number of different power spectra to use
  
***************
*** 51,56 ****
--- 54,109 ----
         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)
+ !dbi params
+        real(dl) gam(nnmax)      
+        real(dl) del1(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     
***************
*** 69,74 ****
--- 122,129 ----
          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
***************
*** 77,82 ****
--- 132,163 ----
           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%del1 = 0.
+          AP%gam = 1.
+         
+ 
+          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)
***************
*** 85,90 ****
--- 166,181 ----
           !Could read spectra from disk here, do other processing, etc.
  
          real(dl) acurv
+ !wig
+         integer :: i        
+ !from Hankel functions
+         real(dl), parameter :: C_const = -0.7296
+ !wkb or uniform approx gives a slightly different constant
+         real(dl), parameter :: D_const = -0.7652
+ !        real(dl), parameter :: C_const = D_const
+ 
+         real(dl), parameter :: pi = 3.14159265
+ !end wig
  
          if (AParamSet%nn > nnmax) then
             write (*,*) 'To use ',AParamSet%nn,'power spectra you need to increase'
***************
*** 95,100 ****
--- 186,347 ----
          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%powosc(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').or.(Power_Name.eq.'power_dbi')) 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: TPL not implemented! '
+                  endif
+                  if ((P%del1(i).ne.0.).or.(P%gam(i).ne.1.)) then
+                     useSecondorder = .false.
+                     write(*,*)'PowerIni: SR second order disabled: DBI slow-roll not implemented!'
+                  endif
+               endif
+ 
+               if ((power_Name.eq.'power_sr').and. &
+                    ((P%del1(i).ne.0.).or.(P%gam(i).ne.1.)) ) then
+                  write(*,*)'power_sr with non-vanishing DBI params!'
+                  write(*,*)'resetting to standard SR'
+                  P%del1(i) = 0.
+                  P%gam(i) = 1.
+               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 &
+                    + (C_const + 2.)*P%del1(i)
+               P%as1(i) = -2.*P%eps1(i) - P%eps2(i) + P%del1(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.
+               P%P_0_tensor(i)=P%P_scalar/(P%as0(i))*P%eps1(i)*16/P%gam(i)
+ 
+            enddo
+ 
+         else
+            stop 'No such initial power spectra'
+         endif
+ !end wig
  
         end subroutine InitializePowers
         
***************
*** 119,126 ****
          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
  
--- 366,405 ----
          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').or.(Power_Name.eq.'power_dbi')) 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,150 ****
     
          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
  
--- 421,463 ----
     
          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))
!  
!  
!          elseif ((Power_Name.eq.'power_sr').or.(Power_Name.eq.'power_dbi')) 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 (curv < 0) TensorPower=TensorPower*tanh(PiByTwo*sqrt(-k**2/curv-3)) 
  
+          if(TensorPower .lt. 0.0) then
+             TensorPower=0.
+          endif
+          
+ !end wig
+ 
         
        end function TensorPower
  
***************
*** 156,161 ****
--- 469,477 ----
           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 ****
--- 496,587 ----
             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%eps3(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)
+ 
+          elseif (Power_Name.eq.'power_dbi') 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%eps3(in)
+             num=num+1 
+             Keys(num) = 'gam'
+             Vals(num) = P%gam(in)
+             num=num+1 
+             Keys(num) = 'del1'
+             Vals(num) = P%del1(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
***************
*** 195,200 ****
--- 597,609 ----
             InitPower%k_0_tensor = Ini_Read_Double_File(Ini,'pivot_tensor',InitPower%k_0_tensor) 
             InitPower%nn = Ini_Read_Int('initial_power_num')
             if (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
+ !end wig
             InitPower%rat(:) = 1
             do i=1, InitPower%nn
  
***************
*** 208,215 ****
--- 617,676 ----
  
                InitPower%ScalarPowerAmp(i) = Ini_Read_Double_Array_File(Ini,'scalar_amp',i,1._dl) 
                !Always need this as may want to set tensor amplitude even if scalars not computed
+ 
+ !wig
+                  InitPower%logomega(i) = Ini_Read_Double_Array_File(Ini,'osc_logomega',i,2.d0)
+                  InitPower%amposc(i) = Ini_Read_Double_Array_File(Ini,'osc_amp',i,0.d0)             
+                  InitPower%omega(i) = 10**(InitPower%logomega(i))
+                  InitPower%phase(i) = Ini_Read_Double_Array_File(Ini,'osc_phase',i,0.d0)
+                  InitPower%logpowosc(i) = Ini_Read_Double_Array_File(Ini,'osc_logpower',i,0.d0)
+                  InitPower%powosc(i) = 10**(InitPower%logpowosc(i))
+ !end wig
+ 
             end do
            
+ !wig
+               InitPower%k_0_osc = Ini_Read_Double('osc_pivot',5d-2)
+ 
+            elseif ((Power_Name.eq.'power_sr').or.(Power_Name.eq.'power_dbi')) then
+            
+               do i=1, InitPower%nn
+               
+ !                 if (P%WantScalars .or. P%WantTransfer) then
+ 
+                  InitPower%eps1(i)= Ini_Read_Double_Array_File(Ini,'eps1',i,0.01d0)
+                  InitPower%eps2(i)=Ini_Read_Double_Array_File(Ini,'eps2',i,0.01d0)
+                  InitPower%eps3(i)=Ini_Read_Double_Array_File(Ini,'eps3',i,0.00d0)
+                  InitPower%gam(i)=Ini_Read_Double_Array_File(Ini,'gam',i,1.00d0)
+                  InitPower%del1(i)=Ini_Read_Double_Array_File(Ini,'del1',i,0.00d0)
+                  InitPower%logeps1invsigma0(i) &
+                       = Ini_Read_Double_Array_File(Ini,'logeps1invsigma0',i,1d0)
+                  InitPower%loginvsigma0(i) = InitPower%logeps1invsigma0(i) &
+                       - log10(InitPower%eps1(i))
+                  InitPower%invsigma0(i) = 10**(InitPower%loginvsigma0(i))
+                  InitPower%Xsigma0(i) = Ini_Read_Double_Array_File(Ini,'Xsigma0',i,0d0)
+                  InitPower%varphi(i) = Ini_Read_Double_Array_file(Ini,'varphi',i,0.d0)
+ 
+ !                 endif
+               
+               enddo
+            
+               InitPower%k_star= Ini_Read_Double('k_star',0.05d0)
+               InitPower%P_scalar= 1e-10*exp(Ini_Read_Double('ln[10^{10} P_scalar]',1d0))
+ 
+            else
+ 
+               stop 'No such initial power spectrum'
+ 
+            endif
+  
+ !end wig
+ 
+ 
+               
+ 
+ 
+           
         end  subroutine InitialPower_ReadParams 
  
  
diff -r -c -b -B cosmomc/distparams.ini cosmomc_wiggles/distparams.ini
*** cosmomc/distparams.ini	2010-05-04 18:31:19.000000000 +0200
--- cosmomc_wiggles/distparams.ini	2010-05-27 15:13:59.873642961 +0200
***************
*** 49,54 ****
--- 49,55 ----
  
  plot_meanlikes = T
  shade_meanlikes = T
+ td_meanlikes = F
  
  # if non-zero, output _thin file, thinned by thin_factor
  thin_factor = 0
***************
*** 95,101 ****
  
  #Parameters to use. If zero use all parameters which have lables.
  plotparams_num = 0
! plotparams = omegabh2 omegadmh2 tau ns
  
  #Get set label to empty to not include a parameter in the parameter_names file
  #lab[asz]=
--- 96,102 ----
  
  #Parameters to use. If zero use all parameters which have lables.
  plotparams_num = 0
! plotparams = omegabh2 omegadmh2 tau ns eps1 eps2 eps3 logomega amp psi logsig
  
  #Get set label to empty to not include a parameter in the parameter_names file
  #lab[asz]=
diff -r -c -b -B cosmomc/params_CMB.paramnames cosmomc_wiggles/params_CMB.paramnames
*** cosmomc/params_CMB.paramnames	2009-09-17 17:47:04.000000000 +0200
--- cosmomc_wiggles/params_CMB.paramnames	2010-06-28 16:33:53.667179421 +0200
***************
*** 5,14 ****
  omegak        \Omega_K
  fnu           f_\nu		#neutrino energy density as fraction of omegadmh2
  w             w			#constant equation of state parameter for scalar field dark energy 
! ns            n_s		#beware that pivot scale can change in .ini file
! nt            n_t
! nrun          n_{run}
! logA          log[10^{10} A_s]
  r             r			#ratio of tensor to scalar primordial amplitudes at pivot scale
  asz           A_{SZ}		#SZ template amplitude, as in WMAP
  omegal*       \Omega_\Lambda
--- 5,19 ----
  omegak        \Omega_K
  fnu           f_\nu		#neutrino energy density as fraction of omegadmh2
  w             w			#constant equation of state parameter for scalar field dark energy 
! eps1	      \epsilon_1-\delta_1	#first slow roll parameter - first sound flow param
! eps2	      \epsilon_2+\delta_1	#second slow roll parameter + first sound flow param
! eps3	      \epsilon_3		#third slow roll parameter
! unused	      log(\epsilon_1)		#unused
! logfreq	      log(\epsilon_1/\sigma_0)	#superimposed oscillation frequency
! amposc	      x\sigma_0			#superimposed oscillation amplitude
! phase	      \psi			#superimposed oscillation phase shift
! logpow	      log(p)			#superimposed power law dependence
! lnA           log[10^{10} P_*]
  r             r			#ratio of tensor to scalar primordial amplitudes at pivot scale
  asz           A_{SZ}		#SZ template amplitude, as in WMAP
  omegal*       \Omega_\Lambda
Only in cosmomc_wiggles: params_CMB.paramnames.powerlaw
diff -r -c -b -B cosmomc/params.ini cosmomc_wiggles/params.ini
*** cosmomc/params.ini	2010-05-14 17:26:20.000000000 +0200
--- cosmomc_wiggles/params.ini	2010-05-27 15:15:42.213060570 +0200
***************
*** 67,84 ****
  
  #use fast-slow parameter distinctions to speed up 
  #(note for basic models WMAP3 code is only ~3x as fast as CAMB)
! use_fast_slow = F
  
  #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
  estimate_propose_matrix = F
  
  #Tolerance on log likelihood to use when estimating best fit point
! delta_loglike = 2
  
  #Scale of proposal relative to covariance; 2.4 is recommended by astro-ph/0405462 for Gaussians
  #If propose_matrix is much broader than the new distribution, make proportionately smaller
--- 67,84 ----
  
  #use fast-slow parameter distinctions to speed up 
  #(note for basic models WMAP3 code is only ~3x as fast as CAMB)
! use_fast_slow = T
  
  #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
  estimate_propose_matrix = F
  
  #Tolerance on log likelihood to use when estimating best fit point
! delta_loglike = 0.1
  
  #Scale of proposal relative to covariance; 2.4 is recommended by astro-ph/0405462 for Gaussians
  #If propose_matrix is much broader than the new distribution, make proportionately smaller
***************
*** 135,145 ****
  #If true, generate checkpoint files and terminated runs can be restarted using exactly the same command
  #and chains continued from where they stopped
  #With checkpoint=T note you must delete all chains/file_root.* files if you want new chains with an old file_root
! checkpoint = F
  
  #CAMB parameters
  #If we are including tensors
! compute_tensors = F
  #Initial power spectrum amplitude point (Mpc^{-1})
  pivot_k = 0.05
  #If using tensors, enforce n_T = -A_T/(8A_s)
--- 135,145 ----
  #If true, generate checkpoint files and terminated runs can be restarted using exactly the same command
  #and chains continued from where they stopped
  #With checkpoint=T note you must delete all chains/file_root.* files if you want new chains with an old file_root
! checkpoint = T
  
  #CAMB parameters
  #If we are including tensors
! compute_tensors = T
  #Initial power spectrum amplitude point (Mpc^{-1})
  pivot_k = 0.05
  #If using tensors, enforce n_T = -A_T/(8A_s)
***************
*** 152,158 ****
  CMB_lensing = T
  #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
--- 152,163 ----
  CMB_lensing = T
  #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 = 1
! k_src_sample_level = 1.
! k_int_sample_level = 1.
! camb_compute_every_l = T
  
  #If action = 1
  redo_likelihoods = T
***************
*** 178,189 ****
  param[fnu] = 0 0 0 0 0
  param[w] = -1 -1 -1 0 0
  
! param[ns] = 0.95 0.5 1.5 0.02 0.01
! param[nt] = 0 0 0 0 0
! param[nrun] = 0 0 0 0 0
! 
! #log[10^10 A_s]
! param[logA] = 3 2.7 4 0.01 0.01
! param[r] = 0 0 0 0 0
  #SZ amplitude, as in WMAP analysis
  param[asz]= 1 0 2 0.4 0.4
--- 183,208 ----
  param[fnu] = 0 0 0 0 0
  param[w] = -1 -1 -1 0 0
  
! ### power law ###
! 
! #param[ns] = 0.95 0.5 1.5 0.02 0.01
! #param[nt] = 0 0 0 0 0
! #param[nrun] = 0 0 0 0 0
! #param[unused] = 0 0 0 0 0
! ### or slow roll (default) ###
! 
! param[eps1] = 0.1 0.00001 0.2 0.01 0.02
! param[eps2] = -0.04 -0.2 0.2  0.01 0.02
! param[eps3] = -0.04 -0.2 0.2  0.01 0.02
! param[unused] = 0 0 0 0 0
! 
! param[logfreq] = 1.5 1 2.5 0.5 0.2
! param[amposc] = 0.05 0 0.45 0.1 0.2
! param[phase] = 0 0 6.28 0.4 0.1
! param[logpow] = 0 0 0 0 0
! 
! #ln[10^10 A_s]
! param[lnA] = 3 2.7 4 0.01 0.01
! param[r] = 1 1 1 0 0
  #SZ amplitude, as in WMAP analysis
  param[asz]= 1 0 2 0.4 0.4
diff -r -c -b -B cosmomc/source/CMB_Cls_simple.f90 cosmomc_wiggles/source/CMB_Cls_simple.f90
*** cosmomc/source/CMB_Cls_simple.f90	2010-05-20 12:13:02.000000000 +0200
--- cosmomc_wiggles/source/CMB_Cls_simple.f90	2010-05-27 15:16:04.453017704 +0200
***************
*** 393,398 ****
--- 393,420 ----
     integer zix
     real redshifts(matter_power_lnzsteps)
  
+ !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)
***************
*** 424,432 ****
              P%Transfer%high_precision=.true.
              P%Transfer%kmax=P%Transfer%kmax + 0.2
             end if
!           AccuracyBoost = AccuracyLevel
!           lAccuracyBoost = AccuracyLevel
!           lSampleBoost = AccuracyLevel
            P%AccurateReionization = .true.
          end if
                  
--- 446,456 ----
              P%Transfer%high_precision=.true.
              P%Transfer%kmax=P%Transfer%kmax + 0.2
             end if
! !wig
! !          AccuracyBoost = AccuracyLevel
! !          lAccuracyBoost = AccuracyLevel
! !          lSampleBoost = AccuracyLevel
! !end wig
            P%AccurateReionization = .true.
          end if
                  
diff -r -c -b -B cosmomc/source/cmbtypes.f90 cosmomc_wiggles/source/cmbtypes.f90
*** cosmomc/source/cmbtypes.f90	2010-05-10 17:28:00.000000000 +0200
--- cosmomc_wiggles/source/cmbtypes.f90	2010-05-27 15:16:25.173754728 +0200
***************
*** 6,12 ****
  
  
  !Number of CMB Cls, 1 for just temperature, 3 (4) for polarization (with B)
!   integer, parameter  :: num_cls  = 3
  
    integer, parameter  :: num_cls_ext=0
     !number of other C_l
--- 6,14 ----
  
  
  !Number of CMB Cls, 1 for just temperature, 3 (4) for polarization (with B)
! !wig
!   integer, parameter  :: num_cls  = 4
! !end wig
  
    integer, parameter  :: num_cls_ext=0
     !number of other C_l
***************
*** 88,94 ****
    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
--- 90,99 ----
    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	2010-05-20 12:13:02.000000000 +0200
--- cosmomc_wiggles/source/driver.F90	2010-05-27 15:16:30.794267702 +0200
***************
*** 57,67 ****
  
  #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')
  
--- 57,76 ----
  
  #ifdef MPI 
  
! !wig
! !        if (instance /= 0) call DoStop('With MPI should not have second parameter')
!         instance_shift=instance
!         write(*,*)'MPI file names shifted by: ',instance_shift
! !end wig
          call mpi_comm_rank(mpi_comm_world,MPIrank,ierror)
  
          instance = MPIrank +1 !start at 1 for chains                                             
! !wig
! !        write (numstr,*) instance
! 	instance_shift = MPIrank + instance_shift
! 	write (numstr,*) instance_shift
! !end wig
! 
          rand_inst = instance
          if (ierror/=MPI_SUCCESS) call DoStop('MPI fail')
  
***************
*** 84,89 ****
--- 93,105 ----
  
          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
  
          checkpoint = Ini_Read_Logical('checkpoint',.false.)
          if (checkpoint) flush_write = .true.
diff -r -c -b -B cosmomc/source/GetDist.f90 cosmomc_wiggles/source/GetDist.f90
*** cosmomc/source/GetDist.f90	2010-05-19 13:23:38.000000000 +0200
--- cosmomc_wiggles/source/GetDist.f90	2010-05-27 16:26:57.773643567 +0200
***************
*** 102,107 ****
--- 102,133 ----
          integer Num_ComparePlots
          logical ::   prob_label = .false.      
              
+ !addon
+ !exact field model
+         logical, parameter :: findMaxValues = .false.
+         logical, parameter :: mapRhoReh = .false.
+         logical, parameter :: mapLFPowerToWreh = .true.
+ 
+ !slow-roll analysis
+         logical, parameter :: adjEps1ToLog = .false.
+         logical, parameter :: extractHinf = .false.
+         logical, parameter :: mapnslogr = .true.
+         logical, parameter :: use2ndSR=.true.        
+         real, parameter :: ln10p5HinfCut=-10
+         real, parameter :: logrCut=-3.7
+ 
+ #ifdef SRREHEAT
+         logical, parameter :: doReheatFromSr = .false.
+         character(len=*), parameter :: reheatModel='smallf'
+ #endif
+ 
+         real, parameter :: wfix = -0.3
+         real, parameter :: pmax = 15, pmin = 0
+         real, parameter :: logmumax=10
+ !end addon
+ 
+ 
+     
  contains
  
    subroutine AdjustPriors
***************
*** 109,116 ****
     !Be careful as this code is parameterisation dependent
     integer i
     real ombh2, chisq 
   
!    stop 'You need to write the AdjustPriors subroutine in GetDist.f90 first!'
    
  !   write (*,*) 'Adjusting priors'
  !   do i=0, nrows-1
--- 135,154 ----
     !Be careful as this code is parameterisation dependent
     integer i
     real ombh2, chisq 
+ !addon
+    integer, parameter :: dp = kind(1._8)
+    real(dp) :: eps1,logeps1,eps2,eps3, lnR,lnRrad
+    real(dp) :: wreh,p,logmu,lnRhoReh,lnRhoEnd   
+ 
+    real :: ln10p5HinfMin=1e30
+    real :: ln10p5HinfMax=-1e30
+    
+    real :: lnRhoRehmax = -1e30, lnRhoEndMax=-1e30, lnRmax = -1e30, lnRradmax = -1e30
+ 
+    real, parameter :: lnRhoNuc = -187.747
   
! !   stop 'You need to write the AdjustPriors subroutine in GetDist.f90 first!'
! !end addon 
    
  !   write (*,*) 'Adjusting priors'
  !   do i=0, nrows-1
***************
*** 122,127 ****
--- 160,373 ----
  !    
  !   end do
  
+ !addon
+    wreh = wfix
+ 
+    
+    if (findMaxValues) then
+       print *,'-->computing maximum values in chains'
+ 
+       do i=0,nrows-1
+          lnR = coldata(14+2,i)
+          lnRhoEnd = coldata(2+30,i)
+          lnRrad = coldata(2+28,i)
+          lnRhoEndMax = max(lnRHoEndMax,lnRhoEnd)
+          lnRmax = max(lnRmax, lnR)
+          lnRradmax = max(lnRradmax, lnRrad)
+       enddo
+ 
+       print *,'lnRhoEndMax= ',lnRhoEndMax
+       print *,'lnRradMax= lnRmax= ',lnRradMax, lnRmax
+    endif
+ 
+    if (mapRhoReh) then
+       print *,'--> rejecting all p <',pmin
+       print *,'--> rejecting all logmu >',logmumax
+       print *,'--> rejecting all lnRhoNuc > lnRhoReh > lnRhoEnd'
+ 
+       do i=0,nrows-1
+          p = coldata(2+10,i)
+          logmu = coldata(11+2,i)
+ 
+          if (mapLFPowerToWreh) wreh = (p-2._dp)/(p+2._dp)
+ 
+          lnR = coldata(14+2,i)
+          lnRrad = coldata(2+28,i)
+          lnRhoReh = coldata(27+2,i)
+          lnRhoEnd = coldata(2+30,i)
+          
+          if ((wreh.lt.-1._dp/3._dp) &
+               .or.(p.lt.pmin).or.(logmu.gt.logmumax) &
+               .or.(lnRhoReh.gt.lnRhoEnd) &
+               .or.(lnRhoReh.lt.lnRhoNuc)) then
+             chisq = 1e30
+             coldata(1,i) = 0.
+             coldata(2,i) = coldata(2,i) + chisq/2
+          else
+             lnRhoRehMax = max(lnRHoRehMax,lnRhoReh)                          
+             lnRhoEndMax = max(lnRHoEndMax,lnRhoEnd)
+             lnRmax = max(lnRmax, lnR)
+             lnRradmax = max(lnRradmax, lnRrad)
+          endif
+ 
+       enddo
+       print *,'Max values after applying prior!'
+       print *,'lnRhoEndMax= ',lnRhoEndMax
+       print *,'lnRhoRehMax= ',lnRhoRehMax      
+       print *,'lnRradMax= lnRmax= ',lnRradMax, lnRmax
+    endif
+       
+    if (adjEps1ToLog) then
+       print *,'----> assumed mapped log(eps1) in col 8'
+       print *,'----> from flat eps1 to flat log(eps1)'
+       read(*,*)
+       do i=0, nrows-1
+          
+          logeps1 = coldata(8+2,i)
+          eps1 = 10.**logeps1
+          
+          coldata(1,i) = coldata(1,i)/(eps1*log(10.))
+          coldata(2,i) = coldata(2,i) 
+ 
+       end do
+    endif
+ 
+ 
+ 
+    if (extractHinf) then
+ 
+       write (*,*)'Computing and adjusting  bounds for ln(10^5 Hinf/mpl) in 17!'
+       
+       do i=0,nrows-1
+          ln10p5HinfMax = max(ln10p5HinfMax,coldata(2+17,i))
+          ln10p5HinfMin = min(ln10p5HinfMin,coldata(2+17,i))         
+ 
+          if (coldata(2+17,i).lt.ln10p5HinfCut) then
+             chisq = 1e30
+             coldata(1,i) = 0.
+             coldata(2,i) = coldata(2,i) + chisq/2
+          endif
+ 
+       enddo
+                
+       write(*,*)'----> ln(10^5 Hinf) Min= Max= ',ln10p5HinfMin,ln10p5HinfMax
+ 
+    endif
+ 
+    if (mapnslogr) then
+ 
+       write (*,*)'Adjusting  bounds for log(r) in 14!'
+ 
+       do i=0,nrows-1
+          
+          if (coldata(2+14,i).lt.logrCut) then
+             chisq = 1e30
+             coldata(1,i) = 0.
+             coldata(2,i) = coldata(2,i) + chisq/2
+          endif
+ 
+       enddo
+    endif
+ 
+ #ifdef SRREHEAT
+    if (doReheatFromSR) then
+       write(*,*)'Adjusting prior for reheat SR analysis'
+ 
+       select case (reheatModel)
+ 
+       case ('largef')
+         
+         
+          do i=0,nrows-1
+             
+             eps1 = 10**(coldata(2+8,i))
+             eps2 = coldata(2+9,i)
+             p = coldata(2+14,i)
+             wreh = (p-2.)/(p+2.)
+ 
+ !large field models cannot produce eps2<=0 (out of prior)
+             if ((eps2.le.0.).or.(p.gt.pmax).or.(p.lt.pmin)) then
+                chisq = 1e30
+                coldata(1,i) = 0.
+                coldata(2,i) = coldata(2,i) + chisq/2
+             endif
+ 
+ !reheating is after end of inflation
+             lnRhoReh = coldata(11+2,i)
+             lnRhoEnd = coldata(12+2,i)
+ 
+ !            print *,'lnrhoreh end=',lnRhoReh,lnRhoEnd,eps1,eps2
+ 
+             if (wreh.lt.-1./3.) then
+                chisq = 1e30
+                coldata(1,i) = 0.
+                coldata(2,i) = coldata(2,i) + chisq/2
+             endif
+             if (lnRhoReh.gt.lnRhoEnd) then
+                chisq = 1e30
+                coldata(1,i) = 0.
+                coldata(2,i) = coldata(2,i) + chisq/2
+             endif
+             if (lnRhoReh.lt.lnRhoNuc) then
+                chisq = 1e30
+                coldata(1,i) = 0.
+                coldata(2,i) = coldata(2,i) + chisq/2            
+             endif
+          enddo
+          
+       case ('smallf')
+ 
+          wreh = wfix
+                  
+          do i=0,nrows-1
+             
+             eps1 = 10**(coldata(2+8,i))
+             eps2 = coldata(2+9,i)
+             eps3 = coldata(2+10,i)
+ 
+             p = coldata(2+14,i)            
+ 
+ 
+ !small field models (out of prior)
+             if (.not.(sr_checkobs_sf(eps1,eps2,eps3)) &
+                  .or.(p.gt.pmax).or.(p.lt.pmin)) then
+                chisq = 1e30
+                coldata(1,i) = 0.
+                coldata(2,i) = coldata(2,i) + chisq/2
+             endif
+ 
+ 
+ !reheating is after end of inflation
+             lnRhoReh = coldata(11+2,i)
+             lnRhoEnd = coldata(12+2,i)
+ 
+             if (wreh.lt.-1./3.) then
+                chisq = 1e30
+                coldata(1,i) = 0.
+                coldata(2,i) = coldata(2,i) + chisq/2
+             endif
+             if (lnRhoReh.gt.lnRhoEnd) then
+                chisq = 1e30
+                coldata(1,i) = 0.
+                coldata(2,i) = coldata(2,i) + chisq/2
+             endif
+             if (lnRhoReh.lt.lnRhoNuc) then
+                chisq = 1e30
+                coldata(1,i) = 0.
+                coldata(2,i) = coldata(2,i) + chisq/2            
+             endif
+          enddo
+          
+       end select
+ 
+       
+    endif
+ #endif
+ !end addon
+ 
+ 
+ 
+ 
    end subroutine AdjustPriors
  
    subroutine MapParameters(invars)
***************
*** 129,135 ****
  ! map parameters in invars: eg. invars(3)=invars(3)*invars(4) 
  
  !    invars(2+13)=invars(17+2)*exp(-invars(2+4)) 
!     stop 'Need to write MapParameters routine first'  
    
    end subroutine MapParameters
  
--- 375,568 ----
  ! map parameters in invars: eg. invars(3)=invars(3)*invars(4) 
  
  !    invars(2+13)=invars(17+2)*exp(-invars(2+4)) 
! 
! 
! !addon
!     integer, save :: counter = 0
!     real :: lnR, lnRhoReh,lnRhoEnd, ns,logr
!     real :: cseps1,eps1mdel1,eps2pdel1,del1
!     real :: as0,ln10p5Hinf
!     real, parameter :: C_const = -0.7296
! 
!     integer, parameter :: dp=kind(1._8)
!     real(dp) :: p, mu, wreh, Pstar,bfold,eps1,eps2,eps3
!     real(dp), dimension(2) :: buffer2D
!     real(dp), dimension(3) :: buffer3D
!         
! 
!     counter = counter + 1
! 
!     if (mapRhoReh) then
! 
!        if (mapLFPowerToWreh) then
!           
!           p = invars(2+10)
!           wreh = (p-2._dp)/(p+2._dp)
!           
!           if (counter.eq.1) print *,'reheat marginalisation on p in col10'
!        else
!           
!           wreh = wfix
!           if (counter.eq.1) print *,'marginalising at fixed wreh= ',wreh
!           
!        endif
!        
!        if (wreh.eq.1./3.) stop 'Treh inversion is singular!'
!        
!        if (counter.eq.1) print *,'mapping lnR -> lnRhoReh in col27'
! 
!        lnR = invars(2+14)
!        
!        lnRhoEnd = invars(2+30)
! 
!        lnRhoReh = 4._dp*((1._dp+wreh)/(1._dp/3._dp-wreh))*lnR &
!             - 4._dp*(1._dp/3._dp + wreh)/(1._dp/3._dp - wreh)*lnRhoEnd/2._dp
!        
!        invars(27+2)=lnRhoReh
! 
!     end if
!     
! 
!     if (extractHinf) then
!        del1 = 0.
! 
!        if (counter.eq.1) print *,'assumiing 10**(col10) is eps1'
! 
!        eps1mdel1 = 10**invars(2+8)
!        eps2pdel1 = invars(2+9)
!        eps1 = eps1mdel1
!        eps2 = eps2pdel1
!        eps3 = invars(2+10)
! 
! 
!        cseps1=eps1mdel1
! !       cseps1 = eps1/invars(2+15)
! !       cseps1= 10.**invars(2+15)
! 
! !       eps1 = 10.**invars(2+8)
! !       eps2 = invars(2+9)
! !       eps3 = invars(2+10)
! 
! 
! 
!        as0 = 1. + eps1mdel1*(-2.*(C_const+1.)) - eps2pdel1*C_const
! 
!        if (use2ndSR) then
!           if (counter.eq.1) print *,'using 2nd order slow-roll!'
!           as0 = as0 + eps1*eps1 * (2.*C_const**2 + 2.*C_const + pi**2/2. - 5.) &
!                + eps1*eps2 * (C_const**2 - C_const + 7.*pi**2/12. - 7.) &
!                + eps2*eps2 * (0.5*C_const**2 + pi**2/8. - 1.) &
!                + eps2*eps3 * (-0.5*C_const**2 + pi**2/24.)
!        endif
!        
!        ln10p5Hinf = + 0.5*invars(2+16) &
!             + 0.5*log(3.141592653589793238*cseps1) - 0.5*log(as0)
!        
! 
!        invars(2+17) = ln10p5Hinf
!     endif
! 
!     if (mapnslogr) then
!        
!        if (counter.eq.1) then          
!           print *,'assumiing 10**(col10) is eps1'
!           print *,'Mapping ns --> col15 & logr --> col14'
!        endif
! 
!        eps1 = 10**invars(2+8)
!        eps2 = invars(2+9)
!        eps3 = invars(2+10)
! 
!        ns = 1 - 2.*eps1 - eps2
!        logr = log10(16*eps1)
! 
!        if (use2ndSR) then
!           if (counter.eq.1) print *,'using 2nd order slow-roll!'
!           ns = ns - 2*eps1**2 - (2*C_const + 3)*eps1*eps2 - C_const*eps2*eps3
!           logr = logr + log10(1 + C_const*eps2)
!        endif
! 
!        invars(2+14) = logr
!        invars(2+15) = ns
!       
!     endif
!        
! 
! #ifdef SRREHEAT
!     if (doReheatFromSR) then
! 
!        if (counter.eq.1) then
!           write(*,*)'************************************************'
!           write(*,*)'Model is: ',reheatModel
!           write(*,*)'Doing Reheating analysis from SR:'
!           write(*,*)'Mapping lnRhoReh -> col11 & lnRhoEnd -> col12'
!           write(*,*)'& -N* -> col13 & p -> col14 & mu -> col15'
!           write(*,*)'************************************************'
!        endif
!           
!        Pstar = exp(invars(2+16))*1.e-10
!              
!        eps1 = 10**(invars(2+8))
!        eps2 = invars(2+9)
!        eps3 = invars(2+10)
!        
!        lnRhoReh = 0.
!        lnRhoEnd = 0.
!        bfold=0.
! 
!        select case(reheatModel)
!           
!        case ('largef')
!                                                    
! !large field models cannot produce eps2<=0 (out of prior)
!           if (eps2.le.0.) return
! 
!           buffer2d = sr_matter_obs_lf(eps1,eps2)
! 
!           p = buffer2d(1)
!           invars(15+2)=p
! 
!           if ((p.gt.pmax).or.(p.lt.pmin)) return
! 
!           lnRhoReh = sr_lnrhoreh_lf(wreh,eps1,eps2,Pstar,bfold)
!           lnRhoEnd = sr_lnrhoend_lf(p,Pstar)
! 
!           invars(11+2)=lnRhoReh
!           invars(12+2)=lnRhoEnd
!           invars(13+2)=bfold
!          
!        case ('smallf')
! 
! !small field models cannot produce any eps1,eps2,eps3<=0 (out of prior)
!           if (.not.sr_checkobs_sf(eps1,eps2,eps3)) return
!           
!           buffer3d = sr_matterovermu_obs_sf(eps1,eps2,eps3)
! 
!           p = buffer3d(1)
!           mu = buffer3d(2)
!           wreh = wfix
! 
!           invars(14+2)=p
!           invars(15+2)=mu
! 
!           if ((p.gt.pmax).or.(p.lt.pmin)) return
! 
!           lnRhoReh = sr_lnrhoreh_sf(wreh,eps1,eps2,eps3,Pstar,bfold)
!           lnRhoEnd = sr_lnrhoend_sf(p,mu,Pstar)
! 
!           invars(11+2) = lnRhoReh
!           invars(12+2) = lnRhoEnd
!           invars(13+2) = bfold
! 
!        case default
! 
!           stop 'model not implemented!'
! 
!        end select
! 
!     endif
! #endif
! !end addon
    
    end subroutine MapParameters
  
diff -r -c -b -B cosmomc/source/Makefile cosmomc_wiggles/source/Makefile
*** cosmomc/source/Makefile	2010-05-04 11:07:37.000000000 +0200
--- cosmomc_wiggles/source/Makefile	2010-08-10 15:27:04.192288768 +0200
***************
*** 1,118 ****
! #You may need to edit the library paths for MKL for Intel
! #Beware of using optmizations that lose accuracy - may give errors when running
  
! ##Uncomment the next line to include dr7 LRG
! EXTDATA = 
! #EXTDATA = LRG
  
! #set WMAP empty not to compile with WMAP
! WMAP = /home/aml1005/WMAP7/likelihood_v4
  
! #Only needed for WMAP
! cfitsio = /usr/local/cfitsio/intel10/64/3.040
  
- #GSL only needed for DR7 LRG
- GSLPATH = /home/aml1005/libs/gsl
  
- IFLAG = -I
- INCLUDE=
- 
- #Intel MPI (assuming mpif90 set to point to ifort)
- #these settings for ifort 11.1 and higher; may need to add explicit link directory otherwise
- F90C     = mpif90
- FFLAGS = -O2 -ip -W0 -WB -openmp -fpp -DMPI -vec_report0 -mkl=parallel
- LAPACKL = -lmkl_lapack 
- 
- #HPCF settings. Use Inteal 9 or 10.1+, not 10.0
- #F90C     = mpif90
- #FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp -DMPI -vec_report0
- #LAPACKL = -L/usr/local/Cluster-Apps/intel/mkl/10.2.2.025/lib/em64t -lmkl_lapack -lmkl -lguide -lpthread
- #GSLPATH = /usr/local/Cluster-Apps/gsl/1.9
- #cfitsio = /usr/local/Cluster-Users/cpac/cmb/2.1.0/cfitsio
- #INCLUDE=
- 
- #COSMOS: use "module load cosmolib latest" 
- #use "runCosmomc" (globally installed) to run, defining required memory usage
- ifeq ($(COSMOHOST),cosmos)
- F90C = ifort
- FFLAGS = -openmp -O3 -w -fpp2 -DMPI
- LAPACKL = -mkl=sequential -lmkl_lapack -lmpi
- cfitsio = $(CFITSIO)
- WMAP = $(COSMOLIB)/WMAP7
- GSLPATH = $(GSL_ROOT)
- endif
- 
- #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
- 
- #G95; make sure LAPACK and MPI libs also compiled with g95
- #F90C     = mpif90
- #FFLAGS = -O2 -DMPI
- #LAPACKL = /LAPACK/lapack_LINUX.a /LAPACK/blas_LINUX.a
- 
- #GFortran: if pre v4.3 add -D__GFORTRAN__
- #F90C     = gfortran
- #FFLAGS = -O2 -ffree-form -x f95-cpp-input 
- #LAPACKL = -Wl,-framework -Wl,accelerate
- 
- #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
- #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
- 
- #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 
- #IFLAG = -p
- 
- #NAGF95, single processor:
- #F90C     = f95
- #FFLAGS = -DNAGF95 -O3
- #LAPACKL = -llapack -lblas -lg2c
- 
- #PGF90
- #F90C = pgf90
- #FFLAGS = -O2 -DESCAPEBACKSLASH
- #LAPACKL = -llapack -lblas
- 
- 
- #Sun, single processor:
- #F90C     = f90
- #FFLAGS = -fast -ftrap=%none
- #LAPACKL = -lsunperf -lfsu
- #LAPACKL = -lsunperf -lfsu -lsocket -lm
- #IFLAG = -M
- 
- #Sun MPI
- #F90C = mpf90
- #FFLAGS =  -O4 -openmp -ftrap=%none -dalign -DMPI
- #LAPACKL = -lsunperf -lfsu -lmpi_mt
- #IFLAG = -M
- 
- #Sun parallel enterprise:
- #F90C     = f95
- #FFLAGS =  -O4 -xarch=native64 -openmp -ftrap=%none
- #LAPACKL = -lsunperf -lfsu
- #IFLAG = -M
- 
- 
- #IBM XL Fortran, multi-processor (run "module load lapack" then run "gmake")
- # See also http://cosmocoffee.info/viewtopic.php?t=326
- #F90C     = xlf90_r $(LAPACK)
- #FFLAGS  = -WF,-DIBMXL -qsmp=omp -qsuffix=f=f90:cpp=F90 -O3 -qstrict -qarch=pwr3 -qtune=pwr3
- #INCLUDE = -lessl
- #LAPACKL =
  
  PROPOSE = propose.o
  CLSFILE = CMB_Cls_simple.o
--- 1,84 ----
! # >>> DESIGNED FOR GMAKE <<<
  
! # Unified Systems makefile for COSMOMC
! # Add FLAGS -DMPI for using MPI
! 
! ext=$(shell uname | cut -c1-3)
! 
! CC=cc
! 
! ifeq ($(ext),IRI)
! F90C= f90
! FFLAGS= -Ofast -mp -n32 -LANG:recursive=ON -lmpi -DMPI
! WMAPFLAGS= $(FFLAGS)
! LAPACKL = -lcomplib.sgimath_mp
! INCLUDE = -I../camb
! CFITSIODIR = /people/ringeval/usr
! GSLDIR =
! endif
! 
! ifeq ($(ext),Lin)
! F90C=gfortran
! FFLAGS= -O -cpp -fopenmp
! WMAPFLAGS= -O
! LAPACKL = -llapack
! INCLUDE = -I../camb 
! CFITSIODIR = /opt/gnu/cfitsio
! GSLDIR = /usr
! endif 
  
! ifeq ($(ext),OSF)
! F90C=f90
! FFLAGS= -omp -O -arch host -math_library fast -tune host -fpe1
! WMAPFLAGS= $(FFLAGS)
! LAPACKL  = -lcxml
! INCLUDE = -I../camb 
! CFITSIODIR = 
! GSLDIR =
! endif
  
! ifeq ($(ext),Sun)
! F90C=f90
! FFLAGS= -O4 -xarch=native64 -openmp -ftrap=%none
! WMAPFLAGS= $(FFLAGS)
! LAPACKL = -lsunperf -lfsu
! INCLUDE = -I../camb -M../camb
! CFITSIODIR = 
! GSLDIR =
! endif
! 
! ifeq ($(ext),AIX)
! F90C    = mpxlf90_r
! FFLAGS  = -O4 -WF,-DIBMXL,-DMPI -qstrict -qsmp=omp -qmaxmem=-1 -qsuffix=f=f90:cpp=F90
! WMAPFLAGS= $(FFLAGS)
! LAPACKL = -lessl
! INCLUDE = -I../camb
! CFITSIODIR = 
! GSLDIR =
! endif
! 
! #EXTDATA = LRG
! EXTDATA = LRG
! EXTINCLUDE = -I$(GSLDIR)/include
! EXTOBJS = bsplinepk.o
! 
! 
! WMAPDIR = ../WMAP
! WMAPINCLUDE = -I$(CFITSIODIR)/include
! WMAPOBJS = read_archive_map.o \
! 	read_fits.o \
! 	healpix_types.o \
! 	br_mod_dist.o \
! 	WMAP_7yr_options.o \
! 	WMAP_7yr_util.o \
! 	WMAP_7yr_gibbs.o \
! 	WMAP_7yr_tt_pixlike.o \
! 	WMAP_7yr_tt_beam_ptsrc_chisq.o \
! 	WMAP_7yr_teeebb_pixlike.o \
! 	WMAP_7yr_tetbeebbeb_pixlike.o \
! 	WMAP_7yr_likelihood.o
  
  
  
  PROPOSE = propose.o
  CLSFILE = CMB_Cls_simple.o
***************
*** 120,132 ****
  #Can use params_H if you prefer more generic parameters
  PARAMETERIZATION = params_CMB.o
  
- F90FLAGS = -DMATRIX_SINGLE $(FFLAGS) $(IFLAG)../camb $(INCLUDE)
  LINKFLAGS = -L../camb -lcamb $(LAPACKL) 
  
  DISTFILES = utils.o ParamNames.o Matrix_utils.o settings.o GetDist.o
  
! OBJFILES= utils.o ParamNames.o Matrix_utils.o settings.o IO.o cmbtypes.o Planck_like.o \
! 	cmbdata.o WeakLen.o bbn.o bao.o lrggettheory.o mpk.o supernovae.o HST.o SDSSLy-a-v3.o \
  	$(CLSFILE) paramdef.o $(PROPOSE) $(PARAMETERIZATION) calclike.o \
  	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
  
--- 86,101 ----
  #Can use params_H if you prefer more generic parameters
  PARAMETERIZATION = params_CMB.o
  
  LINKFLAGS =  -L../camb -lcamb  $(LAPACKL)
  
+ F90FLAGS = -DMATRIX_SINGLE $(FFLAGS) $(INCLUDE)
+ 
  DISTFILES = utils.o ParamNames.o Matrix_utils.o settings.o GetDist.o
  
! 
! 
! OBJFILES = utils.o ParamNames.o Matrix_utils.o settings.o IO.o cmbtypes.o Planck_like.o \
! 	cmbdata.o WeakLen.o bbn.o bao.o lrggettheory.o mpk.o supernovae.o HST.o SDSSLy-a-v3.o\
  	$(CLSFILE) paramdef.o $(PROPOSE) $(PARAMETERIZATION) calclike.o \
  	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
  
***************
*** 130,165 ****
  	$(CLSFILE) paramdef.o $(PROPOSE) $(PARAMETERIZATION) calclike.o \
  	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
  
- 
  ifeq ($(EXTDATA),LRG)
  F90FLAGS += -DDR71RG
! OBJFILES += bsplinepk.o
! GSLINC = -I$(GSLPATH)/include
! LINKFLAGS += -L$(GSLPATH)/lib -lgsl -lgslcblas
! endif
! 
! ifneq ($(WMAP),)
! F90FLAGS += $(IFLAG)$(cfitsio)/include $(IFLAG)$(WMAP)
! LINKFLAGS +=  -L$(cfitsio)/lib -L$(WMAP) -lcfitsio
! 
! OBJFILES += $(WMAP)/read_archive_map.o \
! 	$(WMAP)/read_fits.o \
! 	$(WMAP)/healpix_types.o \
! 	$(WMAP)/WMAP_7yr_options.o \
! 	$(WMAP)/WMAP_7yr_util.o \
! 	$(WMAP)/WMAP_7yr_tt_pixlike.o \
! 	$(WMAP)/WMAP_7yr_teeebb_pixlike.o \
! 	$(WMAP)/WMAP_7yr_likelihood.o \
! 	$(WMAP)/WMAP_7yr_gibbs.o \
! 	$(WMAP)/WMAP_7yr_tt_beam_ptsrc_chisq.o \
! 	$(WMAP)/br_mod_dist.o 
  else
  F90FLAGS += -DNOWMAP
  endif
  
- default: cosmomc
  
! all : cosmomc getdist
  
  utils.o: ../camb/libcamb.a
  settings.o: utils.o
--- 99,124 ----
  	$(CLSFILE) paramdef.o $(PROPOSE) $(PARAMETERIZATION) calclike.o \
  	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
  
  ifeq ($(EXTDATA),LRG)	
  F90FLAGS += -DDR71RG
! EXTINCLUDE = -I$(GSLDIR)/include
! LINKFLAGS += -lgsl -lgslcblas
! OBJFILES += $(EXTOBJS)
! endif	
! 
! ifneq ($(WMAPDIR),)
! F90FLAGS += $(WMAPINCLUDE)
! LINKFLAGS += -L$(CFITSIODIR)/lib -lcfitsio
! OBJFILES += $(WMAPOBJS)
  else
  F90FLAGS += -DNOWMAP
  endif
  
  
! 
! default: cosmomc.$(ext)
! 
! all : cosmomc.$(ext) getdist.$(ext)
  
  utils.o: ../camb/libcamb.a
  settings.o: utils.o
***************
*** 184,197 ****
  MCMC.o: calclike.o
  driver.o: MCMC.o
  
- export FFLAGS
- export F90C
- 
  .f.o:
  	f77 $(F90FLAGS) -c $<
  
  %.o: %.c
! 	$(CC) $(GSLINC) -c $*.c
  
  %.o: %.f90
  	$(F90C) $(F90FLAGS) -c $*.f90
--- 143,160 ----
  MCMC.o: calclike.o
  driver.o: MCMC.o
  
  .f.o:
  	f77 $(F90FLAGS) -c $<
  
  %.o: %.c
! 	$(CC) $(EXTINCLUDE) -c $*.c
! 
! %.o: $(WMAPDIR)/%.f90
! 	$(F90C) $(WMAPFLAGS) $(WMAPINCLUDE) -c $<
! 
! %.o: $(WMAPDIR)/%.F90
! 	$(F90C) $(WMAPFLAGS) $(WMAPINCLUDE) -c $<
! 
  
  %.o: %.f90
  	$(F90C) $(F90FLAGS) -c $*.f90
***************
*** 200,218 ****
  	$(F90C) $(F90FLAGS) -c $*.F90
  
  
! cosmomc: camb $(OBJFILES)
! 	$(F90C) -o ../cosmomc $(OBJFILES) $(LINKFLAGS) $(F90FLAGS)
  
  
! clean: cleancosmomc
! 	rm -f ../camb/*.o ../camb/*.obj ../camb/*.mod
! 
! cleancosmomc:
  	rm -f *.o *.mod *.d *.pc *.obj ../core 
  
! 
! getdist: camb $(DISTFILES)
! 	$(F90C) -o ../getdist $(DISTFILES) $(LINKFLAGS) $(F90FLAGS) 
! 
! camb:
! 	cd ../camb && $(MAKE) --file=Makefile_main libcamb.a
--- 163,174 ----
  	$(F90C) $(F90FLAGS) -c $*.F90
  
  
! cosmomc.$(ext): $(OBJFILES) ../camb/libcamb.a 	
! 	$(F90C) -o ../$@ $(OBJFILES) $(LINKFLAGS) $(F90FLAGS)
  
  
! clean:
  	rm -f *.o *.mod *.d *.pc *.obj ../core
  
! getdist.$(ext): $(DISTFILES)
! 	$(F90C) -o ../$@ $(DISTFILES) $(LINKFLAGS) $(F90FLAGS) 
diff -r -c -b -B cosmomc/source/MCMC.f90 cosmomc_wiggles/source/MCMC.f90
*** cosmomc/source/MCMC.f90	2010-05-10 17:28:00.000000000 +0200
--- cosmomc_wiggles/source/MCMC.f90	2010-05-27 15:18:17.855515009 +0200
***************
*** 215,221 ****
    
      Like = GetLogLike(grid(r))
     
!     if (Feedback > 1) write (*,*) r, 'Likelihood: ', Like, 'Current Like:', CurLike
  
      if ((Like /= logZero) .and. (CurLike > Like .or. randexp1() > Like - CurLike)) then
      !Accept
--- 215,221 ----
    
      Like = GetLogLike(grid(r))
     
!     if (Feedback > 1) write (*,*) r, 'Likelihood: ', Like, 'Current Like:', CurLike, 'Rank', MPIRank
  
      if ((Like /= logZero) .and. (CurLike > Like .or. randexp1() > Like - CurLike)) then
      !Accept
***************
*** 501,507 ****
             output_lines = output_lines +1
             call WriteParams(CurParams, real(mult), CurLike)
            end if
!           if (Feedback > 1) write (*,*) instance, 'Slicing, Current Like:', CurLike
            mult = 1
            if (num_slow /=0) call SliceSampleSlowParam(CurParams, CurLike) 
            if (num_fast /=0) call SliceSampleFastParams(CurParams, CurLike) 
--- 501,507 ----
             output_lines = output_lines +1
             call WriteParams(CurParams, real(mult), CurLike)
            end if
!           if (Feedback > 1) write (*,*) instance, 'Slicing, Current Like:', CurLike, 'Rank', MPIRank
            mult = 1
            if (num_slow /=0) call SliceSampleSlowParam(CurParams, CurLike) 
            if (num_fast /=0) call SliceSampleFastParams(CurParams, CurLike) 
***************
*** 536,542 ****
  
      Like = GetLogLike(Trial) 
  
!     if (Feedback > 1) write (*,*) 'Likelihood: ', Like, 'Current Like:', CurLike
  
  !!!    accpt = (Like /= logZero) .and. (CurLike > Like .or. randexp1() > Like - CurLike) 
               !Include the min() so that compilers not doing optimal compilation don't complain
--- 536,542 ----
  
      Like = GetLogLike(Trial) 
  
!     if (Feedback > 1) write (*,*) 'Likelihood: ', Like, 'Current Like:', CurLike,  'Rank', MPIRank
  
  !!!    accpt = (Like /= logZero) .and. (CurLike > Like .or. randexp1() > Like - CurLike) 
               !Include the min() so that compilers not doing optimal compilation don't complain
diff -r -c -b -B cosmomc/source/params_CMB.f90 cosmomc_wiggles/source/params_CMB.f90
*** cosmomc/source/params_CMB.f90	2010-05-20 12:13:02.000000000 +0200
--- cosmomc_wiggles/source/params_CMB.f90	2010-05-27 15:18:17.855515009 +0200
***************
*** 57,80 ****
         integer, intent(in) :: in
  
  
         if (Power_Name == 'power_tilt') then
  
-        P%InitPower%k_0_scalar = pivot_k
-        P%InitPower%k_0_tensor = pivot_k
- 
         P%InitPower%ScalarPowerAmp(in) = cl_norm*CMB%norm(norm_As)
         P%InitPower%rat(in) = CMB%norm(norm_amp_ratio)
          
         P%InitPower%an(in) = CMB%InitPower(1)
         P%InitPower%ant(in) = CMB%InitPower(2)
         P%InitPower%n_run(in) = CMB%InitPower(3)
!        if (inflation_consistency) then
!          P%InitPower%ant(in) = - CMB%norm(norm_amp_ratio)/8.
!           !note input n_T is ignored, so should be fixed (to anything)
!        end if
         else
           stop 'params_CMB:Wrong initial power spectrum'
         end if
  
      end subroutine SetCAMBInitPower
   
--- 57,118 ----
         integer, intent(in) :: in
  
  
+ !wig
         if (Power_Name == 'power_tilt') then
  
         P%InitPower%ScalarPowerAmp(in) = cl_norm*CMB%norm(norm_As)
         P%InitPower%rat(in) = CMB%norm(norm_amp_ratio)
          
         P%InitPower%an(in) = CMB%InitPower(1)
         P%InitPower%ant(in) = CMB%InitPower(2)
         P%InitPower%n_run(in) = CMB%InitPower(3)
! 
!           P%InitPower%logomega(in) = CMB%InitPower(5)
!           P%InitPower%omega(in) = 10**(P%InitPower%logomega(in))
!           P%InitPower%amposc(in) = CMB%InitPower(6)          
!           P%InitPower%phase(in) = CMB%InitPower(7)
!           P%InitPower%logpowosc(in) = CMB%InitPower(8)
!           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)
!           P%InitPower%eps3(in) =  CMB%InitPower(3)
! 
!           P%InitPower%logeps1invsigma0(in) = CMB%InitPower(5)
!           P%InitPower%loginvsigma0(in) = P%InitPower%logeps1invsigma0(in) &
!                - log10(P%InitPower%eps1(in))
!           P%InitPower%invsigma0(in) = 10**(P%InitPower%loginvsigma0(in))
!           P%InitPower%Xsigma0(in) = CMB%InitPower(6)
!           P%InitPower%varphi(in)  = CMB%InitPower(7)
!           !CMB%InitPower(7) unused.
!           CMB%InitPower(7) = 1./P%InitPower%invsigma0(in)
! 
!        elseif (Power_Name == 'power_dbi') then
! 
!           P%InitPower%P_scalar =  cl_norm*CMB%norm(norm_As)
!           P%InitPower%eps1(in) =  10.**(CMB%InitPower(4))
!           P%InitPower%del1(in) = P%InitPower%eps1(in) - CMB%InitPower(1)
!           P%InitPower%eps2(in) =  CMB%InitPower(2) - P%InitPower%del1(in)
!           P%InitPower%eps3(in) =  CMB%InitPower(3)
! 
!           P%InitPower%logeps1invsigma0(in) = CMB%InitPower(5)
!           P%InitPower%loginvsigma0(in) = P%InitPower%logeps1invsigma0(in) &
!                - log10(P%InitPower%eps1(in))
!           P%InitPower%invsigma0(in) = 10**(P%InitPower%loginvsigma0(in))
!           P%InitPower%Xsigma0(in) = CMB%InitPower(6)
!           P%InitPower%varphi(in)  = CMB%InitPower(7)
!           !CMB%InitPower(7) unused.
!           CMB%InitPower(7) = 1./P%InitPower%invsigma0(in)
! 
!           P%InitPower%gam(in) = 10**(-CMB%InitPower(8))*P%InitPower%eps1(in)
! 
         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	2010-05-19 13:23:38.000000000 +0200
--- cosmomc_wiggles/source/settings.f90	2010-05-27 15:18:17.855515009 +0200
***************
*** 6,11 ****
--- 6,18 ----
    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)
  
***************
*** 19,25 ****
  !                    (e.g. beam uncertainty modes, etc, specific to dataset)
  
    integer, parameter :: num_hard =7 
!   integer, parameter :: num_initpower = 3 
    integer, parameter :: num_freq_params = 1
    integer, parameter :: num_norm = 2 + num_freq_params 
    integer, parameter :: num_nuisance_params= 0
--- 26,35 ----
  !                    (e.g. beam uncertainty modes, etc, specific to dataset)
  
    integer, parameter :: num_hard =7 
! !wig
! !  integer, parameter :: num_initpower = 3 
!   integer, parameter :: num_initpower = 8
! !end wig
    integer, parameter :: num_freq_params = 1
    integer, parameter :: num_norm = 2 + num_freq_params 
    integer, parameter :: num_nuisance_params= 0
***************
*** 63,68 ****
--- 73,81 ----
  
    integer :: num_threads = 0
    integer :: instance = 0
+ !wig
+   integer :: instance_shift = 0
+ !end wig
    integer :: MPIchains = 1, MPIrank = 0
  
    logical :: Use_LSS = .true.
