diff -c -b -B -r cosmomc/camb/cmbmain.f90 cosmomc_wiggles/camb/cmbmain.f90
*** cosmomc/camb/cmbmain.f90	2008-03-17 17:01:34.000000000 +0100
--- cosmomc_wiggles/camb/cmbmain.f90	2008-04-29 11:36:36.000000000 +0200
***************
*** 637,653 ****
  !     polarization will be calculated. For low values of k we
  !     use a logarithmic spacing. closed case dealt with by SetClosedkValues
           if (CP%WantScalars .and. CP%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
--- 637,663 ----
  !     polarization will be calculated. For low values of k we
  !     use a logarithmic spacing. closed case dealt with by SetClosedkValues
           if (CP%WantScalars .and. CP%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
***************
*** 1031,1037 ****
         qmax_int = min(qmax,max_bessels_etak/CP%tau0)
  
  
!        IntSampleBoost=AccuracyBoost
  
  !     Fixing the # of k for the integration. 
        
--- 1041,1049 ----
         qmax_int = min(qmax,max_bessels_etak/CP%tau0)
  
  
! !wig
! !     IntSampleBoost=AccuracyBoost
!       IntSampleBoost=AccuracyBoost/kIntSampleBoost
  
  !     Fixing the # of k for the integration. 
        
diff -c -b -B -r cosmomc/camb/inidriver.F90 cosmomc_wiggles/camb/inidriver.F90
*** cosmomc/camb/inidriver.F90	2008-03-17 17:01:34.000000000 +0100
--- cosmomc_wiggles/camb/inidriver.F90	2008-04-29 11:51:48.000000000 +0200
***************
*** 249,254 ****
--- 249,268 ----
           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'
***************
*** 260,265 ****
--- 274,310 ----
         if (FeedbackLevel > 0) then
           Age = CAMB_GetAge(P) 
           write (*,'("Age of universe/GYr  = ",f7.3)') Age  
+ !wig
+          if (Power_Name.eq.'power_tilt') then
+ 
+              do i=1,P%InitPower%nn
+                 write (*,'("Scalar tilt          = ",f12.8)') P%InitPower%an(i)
+                 write (*,'("Scalar running       = ",f12.8)') P%InitPower%n_run(i)
+                 if (P%WantTensors) then
+                    write (*,'("primordial T/S       = ",f12.8)') P%InitPower%rat(i)
+                    write (*,'("tensor tilt          = ",f12.8)') P%InitPower%ant(i)
+                 endif                
+                 write (*,'("Wiggles amplitude    = ",f12.8)') P%InitPower%amposc(i)
+                 write (*,'("Wiggles frequency    = ",f12.8)') P%InitPower%logomega(i)
+                 write (*,'("Wiggles phase        = ",f12.8)') P%InitPower%phase(i)
+                 write (*,'("Wiggles power        = ",f12.8)') P%InitPower%powosc(i)
+              enddo
+              write (*,'("Wiggles pivot        = ",f12.8)') P%InitPower%k_0_osc
+ 
+           elseif (Power_Name.eq.'power_sr') then
+ 
+              do i=1,P%InitPower%nn
+                 write (*,'("Slow-roll epsilon1   = ",f12.8)') P%InitPower%eps1(i)
+                 write (*,'("Slow-roll epsilon2   = ",f12.8)') P%InitPower%eps2(i)
+ 		write (*,'("Slow-roll epsilon3   = ",f12.8)') P%InitPower%eps3(i)
+                 write (*,'("TPL amplitude        = ",f12.8)') P%InitPower%Xsigma0(i)
+                 write (*,'("TPL scale            = ",f12.8)') P%InitPower%loginvsigma0(i)
+                 write (*,'("TPL phase            = ",f12.8)') P%InitPower%varphi(i)
+              enddo
+           else
+              stop 'Power_Name unknown'
+           endif
+ !end wig
         end if 
  
         call CAMB_GetResults(P)
diff -c -b -B -r cosmomc/camb/Makefile cosmomc_wiggles/camb/Makefile
*** cosmomc/camb/Makefile	2008-03-17 17:03:40.000000000 +0100
--- cosmomc_wiggles/camb/Makefile	2008-04-29 21:53:34.000000000 +0200
***************
*** 1,62 ****
! #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
! #FLAGS = -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
  
  
  #Files containing evolution equations initial power spectrum module
--- 1,31 ----
!  # >>> DESIGNED FOR GMAKE <<<
!  #CAMB System unified Makefile 
! 
! ext=$(shell uname | cut -c1-3)
! ifeq ($(ext),IRI)
! F90C=f90	
! FFLAGS = -Ofast -mp -n32 -LANG:recursive=ON
! endif
! 
! ifeq ($(ext),Lin)
! F90C=ifort
! FFLAGS= -O -openmp #-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
***************
*** 89,97 ****
  CAMBOBJ       = utils.o subroutines.o inifile.o $(POWERSPECTRUM).o recfast.o $(REIONIZATION).o modules.o \
  	bessels.o $(EQUATIONS).o $(NONLINEAR).o lensing.o cmbmain.o camb.o
  
! default: camb
  
! all: camb $(CAMBLIB)
  
  subroutines: utils.o
  $(POWERSPECTRUM): subroutines.o
--- 58,66 ----
  CAMBOBJ       = utils.o subroutines.o inifile.o $(POWERSPECTRUM).o recfast.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: utils.o
  $(POWERSPECTRUM): subroutines.o
***************
*** 106,112 ****
  camb.o: cmbmain.o
  
  
! camb: $(CAMBOBJ) $(DRIVER)
  	$(F90C) $(F90FLAGS) $(CAMBOBJ) $(DRIVER) -o $@
  
  $(CAMBLIB): $(CAMBOBJ)
--- 75,81 ----
  camb.o: cmbmain.o
  
  
! camb.$(ext): $(CAMBOBJ) $(DRIVER)
  	$(F90C) $(F90FLAGS) $(CAMBOBJ) $(DRIVER) -o $@
  
  $(CAMBLIB): $(CAMBOBJ)
diff -c -b -B -r cosmomc/camb/modules.f90 cosmomc_wiggles/camb/modules.f90
*** cosmomc/camb/modules.f90	2008-03-17 17:01:34.000000000 +0100
--- cosmomc_wiggles/camb/modules.f90	2008-04-29 17:55:06.000000000 +0200
***************
*** 46,53 ****
          character(LEN=*), parameter :: version = 'Mar_08'
          
          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
  
--- 46,55 ----
          character(LEN=*), parameter :: version = 'Mar_08'
          
          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
  
***************
*** 181,196 ****
  
        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
--- 183,200 ----
  
        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
***************
*** 206,212 ****
        integer, parameter:: l0max=4000
  
  !     lmax is max possible number of l's evaluated
!       integer, parameter :: lmax_arr = 100+l0max/7
   
          contains
        
--- 210,223 ----
        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
        
***************
*** 552,557 ****
--- 563,587 ----
          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
***************
*** 668,673 ****
--- 698,706 ----
          end if
          lSet%l0=lind
          lSet%l(1:lind) = ls(1:lind)
+ !wig        
+         endif
+ !end wig
          
        end subroutine initlval
  
***************
*** 687,693 ****
  
        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
--- 720,727 ----
  
        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
***************
*** 703,708 ****
--- 737,756 ----
                         +(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
  
***************
*** 2009,2014 ****
--- 2057,2066 ----
               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 -c -b -B -r cosmomc/camb/params.ini cosmomc_wiggles/camb/params.ini
*** cosmomc/camb/params.ini	2008-03-26 15:24:18.000000000 +0100
--- cosmomc_wiggles/camb/params.ini	2008-04-29 12:26:20.000000000 +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}
***************
*** 32,42 ****
  #Main cosmological parameters, neutrino masses are assumed degenerate
  # If use_phyical set phyiscal densities in baryone, CDM and neutrinos + Omega_k
  use_physical   = T
! ombh2          = 0.022
! omch2          = 0.12
  omnuh2         = 0
  omk            = 0
! hubble         = 70
  #effective equation of state parameter for dark energy, assumed constant
  w              = -1
  #constant comoving sound speed of the dark energy (1=quintessence)
--- 32,42 ----
  #Main cosmological parameters, neutrino masses are assumed degenerate
  # If use_phyical set phyiscal densities in baryone, CDM and neutrinos + Omega_k
  use_physical   = T
! ombh2          = 0.0221
! omch2          = 0.1059
  omnuh2         = 0
  omk            = 0
! hubble         = 72.46
  #effective equation of state parameter for dark energy, assumed constant
  w              = -1
  #constant comoving sound speed of the dark energy (1=quintessence)
***************
*** 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,136 ----
  RECFAST_Heswitch = 6
  
  
+ #primordial power_tilted wiggle
+ osc_logomega(1)           = 1.
+ osc_amp(1)                = 0.1
+ osc_phase(1)              = 0.
+ osc_logpower(1)           = 0.
+ osc_pivot                 = 0.05
+ 
+ #####################################################################################
+ 
+ #for power_sr power spectra
+ 
+ #slow-roll parameters
+ ln[10^{10} P_scalar] = 0.303161E+01
+ eps1(1) = 0.210391E-02
+ eps2(1) = 0.421584E-01
+ eps3(1) = 0
+ 
+ #wiggles parameters
+ logeps1invsigma0(1) = 0.223264E+01
+ Xsigma0(1) = 0.268451E+00
+ varphi(1) = 0.322895E+01
+ 
+ #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
***************
*** 176,188 ****
  
  #Increase accuracy_boost to decrease time steps, use more k values,  etc.
  #Decrease to speed up at cost of worse accuracy. Suggest 0.8 to 3.
! accuracy_boost          = 1
  
  #Larger to keep more terms in the hierarchy evolution. 
! l_accuracy_boost        = 1
  
  #Increase to use more C_l values for interpolation.
  #Increasing a bit will improve the polarization accuracy at l up to 200 -
  #interpolation errors may be up to 3%
  #Decrease to speed up non-flat models a bit
  l_sample_boost          = 1
--- 211,235 ----
  
  #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          = 8
  
  #Larger to keep more terms in the hierarchy evolution. 
! l_accuracy_boost        = 2
  
  #Increase to use more C_l values for interpolation.
  #Increasing a bit will improve the polarization accuracy at l up to 200 -
  #interpolation errors may be up to 3%
  #Decrease to speed up non-flat models a bit
  l_sample_boost          = 1
+ 
+ # Boost the sampling on k sources only, in unit of AccuracyBoost
+ k_src_sample_boost      = 1.
+ 
+ # Boost the sampling on k line of sight integral only, in unit of AccuracyBoost
+ k_int_sample_boost      = 1.
+ 
+ #Use spline interpolation for the Cls
+ do_cls_interpol         = F
+  
+ # Remove undersampling of the Cls: if true lSampleBoost has no effect
+ compute_every_l         = T
diff -c -b -B -r 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	2008-04-29 12:28:29.000000000 +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,106 ----
         real(dl) k_0_scalar, k_0_tensor
         real(dl) ScalarPowerAmp(nnmax)
   
+ !wig
+ ! slow-roll power
+        real(dl) P_scalar
+        real(dl) k_star
+        real(dl) eps1(nnmax)
+        real(dl) eps2(nnmax)
+        real(dl) eps3(nnmax)
+ ! tpl power
+        real(dl) invsigma0(nnmax)
+        real(dl) loginvsigma0(nnmax)
+        real(dl) logeps1invsigma0(nnmax)
+        real(dl) Xsigma0(nnmax)
+        real(dl) varphi(nnmax) 
+ ! osc power
+        real(dl) amposc(nnmax)
+        real(dl) omega(nnmax)
+        real(dl) logomega(nnmax)
+        real(dl) phase(nnmax)
+        real(dl) powosc(nnmax)
+        real(dl) logpowosc(nnmax)
+        real(dl) k_0_osc
+ !Internal parameters              
+        real(dl) ks_pivot          !scalar pivot wave number
+        real(dl) kt_pivot          !tensor pivot wave number
+        real(dl) ko_pivot          !wiggles pivot wave number
+        real(dl) P_0_scalar(nnmax) !scalar normalisation during inflation
+        real(dl) P_0_tensor(nnmax) !tensor normalisation during inflation
+        real(dl) as0(nnmax)        !scalar power spectrum coefficients: order 0
+        real(dl) as1(nnmax)        !order 1
+        real(dl) as2(nnmax)        !order 2
+        real(dl) amps0(nnmax)      !amplitudes of tpl effect in scalar power: order 0
+        real(dl) amps1(nnmax)      !order 1
+        real(dl) amps2(nnmax)      !order 2
+        real(dl) args0(nnmax)      !argument of the modulations: order 0
+        real(dl) args1(nnmax)      !order 1
+        real(dl) args2(nnmax)      !order 2
+        real(dl) at0(nnmax)        !Tensor power spectrum coefficients: same as scalars
+        real(dl) at1(nnmax)
+        real(dl) at2(nnmax)
+        real(dl) ampt0(nnmax) 
+        real(dl) ampt1(nnmax)
+        real(dl) ampt2(nnmax)
+        real(dl) argt0(nnmax)
+        real(dl) argt1(nnmax) 
+        real(dl) argt2(nnmax)
+ !end wig
+         
        end Type InitialPowerParams
  
        real(dl) curv  !Curvature contant, set in InitializePowers     
***************
*** 69,74 ****
--- 119,126 ----
          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 ****
--- 129,156 ----
           AP%k_0_tensor = 0.05
           AP%ScalarPowerAmp = 1
  
+          AP%amposc = 0.
+          AP%logomega = 2.
+          AP%omega = 10**(AP%logomega)
+          AP%phase = 0.
+          AP%logpowosc = 0.
+          AP%powosc = 10**(AP%logpowosc)
+          AP%k_0_osc = 0.05
+          
+          AP%eps1 = 0.01
+          AP%eps2 = 0.01 
+          AP%eps3 = 0.
+          AP%P_scalar=1.d-9
+          AP%k_star = 0.05            
+ 
+          AP%Xsigma0 = 0.                
+          AP%logeps1invsigma0 = 0. 
+          AP%loginvsigma0 = AP%logeps1invsigma0 - log10(AP%eps1) 
+          AP%invsigma0=10**(AP%loginvsigma0)
+          AP%varphi=0.
+         
+ !end wig
+ 
         end subroutine SetDefPowerParams
  
         subroutine InitializePowers(AParamSet,acurv)
***************
*** 85,90 ****
--- 159,169 ----
           !Could read spectra from disk here, do other processing, etc.
  
          real(dl) acurv
+ !wig
+         integer :: i
+         real(dl), parameter :: C_const = -0.7296
+         real(dl), parameter :: pi = 3.14159265
+ !end wig
  
          if (AParamSet%nn > nnmax) then
             write (*,*) 'To use ',AParamSet%nn,'power spectra you need to increase'
***************
*** 95,100 ****
--- 174,318 ----
          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') then
+            if(abs(curv) .gt. 1.d-3) then
+               write(*,*)'Non flat cosmology, curv = ',curv
+               write(*,*)'Cannot use SR or TPL spectrum.'
+               stop
+            endif
+ 
+            P%ks_pivot = P%k_star
+            P%kt_pivot = P%k_star
+            P%ko_pivot = P%k_star
+ 
+            do i=1, P%nn        
+ 
+ !TODO: derive second order tpl corrections?
+               if (useSecondOrder) then
+                  if (P%Xsigma0(i).ne.0.) then
+                     useSecondOrder=.false.
+                     write(*,*)'PowerIni: SR second order disabled!'
+                  endif
+               endif
+ 
+ !first order tpl corrections
+               
+               P%loginvsigma0(i) = P%logeps1invsigma0(i) - log10(P%eps1(i))
+               P%invsigma0(i)=10**(P%loginvsigma0(i))
+               
+               P%amps0(i) = 2.*P%Xsigma0(i)
+               P%amps1(i) = pi*P%Xsigma0(i)*(2.*P%eps1(i) + P%eps2(i))
+               P%amps2(i) = 0.
+               if (useTplPhase) then
+                  P%args0(i) = 2.*P%invsigma0(i)*(1. + P%eps1(i)) + P%varphi(i)
+               else
+                  P%args0(i) = P%varphi(i)
+               endif
+               P%args1(i) = 2.*P%invsigma0(i)*P%eps1(i)
+               P%args2(i) = 0.
+ 
+               P%ampt0(i) = P%amps0(i)
+               P%ampt1(i) = pi*P%Xsigma0(i)*2.*P%eps1(i)
+               P%ampt2(i) = 0.
+               P%argt0(i) = P%args0(i)
+               P%argt1(i) = P%args1(i)
+               P%argt2(i) = 0.
+ 
+ 
+ !first and second order slowroll
+ 
+               P%as0(i) = 1. + P%eps1(i)*(-2.*(C_const+1.))-P%eps2(i)*C_const
+               P%as1(i) = -2.*P%eps1(i) - P%eps2(i)
+               P%as2(i) = 0.
+                             
+               P%at0(i) = 1. + P%eps1(i)*(-2.*(C_const+1.))
+               P%at1(i) = -2.*P%eps1(i)
+               P%at2(i) = 0.
+ 
+               if (useSecondOrder) then
+                                   
+                  P%as0(i) = P%as0(i) &
+                       + P%eps1(i)*P%eps1(i) * (2.*C_const**2 + 2.*C_const + pi**2/2. - 5.) &
+                       + P%eps1(i)*P%eps2(i) * (C_const**2 - C_const + 7.*pi**2/12. - 7.) &
+                       + P%eps2(i)*P%eps2(i) * (0.5*C_const**2 + pi**2/8. - 1.) &
+                       + P%eps2(i)*P%eps3(i) * (-0.5*C_const**2 + pi**2/24.)
+ 
+                  P%as1(i) = P%as1(i) &
+                       + P%eps1(i)*P%eps1(i) * 2.*(2.*C_const + 1.) &
+                       + P%eps1(i)*P%eps2(i) * (2.*C_const - 1.) &
+                       + P%eps2(i)*P%eps2(i) * C_const &
+                       + P%eps2(i)*P%eps3(i) * (-C_const)
+ 
+                  P%as2(i) = P%as2(i) &
+                       + P%eps1(i)*P%eps1(i) * 4. &
+                       + P%eps1(i)*P%eps2(i) * 2. &
+                       + P%eps2(i)*P%eps2(i) * 1. &
+                       + P%eps2(i)*P%eps3(i) * (-1.)
+ 
+                  
+                  P%at0(i) = P%at0(i) &
+                       +  P%eps1(i)*P%eps1(i) * (2.*C_const**2 + 2.*C_const + pi**2/2. - 5.) &
+                       +  P%eps1(i)*P%eps2(i) * (-C_const**2 - 2.*C_const + pi**2/12. - 2.)
+                  
+                  P%at1(i) = P%at1(i) &
+                       + P%eps1(i)*P%eps1(i) * 2.*(2.*C_const + 1.) &
+                       + P%eps1(i)*P%eps2(i) * 2.*(-C_const - 1.)
+ 
+                  P%at2(i) = P%at2(i)  &
+                       + P%eps1(i)*P%eps1(i) * 4. &
+                       + P%eps1(i)*P%eps2(i) * (-2.)
+                       
+               endif
+ 
+               
+               !Relative normalisation of scalar and tensor spectra
+                  
+               P%P_0_scalar(i)=P%P_scalar/(P%as0(i))
+               P%P_0_tensor(i)=P%P_scalar/(P%as0(i))*P%eps1(i)*16.
+ 
+            enddo
+ 
+         else
+            stop 'No such initial power spectra'
+         endif
+ !end wig
  
         end subroutine InitializePowers
         
***************
*** 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
  
--- 337,376 ----
          real(dl) ScalarPower,k, lnrat
          integer in
  
! !wig
!         real :: angle
!         
! !          lnrat = log(k/P%k_0_scalar)
! !          ScalarPower=P%ScalarPowerAmp(in)*exp((P%an(in)-1)*lnrat + P%n_run(in)/2*lnrat**2)   
!                  
!         lnrat = log(k/P%ks_pivot)
! 
!         if (Power_Name.eq.'power_tilt') then
! 
!            angle = P%args0(in) + P%args1(in) * ((k/P%ko_pivot)**P%powosc(in)-1.)
! 
!            ScalarPower = exp(P%as1(in)*lnrat + P%as2(in)*lnrat**2) &
!                 * (1. - P%amps0(in)*cos(angle))
! 
!         elseif (Power_Name.eq.'power_sr') then
! 
!            angle = P%args0(in) + P%args1(in) * log(k/P%ko_pivot) 
! 
!            ScalarPower = (P%as0(in) + P%as1(in)*lnrat + 0.5*P%as2(in)*lnrat**2) &
!                 * (1. - P%amps0(in)*cos(angle)) - P%amps1(in)*sin(angle)
! 
!         else
!            stop 'No such initial power spectra'
!         endif
! 
!         ScalarPower = ScalarPower*P%P_0_scalar(in)
! 
!         if(ScalarPower .lt. 0.0) then
!            ScalarPower =0.
!         end if
! 
! 
! !end wig
   
        end function ScalarPower
  
***************
*** 142,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
  
--- 392,434 ----
     
          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') 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 ****
--- 440,448 ----
           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 ****
--- 467,522 ----
             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)
+ 
+          else
+             stop 'No such initial power spectra'
+          endif
+ !end wig
           Power_Descript = num
  
         end  function Power_Descript
***************
*** 195,200 ****
--- 532,544 ----
             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 ****
--- 552,609 ----
                   
                   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') 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%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 -c -b -B -r cosmomc/distparams.ini cosmomc_wiggles/distparams.ini
*** cosmomc/distparams.ini	2008-03-07 15:48:00.000000000 +0100
--- cosmomc_wiggles/distparams.ini	2008-04-29 14:14:50.000000000 +0200
***************
*** 1,7 ****
  #Params for "getdist" - for processing .txt chain information
  
  #if zero, columnnum calculated automatically as total number of columns
! columnnum = 22
  file_root = chains/test
  out_root = 
  out_dir = 
--- 1,7 ----
  #Params for "getdist" - for processing .txt chain information
  
  #if zero, columnnum calculated automatically as total number of columns
! columnnum = 0
  file_root = chains/test
  out_root = 
  out_dir = 
***************
*** 42,47 ****
--- 42,48 ----
  
  plot_meanlikes = T
  shade_meanlikes = T
+ td_meanlikes = F
  
  # if non-zero, output _thin file, thinned by thin_factor
  thin_factor = 0
***************
*** 82,88 ****
  #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 =
--- 83,89 ----
  #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')
  matlab_colscheme =
***************
*** 103,135 ****
  #w
  lab7 = 
  
! lab8 = n_s
! #n_t
! lab9 =  
! lab10 = n_{run}
  #markerx adds vertical line to MatLab 1D plot
  #marker10 = 0
  
! lab11 = log[10^{10} A_s]
! #amp ratio
! lab12 = r
! lab13 =  A_{SZ}
! 
! lab14 = \Omega_\Lambda
! lab15 = Age/GYr
! lab16 = \Omega_m
! lab17 = \sigma_8
! lab18 = z_{re}
  #r_{10}
! lab19 = 
! lab20 = H_0
  
  #Need to give limits if prior cuts off distribution where not very small
  limits4 = 0.01 N
  limits6 = 0 N
  limits7 = -1 N
  limits12 = 0 N
! limits13 = 0 2
  
  #all_limits sets all limitsxx for all variables to the same; can be useful for bins
  all_limits = 
--- 104,148 ----
  #w
  lab7 = 
  
! #n_s eps1
! lab8 = \epsilon_1
! #n_t eps2
! lab9 = \epsilon_2
! #n_{run}
! lab10 = \epsilon_3
  #markerx adds vertical line to MatLab 1D plot
  #marker10 = 0
  
! #log(\omega) or log(\epsilon_1 / \sigma_0)
! lab11 = log(\epsilon_1 / \sigma_0)
! #amp osc A_\omega or |x| \sigma_0
! lab12 = |x| \sigma_0
! #phase
! lab13 = \psi
! #log(p) or log(\sigma_0)
! lab14 = log(\sigma_0)
! #ln[10^{10} A_s] or ln[10^{10} P_scalar]
! lab15 = ln[10^{10} P_*]
! 
! lab16 =
! lab17 = A_{SZ}
! 
! lab18 = \Omega_\Lambda
! lab19 = Age/GYr
! lab20 = \Omega_m
! lab21 = \sigma_8
! lab22 = z_{re}
  #r_{10}
! lab23 = 
! lab24 = H_0
  
  #Need to give limits if prior cuts off distribution where not very small
  limits4 = 0.01 N
  limits6 = 0 N
  limits7 = -1 N
+ limits8 = 0 N
  limits12 = 0 N
! #limits13 = 0 2
  
  #all_limits sets all limitsxx for all variables to the same; can be useful for bins
  all_limits = 
diff -c -b -B -r cosmomc/params.ini cosmomc_wiggles/params.ini
*** cosmomc/params.ini	2008-03-07 19:11:58.000000000 +0100
--- cosmomc_wiggles/params.ini	2008-04-29 14:20:07.000000000 +0200
***************
*** 58,75 ****
  
  #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
--- 58,75 ----
  
  #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
***************
*** 126,137 ****
  #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)
--- 126,137 ----
  #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)
***************
*** 141,147 ****
  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
--- 141,152 ----
  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
***************
*** 174,189 ****
  #w
  param7 = -1 -1 -1 0 0
  
! #n_s
! param8 = 0.95 0.5 1.5 0.02 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 2.7 4 0.01 0.01
! #amp_ratio
! param12 = 0 0 0 0 0
  #SZ amplitude, as in WMAP analysis
! param13 = 1 0 2 0.4 0.4
--- 179,242 ----
  #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 = 0.95 0.5 1.5 0.02 0.01
! ##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
! 
! ##log[10^10 A_s]
! #param15 = 3 2.7 4 0.01 0.01
! ##amp_ratio
! #param16 = 0 0 0 0 0
! ##SZ amplitude, as in WMAP analysis
! #param17 = 1 0 2 0.4 0.4
! 
!  
! ########################################################################################
! #		        	For Power_Name = power_sr                              #
! ########################################################################################
!  
! #eps1
! param8 = 0.01 0.00001 0.2 0.01 0.02
! #eps2
! param9 = -0.04 -0.2 0.2  0.01 0.02
! #eps3
! param10 = -0.04 -0.2 0.2  0.01 0.02
!  
! #log(eps1*invsigma0)
! param11 = 1.5 1 2.5 0.5 0.2
! #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
  #SZ amplitude, as in WMAP analysis
! param17 = 1 0 2 0.4 0.4
! 
diff -c -b -B -r cosmomc/source/CMB_Cls_simple.f90 cosmomc_wiggles/source/CMB_Cls_simple.f90
*** cosmomc/source/CMB_Cls_simple.f90	2008-03-13 20:19:20.000000000 +0100
--- cosmomc_wiggles/source/CMB_Cls_simple.f90	2008-04-29 14:20:47.000000000 +0200
***************
*** 288,293 ****
--- 288,315 ----
     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)
