diff -r -c -b -B -N cosmomc/camb/cmbmain.f90 cosmomc_wiggles/camb/cmbmain.f90
*** cosmomc/camb/cmbmain.f90	2006-08-28 18:54:02.000000000 +0100
--- cosmomc_wiggles/camb/cmbmain.f90	2006-09-13 12:19:49.000000000 +0100
***************
*** 638,654 ****
  !     polarization will be calculated. For low values of k we
  !     use a logarithmic spacing. closed case dealt with by SetClosedkValues
           if (CP%WantScalars .and. CP%Reionization .and. CP%AccuratePolarization) then
!             dlnk0=2._dl/10/AccuracyBoost
              !Need this to get accurate low l polarization
           else
!             dlnk0=5._dl/10/AccuracyBoost
              if (CP%closed) dlnk0=dlnk0/2
           end if
  
           if (CP%AccurateReionization) dlnk0 = dlnk0/2
! 
!          dkn1=0.6_dl/taurst/AccuracyBoost 
!          dkn2=1.1_dl/taurst/AccuracyBoost 
            if (CP%WantTensors .or. CP%WantVectors) then
                dkn1=dkn1  *0.8_dl
                dlnk0=dlnk0/2 !*0.3_dl
--- 638,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%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
***************
*** 1036,1042 ****
        integer lognum
        real(dl)  IntSampleBoost
  
!       IntSampleBoost=AccuracyBoost
  
  !     Fixing the # of k for the integration. 
        
--- 1045,1053 ----
        integer lognum
        real(dl)  IntSampleBoost
  
! !wig
! !     IntSampleBoost=AccuracyBoost
!       IntSampleBoost=AccuracyBoost/kIntSampleBoost
  
  !     Fixing the # of k for the integration. 
        
diff -r -c -b -B -N cosmomc/camb/inidriver.F90 cosmomc_wiggles/camb/inidriver.F90
*** cosmomc/camb/inidriver.F90	2006-07-26 17:02:29.000000000 +0100
--- cosmomc_wiggles/camb/inidriver.F90	2006-09-13 12:19:49.000000000 +0100
***************
*** 194,200 ****
--- 194,207 ----
  
             P%InitPower%nn = Ini_Read_Int('initial_power_num')
             if (P%InitPower%nn>nnmax) stop 'Too many initial power spectra - increase nnmax in InitialPower'
+ !wig
+ 
+        Power_Name = Ini_Read_String('power_name')
+        
+        if (Power_Name.eq.'power_tilt') then
+ 
             P%InitPower%rat(:) = 1
+           
             do i=1, P%InitPower%nn
                write (numstr,*) i 
                numstr=adjustl(numstr)
***************
*** 198,208 ****
--- 205,218 ----
             do i=1, P%InitPower%nn
                write (numstr,*) i 
                numstr=adjustl(numstr)