***************
*** 319,327 ****
              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
                  
--- 341,351 ----
              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 -c -b -B -r cosmomc/source/cmbtypes.f90 cosmomc_wiggles/source/cmbtypes.f90
*** cosmomc/source/cmbtypes.f90	2008-03-06 19:01:34.000000000 +0100
--- cosmomc_wiggles/source/cmbtypes.f90	2008-04-29 14:21:03.000000000 +0200
***************
*** 6,13 ****
  
  
  !Number of Cls, 1 for just temperature, 3 (4) for polarization (with B)
!   integer, parameter  :: num_cls  = 3
! 
  !l_max. Tensors are not computed unless compute_tensors = T in input file
  !Make these multiples of 50, should be 50 more than you need accurately
    integer, parameter :: lmax = 2100, lmax_tensor = 400
--- 6,14 ----
  
  
  !Number of Cls, 1 for just temperature, 3 (4) for polarization (with B)
! !wig
!   integer, parameter  :: num_cls  = 4
! !end wig
  !l_max. Tensors are not computed unless compute_tensors = T in input file
  !Make these multiples of 50, should be 50 more than you need accurately
    integer, parameter :: lmax = 2100, lmax_tensor = 400
***************
*** 62,68 ****
    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
--- 63,72 ----
    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 -c -b -B -r cosmomc/source/driver.F90 cosmomc_wiggles/source/driver.F90