+              
+               if (P%WantScalars .or. P%WantTransfer) then
                P%InitPower%an(i) = &
                     Ini_Read_Double('scalar_spectral_index('//trim(numstr)//')')
  
                P%InitPower%n_run(i) = &
                     Ini_Read_Double('scalar_nrun('//trim(numstr)//')',0._dl)
+               end if
      
                if (P%WantTensors) then
                   P%InitPower%ant(i) = Ini_Read_Double('tensor_spectral_index('//trim(numstr)//')')
***************
*** 211,218 ****
--- 221,266 ----
  
                P%InitPower%ScalarPowerAmp(i) = Ini_Read_Double('scalar_amp('//trim(numstr)//')',1.d0) 
                !Always need this as may want to set tensor amplitude even if scalars not computed
+              
+              P%InitPower%logomega(i) = Ini_Read_Double('osc_logomega('//trim(numstr)//')',2.d0)
+              P%InitPower%amposc(i) = Ini_Read_Double('osc_amp('//trim(numstr)//')',0.d0)             
+              P%InitPower%omega(i) = 10**(P%InitPower%logomega(i))
+              P%InitPower%phase(i) = Ini_Read_Double('osc_phase('//trim(numstr)//')',0.d0)
+              P%InitPower%logpowosc(i) = Ini_Read_Double('osc_logpower('//trim(numstr)//')',0.d0)
+              P%InitPower%powosc(i) = 10**(P%InitPower%logpowosc(i))
             end do
        
+           P%InitPower%k_0_osc = Ini_Read_Double('osc_pivot',1d-2)          
+ 
+        elseif (Power_Name.eq.'power_sr') then
+           
+           do i=1, P%InitPower%nn
+              write (numstr,*) i 
+              numstr=adjustl(numstr)
+              
+              if (P%WantScalars .or. P%WantTransfer) then
+                 P%InitPower%eps1(i)= Ini_Read_Double('eps1('//trim(numstr)//')',0.01d0)
+                 P%InitPower%eps2(i)=Ini_Read_Double('eps2('//trim(numstr)//')',0.01d0)
+ 		P%InitPower%eps3(i)=Ini_Read_Double('eps3('//trim(numstr)//')',0.00d0)
+                 P%InitPower%logeps1invsigma0(i) &
+ 		= Ini_Read_Double('logeps1invsigma0('//trim(numstr)//')',1d0)
+                 P%InitPower%loginvsigma0(i) = P%InitPower%logeps1invsigma0(i) &
+                      - log10(P%InitPower%eps1(i))
+                 P%InitPower%invsigma0(i) = 10**(P%InitPower%loginvsigma0(i))
+                 P%InitPower%Xsigma0(i) = Ini_Read_Double('Xsigma0('//trim(numstr)//')',0d0)
+                 P%InitPower%varphi(i) = Ini_Read_Double('varphi('//trim(numstr)//')',0.d0)
+                 P%InitPower%k_star= Ini_Read_Double('k_star',0.05d0)
+                 P%InitPower%P_scalar= 1e-10*exp(Ini_Read_Double('ln[10^{10} P_scalar]',1d0))
+              endif
+              
+           enddo
+           
+        else
+           stop 'No such initial power spectrum'
+        endif
+ 
+ !end wig
+ 
              if (P%WantScalars .or. P%WantTransfer) then
              P%Scalar_initial_condition = Ini_Read_Int('initial_condition',initial_adiabatic)
              if (P%Scalar_initial_condition == initial_vector) then
***************
*** 278,283 ****
--- 326,344 ----
           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'
***************
*** 289,294 ****
--- 350,386 ----
         if (FeedbackLevel > 0) then
           Age = CAMB_GetAge(P) 
           write (*,'("Age of universe/GYr  = ",f7.3)') Age  
+ !wig
+          if (Power_Name.eq.'power_tilt') then
+ 
+              do i=1,P%InitPower%nn
+                 write (*,'("Scalar tilt          = ",f12.8)') P%InitPower%an(i)
+                 write (*,'("Scalar running       = ",f12.8)') P%InitPower%n_run(i)
+                 if (P%WantTensors) then
+                    write (*,'("primordial T/S       = ",f12.8)') P%InitPower%rat(i)
+                    write (*,'("tensor tilt          = ",f12.8)') P%InitPower%ant(i)
+                 endif                
+                 write (*,'("Wiggles amplitude    = ",f12.8)') P%InitPower%amposc(i)
+                 write (*,'("Wiggles frequency    = ",f12.8)') P%InitPower%logomega(i)
+                 write (*,'("Wiggles phase        = ",f12.8)') P%InitPower%phase(i)
+                 write (*,'("Wiggles power        = ",f12.8)') P%InitPower%powosc(i)
+              enddo
+              write (*,'("Wiggles pivot        = ",f12.8)') P%InitPower%k_0_osc
+ 
+           elseif (Power_Name.eq.'power_sr') then
+ 
+              do i=1,P%InitPower%nn
+                 write (*,'("Slow-roll epsilon1   = ",f12.8)') P%InitPower%eps1(i)
+                 write (*,'("Slow-roll epsilon2   = ",f12.8)') P%InitPower%eps2(i)
+ 		write (*,'("Slow-roll epsilon3   = ",f12.8)') P%InitPower%eps3(i)
+                 write (*,'("TPL amplitude        = ",f12.8)') P%InitPower%Xsigma0(i)
+                 write (*,'("TPL scale            = ",f12.8)') P%InitPower%loginvsigma0(i)
+                 write (*,'("TPL phase            = ",f12.8)') P%InitPower%varphi(i)
+              enddo
+           else
+              stop 'Power_Name unknown'
+           endif
+ !end wig
         end if 
  
         call CAMB_GetResults(P)
diff -r -c -b -B -N cosmomc/camb/Makefile cosmomc_wiggles/camb/Makefile
*** cosmomc/camb/Makefile	2006-03-14 11:07:24.000000000 +0000
--- cosmomc_wiggles/camb/Makefile	2006-09-13 12:19:49.000000000 +0100
***************
*** 1,57 ****
! #CAMB Makefile
! #Edit for your compiler
! #Note there are many ifc versions, some of which behave oddly
! 
! #Intel 8 , -openmp toggles mutli-processor:
! F90C     = ifort
! FFLAGS = -openmp -ip -O2 -Vaxlib -W0 -WB -fpp2
! 
! # Intel 9 on IA-64 (eg. COSMOS)
! # (do "module load icomp90" before compiling)
! #F90C = ifort
! #FFLAGS = -openmp -fpp2 -w -O3 -ip -mP2OPT_hlo_prefetch=F
! 
! #Intel ifc, add -openmp for multi-processor (some have bugs):
! #F90C     = ifc
! #FFLAGS = -O2 -Vaxlib -ip -W0 -WB -quiet -fpp2
! #some systems can can also add e.g. -tpp7 -xW
! 
! #SGI, -mp toggles multi-processor. Use -O2 if -Ofast gives problems.
! #F90C     = f90
! #FFLAGS  = -Ofast -mp
! 
! #Digital/Compaq fortran, -omp toggles multi-processor
! #F90C    = f90
! #FFLAGS  = -omp -O4 -arch host -math_library fast -tune host -fpe1
! 
! #Absoft ProFortran, single processor:
! #F90C     = f95
! #FFLAGS = -O2 -cpu:athlon -s -lU77 -w -YEXT_NAMES="LCS" -YEXT_SFX="_"
! 
! #NAGF95, single processor:
! #F90C     = f95
! #FFLAGS = -DNAGF95 -O3
! 
! #G95 compiler
! #F90C   = g95
! #FFLAGS = -O2
! 
! #Sun, single processor:
! #F90C     = f90
! #FFLAGS = -O5
! 
! #Sun V880
! #F90C = mpf90
! #FFLAGS =  -O4 -openmp -ftrap=%none -dalign -DMPI
! 
! #Sun parallel enterprise:
! #F90C     = f95
! #FFLAGS =  -O2 -xarch=native64 -openmp -ftrap=%none
! #try removing -openmp if get bus errors. -03, -04 etc are dodgy.
! 
! #IBM XL Fortran, multi-processor (run gmake)
! #F90C     = xlf90_r
! #FFLAGS  = -qsmp=omp -qsuffix=f=f90:cpp=F90 -O3 -qstrict -qarch=pwr3 -qtune=pwr3
  
  
  #Files containing evolution equations initial power spectrum module
--- 1,33 ----
! # >>> DESIGNED FOR GMAKE <<<
! 
! #CAMB System unified Makefile
! 
! ext=$(shell uname | cut -c1-3)
! 
! ifeq ($(ext),IRI)
! F90C=f90
! FFLAGS = -Ofast -mp -n32 -LANG:recursive=ON
! endif
! 
! ifeq ($(ext),Lin)
! F90C=ifort
! FFLAGS= -O -openmp -fpp
! endif 
! 
! ifeq ($(ext),OSF)
! F90C=f90
! FFLAGS= -omp -O -arch host -math_library fast -tune host -fpe1
! endif
! 
! ifeq ($(ext),Sun)
! F90C=f90
! FFLAGS= -O4 -xarch=native64 -openmp -ftrap=%none
! endif
! 
! ifeq ($(ext),AIX)
! F90C=xlf90_r
! FFLAGS= -O4 -qsmp=omp -qmaxmem=-1 -qstrict -qfree=f90 -qsuffix=f=f90:cpp=F90
! endif
  
  
  #Files containing evolution equations initial power spectrum module
***************
*** 83,99 ****
  CAMBOBJ       = subroutines.o inifile.o $(POWERSPECTRUM).o recfast.o modules.o \
  	bessels.o $(EQUATIONS).o $(NONLINEAR).o lensing.o cmbmain.o camb.o
  
! default: camb
  
! all: camb $(CAMBLIB)
  
! camb: $(CAMBOBJ) $(DRIVER)
  	$(F90C) $(F90FLAGS) $(CAMBOBJ) $(DRIVER) -o $@
  
  $(CAMBLIB): $(CAMBOBJ)
  	ar -r $@ $?
  
! camb_fits: writefits.f90 $(CAMBOBJ) $(DRIVER)
  	$(F90C) $(F90FLAGS) -I$(HEALPIXDIR)/include $(CAMBOBJ) writefits.f90 $(DRIVER) $(HEALPIXLD) -DWRITE_FITS -o $@
  
  %.o: %.f90
--- 59,75 ----
  CAMBOBJ       = subroutines.o inifile.o $(POWERSPECTRUM).o recfast.o modules.o \
  	bessels.o $(EQUATIONS).o $(NONLINEAR).o lensing.o cmbmain.o camb.o
  
! default: camb.$(ext)
  
! all: camb.$(ext) $(CAMBLIB)
  
! camb.$(ext): $(CAMBOBJ) $(DRIVER)
  	$(F90C) $(F90FLAGS) $(CAMBOBJ) $(DRIVER) -o $@
  
  $(CAMBLIB): $(CAMBOBJ)
  	ar -r $@ $?
  
! camb_fits.$(ext): writefits.f90 $(CAMBOBJ) $(DRIVER)
  	$(F90C) $(F90FLAGS) -I$(HEALPIXDIR)/include $(CAMBOBJ) writefits.f90 $(DRIVER) $(HEALPIXLD) -DWRITE_FITS -o $@
  
  %.o: %.f90
diff -r -c -b -B -N cosmomc/camb/modules.f90 cosmomc_wiggles/camb/modules.f90
*** cosmomc/camb/modules.f90	2006-08-23 19:28:12.000000000 +0100
--- cosmomc_wiggles/camb/modules.f90	2006-09-13 12:19:49.000000000 +0100
***************
*** 33,40 ****
          character(LEN=*), parameter :: version = 'Aug_06'
          
          integer :: FeedbackLevel = 0 !if >0 print out useful information about the model
! 
!         logical, parameter :: DebugMsgs=.false. !Set to true to view progress and timing
  
          logical, parameter :: DebugEvolution = .false. !Set to true to do all the evolution for all k
  
--- 33,42 ----
          character(LEN=*), parameter :: version = 'Aug_06'
          
          integer :: FeedbackLevel = 0 !if >0 print out useful information about the model
! !wig
! !        logical, parameter :: DebugMsgs=.false. !Set to true to view progress and timing
!         logical, parameter :: DebugMsgs=.true.
! !end wig
  
          logical, parameter :: DebugEvolution = .false. !Set to true to do all the evolution for all k
  
***************
*** 176,189 ****
        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
--- 178,197 ----
        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
***************
*** 199,205 ****
        integer, parameter:: l0max=4000
  
  !     lmax is max possible number of l's evaluated
!       integer, parameter :: lmax_arr = 100+l0max/7
   
          contains
        
--- 207,220 ----
        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
        
***************
*** 594,599 ****
--- 609,633 ----
          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
***************
*** 695,700 ****
--- 729,737 ----
          end if
          lSet%l0=lind
          lSet%l(1:lind) = ls(1:lind)
+ !wig        
+         endif
+ !end wig
          
        end subroutine initlval
  
***************
*** 714,720 ****
  
        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
--- 751,758 ----
  
        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
***************
*** 730,735 ****
--- 768,787 ----
                         +(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
  
***************
*** 1975,1980 ****
--- 2027,2036 ----
               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 -N cosmomc/camb/params.ini cosmomc_wiggles/camb/params.ini
*** cosmomc/camb/params.ini	2006-08-12 16:34:12.000000000 +0100
--- cosmomc_wiggles/camb/params.ini	2006-09-13 12:19:49.000000000 +0100
***************
*** 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)
***************
*** 66,79 ****
  
  #Reionization (assumed sharp), ignored unless reionization = T
  reionization         = T
! re_use_optical_depth = F
! re_optical_depth     = 0.04
  #If re_use_optical_depth = F then use following, otherwise ignored
  re_redshift          = 12
  re_ionization_frac   = 1
  
  #Initial power spectrum, amplitude, spectral index and running
  initial_power_num         = 1
  scalar_amp(1)             = 2.3e-9
  scalar_spectral_index(1)  = 1
  scalar_nrun(1)            = 0
--- 66,87 ----
  
  #Reionization (assumed sharp), ignored unless reionization = T
  reionization         = T
! re_use_optical_depth = T
! re_optical_depth     = 0.094
  #If re_use_optical_depth = F then use following, otherwise ignored
  re_redshift          = 12
  re_ionization_frac   = 1
  
  #Initial power spectrum, amplitude, spectral index and running
  initial_power_num         = 1
+ 
+ #Type of initial power spectrum: power_tilt or power_sr
+ power_name = power_sr
+ 
+ ######################################################################################
+ 
+ #for power_tilt power spectra
+ 
  scalar_amp(1)             = 2.3e-9
  scalar_spectral_index(1)  = 1
  scalar_nrun(1)            = 0
***************
*** 82,87 ****
--- 90,122 ----
  initial_ratio(1)          = 1
  #note vector modes use the scalar settings above
  
+ #primordial power_tilted wiggle
+ osc_logomega(1)           = 1.
+ osc_amp(1)                = 0.1
+ osc_phase(1)              = 0.
+ osc_logpower(1)           = 0.
+ osc_pivot                 = 0.05
+ 
+ #####################################################################################
+ 
+ #for power_sr power spectra
+ 
+ #slow-roll parameters
+ ln[10^{10} P_scalar] = 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
***************
*** 114,124 ****
  
  
  #Output files not produced if blank. make camb_fits to use use the FITS setting.
! scalar_output_file = scalCls.dat
! vector_output_file = vecCls.dat
! tensor_output_file = tensCls.dat
! total_output_file  = totCls.dat
! lensed_output_file = lensedCls.dat
  FITS_filename      = scalCls.fits
  
  ##Optional parameters to control the computation speed,accuracy and feedback
--- 149,159 ----
  
  
  #Output files not produced if blank. make camb_fits to use use the FITS setting.
! scalar_output_file = scalar_cls.dat
! vector_output_file = vector_cls.dat
! tensor_output_file = tensor_cls.dat
! total_output_file  = total_cls.dat
! lensed_output_file = lensed_cls.dat
  FITS_filename      = scalCls.fits
  
  ##Optional parameters to control the computation speed,accuracy and feedback
***************
*** 160,172 ****
  
  #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
--- 195,219 ----
  
  #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 -r -c -b -B -N cosmomc/camb/power_tilt.f90 cosmomc_wiggles/camb/power_tilt.f90
*** cosmomc/camb/power_tilt.f90	2005-03-31 00:37:13.000000000 +0100
--- cosmomc_wiggles/camb/power_tilt.f90	2006-09-13 12:19:49.000000000 +0100
***************
*** 33,40 ****
  
       private
      
!       character(LEN=*), parameter :: Power_Name = 'power_tilt'
! 
        integer, parameter :: nnmax= 5 
         !Maximum possible number of different power spectra to use
  
--- 33,43 ----
  
       private
      
! !wig      character(LEN=*), parameter :: Power_Name = 'power_tilt'
!       character(80) :: Power_Name = 'power_sr'
!       logical, parameter :: useTplPhase = .false.      
!       logical, 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     
***************
*** 70,75 ****
--- 120,127 ----
          Type (InitialPowerParams) :: AP
  
           AP%nn     = 1 !number of initial power spectra
+ !wig
+            
           AP%an     = 1 !scalar spectral index
           AP%n_run   = 0 !running of scalar spectral index
           AP%ant    = 0 !tensor spectra index
***************
*** 78,83 ****
--- 130,157 ----
           AP%k_0_tensor = 0.05
           AP%ScalarPowerAmp = 1
  
+          AP%amposc = 0.
+          AP%logomega = 2.
+          AP%omega = 10**(AP%logomega)
+          AP%phase = 0.
+          AP%logpowosc = 0.
+          AP%powosc = 10**(AP%logpowosc)
+          AP%k_0_osc = 0.05
+          
+          AP%eps1 = 0.01
+          AP%eps2 = 0.01 
+          AP%eps3 = 0.
+          AP%P_scalar=1.d-9
+          AP%k_star = 0.05            
+ 
+          AP%Xsigma0 = 0.                
+          AP%logeps1invsigma0 = 0. 
+          AP%loginvsigma0 = AP%logeps1invsigma0 - log10(AP%eps1) 
+          AP%invsigma0=10**(AP%loginvsigma0)
+          AP%varphi=0.
+         
+ !end wig
+ 
         end subroutine SetDefPowerParams
  
         subroutine InitializePowers(AParamSet,acurv)
***************
*** 86,91 ****
--- 160,170 ----
           !Could read spectra from disk here, do other processing, etc.
  
          real(dl) acurv
+ !wig
+         integer :: i
+         real(dl), parameter :: C_const = -0.7296
+         real(dl), parameter :: pi = 3.14159265
+ !end wig
  
          P = AParamSet
          if (P%nn > nnmax) then
***************
*** 96,101 ****
--- 175,319 ----
          curv=acurv         
  
  !Write implementation specific code here...
+ !wig
+         if (Power_Name.eq.'power_tilt') then
+            
+            P%ks_pivot = P%k_0_scalar
+            P%kt_pivot = P%k_0_tensor
+            P%ko_pivot = P%k_0_osc
+ 
+            do i=1, P%nn        
+                                                 
+               P%omega(i) = 10**(P%logomega(i))
+               P%powosc(i) = 10**(P%logpowosc(i))
+ 
+               P%as0(i) = 0.
+               P%as1(i) = P%an(i)-1
+               P%as2(i) = P%n_run(i)/2.
+               P%amps0(i) = P%amposc(i)
+               P%amps1(i) = 0.
+               P%amps2(i) = 0.
+               P%args0(i) = P%phase(i)
+               P%args1(i) = P%omega(i)/P%powosc(i)
+               P%args2(i) = 0.
+               
+               P%at0(i) = 0.
+               P%at1(i) = P%ant(i)
+               P%at2(i) = 0.
+               P%ampt0(i) = P%amps0(i)
+               P%ampt1(i) = 0.
+               P%ampt2(i) = 0.
+               P%argt0(i) = P%phase(i)
+               P%argt1(i) = P%omega(i)/P%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
         
***************
*** 120,127 ****
          real(dl) ScalarPower,k, lnrat
          integer in
  
!           lnrat = log(k/P%k_0_scalar)
!           ScalarPower=P%ScalarPowerAmp(in)*exp((P%an(in)-1)*lnrat + P%n_run(in)/2*lnrat**2)   
   
        end function ScalarPower
  
--- 338,377 ----
          real(dl) ScalarPower,k, lnrat
          integer in
  
! !wig
!         real :: angle
!         
! !          lnrat = log(k/P%k_0_scalar)
! !          ScalarPower=P%ScalarPowerAmp(in)*exp((P%an(in)-1)*lnrat + P%n_run(in)/2*lnrat**2)   
!                  
!         lnrat = log(k/P%ks_pivot)
! 
!         if (Power_Name.eq.'power_tilt') then
! 
!            angle = P%args0(in) + P%args1(in) * ((k/P%ko_pivot)**P%powosc(in)-1.)
! 
!            ScalarPower = exp(P%as1(in)*lnrat + P%as2(in)*lnrat**2) &
!                 * (1. - P%amps0(in)*cos(angle))
! 
!         elseif (Power_Name.eq.'power_sr') then
! 
!            angle = P%args0(in) + P%args1(in) * log(k/P%ko_pivot) 
! 
!            ScalarPower = (P%as0(in) + P%as1(in)*lnrat + 0.5*P%as2(in)*lnrat**2) &
!                 * (1. - P%amps0(in)*cos(angle)) - P%amps1(in)*sin(angle)
! 
!         else
!            stop 'No such initial power spectra'
!         endif
! 
!         ScalarPower = ScalarPower*P%P_0_scalar(in)
! 
!         if(ScalarPower .lt. 0.0) then
!            ScalarPower =0.
!         end if
! 
! 
! !end wig
   
        end function ScalarPower
  
***************
*** 142,151 ****
          real(dl), parameter :: PiByTwo=3.14159265d0/2._dl
     
          integer in
  
-         TensorPower=P%rat(in)*P%ScalarPowerAmp(in)*exp(P%ant(in)*log(k/P%k_0_tensor))
          if (curv < 0) TensorPower=TensorPower*tanh(PiByTwo*sqrt(-k**2/curv-3)) 
  
        end function TensorPower
  
        !Get parameters describing parameterisation (for FITS file)
--- 392,434 ----
          real(dl), parameter :: PiByTwo=3.14159265d0/2._dl
     
          integer in
+ !wig
+         real :: lnrat,angle
+ 
+ !        TensorPower=P%rat(in)*P%ScalarPowerAmp(in)*exp(P%ant(in)*log(k/P%k_0_tensor))
+ !        if (curv < 0) TensorPower=TensorPower*tanh(PiByTwo*sqrt(-k**2/curv-3)) 
+         
+         lnrat = log(k/P%kt_pivot)
+ 
+         if (Power_Name.eq.'power_tilt') then
+ 
+            angle = P%argt0(in) + P%argt1(in) * ((k/P%ko_pivot)**P%powosc(in)-1.)
+ 
+            TensorPower = exp(P%at1(in)*lnrat) &
+                 * (1. - P%ampt0(in)*cos(angle))
+ 
+ 
+         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
  
        !Get parameters describing parameterisation (for FITS file)
***************
*** 156,161 ****
--- 439,447 ----
           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 ****
--- 466,521 ----
             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
diff -r -c -b -B -N cosmomc/distparams.ini cosmomc_wiggles/distparams.ini
*** cosmomc/distparams.ini	2006-08-17 17:38:45.000000000 +0100
--- cosmomc_wiggles/distparams.ini	2006-09-13 12:19:49.000000000 +0100
***************
*** 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 = 
  
--- 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 = 
  
***************
*** 40,45 ****
--- 40,46 ----
  
  plot_meanlikes = T
  shade_meanlikes = T
+ td_meanlikes = F
  
  # if non-zero, output _thin file, thinned by thin_factor
  thin_factor = 0
***************
*** 80,86 ****
  #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 =
--- 81,87 ----
  #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 =
***************
*** 101,129 ****
  #w
  lab7 = 
  
! lab8 = n_s
! #n_t
! lab9 =  
! lab10 = n_{run}
! 
! lab11 = log[10^{10} A_s]
! #amp ratio
! lab12 = r
! lab13 =  
! 
! 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
  #all_limits sets all limitsxx for all variables to the same; can be useful for bins
  all_limits = 
--- 102,142 ----
  #w
  lab7 = 
  
! #n_s or \epsilon_1
! lab8 = \epsilon_1
! #n_t or \epsilon_2
! lab9 =  \epsilon_2
! #n_{run} or \epsilon_3
! lab10 = \epsilon_3
! 
! #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 =
! 
! 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
  #all_limits sets all limitsxx for all variables to the same; can be useful for bins
  all_limits = 
diff -r -c -b -B -N cosmomc/params.ini cosmomc_wiggles/params.ini
*** cosmomc/params.ini	2006-08-28 19:58:08.000000000 +0100
--- cosmomc_wiggles/params.ini	2006-09-13 12:19:49.000000000 +0100
***************
*** 54,71 ****
  
  #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
--- 54,71 ----
  
  #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
***************
*** 122,133 ****
  #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)
--- 122,133 ----
  #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)
***************
*** 137,143 ****
  CMB_lensing = F
  #increase accuracy_level to run CAMB on higher accuracy
  #(default is about 0.3%, accuracy_level=2 around 0.1% at high l)
! accuracy_level = 1
  
  #If action = 1
  redo_likelihoods = T
--- 137,148 ----
  CMB_lensing = F
  #increase accuracy_level to run CAMB on higher accuracy
  #(default is about 0.3%, accuracy_level=2 around 0.1% at high l)
! accuracy_level = 6
! l_accuracy_level = 2
! l_sample_level = 1
! k_src_sample_level = 1.
! k_int_sample_level = 1.
! camb_compute_every_l = T
  
  #If action = 1
  redo_likelihoods = T
***************
*** 170,186 ****
  #w
  param7 = -1 -1 -1 0 0
  
- #n_s
- param8 = 0.95 0.5 1.5 0.03 0.01
- #n_t
- param9 = 0 0 0 0 0
- #n_run
- param10 = 0 0 0 0 0
- 
- #log[10^10 A_s]
- param11 = 3 2.7 4 0.01 0.01
- #amp_ratio
- param12 = 0 0 0 0 0
- #unused fast parameter
- param13 = 0 0 0 0 0
  
--- 175,236 ----
  #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.03 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
+ ##unused fast parameter
+ #param17 = 0 0 0 0 0
+ 
+ ########################################################################################
+ #		        	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
+ #unused
+ param17 = 0 0 0 0 0
diff -r -c -b -B -N cosmomc/source/CMB_Cls_simple.f90 cosmomc_wiggles/source/CMB_Cls_simple.f90
*** cosmomc/source/CMB_Cls_simple.f90	2006-08-17 13:16:15.000000000 +0100
--- cosmomc_wiggles/source/CMB_Cls_simple.f90	2006-09-13 12:19:49.000000000 +0100
***************
*** 266,271 ****
--- 266,293 ----
     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)
***************
*** 297,305 ****
              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
                  
--- 319,329 ----
              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 -N cosmomc/source/cmbtypes.f90 cosmomc_wiggles/source/cmbtypes.f90
*** cosmomc/source/cmbtypes.f90	2006-08-11 17:16:55.000000000 +0100
--- cosmomc_wiggles/source/cmbtypes.f90	2006-09-13 12:19:49.000000000 +0100
***************
*** 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 -r -c -b -B -N cosmomc/source/distparams.ini cosmomc_wiggles/source/distparams.ini
*** cosmomc/source/distparams.ini	1970-01-01 01:00:00.000000000 +0100
--- cosmomc_wiggles/source/distparams.ini	2006-09-13 12:19:49.000000000 +0100
***************
*** 0 ****
--- 1,155 ----
+ #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 = 
+ 
+ #If generated chain at higher temperature can cool for computing results
+ cool = 1
+ 
+ #If 0 assume 1 and no chain filename prefixes
+ chain_num = 4
+ first_chain = 
+ exclude_chain = 
+ 
+ #width of Gaussian smoothing - Should check plots are robust to changes in 
+ #this parameter. Narrow diagonal distributions need larger number
+ #Can also check plots by comparing with setting smoothing=F below
+ num_bins = 20
+ 
+ #For disgarding burn-in if using raw chains
+ #if < 1 interpreted as a fraction of the total number of rows (0.3 ignores first 30% of lines)
+ ignore_rows = 0.5
+ 
+ #if T produced B&W printer friendly output
+ B&W = F
+ #version of MatLab you are using
+ matlab_version = 7
+ 
+ #Switches; auto_label labels all parameters by their number
+ no_plots = F
+ no_tests = F
+ auto_label = F
+ #samples_are_chains = F can be useful for other samples when first two columns not present
+ samples_are_chains = T
+ 
+ #Include these in 1D plots for comparison - must have same parameters
+ compare_num = 0
+ compare1 = basic6_cmb
+ 
+ plot_meanlikes = T
+ shade_meanlikes = T
+ td_meanlikes = F
+ 
+ # if non-zero, output _thin file, thinned by thin_factor
+ thin_factor = 0
+ #Do probabilistic importance sampling to single samples
+ make_single_samples = F
+ single_thin = 4
+ 
+ #Do simple importance sampling
+ adjust_priors = F
+ #Map parameters to other derived quantities
+ map_params = F
+ 
+ #Use a Gaussian smoothing with width of the bin size
+ #Otherwise use top hat bins
+ smoothing = T
+ 
+ contour1 = 0.68
+ contour2 = 0.95
+ 
+ #if we only want 2D plots agains a particular variable
+ plot_2D_param = 0
+ 
+ #if above zero, instead plot just these combinations:
+ #if both zero it will plot most correlated variables
+ plot_2D_num = 0
+ plot1 = 8 1
+ plot2 =
+ 
+ #number of sample plots, colored by third parameter
+ #if last parameter is 0 or -1 colored by the parameter most correlated
+ #with one of the eigenvector directions (e.g. parallel or orthogonal to degeneracy)
+ num_3D_plots = 1
+ 3D_plot1 = 1 2 3
+ 
+ #Output 2D plots for param combos with 1D marginalized plots along the diagonal
+ triangle_plot = T
+ 
+ #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 =
+ 
+ #Parameters to use. If zero use all parameters which have lables.
+ plotparams_num = 0
+ plotparams = 1 2 3 11
+ 
+ #labels for the parameters
+ lab1 = \Omega_b h^2
+ lab2 = \Omega_c h^2
+ lab3 = \theta
+ lab4 = \tau
+ #Omk
+ lab5 =
+ #nufrac 
+ lab6 = 
+ #w
+ lab7 = 
+ 
+ #n_s or \epsilon_1
+ lab8 = \epsilon_1
+ #n_t or \epsilon_2
+ lab9 =  \epsilon_2
+ #n_{run} or \epsilon_3
+ lab10 = \epsilon_3
+ 
+ #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 =
+ 
+ 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
+ #all_limits sets all limitsxx for all variables to the same; can be useful for bins
+ all_limits = 
+ 
+ #compute two-tail marginalized limits irrespective of limits settings above
+ #(otherwise limits are two-tail only for those parameters without limits)
+ force_twotail = F
+ 
+ #PCA - analysis output in file file_root.PCA
+ #number of parameter to do PCA for
+ PCA_num = 0
+ PCA_normparam = 13
+ #The parameters to use
+ PCA_params = 3 13 16
+ #L for log(x), M for log(-x), N for no log
+ PCA_func   = LLLL
diff -r -c -b -B -N cosmomc/source/driver.F90 cosmomc_wiggles/source/driver.F90
*** cosmomc/source/driver.F90	2006-08-17 13:16:15.000000000 +0100
--- cosmomc_wiggles/source/driver.F90	2006-09-13 12:19:49.000000000 +0100
***************
*** 51,62 ****
  
  #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')
  
--- 51,70 ----
  
  #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')
  
***************
*** 78,83 ****
--- 86,98 ----
  
          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 -N cosmomc/source/GetDist.f90 cosmomc_wiggles/source/GetDist.f90
*** cosmomc/source/GetDist.f90	2006-08-18 08:50:43.000000000 +0100
--- cosmomc_wiggles/source/GetDist.f90	2006-09-13 12:19:49.000000000 +0100
***************
*** 68,73 ****
--- 68,76 ----
          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)   
***************
*** 1201,1215 ****
   !             write (aunit,fmt) (/(I, I=ix_min(j), ix_max(j), 1)/)*width(j)+center(j)
                if (DoShade) then
                if (shade_meanlikes) then
                  write (aunit,'(a)') 'load (fullfile(''plot_data'',''' // trim(plotfile) //'_likes''));'
                  write (aunit,*) 'contourf(x1,x2,'//trim(plotfile)//'_likes,64);'
                else     
                  write (aunit,*) 'contourf(x1,x2,pts,64);'
                end if
! 
                write (aunit,*) 'set(gca,''climmode'',''manual''); shading flat; hold on;'
!               end if
!               
             
                if (num_contours /= 0) then
  
--- 1204,1231 ----
   !             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(''plot_data'',''' // trim(plotfile) //'_likes''));'
+                  if (td_meanlikes) then
+                     write (aunit,*) 'surfl(x1,x2,'//trim(plotfile)//'_likes);'
+                  else
                  write (aunit,*) 'contourf(x1,x2,'//trim(plotfile)//'_likes,64);'
+                  endif
+ !end 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
!               endif
             
                if (num_contours /= 0) then
  
***************
*** 1516,1521 ****
--- 1532,1540 ----
            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.)
***************
*** 2052,2058 ****
--- 2071,2084 ----
                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 -r -c -b -B -N cosmomc/source/Makefile cosmomc_wiggles/source/Makefile
*** cosmomc/source/Makefile	2006-04-18 22:27:46.000000000 +0100
--- cosmomc_wiggles/source/Makefile	2006-09-13 12:19:49.000000000 +0100
***************
*** 1,154 ****
! #Settings here are for Intel ifort on CITA cluster; uncomment others below
! #You may need to edit the library paths for MKL for Intel
! #Beware of using optmizations which lose accuracy - may give errors when running
! #With Visual Fortan open cosmomc.dsw workspace file instead.
! 
! #Edit for CFITSIO, HEALPIX directories
! cfitsio = /home/pascal/include/cfitsio
! healpix = /home/pascal/include/Healpix_2.00
! WMAP3 = /home/antlewis/WMAP3
! 
! IFLAG = -I
! 
! #Intel MPI (assuming mpif77 set to point to ifc)
! #change -lmkl_ia32 to -lmkl_p3 for MKL versions earlier than 6 (6 needed for ifc 8+) 
! F90C     = mpif77
! FFLAGS = -O2 -Vaxlib -ip -W0 -WB -openmp -fpp -DMPI
! LAPACKL = -L/opt/intel/mkl61/lib/32 -lmkl_lapack -lmkl_ia32 -lguide -lpthread
! 
! # Use with Intel 9 on the Altix (eg. COSMOS)
! # (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
! 
! #Intel Itanium efc, -openmp toggles multi-processor:
! #F90C     = efc
! #FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp
! #LAPACKL = -L/opt/intel/mkl/lib/64 -lmkl_lapack -lmkl_ipf -lguide -lpthread -static
! 
! #Intel Itanium efc MPI:
! # mpirun -np 4 ./cosmomc params.ini
! #F90C     = efc
! #FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp -lmpi -lscs  -ldl -DMPI
! #LAPACKL = -L/opt/intel/mkl/lib/64 -lmkl_lapack -lmkl_ipf -limf
! 
! #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
! 
! #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
  
  #Can use params_H if you prefer more generic parameters
  PARAMETERIZATION = params_CMB.o
  
! WMAP3_inc = $(IFLAG)$(cfitsio)/include $(IFLAG)$(healpix)/include $(IFLAG)$(WMAP3)
! WMAP3_lib = -L$(cfitsio)/lib -L$(healpix)/lib -L$(WMAP3) \
! 	-lhealpix -lcfitsio
  
! CLSLIB  = $(WMAP3_lib) -L../camb -lcamb $(LAPACKL) 
  
! F90FLAGS = $(FFLAGS) $(WMAP3_inc) $(IFLAG)../camb $(INCLUDE)
  
  DISTFILES = utils.o settings.o GetDist.o
  
! WMAPOBJS = $(WMAP3)/read_archive_map.o \
! 	$(WMAP3)/read_fits.o \
! 	$(WMAP3)/WMAP_3yr_options.o \
! 	$(WMAP3)/WMAP_3yr_util.o \
! 	$(WMAP3)/WMAP_3yr_tt_pixlike.o \
! 	$(WMAP3)/WMAP_3yr_tt_beam_and_ptsrc_corr.o \
! 	$(WMAP3)/WMAP_3yr_teeebb_pixlike.o \
! 	$(WMAP3)/WMAP_3yr_likelihood.o
  
  OBJFILES= $(WMAPOBJS) utils.o settings.o cmbtypes.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
  
  .f.o:
  	f77 $(F90FLAGS) -c $<
  
  %.o: %.f90
! 	$(F90C) $(F90FLAGS) -c $*.f90
  
  %.o: %.F90
! 	$(F90C) $(F90FLAGS) -c $*.F90
! 
  
- cosmomc: $(OBJFILES) 	
- 	$(F90C) -o ../cosmomc $(OBJFILES) $(CLSLIB) $(F90FLAGS)
  
  
  clean:
  	rm -f *.o *.mod *.d *.pc *.obj ../core
  
! getdist: $(DISTFILES)
! 	$(F90C) -o ../getdist $(DISTFILES) $(CLSLIB) $(F90FLAGS) 
--- 1,108 ----
! # >>> 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
! LAPACKL = -lcomplib.sgimath_mp
! INCLUDE = -I../camb
! CFITSIODIR = /people/ringeval/usr
! HEALPIXDIR = /people/ringeval/usr
! endif
! 
! ifeq ($(ext),Lin)
! F90C=ifort
! FFLAGS= -O -openmp -fpp
! WMAPFLAGS= -O -save -fpp
! LAPACKL = -lmkl -lmkl_lapack
! INCLUDE = -I../camb 
! CFITSIODIR = /opt/intel/cfitsio-3.006
! HEALPIXDIR = /opt/intel/healpix-2.0.1
! 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 = 
! HEALPIXDIR = 
! endif
! 
! ifeq ($(ext),Sun)
! F90C=f90
! FFLAGS= -O4 -xarch=native64 -openmp -ftrap=%none
! WMAPFLAGS= $(FFLAGS)
! LAPACKL = -lsunperf -lfsu
! INCLUDE = -I../camb -M../camb
! CFITSIODIR = 
! HEALPIXDIR = 
! endif
! 
! ifeq ($(ext),AIX)
! F90C    = mpxlf90_r
! FFLAGS  = -O4 -q64 -WF,-DIBMXL,-DMPI -qstrict -qsmp=omp -qmaxmem=-1 -qsuffix=f=f90:cpp=F90
! WMAPFLAGS= $(FFLAGS)
! INCLUDE = -I../camb
! LAPACKL = -lessl -llapack
! CFITSIODIR = 
! HEALPIXDIR = 
! endif
  
+ WMAPDIR = ../WMAP
  PROPOSE = propose.o
  CLSFILE = CMB_Cls_simple.o
  
  #Can use params_H if you prefer more generic parameters
  PARAMETERIZATION = params_CMB.o
  
! WMAPINCLUDE =  -I$(CFITSIODIR)/include -I$(HEALPIXDIR)/include
  
! CLSLIB  = -L$(HEALPIXDIR)/lib -L$(CFITSIODIR)/lib \
! 	-L../camb -lcamb -lhealpix -lcfitsio $(LAPACKL)
  
! F90FLAGS = $(FFLAGS) $(INCLUDE)
  
  DISTFILES = utils.o settings.o GetDist.o
  
! WMAPOBJS = read_archive_map.o read_fits.o WMAP_3yr_options.o WMAP_3yr_util.o \
! 	WMAP_3yr_tt_pixlike.o WMAP_3yr_tt_beam_and_ptsrc_corr.o \
! 	WMAP_3yr_teeebb_pixlike.o WMAP_3yr_likelihood.o
  
  OBJFILES= $(WMAPOBJS) utils.o settings.o cmbtypes.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)
  
  .f.o:
  	f77 $(F90FLAGS) -c $<
  
+ %.o: $(WMAPDIR)/%.f90
+ 	$(F90C) $(WMAPFLAGS) $(WMAPINCLUDE) -c $<
+ 
  %.o: %.f90
! 	$(F90C) $(F90FLAGS) -c $<
  
  %.o: %.F90
! 	$(F90C) $(F90FLAGS) -c $<
  
  
+ cosmomc.$(ext): $(OBJFILES) 	
+ 	$(F90C) -o ../$@ $(OBJFILES) $(CLSLIB) $(F90FLAGS)
  
  clean:
  	rm -f *.o *.mod *.d *.pc *.obj ../core
  
! getdist.$(ext): $(DISTFILES)
! 	$(F90C) -o ../$@ $(DISTFILES) $(CLSLIB) $(F90FLAGS) 
diff -r -c -b -B -N cosmomc/source/MCMC.f90 cosmomc_wiggles/source/MCMC.f90
*** cosmomc/source/MCMC.f90	2006-09-04 07:37:02.000000000 +0100
--- cosmomc_wiggles/source/MCMC.f90	2006-09-13 12:19:49.000000000 +0100
***************
*** 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 -r -c -b -B -N cosmomc/source/params_CMB.f90 cosmomc_wiggles/source/params_CMB.f90
*** cosmomc/source/params_CMB.f90	2006-08-17 13:23:25.000000000 +0100
--- cosmomc_wiggles/source/params_CMB.f90	2006-09-13 12:19:49.000000000 +0100
***************
*** 82,105 ****
         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
   
--- 82,124 ----
         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 -r -c -b -B -N cosmomc/source/params.ini cosmomc_wiggles/source/params.ini
*** cosmomc/source/params.ini	1970-01-01 01:00:00.000000000 +0100
--- cosmomc_wiggles/source/params.ini	2006-09-13 12:19:49.000000000 +0100
***************
*** 0 ****
--- 1,236 ----
+ #Sample parameters for cosmomc in default parameterization
+ 
+ #Root name for files produced
+ file_root = chains/test
+ 
+ #action = 0:  MCMC, action=1: postprocess .data file, action=2: find best fit point only
+ action = 0
+ 
+ #Maximum number of chain steps
+ samples = 200000
+ 
+ #Feedback level ( 2=lots,1=chatty,0=none)
+ feedback = 1
+ 
+ #Temperature at which to Monte-Carlo
+ temperature = 1
+ 
+ #filenames for CMB datasets
+ #Note you may need to change lmax in cmbtypes.f90 to use small scales (e.g. lmax=2100)
+ cmb_numdatasets = 1
+ cmb_dataset1 = WMAP
+ 
+ cmb_dataset2 = data/ACBAR_lge800_.dataset
+ cmb_dataset3 = data/CBI2.0_2_to_14.newdat
+ cmb_dataset4 = data/B03_NA_21July05.newdat
+ 
+ #filenames for matter power spectrum datasets, incl twodf
+ mpk_numdatasets = 1
+ mpk_dataset1 = data/2df_2005.dataset
+ #mpk_dataset1 = data/sdss.dataset
+ 
+ #if true, use HALOFIT for non-linear corrections (astro-ph/0207664).
+ #note lyman-alpha (lya) code assumes linear spectrum
+ nonlinear_pk = F
+ 
+ use_CMB = T
+ use_HST = F
+ use_mpk = F
+ use_clusters = F
+ use_BBN = F
+ use_Age_Tophat_Prior = T
+ use_SN = F
+ use_lya = F
+ use_min_zre = 0
+ 
+ #Force computation of sigma_8 even if use_mpk = F
+ get_sigma8 = F
+ 
+ #1: Simple Metropolis, 2: slice sampling, 3: slice sampling fast parameters, 4: directional gridding
+ sampling_method = 1
+ 
+ #if sampling_method =4, iterations per gridded direction
+ directional_grid_steps = 20
+ 
+ #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
+ #Generally make smaller if your acceptance rate is too low
+ propose_scale = 2.4
+ 
+ #Increase to oversample fast parameters more, e.g. if space is odd shape
+ oversample_fast = 1
+ 
+ #if non-zero number of steps between sample info dumped to file file_root.data
+ indep_sample = 0
+ 
+ #number of samples to disgard at start; usually set to zero and remove later
+ burn_in = 0
+ 
+ #If zero set automatically
+ num_threads = 0
+ 
+ #MPI mode multi-chain options (recommended)
+ #MPI_Converge_Stop is a (variance of chain means)/(mean of variances) parameter that can be used to stop the chains
+ #Set to a negative number not to use this feature. Does not guarantee good accuracy of confidence limits.
+ MPI_Converge_Stop = 0.03
+ 
+ #Do initial period of slice sampling; may be good idea if 
+ #cov matrix or widths are likely to be very poor estimates
+ MPI_StartSliceSampling  = F
+ 
+ #Can optionally also check for convergence of confidence limits (after MPI_Converge_Stop reached)
+ #Can be good idea as small value of MPI_Converge_Stop does not (necessarily) imply good exploration of tails
+ MPI_Check_Limit_Converge = F
+ 
+ #if MPI_Check_Limit_Converge = T, give tail fraction to check (checks both tails):
+ MPI_Limit_Converge = 0.025
+ #permitted quantile chain variance in units of the standard deviation (small values v slow):
+ MPI_Limit_Converge_Err = 0.2
+ #which parameter's tails to check. If zero, check all parameters:
+ MPI_Limit_Param = 0
+ 
+ #if MPI_LearnPropose = T, the proposal density is continally updated from the covariance of samples so far (since burn in)
+ MPI_LearnPropose = T
+ #can set a value of converge at which to stop updating covariance (so that it becomes rigorously Markovian)
+ #e.g. MPI_R_StopProposeUpdate = 0.4 will stop updating when (variance of chain means)/(mean of variances) < 0.4
+ MPI_R_StopProposeUpdate = 0
+ 
+ #If have covmat, R to reach before updating proposal density (increase if covmat likely to be poor)
+ #Only used if not varying new parameters that are fixed in covmat
+ MPI_Max_R_ProposeUpdate = 2
+ #As above, but used if varying new parameters that were fixed in covmat
+ MPI_Max_R_ProposeUpdateNew = 30
+ 
+ #if blank this is set from system clock
+ rand_seed = 
+ 
+ #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)
+ inflation_consistency = F
+ 
+ #Whether the CMB should be lensed (slows a lot unless also computing matter power)
+ CMB_lensing = F
+ #increase accuracy_level to run CAMB on higher accuracy
+ #(default is about 0.3%, accuracy_level=2 around 0.1% at high l)
+ accuracy_level = 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
+ redo_theory = F
+ redo_cls = F
+ redo_pk = F
+ redo_skip = 0
+ redo_outroot = 
+ redo_thin = 1
+ redo_add = F
+ redo_from_text = F
+ #If large difference in log likelihoods may need to offset to give sensible weights
+ #for exp(difference in likelihoods)
+ redo_likeoffset =  0
+ 
+ #parameter start center, min, max, start width, st. dev. estimate
+ #ombh2
+ param1 = 0.0223 0.005 0.1 0.001 0.001
+ #omdmh2
+ param2 = 0.105 0.01 0.99 0.01 0.01
+ #theta
+ param3 = 1.04 0.5 10 0.002 0.002
+ #tau
+ param4 = 0.09 0.01 0.8 0.03 0.03
+ 
+ #omk
+ param5 = 0 0 0 0 0
+ #nufrac
+ param6 = 0 0 0 0 0
+ #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.03 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
+ ##unused fast parameter
+ #param17 = 0 0 0 0 0
+ 
+ ########################################################################################
+ #		        	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
+ #unused
+ param17 = 0 0 0 0 0
diff -r -c -b -B -N cosmomc/source/settings.f90 cosmomc_wiggles/source/settings.f90
*** cosmomc/source/settings.f90	2006-08-11 17:16:55.000000000 +0100
--- cosmomc_wiggles/source/settings.f90	2006-09-13 12:19:49.000000000 +0100
***************
*** 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
  
***************
*** 44,49 ****
--- 54,62 ----
  
    integer :: num_threads = 0
    integer :: instance = 0
+ !wig
+   integer :: instance_shift = 0
+ !end wig
    integer :: MPIchains = 1, MPIrank = 0
  
    logical :: Use_LSS = .true.