*** cosmomc/source/driver.F90	2008-04-04 10:55:33.000000000 +0200
--- cosmomc_wiggles/source/driver.F90	2008-04-29 14:21:18.000000000 +0200
***************
*** 53,64 ****
  
  #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')
  
--- 53,72 ----
  
  #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')
  
***************
*** 80,85 ****
--- 88,100 ----
  
          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 -c -b -B -r cosmomc/source/GetDist.f90 cosmomc_wiggles/source/GetDist.f90
*** cosmomc/source/GetDist.f90	2008-04-04 10:55:33.000000000 +0200
--- cosmomc_wiggles/source/GetDist.f90	2008-04-29 14:21:29.000000000 +0200
***************
*** 78,83 ****
--- 78,86 ----
          logical force_twotail 
          logical smoothing, plot_meanlikes, plot_NDcontours
          logical shade_meanlikes, make_single_samples
+ !addon
+         logical td_meanlikes
+ !end addon
          integer single_thin,cust2DPlots(max_cols**2)
          integer ix_min(max_cols),ix_max(max_cols)
          real center(max_cols),width(max_cols)   
***************
*** 1209,1221 ****
   !             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('"//trim(plot_data_dir)//"','" // 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
                
             
--- 1212,1240 ----
   !             write (aunit,fmt) (/(I, I=ix_min(j), ix_max(j), 1)/)*width(j)+center(j)
                if (DoShade) then
                if (shade_meanlikes) then
+ !addon
                  write (aunit,'(a)') "load (fullfile('"//trim(plot_data_dir)//"','" // 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 addon
+ 
                else     
                  write (aunit,*) 'contourf(x1,x2,pts,64);'
                end if
! !addon
! !              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(''perspective'');'
!                  write (aunit,*) 'set(gca,''climmode'',''manual'');alpha(0.9);axis tight'
!               else
                write (aunit,*) 'set(gca,''climmode'',''manual''); shading flat; hold on;'
+               endif
+ !end addon
+ 
                end if
                
             
***************
*** 1537,1542 ****
--- 1556,1564 ----
            write (*,*) 'WARNING: Mapping params - .covmat file is new params.'
  
      shade_meanlikes = Ini_Read_Logical('shade_meanlikes',.false.)   
+ !addon
+     td_meanlikes = Ini_Read_Logical('td_meanlikes',.false.)  
+ !end addon
  
      matlab_col = Ini_Read_String('matlab_colscheme')
      smoothing = Ini_Read_Logical('smoothing',.true.)
***************
*** 2089,2095 ****
--- 2111,2124 ----
                write (51,*)  'set(gcf, ''PaperUnits'',''inches'');'
       
                if (plot_row < plot_col .and. plot_row*plot_col>9) then
+ !addon
+ !            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 addon           
                else
                  write (51,*) 'set(gcf, ''PaperPosition'',[ 0 0 8 8]);';
                end if
diff -c -b -B -r cosmomc/source/Makefile cosmomc_wiggles/source/Makefile
*** cosmomc/source/Makefile	2008-03-13 20:21:49.000000000 +0100
--- cosmomc_wiggles/source/Makefile	2008-04-29 20:34:09.000000000 +0200
***************
*** 1,103 ****
! #You may need to edit the library paths for MKL for Intel
! #Beware of using optmizations that lose accuracy - may give errors when running
  
! #Edit for CFITSIO directories
! cfitsio = /usr/local/cfitsio/intel10/64/3.040/lib
! #/usr/local/Cluster-Users/cpac/cmb/2.1.0/cfitsio
! 
! WMAP = /home/aml1005/WMAP5/likelihood_v3
! 
! IFLAG = -I
! INCLUDE=
! 
! #Intel MPI (assuming mpif77 set to point to ifort)
! #change -lmkl_ia32 to -lmkl_p3 for MKL versions earlier than 6 (6 needed for ifc 8+) 
! F90C     = mpif90
! FFLAGS = -O2 -ip -W0 -WB -openmp -fpp -DMPI -vec_report0
! LAPACKL = -L/usr/local/intel/mkl/9.1/lib/em64t -lmkl_lapack -lmkl_em64t -lguide -lpthread
! 
! #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/9.0.018/lib/em64t -lmkl_lapack -lmkl -lguide -lpthread
! #INCLUDE=
! 
! # Use with Intel 9 on the Altix (eg. COSMOS) [this is old]
! # (do "module load icomp90" before compiling, or add to .bashrc)
! #F90C = ifort
! #FFLAGS = -O3 -w -fpp2 -DMPI 
! #LAPACKL = -lscs -lguide -lmpi -ldl
! 
! #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
! 
! #Intel fortran 8 on COSMOS
! #F90C     = ifort
! #FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp2 -lmpi -lscs  -ldl -DMPI
! #LAPACKL = -L/opt/intel/mkl80/lib/64 -lmkl_lapack -lmkl_ipf -limf
! 
! #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
! 
! #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     = gfc
! #FFLAGS = -O2 -ffree-form 
! #LAPACKL = -Wl,-framework -Wl,accelerate
! 
! #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
! 
! #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,57 ----
! # >>> 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= -Ofast -mp -n32 -LANG:recursive=ON -lmpi -DMPI
! WMAPFLAGS= $(FFLAGS)
! LAPACKL = -lcomplib.sgimath_mp
! INCLUDE = -I../camb
! CFITSIODIR = /people/ringeval/usr
! endif
! 
! ifeq ($(ext),Lin)
! F90C=ifort
! FFLAGS= -O -openmp -fpp
! WMAPFLAGS= -O -fpp
! LAPACKL = -L/opt/intel/mkl/9.1/lib/64 -lmkl -lmkl_lapack
! INCLUDE = -I../camb 
! CFITSIODIR = /opt/intel/cfitsio
! 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 = 
! endif
! 
! ifeq ($(ext),Sun)
! F90C=f90
! FFLAGS= -O4 -xarch=native64 -openmp -ftrap=%none
! WMAPFLAGS= $(FFLAGS)
! LAPACKL = -lsunperf -lfsu
! INCLUDE = -I../camb -M../camb
! CFITSIODIR = 
! 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 = 
! endif
! 
! 
! WMAPDIR = ../WMAP
  
  PROPOSE = propose.o
  CLSFILE = CMB_Cls_simple.o
***************
*** 105,138 ****
  #Can use params_H if you prefer more generic parameters
  PARAMETERIZATION = params_CMB.o
  
! WMAP_inc = $(IFLAG)$(cfitsio)/include $(IFLAG)$(WMAP)
! WMAP_lib = -L$(cfitsio)/lib -L$(WMAP) -lcfitsio
  
! CLSLIB  = $(WMAP_lib) -L../camb -lcamb $(LAPACKL) 
  
! F90FLAGS = -DMATRIX_SINGLE $(FFLAGS) $(WMAP_inc) $(IFLAG)../camb $(INCLUDE)
  
  DISTFILES = utils.o Matrix_utils.o settings.o GetDist.o
  
! WMAPOBJS = $(WMAP)/read_archive_map.o \
! 	$(WMAP)/read_fits.o \
! 	$(WMAP)/WMAP_5yr_options.o \
! 	$(WMAP)/WMAP_5yr_util.o \
! 	$(WMAP)/WMAP_5yr_tt_pixlike.o \
! 	$(WMAP)/WMAP_5yr_teeebb_pixlike.o \
! 	$(WMAP)/WMAP_5yr_likelihood.o \
! 	$(WMAP)/WMAP_5yr_gibbs.o \
! 	$(WMAP)/WMAP_5yr_tt_beam_ptsrc_chisq.o \
! 	$(WMAP)/br_mod_dist.o 
  
  OBJFILES= $(WMAPOBJS) utils.o Matrix_utils.o settings.o cmbtypes.o Planck_like.o \
  	cmbdata.o WeakLen.o mpk.o supernovae.o SDSSLy-a-v3.o\
  	$(CLSFILE) paramdef.o $(PROPOSE) $(PARAMETERIZATION) calclike.o \
  	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
  
! default: cosmomc
  
! all : cosmomc getdist
  
  settings.o: utils.o
  cmbtypes.o: settings.o
--- 59,94 ----
  #Can use params_H if you prefer more generic parameters
  PARAMETERIZATION = params_CMB.o
  
! WMAPINCLUDE = -I$(CFITSIODIR)/include
  
! CLSLIB = -L$(CFITSIODIR)/lib -L../camb -lcamb -lcfitsio $(LAPACKL)
  
! F90FLAGS = -DMATRIX_SINGLE $(FFLAGS) $(INCLUDE)
  
  DISTFILES = utils.o Matrix_utils.o settings.o GetDist.o
  
! WMAPOBJS = read_archive_map.o \
! 	read_fits.o \
! 	healpix_types.o \
! 	br_mod_dist.o \
! 	WMAP_5yr_options.o \
! 	WMAP_5yr_util.o \
! 	WMAP_5yr_gibbs.o \
! 	WMAP_5yr_tt_pixlike.o \
! 	WMAP_5yr_tt_beam_ptsrc_chisq.o \
! 	WMAP_5yr_teeebb_pixlike.o \
! 	WMAP_5yr_tetbeebbeb_pixlike.o \
! 	WMAP_5yr_likelihood.o
! 
  
  OBJFILES= $(WMAPOBJS) utils.o Matrix_utils.o settings.o cmbtypes.o Planck_like.o \
  	cmbdata.o WeakLen.o mpk.o supernovae.o SDSSLy-a-v3.o\
  	$(CLSFILE) paramdef.o $(PROPOSE) $(PARAMETERIZATION) calclike.o \
  	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
  
! default: cosmomc.$(ext)
  
! all : cosmomc.$(ext) getdist.$(ext)
  
  settings.o: utils.o
  cmbtypes.o: settings.o
***************
*** 156,161 ****
--- 112,124 ----
  .f.o:
  	f77 $(F90FLAGS) -c $<
  
+ %.o: $(WMAPDIR)/%.f90
+ 	$(F90C) $(WMAPFLAGS) $(WMAPINCLUDE) -c $<
+ 
+ %.o: $(WMAPDIR)/%.F90
+ 	$(F90C) $(WMAPFLAGS) $(WMAPINCLUDE) -c $<
+ 
+ 
  %.o: %.f90
  	$(F90C) $(F90FLAGS) -c $*.f90
  
***************
*** 163,174 ****
  	$(F90C) $(F90FLAGS) -c $*.F90
  
  
! cosmomc: $(OBJFILES) ../camb/libcamb.a 	
! 	$(F90C) -o ../cosmomc $(OBJFILES) $(CLSLIB) $(F90FLAGS)
  
  
  clean:
  	rm -f *.o *.mod *.d *.pc *.obj ../core
  
! getdist: $(DISTFILES)
! 	$(F90C) -o ../getdist $(DISTFILES) $(CLSLIB) $(F90FLAGS) 
--- 126,137 ----
  	$(F90C) $(F90FLAGS) -c $*.F90
  
  
! cosmomc.$(ext): $(OBJFILES) ../camb/libcamb.a 	
! 	$(F90C) -o ../$@ $(OBJFILES) $(CLSLIB) $(F90FLAGS)
  
  
  clean:
  	rm -f *.o *.mod *.d *.pc *.obj ../core
  
! getdist.$(ext): $(DISTFILES)
! 	$(F90C) -o ../$@ $(DISTFILES) $(CLSLIB) $(F90FLAGS) 
diff -c -b -B -r cosmomc/source/MCMC.f90 cosmomc_wiggles/source/MCMC.f90
*** cosmomc/source/MCMC.f90	2008-03-05 11:39:29.000000000 +0100
--- cosmomc_wiggles/source/MCMC.f90	2008-04-29 14:37:28.000000000 +0200
***************
*** 214,220 ****
    
      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
--- 214,220 ----
    
      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
***************
*** 499,505 ****
             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) 
--- 499,505 ----
             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) 
***************
*** 534,540 ****
  
      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
--- 534,540 ----
  
      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 -c -b -B -r cosmomc/source/params_CMB.f90 cosmomc_wiggles/source/params_CMB.f90
*** cosmomc/source/params_CMB.f90	2008-04-04 10:55:33.000000000 +0200
--- cosmomc_wiggles/source/params_CMB.f90	2008-04-29 14:38:24.000000000 +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,99 ----
         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(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)
!           P%InitPower%eps3(in) =  CMB%InitPower(3)
! 
!           P%InitPower%logeps1invsigma0(in) = CMB%InitPower(4)
!           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(5)
!           P%InitPower%varphi(in)  = CMB%InitPower(6)
!           !CMB%InitPower(7) unused.
!           CMB%InitPower(7) = 1./P%InitPower%invsigma0(in)
! 
         else
           stop 'params_CMB:Wrong initial power spectrum'
         end if
+ !end wig
  
      end subroutine SetCAMBInitPower
   
diff -c -b -B -r cosmomc/source/settings.f90 cosmomc_wiggles/source/settings.f90
*** cosmomc/source/settings.f90	2008-04-04 10:55:34.000000000 +0200
--- cosmomc_wiggles/source/settings.f90	2008-04-29 14:38:44.000000000 +0200
***************
*** 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
  
***************
*** 48,53 ****
--- 58,66 ----
  
    integer :: num_threads = 0
    integer :: instance = 0
+ !wig
+   integer :: instance_shift = 0
+ !end wig
    integer :: MPIchains = 1, MPIrank = 0
  
    logical :: Use_LSS = .true.
