diff --git a/README b/README
new file mode 100644
index 0000000..a220f52
--- /dev/null
+++ b/README
@@ -0,0 +1,76 @@
+**********************************************************************
+Adds trans-Planckian like modulation to primordial power spectra
+expanded at first order in Hubble and sound flow parameters (valid for
+k-inflation). A switch allows the inclusion of second order
+corrections but in the Hubble flow functions only (valid for
+canonically normalised fields).
+**********************************************************************
+
+-These patches modify the cosmomc or camb codes and therefore may have
+bugs that are not present in the originals. The authors of the
+original versions should not be bothered with problems resulting in
+the use of these patches. Feel free to insult me, or at least contact
+me, in case of bugs!
+
+-Apply the patch in the cosmomc/ directory:
+patch -p1 < purpose_MMYY.patch
+
+*Modifications traced by the comment "wig" in the source files.
+
+*Please note that the genuine Makefiles have been heavily changed and
+may require edition to work on your system. They assume the existence
+of a directory "cosmomc/WMAP" which points to the WMAP team likelihood
+code directory. The WMAP data are also assumed to reside in the
+"cosmomc/data" directory.
+
+
+*Modified files and settings 
+
+- camb/cmbmain.f90
+New switches "kSrcSampleBoost" and "kIntSampleBoost" to check accuracy.
+
+- camb/modules.f90 
+New switches "DoClsInterpol" and "ComputeEveryl". AccuracyBoost and
+lAccuracyBoost default values modified. points=10000 to accurately
+sample oscillatory matter power spectra.
+
+- camb/power_sr.f90
+Implementation of Slow-Roll/k-inflation power spectrum and its
+trans-Planckian extension.
+
+- camb/inidriver.F90
+New keys for the above material.
+
+- camb/params.ini
+New keys for the above material.
+
+- cosmomc/CMB_Cls_simple.f90
+Added screen output to check camb accuracy default values. New
+switches kSrcSampleLevel, kIntSampleLevel, CambComputeEveryl to
+control the camb accuracy from cosmomc.
+
+- cosmomc/cmbtypes.f90
+Default value for "write_all_cls" to true.
+
+- cosmomc/params_CMB.f90
+Take care of the new power spectra parameters
+
+- cosmomc/settings.f90
+num_init_power modified
+
+- cosmomc/MCMC.f90
+Add MpiRank output in the logs
+
+- cosmomc/GetDist.f90 
+3d transparent plots of the mean likelihoods. Try colormap('copper').
+
+- cosmomc/driver.F90
+2nd instance in mpirun shifts the file names.
+
+- cosmomc/params.ini
+New keys for the above material.
+
+- cosmomc/distparams.ini
+New keys for the above material.
+
+
diff --git a/camb/Makefile b/camb/Makefile
index 5295384..c9140fd 100644
--- a/camb/Makefile
+++ b/camb/Makefile
@@ -1,77 +1,116 @@
-#CAMB Makefile
-
-#Set FISHER=Y to compile bispectrum fisher matrix code
-FISHER=
-
-#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
-ifneq ($(FISHER),)
-FFLAGS += -mkl
+# >>> 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
+LFLAGS= 
 endif
 
-#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
+ifeq ($(ext),Lin)
+F90C=gfortran
+FFLAGS= -O -fopenmp
+LFLAGS= 
+endif 
 
-#G95 compiler
-#F90C   = g95
-#FFLAGS = -O2
-
-#Gfortran compiler: if pre v4.3 add -D__GFORTRAN__
-#F90C     = gfortran
-#FFLAGS =  -O2  
-
-#SGI, -mp toggles multi-processor. Use -O2 if -Ofast gives problems.
-#F90C     = f90
-#FFLAGS  = -Ofast -mp
+ifeq ($(ext),OSF)
+F90C=f90
+FFLAGS= -omp -O -arch host -math_library fast -tune host -fpe1
+LFLAGS= 
+endif
 
-#Digital/Compaq fortran, -omp toggles multi-processor
-#F90C    = f90
-#FFLAGS  = -omp -O4 -arch host -math_library fast -tune host -fpe1
+ifeq ($(ext),Sun)
+F90C=f90
+FFLAGS= -O4 -xarch=native64 -openmp -ftrap=%none
+LFLAGS= 
+endif
 
-#Absoft ProFortran, single processor:
-#F90C     = f95
-#FFLAGS = -O2 -cpu:athlon -s -lU77 -w -YEXT_NAMES="LCS" -YEXT_SFX="_"
+ifeq ($(ext),AIX)
+F90C=xlf90_r
+FFLAGS= -O4 -qsmp=omp -qmaxmem=-1 -qstrict -qfree=f90 -qsuffix=f=f90:cpp=F90
+LFLAGS= 
+endif
 
-#NAGF95, single processor:
-#F90C     = f95
-#FFLAGS = -DNAGF95 -O3
 
-#PGF90
-#F90C = pgf90
-#FFLAGS = -O2 -DESCAPEBACKSLASH
+#Files containing evolution equations initial power spectrum module
+EQUATIONS     = equations
+POWERSPECTRUM = power_sr
+REIONIZATION = reionization
+RECOMBINATION = recfast
+BISPECTRUM = SeparableBispectrum
+DENABLE_FISHER=
 
-#Sun V880
-#F90C = mpf90
-#FFLAGS =  -O4 -openmp -ftrap=%none -dalign -DMPI
+#Module doing non-linear scaling
+NONLINEAR     = halofit
 
-#Sun parallel enterprise:
-#F90C     = f95
-#FFLAGS =  -O2 -xarch=native64 -openmp -ftrap=%none
-#try removing -openmp if get bus errors. -03, -04 etc are dodgy.
+#Driver program
+DRIVER        = inidriver.F90
+#DRIVER        = sigma8.f90
+#DRIVER        = tester.f90
 
-#IBM XL Fortran, multi-processor (run gmake)
-#F90C     = xlf90_r
-#FFLAGS  = -DESCAPEBACKSLASH -DIBMXL -qsmp=omp -qsuffix=f=f90:cpp=F90 -O3 -qstrict -qarch=pwr3 -qtune=pwr3
 
-#Settings for building camb_fits
-#Location of FITSIO and name of library
-FITSDIR       = /home/cpac/cpac-tools/lib
-FITSLIB       = cfitsio
-#Location of HEALPIX for building camb_fits
-HEALPIXDIR    = /home/cpac/cpac-tools/healpix
+CAMBLIB       = libcamb.a
 
-ifneq ($(FISHER),)
+ifneq ($(DENABLE_FISHER),)
 FFLAGS += -DFISHER
+LFLAGS += -llapack -lblas
 EXTCAMBFILES = Matrix_utils.o
 else
 EXTCAMBFILES = 
 endif
 
-include ./Makefile_main
\ No newline at end of file
+HEALPIXDIR=
+
+#Shouldn't need to change anything else...
+
+F90FLAGS      = $(FFLAGS)
+HEALPIXLD     = -L$(HEALPIXDIR)/lib -lhealpix
+
+CAMBOBJ       = constants.o utils.o $(EXTCAMBFILES) subroutines.o inifile.o $(POWERSPECTRUM).o $(RECOMBINATION).o $(REIONIZATION).o modules.o \
+        bessels.o $(EQUATIONS).o $(NONLINEAR).o lensing.o $(BISPECTRUM).o \
+	cmbmain.o camb.o
+
+default: camb.$(ext)
+
+all: camb.$(ext) $(CAMBLIB)
+
+subroutines.o: constants.o utils.o
+$(POWERSPECTRUM): subroutines.o  inifile.o
+$(RECOMBINATION).o: subroutines.o inifile.o
+$(REIONIZATION).o: constants.o inifile.o
+modules.o: $(REIONIZATION).o $(POWERSPECTRUM).o $(RECOMBINATION).o
+bessels.o: modules.o
+$(EQUATIONS): bessels.o
+$(NONLINEAR).o:  modules.o
+lensing.o: bessels.o
+$(BISPECTRUM).o: lensing.o modules.o
+cmbmain.o: lensing.o $(NONLINEAR).o $(EQUATIONS).o $(BISPECTRUM).o
+camb.o: cmbmain.o
+Matrix_utils.o: utils.o
+
+camb.$(ext): $(CAMBOBJ) $(DRIVER)
+	$(F90C) $(F90FLAGS) $(CAMBOBJ) $(DRIVER)  $(LFLAGS) -o $@
+
+$(CAMBLIB): $(CAMBOBJ)
+	ar -r $@ $?
+
+camb_fits: writefits.f90 $(CAMBOBJ) $(DRIVER)
+	$(F90C) $(F90FLAGS) -I$(HEALPIXDIR)/include tils.o $(CAMBOBJ) writefits.f90 $(DRIVER) $(HEALPIXLD)  $(LFLAGS) -DWRITE_FITS -o $@
+
+%.o: %.f90
+	$(F90C) $(F90FLAGS) -c $*.f90
+
+utils.o:
+	$(F90C) $(F90FLAGS) -c utils.F90	
+
+$(BISPECTRUM).o:
+	$(F90C) $(F90FLAGS) -c $(BISPECTRUM).F90	
+
+Matrix_utils.o:
+	$(F90C) $(F90FLAGS) -c Matrix_utils.F90	
+
+clean:
+	-rm -f *.o *.a *.d core *.mod
+
+
diff --git a/camb/cmbmain.f90 b/camb/cmbmain.f90
index f290c72..3ca0794 100644
--- a/camb/cmbmain.f90
+++ b/camb/cmbmain.f90
@@ -645,17 +645,24 @@ end if
 
          SourceAccuracyBoost = AccuracyBoost     
          if (CP%WantScalars .and. CP%Reion%Reionization .and. CP%AccuratePolarization) then
-            dlnk0=2._dl/10/SourceAccuracyBoost
+!wig
+!            dlnk0=2._dl/10/SourceAccuracyBoost
+            dlnk0=2._dl/10/SourceAccuracyBoost/kSrcSampleBoost
             !Need this to get accurate low l polarization
          else
-            dlnk0=5._dl/10/SourceAccuracyBoost
+!            dlnk0=5._dl/10/SourceAccuracyBoost
+             dlnk0=5._dl/10/SourceAccuracyBoost/kSrcSampleBoost
+!end wig
             if (CP%closed) dlnk0=dlnk0/2
          end if
 
          if (CP%AccurateReionization) dlnk0 = dlnk0/2
-
-         dkn1=0.6_dl/taurst/SourceAccuracyBoost   
-         dkn2=0.9_dl/taurst/SourceAccuracyBoost 
+!wig
+!         dkn1=0.6_dl/taurst/SourceAccuracyBoost   
+!         dkn2=0.9_dl/taurst/SourceAccuracyBoost 
+         dkn1=0.6_dl/taurst/SourceAccuracyBoost/kSrcSampleBoost
+         dkn2=0.9_dl/taurst/SourceAccuracyBoost/kSrcSampleBoost
+!end wig
          if (HighAccuracyDefault) dkn2=dkn2/1.2
          if (CP%WantTensors .or. CP%WantVectors) then
               dkn1=dkn1  *0.8_dl
@@ -1024,8 +1031,10 @@ end if
       
 
        qmax_int = min(qmax,max_bessels_etak/CP%tau0)
-
-       IntSampleBoost=AccuracyBoost  
+!wig
+!       IntSampleBoost=AccuracyBoost  
+       IntSampleBoost=AccuracyBoost*kIntSampleBoost
+!end wig
        if (do_bispectrum) then
         IntSampleBoost = IntSampleBoost * 2  
         if (hard_bispectrum) IntSampleBoost = IntSampleBoost * 2  
diff --git a/camb/inidriver.F90 b/camb/inidriver.F90
index b67ce83..476361b 100644
--- a/camb/inidriver.F90
+++ b/camb/inidriver.F90
@@ -270,6 +270,20 @@
          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'
@@ -281,6 +295,23 @@
        if (FeedbackLevel > 0) then
          Age = CAMB_GetAge(P) 
          write (*,'("Age of universe/GYr  = ",f7.3)') Age  
+!wig
+         if ((Power_Name.eq.'power_sr').or.(Power_Name.eq.'power_dbi')) then
+
+             do i=1,P%InitPower%nn
+                write (*,'("Slow-roll epsilon1   = ",f12.8)') P%InitPower%eps1(i)
+                write (*,'("Slow-roll epsilon2   = ",f12.8)') P%InitPower%eps2(i)
+		write (*,'("Slow-roll epsilon3   = ",f12.8)') P%InitPower%eps3(i)
+                write (*,'("DBI       delta1     = ",f12.8)') P%InitPower%del1(i)
+                write (*,'("DBI       gamma      = ",f12.8)') P%InitPower%gam(i)
+                write (*,'("TPL amplitude        = ",f12.8)') P%InitPower%Xsigma0(i)
+                write (*,'("TPL scale            = ",f12.8)') P%InitPower%loginvsigma0(i)
+                write (*,'("TPL phase            = ",f12.8)') P%InitPower%varphi(i)
+             enddo
+          else
+             stop 'Power_Name unknown'
+          endif
+!end wig
        end if 
 
        call CAMB_GetResults(P)
diff --git a/camb/modules.f90 b/camb/modules.f90
index d3a041e..63c3e25 100644
--- a/camb/modules.f90
+++ b/camb/modules.f90
@@ -37,8 +37,10 @@
         character(LEN=*), parameter :: version = 'July_11'
         
         integer :: FeedbackLevel = 0 !if >0 print out useful information about the model
-
-        logical, parameter :: DebugMsgs=.false. !Set to true to view progress and timing
+!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
  
@@ -186,15 +188,17 @@
 
       real(dl) :: lSampleBoost=1._dl
           !Increase lSampleBoost to increase sampling in lSamp%l for Cl interpolation
-          
-      real(dl) :: AccuracyBoost =1._dl  
 
+!wig          
+      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
-
-      real(sp) :: lAccuracyBoost=1. 
+!wig
+      real(sp) :: lAccuracyBoost=2._dl
+!end wig
           !Boost number of multipoles integrated in Boltzman heirarchy
 
       integer, parameter :: lmin = 2  
@@ -211,7 +215,13 @@
 
 !     lmax is max possible number of l's evaluated
       integer, parameter :: lmax_arr = l0max 
- 
+
+!wig
+      real(dl) :: kSrcSampleBoost=1._dl
+      real(dl) :: kIntSampleBoost=1._dl
+      logical :: DoClsInterpol=.true.
+      logical :: ComputeEveryl=.true.
+!end wig 
         contains
       
 
@@ -573,6 +583,25 @@
         integer lind, lvar, step,top,bot,ls(lmax_arr)
         real(dl) AScale
     
+!wig        
+        if (ComputeEveryl) then
+           write(*,*) 'Warning: computing every l'           
+           if (max_l.gt.lmax_arr) then
+              write(*,*)'lmax_arr maxl',lmax_arr,max_l
+              stop
+           endif
+           lind=0
+           do lvar=2, max_l
+              lind=lind+1
+              ls(lind)=lvar 
+           end do
+
+           lSet%l0=lind
+           lSet%l = ls(1:lind)
+
+        else
+!end wig
+
         Ascale=scale/lSampleBoost       
        
         if (lSampleBoost >=50) then
@@ -708,6 +737,9 @@
         end if
         lSet%l0=lind
         lSet%l(1:lind) = ls(1:lind)
+!wig        
+        endif
+!end wig
         
       end subroutine initlval
 
@@ -727,7 +759,8 @@
 
       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
@@ -743,6 +776,20 @@
                        +(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
 
@@ -1949,6 +1996,10 @@
              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 --git a/camb/params.ini b/camb/params.ini
index 95893cc..a071b35 100644
--- a/camb/params.ini
+++ b/camb/params.ini
@@ -1,12 +1,12 @@
 #Parameters for CAMB
 
 #output_root is prefixed to output file names
-output_root = test
+output_root = tpl
 
 #What to do
 get_scalar_cls = T
 get_vector_cls = F
-get_tensor_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}
@@ -65,18 +65,27 @@ nu_mass_degeneracies = 0
 #Fraction of total omega_nu h^2 accounted for by each eigenstate, eg. 0.5 0.5
 nu_mass_fractions = 1
 
-#Initial power spectrum, amplitude, spectral index and running. Pivot k in Mpc^{-1}.
-initial_power_num         = 1
-pivot_scalar              = 0.05
-pivot_tensor              = 0.05
-scalar_amp(1)             = 2.1e-9
-scalar_spectral_index(1)  = 0.96
-scalar_nrun(1)            = 0
-tensor_spectral_index(1)  = 0
-#ratio is that of the initial tens/scal power spectrum amplitudes
-initial_ratio(1)          = 1
-#note vector modes use the scalar settings above
+#Type of initial power spectrum: power_sr or power_dbi
+initial_power_num = 1
+power_name = power_sr
+ 
+#slow-roll parameters
+ln[10^{10} P_scalar] = 0.303161E+01
+eps1(1) = 0.210391E-02
+eps2(1) = 0.421584E-01
+eps3(1) = 0
 
+#dbi parameter
+gam(1) = 1.
+del1(1) = 0.
+
+#wiggles parameters
+logeps1invsigma0(1) = 1.5
+Xsigma0(1) = 0.1
+varphi(1) = 0
+
+#tpl pivot
+k_star = 0.05
 
 #Reionization, ignored unless reionization = T, re_redshift measures where x_e=0.5
 reionization         = T
@@ -98,6 +107,7 @@ RECFAST_fudge_He = 0.86
 RECFAST_Heswitch = 6
 RECFAST_Hswitch  = T
 
+
 #Initial scalar perturbation mode (adiabatic=1, CDM iso=2, Baryon iso=3, 
 # neutrino density iso =4, neutrino velocity iso = 5) 
 initial_condition   = 1
@@ -206,7 +216,7 @@ high_accuracy_default=F
 
 #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
+accuracy_boost          = 4
 
 #Larger to keep more terms in the hierarchy evolution. 
 l_accuracy_boost        = 1
@@ -216,3 +226,15 @@ l_accuracy_boost        = 1
 #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 --git a/camb/power_sr.f90 b/camb/power_sr.f90
new file mode 100644
index 0000000..5a518d2
--- /dev/null
+++ b/camb/power_sr.f90
@@ -0,0 +1,544 @@
+!This module provides the initial power spectra, parameterized as an expansion in ln k
+!
+! ln P = ln A_s + (n_s -1)*ln(k/k_0) + n_{run}/2 * ln(k/k_0)^2 
+!
+! so if n_run = 0
+!
+! P = A_s (k/k_0_scalar)^(n_s-1)
+!
+!for the scalar spectrum, when an(in) is the in'th spectral index. k_0_scalar
+!is a pivot scale, fixed here to 0.05/Mpc (change it below as desired).
+!
+!This module uses the same inputs an(in), ant(in) and rat(in) as CMBFAST, however here
+!rat(in) is used to set the ratio of the initial power spectra, so here
+!
+!** rat(in) is not the Cl quadrupole ratio ***
+!
+!in general models the quadrupole ratio depends in a complicated way on the ratio of the initial
+!power spectra
+
+!The absolute normalization of the Cls is unimportant here, but the relative ratio
+!of the tensor and scalar Cls generated with this module will be correct for general models
+
+
+!The OutputNormalization parameter controls the final output
+!Absolute Cls can be obtained by setting OuputNormalization=outNone, otherwise the overall normalization
+!of the power spectra doesn't matter
+
+!This version December 2003 - changed default tensor pivot to 0.05 (consistent with CMBFAST 4.5)
+
+     module InitialPower   
+     use Precision
+     implicit none   
+
+     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
+
+      Type InitialPowerParams
+
+       integer nn  !Must have this variable
+      !The actual number of power spectra to use
+  
+      !For the default implementation return power spectra based on spectral indices
+       real(dl) an(nnmax) !scalar spectral indices
+       real(dl) n_run(nnmax) !running of spectral index 
+       real(dl) ant(nnmax) !tensor spectral indices
+       real(dl) rat(nnmax) !ratio of scalar to tensor initial power spectrum amplitudes
+       real(dl) k_0_scalar, k_0_tensor
+       real(dl) ScalarPowerAmp(nnmax)
+ 
+!wig
+! slow-roll power
+       real(dl) P_scalar
+       real(dl) k_star
+       real(dl) eps1(nnmax)
+       real(dl) eps2(nnmax)
+       real(dl) eps3(nnmax)
+!dbi params
+       real(dl) gam(nnmax)      
+       real(dl) del1(nnmax)
+! tpl power
+       real(dl) invsigma0(nnmax)
+       real(dl) loginvsigma0(nnmax)
+       real(dl) logeps1invsigma0(nnmax)
+       real(dl) Xsigma0(nnmax)
+       real(dl) varphi(nnmax) 
+!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     
+      
+      Type(InitialPowerParams) :: P
+  
+!Make things visible as neccessary...
+ 
+      public InitialPowerParams, InitialPower_ReadParams, InitializePowers, ScalarPower, TensorPower
+      public nnmax,Power_Descript, Power_Name, SetDefPowerParams
+!      public 
+    contains
+       
+  
+       subroutine SetDefPowerParams(AP)
+        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
+         AP%rat    = 1
+         AP%k_0_scalar = 0.05
+         AP%k_0_tensor = 0.05
+         AP%ScalarPowerAmp = 1      
+         
+         AP%eps1 = 0.01
+         AP%eps2 = 0.01 
+         AP%eps3 = 0.
+         AP%P_scalar=1.d-9
+         AP%k_star = 0.05            
+
+         AP%del1 = 0.
+         AP%gam = 1.
+        
+
+         AP%Xsigma0 = 0.                
+         AP%logeps1invsigma0 = 0. 
+         AP%loginvsigma0 = AP%logeps1invsigma0 - log10(AP%eps1) 
+         AP%invsigma0=10**(AP%loginvsigma0)
+         AP%varphi=0.
+        
+!end wig
+
+       end subroutine SetDefPowerParams
+
+       subroutine InitializePowers(AParamSet,acurv)
+         Type (InitialPowerParams) :: AParamSet
+         !Called before computing final Cls in cmbmain.f90
+         !Could read spectra from disk here, do other processing, etc.
+
+        real(dl) acurv
+!wig
+        integer :: i        
+!from Hankel functions
+        real(dl), parameter :: C_const = -0.7296
+!wkb or uniform approx gives a slightly different constant
+        real(dl), parameter :: D_const = -0.7652
+!        real(dl), parameter :: C_const = D_const
+
+        real(dl), parameter :: pi = 3.14159265
+!end wig
+
+        if (AParamSet%nn > nnmax) then
+           write (*,*) 'To use ',AParamSet%nn,'power spectra you need to increase'
+           write (*,*) 'nnmax in power_tilt.f90, currently ',nnmax
+        end if
+        P = AParamSet
+
+        curv=acurv         
+
+!Write implementation specific code here...
+!wig
+  
+        if ((Power_Name.eq.'power_sr').or.(Power_Name.eq.'power_dbi')) then
+           if(abs(curv) .gt. 1.d-3) then
+              write(*,*)'Non flat cosmology, curv = ',curv
+              write(*,*)'Cannot use SR or TPL spectrum.'
+              stop
+           endif
+
+           P%ks_pivot = P%k_star
+           P%kt_pivot = P%k_star
+           P%ko_pivot = P%k_star
+
+
+           do i=1, P%nn        
+
+!TODO: derive second order tpl corrections?
+              if (useSecondOrder) then
+                 if (P%Xsigma0(i).ne.0.) then
+                    useSecondOrder=.false.
+                    write(*,*)'PowerIni: SR second order disabled: TPL not implemented! '
+                 endif
+                 if ((P%del1(i).ne.0.).or.(P%gam(i).ne.1.)) then
+                    useSecondorder = .false.
+                    write(*,*)'PowerIni: SR second order disabled: DBI slow-roll not implemented!'
+                 endif
+              endif
+
+              if ((power_Name.eq.'power_sr').and. &
+                   ((P%del1(i).ne.0.).or.(P%gam(i).ne.1.)) ) then
+                 write(*,*)'power_sr with non-vanishing DBI params!'
+                 write(*,*)'resetting to standard SR'
+                 P%del1(i) = 0.
+                 P%gam(i) = 1.
+              endif
+
+
+!first order tpl corrections
+              
+              P%loginvsigma0(i) = P%logeps1invsigma0(i) - log10(P%eps1(i))
+              P%invsigma0(i)=10**(P%loginvsigma0(i))
+              
+              P%amps0(i) = 2.*P%Xsigma0(i)
+              P%amps1(i) = pi*P%Xsigma0(i)*(2.*P%eps1(i) + P%eps2(i))
+              P%amps2(i) = 0.
+              if (useTplPhase) then
+                 P%args0(i) = 2.*P%invsigma0(i)*(1. + P%eps1(i)) + P%varphi(i)
+              else
+                 P%args0(i) = P%varphi(i)
+              endif
+              P%args1(i) = 2.*P%invsigma0(i)*P%eps1(i)
+              P%args2(i) = 0.
+
+              P%ampt0(i) = P%amps0(i)
+              P%ampt1(i) = pi*P%Xsigma0(i)*2.*P%eps1(i)
+              P%ampt2(i) = 0.
+              P%argt0(i) = P%args0(i)
+              P%argt1(i) = P%args1(i)
+              P%argt2(i) = 0.
+
+
+!first and second order slowroll
+
+              P%as0(i) = 1. + P%eps1(i)*(-2.*(C_const+1.))-P%eps2(i)*C_const &
+                   + (C_const + 2.)*P%del1(i)
+              P%as1(i) = -2.*P%eps1(i) - P%eps2(i) + P%del1(i)
+              P%as2(i) = 0.
+                            
+              P%at0(i) = 1. + P%eps1(i)*(-2.*(C_const+1.))
+              P%at1(i) = -2.*P%eps1(i)
+              P%at2(i) = 0.
+
+              if (useSecondOrder) then
+                                                   
+
+                 P%as0(i) = P%as0(i) &
+                      + P%eps1(i)*P%eps1(i) * (2.*C_const**2 + 2.*C_const + pi**2/2. - 5.) &
+                      + P%eps1(i)*P%eps2(i) * (C_const**2 - C_const + 7.*pi**2/12. - 7.) &
+                      + P%eps2(i)*P%eps2(i) * (0.5*C_const**2 + pi**2/8. - 1.) &
+                      + P%eps2(i)*P%eps3(i) * (-0.5*C_const**2 + pi**2/24.)
+
+                 P%as1(i) = P%as1(i) &
+                      + P%eps1(i)*P%eps1(i) * 2.*(2.*C_const + 1.) &
+                      + P%eps1(i)*P%eps2(i) * (2.*C_const - 1.) &
+                      + P%eps2(i)*P%eps2(i) * C_const &
+                      + P%eps2(i)*P%eps3(i) * (-C_const)
+
+                 P%as2(i) = P%as2(i) &
+                      + P%eps1(i)*P%eps1(i) * 4. &
+                      + P%eps1(i)*P%eps2(i) * 2. &
+                      + P%eps2(i)*P%eps2(i) * 1. &
+                      + P%eps2(i)*P%eps3(i) * (-1.)
+
+                 
+                 P%at0(i) = P%at0(i) &
+                      +  P%eps1(i)*P%eps1(i) * (2.*C_const**2 + 2.*C_const + pi**2/2. - 5.) &
+                      +  P%eps1(i)*P%eps2(i) * (-C_const**2 - 2.*C_const + pi**2/12. - 2.)
+                 
+                 P%at1(i) = P%at1(i) &
+                      + P%eps1(i)*P%eps1(i) * 2.*(2.*C_const + 1.) &
+                      + P%eps1(i)*P%eps2(i) * 2.*(-C_const - 1.)
+
+                 P%at2(i) = P%at2(i)  &
+                      + P%eps1(i)*P%eps1(i) * 4. &
+                      + P%eps1(i)*P%eps2(i) * (-2.)
+                      
+              endif
+
+              
+              !Relative normalisation of scalar and tensor spectra
+                 
+              P%P_0_scalar(i)=P%P_scalar/(P%as0(i))
+!              P%P_0_tensor(i)=P%P_scalar/(P%as0(i))*P%eps1(i)*16.
+              P%P_0_tensor(i)=P%P_scalar/(P%as0(i))*P%eps1(i)*16/P%gam(i)
+
+           enddo
+
+        else
+           stop 'No such initial power spectra'
+        endif
+!end wig
+
+       end subroutine InitializePowers
+       
+
+      function ScalarPower(k,in)
+
+       !"in" gives the index of the power to return for this k
+       !ScalarPower = const for scale invariant spectrum
+       !The normalization is defined so that for adiabatic perturbations the gradient of the 3-Ricci 
+       !scalar on co-moving hypersurfaces receives power
+       ! < |D_a R^{(3)}|^2 > = int dk/k 16 k^6/S^6 (1-3K/k^2)^2 ScalarPower(k) 
+       !In other words ScalarPower is the power spectrum of the conserved curvature perturbation given by
+       !-chi = \Phi + 2/3*\Omega^{-1} \frac{H^{-1}\Phi' - \Psi}{1+w}
+       !(w=p/\rho), so < |\chi(x)|^2 > = \int dk/k ScalarPower(k).
+       !Near the end of inflation chi is equal to 3/2 Psi.
+       !Here nu^2 = (k^2 + curv)/|curv| 
+
+       !This power spectrum is also used for isocurvature modes where 
+       !< |\Delta(x)|^2 > = \int dk/k ScalarPower(k)
+       !For the isocurvture velocity mode ScalarPower is the power in the neutrino heat flux.
+
+        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_sr').or.(Power_Name.eq.'power_dbi')) then
+
+           angle = P%args0(in) + P%args1(in) * log(k/P%ko_pivot) 
+
+           ScalarPower = (P%as0(in) + P%as1(in)*lnrat + 0.5*P%as2(in)*lnrat**2) &
+                * (1. - P%amps0(in)*cos(angle)) - P%amps1(in)*sin(angle)
+
+        else
+           stop 'No such initial power spectra'
+        endif
+
+        ScalarPower = ScalarPower*P%P_0_scalar(in)
+
+        if(ScalarPower .lt. 0.0) then
+           ScalarPower =0.
+        end if
+
+
+!end wig
+     
+!         ScalarPower = ScalarPower * (1 + 0.1*cos( lnrat*30 ) )
+ 
+      end function ScalarPower
+
+      
+      function TensorPower(k,in)
+      
+       !TensorPower= const for scale invariant spectrum
+       !The normalization is defined so that
+       ! < h_{ij}(x) h^{ij}(x) > = \sum_nu nu /(nu^2-1) (nu^2-4)/nu^2 TensorPower(k)
+       !for a closed model
+       ! < h_{ij}(x) h^{ij}(x) > = int d nu /(nu^2+1) (nu^2+4)/nu^2 TensorPower(k)
+       !for an open model
+       !"in" gives the index of the power spectrum to return 
+       !Here nu^2 = (k^2 + 3*curv)/|curv| 
+
+
+        real(dl) TensorPower,k   
+        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_sr').or.(Power_Name.eq.'power_dbi')) then
+ 
+            angle = P%argt0(in) + P%argt1(in) * log(k/P%ko_pivot)
+ 
+            TensorPower = (P%at0(in) + P%at1(in)*lnrat + 0.5*P%at2(in)*lnrat**2) &
+                 * (1. - P%ampt0(in)*cos(angle)) - P%ampt1(in)*sin(angle)
+ 
+         else
+            stop 'No such initial power spectra'
+         endif
+ 
+         TensorPower = TensorPower*P%P_0_tensor(in)
+  
+        if (curv < 0) TensorPower=TensorPower*tanh(PiByTwo*sqrt(-k**2/curv-3)) 
+
+         if(TensorPower .lt. 0.0) then
+            TensorPower=0.
+         endif
+         
+!end wig
+
+       
+      end function TensorPower
+
+      !Get parameters describing parameterisation (for FITS file)
+     function Power_Descript(in, Scal, Tens, Keys, Vals)
+         character(LEN=8), intent(out) :: Keys(*)
+         real(dl), intent(out) :: Vals(*)
+         integer, intent(IN) :: in
+         logical, intent(IN) :: Scal, Tens
+         integer num, Power_Descript
+         num=0
+!wig
+        
+
+         if (Power_Name.eq.'power_sr') then
+
+            num=num+1
+            Keys(num) = 'k_star'
+            Vals(num) = P%k_star
+            num=num+1 
+            Keys(num) = 'P_scalar'
+            Vals(num) = P%P_scalar     
+            num=num+1 
+            Keys(num) = 'eps1'
+            Vals(num) = P%eps1(in)
+            num=num+1 
+            Keys(num) = 'eps2'
+            Vals(num) = P%eps2(in)
+            num=num+1 
+            Keys(num) = 'eps3'
+            Vals(num) = P%eps3(in)
+            num=num+1 
+            Keys(num) = 'isigma0'
+            Vals(num) = P%invsigma0(in)
+            num=num+1
+            Keys(num) = 'lisigma0'
+            Vals(num) = P%loginvsigma0(in)
+            num=num+1 
+            Keys(num) = 'Xsigma0'
+            Vals(num) = P%Xsigma0(in)
+            num=num+1 
+            Keys(num) = 'varphi'
+            Vals(num) = P%varphi(in)
+
+         elseif (Power_Name.eq.'power_dbi') then
+
+            num=num+1
+            Keys(num) = 'k_star'
+            Vals(num) = P%k_star
+            num=num+1 
+            Keys(num) = 'P_scalar'
+            Vals(num) = P%P_scalar     
+            num=num+1 
+            Keys(num) = 'eps1'
+            Vals(num) = P%eps1(in)
+            num=num+1 
+            Keys(num) = 'eps2'
+            Vals(num) = P%eps2(in)
+            num=num+1 
+            Keys(num) = 'eps3'
+            Vals(num) = P%eps3(in)
+            num=num+1 
+            Keys(num) = 'gam'
+            Vals(num) = P%gam(in)
+            num=num+1 
+            Keys(num) = 'del1'
+            Vals(num) = P%del1(in)
+            num=num+1 
+            Keys(num) = 'isigma0'
+            Vals(num) = P%invsigma0(in)
+            num=num+1
+            Keys(num) = 'lisigma0'
+            Vals(num) = P%loginvsigma0(in)
+            num=num+1 
+            Keys(num) = 'Xsigma0'
+            Vals(num) = P%Xsigma0(in)
+            num=num+1 
+            Keys(num) = 'varphi'
+            Vals(num) = P%varphi(in)
+
+         else
+            stop 'No such initial power spectra'
+         endif
+!end wig
+         Power_Descript = num
+
+       end  function Power_Descript
+        
+       subroutine InitialPower_ReadParams(InitPower, Ini, WantTensors)
+          use IniFile
+          Type(InitialPowerParams) :: InitPower
+          Type(TIniFile) :: Ini
+          logical, intent(in) :: WantTensors
+          integer i
+          
+           InitPower%k_0_scalar = Ini_Read_Double_File(Ini,'pivot_scalar',InitPower%k_0_scalar)
+           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_sr').or.(Power_Name.eq.'power_dbi')) then
+           
+              do i=1, InitPower%nn
+              
+!                 if (P%WantScalars .or. P%WantTransfer) then
+
+                 InitPower%eps1(i)= Ini_Read_Double_Array_File(Ini,'eps1',i,0.01d0)
+                 InitPower%eps2(i)=Ini_Read_Double_Array_File(Ini,'eps2',i,0.01d0)
+                 InitPower%eps3(i)=Ini_Read_Double_Array_File(Ini,'eps3',i,0.00d0)
+                 InitPower%gam(i)=Ini_Read_Double_Array_File(Ini,'gam',i,1.00d0)
+                 InitPower%del1(i)=Ini_Read_Double_Array_File(Ini,'del1',i,0.00d0)
+                 InitPower%logeps1invsigma0(i) &
+                      = Ini_Read_Double_Array_File(Ini,'logeps1invsigma0',i,1d0)
+                 InitPower%loginvsigma0(i) = InitPower%logeps1invsigma0(i) &
+                      - log10(InitPower%eps1(i))
+                 InitPower%invsigma0(i) = 10**(InitPower%loginvsigma0(i))
+                 InitPower%Xsigma0(i) = Ini_Read_Double_Array_File(Ini,'Xsigma0',i,0d0)
+                 InitPower%varphi(i) = Ini_Read_Double_Array_file(Ini,'varphi',i,0.d0)
+
+!                 endif
+              
+              enddo
+           
+              InitPower%k_star= Ini_Read_Double('k_star',0.05d0)
+              InitPower%P_scalar= 1e-10*exp(Ini_Read_Double('ln[10^{10} P_scalar]',1d0))
+
+           else
+
+              stop 'No such initial power spectrum'
+
+           endif
+ 
+!end wig
+
+
+              
+
+
+          
+       end  subroutine InitialPower_ReadParams 
+
+
+     end module InitialPower
diff --git a/distparams.ini b/distparams.ini
index 95e67d5..be37ce6 100644
--- a/distparams.ini
+++ b/distparams.ini
@@ -49,6 +49,7 @@ 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
@@ -94,7 +95,7 @@ matlab_colscheme =
 
 #Parameters to use. If zero use all parameters which have lables.
 plotparams_num = 0
-plotparams = omegabh2 omegadmh2 tau ns
+plotparams = omegabh2 omegadmh2 tau ns eps1 eps2 eps3 logomega amp psi logsig
 
 #Get set label to empty to not include a parameter in the parameter_names file
 #lab[asz]=
diff --git a/params.ini b/params.ini
index c30e575..a36ad2e 100644
--- a/params.ini
+++ b/params.ini
@@ -67,18 +67,18 @@ 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 = F
+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
+#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
+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
@@ -142,7 +142,7 @@ highL_unlensed_cl_template = ./camb/HighLExtrapTemplate_lenspotentialCls.dat
 
 #CAMB parameters
 #If we are including tensors
-compute_tensors = F
+compute_tensors = T
 #Initial power spectrum amplitude point (Mpc^{-1})
 pivot_k = 0.05
 #If using tensors, enforce n_T = -A_T/(8A_s)
@@ -157,7 +157,12 @@ high_accuracy_default = F
 #increase accuracy_level to run CAMB on higher accuracy
 #(default is about 0.3%, 0.1% if high_accuracy_default =T)
 #accuracy_level can be increased to check stability/higher accuracy, but inefficient
-accuracy_level = 1
+accuracy_level = 4
+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
@@ -183,12 +188,19 @@ param[omegak] = 0 0 0 0 0
 param[fnu] = 0 0 0 0 0
 param[w] = -1 -1 -1 0 0
 
-param[ns] = 0.95 0.5 1.5 0.02 0.01
-param[nt] = 0 0 0 0 0
-param[nrun] = 0 0 0 0 0
+param[sr1] = 0.1 0.00001 0.2 0.01 0.02
+param[sr2] = -0.04 -0.2 0.2  0.01 0.02
+param[sr3] = -0.04 -0.2 0.2  0.01 0.02
+param[unused] = 0 0 0 0 0
 
-#log[10^10 A_s]
-param[logA] = 3 2.7 4 0.01 0.01
-param[r] = 0 0 0 0 0
+param[logfreq] = 1.5 1 2.5 0.5 0.2
+param[amposc] = 0.05 0 0.45 0.1 0.2
+param[phase] = 0 0 6.28 0.4 0.1
+
+param[logsnd] = 0 0 0 0 0
+
+#ln[10^10 A_s]
+param[lnA] = 3 2.7 4 0.01 0.01
+param[r] = 1 1 1 0 0
 #SZ amplitude, as in WMAP analysis
 param[asz]= 1 0 2 0.4 0.4
diff --git a/params_CMB.paramnames b/params_CMB.paramnames
index 0c08e67..8fa3488 100644
--- a/params_CMB.paramnames
+++ b/params_CMB.paramnames
@@ -1,20 +1,25 @@
-omegabh2      \Omega_b h^2	#physical baryon density
-omegadmh2     \Omega_{DM} h^2	#physical dark matter density, including CDM and massive neutrinos
-theta         \theta		#100 times the ratio of the angular diameter distance to the LSS sound horizon (approx)
+omegabh2      \Omega_b h^2		#physical baryon density
+omegadmh2     \Omega_{DM} h^2		#physical dark matter density, including CDM and massive neutrinos
+theta         \theta			#100 times the ratio of the angular diameter distance to the LSS sound horizon (approx)
 tau           \tau
 omegak        \Omega_K
-fnu           f_\nu		#neutrino energy density as fraction of omegadmh2
-w             w			#constant equation of state parameter for scalar field dark energy 
-ns            n_s		#beware that pivot scale can change in .ini file
-nt            n_t
-nrun          n_{run}
-logA          log[10^{10} A_s]
-r             r			#ratio of tensor to scalar primordial amplitudes at pivot scale
-asz           A_{SZ}		#SZ template amplitude, as in WMAP
+fnu           f_\nu			#neutrino energy density as fraction of omegadmh2
+w             w				#constant equation of state parameter for scalar field dark energy 
+sr1	      \epsilon_1-\delta_1	#first slow roll parameter - first sound flow param
+sr2	      \epsilon_2+\delta_1	#second slow roll parameter + first sound flow param
+sr3	      \epsilon_3		#third slow roll parameter
+unused	      log(\epsilon_1)		#unused
+logfreq	      log(\epsilon_1/\sigma_0)	#superimposed oscillation frequency
+amposc	      x\sigma_0			#superimposed oscillation amplitude
+phase	      \psi			#superimposed oscillation phase shift
+logsnd	      log(\epsilon_1 c_s)	#speed of sound
+lnA           ln[10^{10} P_*]
+r             r				#ratio of tensor to scalar primordial amplitudes at pivot scale
+asz           A_{SZ}			#SZ template amplitude, as in WMAP
 omegal*       \Omega_\Lambda
 age*          Age/Gyr
 omegam*       \Omega_m
 sigma8*       \sigma_8
 zrei*         z_{re}
-r10*          r_{10}		#tensor-scalar C_l amplitude at l=10
-H0*           H_0		#hubble parameter is H0 km/s/Mpc
\ No newline at end of file
+r10*          r_{10}			#tensor-scalar C_l amplitude at l=10
+H0*           H_0			#hubble parameter is H0 km/s/Mpc
diff --git a/source/CMB_Cls_simple.f90 b/source/CMB_Cls_simple.f90
index bb64245..d28e789 100644
--- a/source/CMB_Cls_simple.f90
+++ b/source/CMB_Cls_simple.f90
@@ -394,6 +394,28 @@ contains
    integer zix
    real redshifts(matter_power_lnzsteps)
 
+!wig
+   AccuracyBoost = AccuracyLevel
+   lAccuracyBoost = lAccuracyLevel
+   lSampleBoost = lSampleLevel
+   kSrcSampleBoost = kSrcSampleLevel
+   kIntSampleBoost = kIntSampleLevel
+   ComputeEveryl = CambComputeEveryl
+
+   write(*,*)
+   write(*,*) '>>>>>>>>>>>>>> accuracy settings <<<<<<<<<<<<<<'
+   write(*,*) 'Power_Name     =    ',Power_Name
+   write(*,*) 'write_all_cls  =   ',write_all_cls
+   write(*,*) 'AccuracyBoost  = ',AccuracyBoost
+   write(*,*) 'lAccuracyBoost = ',lAccuracyBoost          
+   write(*,*) 'DoClsInterpol  =   ',DoClsInterpol
+   write(*,*) 'ComputeEveryl  =   ',ComputeEveryl
+   write(*,*) '>>>>>>>>>>>>>> accuracy settings <<<<<<<<<<<<<<'
+   write(*,*)
+
+   if ((Power_Name == 'power_sr').and.(.not.compute_tensors)) stop 'Slow-Roll requires tensors'
+!end wig
+
         Threadnum =num_threads
         w_lam = -1
         call CAMB_SetDefParams(P)
@@ -426,9 +448,11 @@ contains
             P%Transfer%high_precision=.true.
             P%Transfer%kmax=P%Transfer%kmax + 0.2
            end if
-          AccuracyBoost = AccuracyLevel
-          lAccuracyBoost = AccuracyLevel
-          lSampleBoost = AccuracyLevel
+!wig
+!          AccuracyBoost = AccuracyLevel
+!          lAccuracyBoost = AccuracyLevel
+!          lSampleBoost = AccuracyLevel
+!end wig
           P%AccurateReionization = .true.
         end if
                 
diff --git a/source/GetDist.f90 b/source/GetDist.f90
index b4094ce..11fb2bf 100644
--- a/source/GetDist.f90
+++ b/source/GetDist.f90
@@ -104,6 +104,32 @@ module MCSamples
         integer Num_ComparePlots
         logical ::   prob_label = .false.      
             
+!addon
+!exact field model
+        logical, parameter :: findMaxValues = .false.
+        logical, parameter :: mapRhoReh = .false.
+        logical, parameter :: mapLFPowerToWreh = .true.
+
+!slow-roll analysis
+        logical, parameter :: adjEps1ToLog = .false.
+        logical, parameter :: extractHinf = .false.
+        logical, parameter :: mapnslogr = .true.
+        logical, parameter :: use2ndSR=.true.        
+        real, parameter :: ln10p5HinfCut=-10
+        real, parameter :: logrCut=-3.7
+
+#ifdef SRREHEAT
+        logical, parameter :: doReheatFromSr = .false.
+        character(len=*), parameter :: reheatModel='smallf'
+#endif
+
+        real, parameter :: wfix = -0.3
+        real, parameter :: pmax = 15, pmin = 0
+        real, parameter :: logmumax=10
+!end addon
+
+
+    
 contains
 
   subroutine AdjustPriors
@@ -111,8 +137,20 @@ contains
    !Be careful as this code is parameterisation dependent
    integer i
    real ombh2, chisq 
+!addon
+   integer, parameter :: dp = kind(1._8)
+   real(dp) :: eps1,logeps1,eps2,eps3, lnR,lnRrad
+   real(dp) :: wreh,p,logmu,lnRhoReh,lnRhoEnd   
+
+   real :: ln10p5HinfMin=1e30
+   real :: ln10p5HinfMax=-1e30
+   
+   real :: lnRhoRehmax = -1e30, lnRhoEndMax=-1e30, lnRmax = -1e30, lnRradmax = -1e30
+
+   real, parameter :: lnRhoNuc = -187.747
  
-   stop 'You need to write the AdjustPriors subroutine in GetDist.f90 first!'
+!   stop 'You need to write the AdjustPriors subroutine in GetDist.f90 first!'
+!end addon 
   
 !   write (*,*) 'Adjusting priors'
 !   do i=0, nrows-1
@@ -124,6 +162,214 @@ contains
 !    
 !   end do
 
+!addon
+   wreh = wfix
+
+   
+   if (findMaxValues) then
+      print *,'-->computing maximum values in chains'
+
+      do i=0,nrows-1
+         lnR = coldata(14+2,i)
+         lnRhoEnd = coldata(2+30,i)
+         lnRrad = coldata(2+28,i)
+         lnRhoEndMax = max(lnRHoEndMax,lnRhoEnd)
+         lnRmax = max(lnRmax, lnR)
+         lnRradmax = max(lnRradmax, lnRrad)
+      enddo
+
+      print *,'lnRhoEndMax= ',lnRhoEndMax
+      print *,'lnRradMax= lnRmax= ',lnRradMax, lnRmax
+   endif
+
+   if (mapRhoReh) then
+      print *,'--> rejecting all p <',pmin
+      print *,'--> rejecting all logmu >',logmumax
+      print *,'--> rejecting all lnRhoNuc > lnRhoReh > lnRhoEnd'
+
+      do i=0,nrows-1
+         p = coldata(2+10,i)
+         logmu = coldata(11+2,i)
+
+         if (mapLFPowerToWreh) wreh = (p-2._dp)/(p+2._dp)
+
+         lnR = coldata(14+2,i)
+         lnRrad = coldata(2+28,i)
+         lnRhoReh = coldata(27+2,i)
+         lnRhoEnd = coldata(2+30,i)
+         
+         if ((wreh.lt.-1._dp/3._dp) &
+              .or.(p.lt.pmin).or.(logmu.gt.logmumax) &
+              .or.(lnRhoReh.gt.lnRhoEnd) &
+              .or.(lnRhoReh.lt.lnRhoNuc)) then
+            chisq = 1e30
+            coldata(1,i) = 0.
+            coldata(2,i) = coldata(2,i) + chisq/2
+         else
+            lnRhoRehMax = max(lnRHoRehMax,lnRhoReh)                          
+            lnRhoEndMax = max(lnRHoEndMax,lnRhoEnd)
+            lnRmax = max(lnRmax, lnR)
+            lnRradmax = max(lnRradmax, lnRrad)
+         endif
+
+      enddo
+      print *,'Max values after applying prior!'
+      print *,'lnRhoEndMax= ',lnRhoEndMax
+      print *,'lnRhoRehMax= ',lnRhoRehMax      
+      print *,'lnRradMax= lnRmax= ',lnRradMax, lnRmax
+   endif
+      
+   if (adjEps1ToLog) then
+      print *,'----> assumed mapped log(eps1) in col 8'
+      print *,'----> from flat eps1 to flat log(eps1)'
+      read(*,*)
+      do i=0, nrows-1
+         
+         logeps1 = coldata(8+2,i)
+         eps1 = 10.**logeps1
+         
+         coldata(1,i) = coldata(1,i)/(eps1*log(10.))
+         coldata(2,i) = coldata(2,i) 
+
+      end do
+   endif
+
+
+
+   if (extractHinf) then
+
+      write (*,*)'Computing and adjusting  bounds for ln(10^5 Hinf/mpl) in 17!'
+      
+      do i=0,nrows-1
+         ln10p5HinfMax = max(ln10p5HinfMax,coldata(2+17,i))
+         ln10p5HinfMin = min(ln10p5HinfMin,coldata(2+17,i))         
+
+         if (coldata(2+17,i).lt.ln10p5HinfCut) then
+            chisq = 1e30
+            coldata(1,i) = 0.
+            coldata(2,i) = coldata(2,i) + chisq/2
+         endif
+
+      enddo
+               
+      write(*,*)'----> ln(10^5 Hinf) Min= Max= ',ln10p5HinfMin,ln10p5HinfMax
+
+   endif
+
+   if (mapnslogr) then
+
+      write (*,*)'Adjusting  bounds for log(r) in 14!'
+
+      do i=0,nrows-1
+         
+         if (coldata(2+14,i).lt.logrCut) then
+            chisq = 1e30
+            coldata(1,i) = 0.
+            coldata(2,i) = coldata(2,i) + chisq/2
+         endif
+
+      enddo
+   endif
+
+#ifdef SRREHEAT
+   if (doReheatFromSR) then
+      write(*,*)'Adjusting prior for reheat SR analysis'
+
+      select case (reheatModel)
+
+      case ('largef')
+        
+        
+         do i=0,nrows-1
+            
+            eps1 = 10**(coldata(2+8,i))
+            eps2 = coldata(2+9,i)
+            p = coldata(2+14,i)
+            wreh = (p-2.)/(p+2.)
+
+!large field models cannot produce eps2<=0 (out of prior)
+            if ((eps2.le.0.).or.(p.gt.pmax).or.(p.lt.pmin)) then
+               chisq = 1e30
+               coldata(1,i) = 0.
+               coldata(2,i) = coldata(2,i) + chisq/2
+            endif
+
+!reheating is after end of inflation
+            lnRhoReh = coldata(11+2,i)
+            lnRhoEnd = coldata(12+2,i)
+
+!            print *,'lnrhoreh end=',lnRhoReh,lnRhoEnd,eps1,eps2
+
+            if (wreh.lt.-1./3.) then
+               chisq = 1e30
+               coldata(1,i) = 0.
+               coldata(2,i) = coldata(2,i) + chisq/2
+            endif
+            if (lnRhoReh.gt.lnRhoEnd) then
+               chisq = 1e30
+               coldata(1,i) = 0.
+               coldata(2,i) = coldata(2,i) + chisq/2
+            endif
+            if (lnRhoReh.lt.lnRhoNuc) then
+               chisq = 1e30
+               coldata(1,i) = 0.
+               coldata(2,i) = coldata(2,i) + chisq/2            
+            endif
+         enddo
+         
+      case ('smallf')
+
+         wreh = wfix
+                 
+         do i=0,nrows-1
+            
+            eps1 = 10**(coldata(2+8,i))
+            eps2 = coldata(2+9,i)
+            eps3 = coldata(2+10,i)
+
+            p = coldata(2+14,i)            
+
+
+!small field models (out of prior)
+            if (.not.(sr_checkobs_sf(eps1,eps2,eps3)) &
+                 .or.(p.gt.pmax).or.(p.lt.pmin)) then
+               chisq = 1e30
+               coldata(1,i) = 0.
+               coldata(2,i) = coldata(2,i) + chisq/2
+            endif
+
+
+!reheating is after end of inflation
+            lnRhoReh = coldata(11+2,i)
+            lnRhoEnd = coldata(12+2,i)
+
+            if (wreh.lt.-1./3.) then
+               chisq = 1e30
+               coldata(1,i) = 0.
+               coldata(2,i) = coldata(2,i) + chisq/2
+            endif
+            if (lnRhoReh.gt.lnRhoEnd) then
+               chisq = 1e30
+               coldata(1,i) = 0.
+               coldata(2,i) = coldata(2,i) + chisq/2
+            endif
+            if (lnRhoReh.lt.lnRhoNuc) then
+               chisq = 1e30
+               coldata(1,i) = 0.
+               coldata(2,i) = coldata(2,i) + chisq/2            
+            endif
+         enddo
+         
+      end select
+
+      
+   endif
+#endif
+!end addon
+
+
+
+
   end subroutine AdjustPriors
 
   subroutine MapParameters(invars)
@@ -131,7 +377,194 @@ contains
 ! map parameters in invars: eg. invars(3)=invars(3)*invars(4) 
 
 !    invars(2+13)=invars(17+2)*exp(-invars(2+4)) 
-    stop 'Need to write MapParameters routine first'  
+
+
+!addon
+    integer, save :: counter = 0
+    real :: lnR, lnRhoReh,lnRhoEnd, ns,logr
+    real :: cseps1,eps1mdel1,eps2pdel1,del1
+    real :: as0,ln10p5Hinf
+    real, parameter :: C_const = -0.7296
+
+    integer, parameter :: dp=kind(1._8)
+    real(dp) :: p, mu, wreh, Pstar,bfold,eps1,eps2,eps3
+    real(dp), dimension(2) :: buffer2D
+    real(dp), dimension(3) :: buffer3D
+        
+
+    counter = counter + 1
+
+    if (mapRhoReh) then
+
+       if (mapLFPowerToWreh) then
+          
+          p = invars(2+10)
+          wreh = (p-2._dp)/(p+2._dp)
+          
+          if (counter.eq.1) print *,'reheat marginalisation on p in col10'
+       else
+          
+          wreh = wfix
+          if (counter.eq.1) print *,'marginalising at fixed wreh= ',wreh
+          
+       endif
+       
+       if (wreh.eq.1./3.) stop 'Treh inversion is singular!'
+       
+       if (counter.eq.1) print *,'mapping lnR -> lnRhoReh in col27'
+
+       lnR = invars(2+14)
+       
+       lnRhoEnd = invars(2+30)
+
+       lnRhoReh = 4._dp*((1._dp+wreh)/(1._dp/3._dp-wreh))*lnR &
+            - 4._dp*(1._dp/3._dp + wreh)/(1._dp/3._dp - wreh)*lnRhoEnd/2._dp
+       
+       invars(27+2)=lnRhoReh
+
+    end if
+    
+
+    if (extractHinf) then
+       del1 = 0.
+
+       if (counter.eq.1) print *,'assumiing 10**(col10) is eps1'
+
+       eps1mdel1 = 10**invars(2+8)
+       eps2pdel1 = invars(2+9)
+       eps1 = eps1mdel1
+       eps2 = eps2pdel1
+       eps3 = invars(2+10)
+
+
+       cseps1=eps1mdel1
+!       cseps1 = eps1/invars(2+15)
+!       cseps1= 10.**invars(2+15)
+
+!       eps1 = 10.**invars(2+8)
+!       eps2 = invars(2+9)
+!       eps3 = invars(2+10)
+
+
+
+       as0 = 1. + eps1mdel1*(-2.*(C_const+1.)) - eps2pdel1*C_const
+
+       if (use2ndSR) then
+          if (counter.eq.1) print *,'using 2nd order slow-roll!'
+          as0 = as0 + eps1*eps1 * (2.*C_const**2 + 2.*C_const + pi**2/2. - 5.) &
+               + eps1*eps2 * (C_const**2 - C_const + 7.*pi**2/12. - 7.) &
+               + eps2*eps2 * (0.5*C_const**2 + pi**2/8. - 1.) &
+               + eps2*eps3 * (-0.5*C_const**2 + pi**2/24.)
+       endif
+       
+       ln10p5Hinf = + 0.5*invars(2+16) &
+            + 0.5*log(3.141592653589793238*cseps1) - 0.5*log(as0)
+       
+
+       invars(2+17) = ln10p5Hinf
+    endif
+
+    if (mapnslogr) then
+       
+       if (counter.eq.1) then          
+          print *,'assumiing 10**(col10) is eps1'
+          print *,'Mapping ns --> col15 & logr --> col14'
+       endif
+
+       eps1 = 10**invars(2+8)
+       eps2 = invars(2+9)
+       eps3 = invars(2+10)
+
+       ns = 1 - 2.*eps1 - eps2
+       logr = log10(16*eps1)
+
+       if (use2ndSR) then
+          if (counter.eq.1) print *,'using 2nd order slow-roll!'
+          ns = ns - 2*eps1**2 - (2*C_const + 3)*eps1*eps2 - C_const*eps2*eps3
+          logr = logr + log10(1 + C_const*eps2)
+       endif
+
+       invars(2+14) = logr
+       invars(2+15) = ns
+      
+    endif
+       
+
+#ifdef SRREHEAT
+    if (doReheatFromSR) then
+
+       if (counter.eq.1) then
+          write(*,*)'************************************************'
+          write(*,*)'Model is: ',reheatModel
+          write(*,*)'Doing Reheating analysis from SR:'
+          write(*,*)'Mapping lnRhoReh -> col11 & lnRhoEnd -> col12'
+          write(*,*)'& -N* -> col13 & p -> col14 & mu -> col15'
+          write(*,*)'************************************************'
+       endif
+          
+       Pstar = exp(invars(2+16))*1.e-10
+             
+       eps1 = 10**(invars(2+8))
+       eps2 = invars(2+9)
+       eps3 = invars(2+10)
+       
+       lnRhoReh = 0.
+       lnRhoEnd = 0.
+       bfold=0.
+
+       select case(reheatModel)
+          
+       case ('largef')
+                                                   
+!large field models cannot produce eps2<=0 (out of prior)
+          if (eps2.le.0.) return
+
+          buffer2d = sr_matter_obs_lf(eps1,eps2)
+
+          p = buffer2d(1)
+          invars(15+2)=p
+
+          if ((p.gt.pmax).or.(p.lt.pmin)) return
+
+          lnRhoReh = sr_lnrhoreh_lf(wreh,eps1,eps2,Pstar,bfold)
+          lnRhoEnd = sr_lnrhoend_lf(p,Pstar)
+
+          invars(11+2)=lnRhoReh
+          invars(12+2)=lnRhoEnd
+          invars(13+2)=bfold
+         
+       case ('smallf')
+
+!small field models cannot produce any eps1,eps2,eps3<=0 (out of prior)
+          if (.not.sr_checkobs_sf(eps1,eps2,eps3)) return
+          
+          buffer3d = sr_matterovermu_obs_sf(eps1,eps2,eps3)
+
+          p = buffer3d(1)
+          mu = buffer3d(2)
+          wreh = wfix
+
+          invars(14+2)=p
+          invars(15+2)=mu
+
+          if ((p.gt.pmax).or.(p.lt.pmin)) return
+
+          lnRhoReh = sr_lnrhoreh_sf(wreh,eps1,eps2,eps3,Pstar,bfold)
+          lnRhoEnd = sr_lnrhoend_sf(p,mu,Pstar)
+
+          invars(11+2) = lnRhoReh
+          invars(12+2) = lnRhoEnd
+          invars(13+2) = bfold
+
+       case default
+
+          stop 'model not implemented!'
+
+       end select
+
+    endif
+#endif
+!end addon
   
   end subroutine MapParameters
 
diff --git a/source/MCMC.f90 b/source/MCMC.f90
index 9c3285d..c731cb4 100644
--- a/source/MCMC.f90
+++ b/source/MCMC.f90
@@ -215,7 +215,7 @@ contains
   
     Like = GetLogLike(grid(r))
    
-    if (Feedback > 1) write (*,*) r, 'Likelihood: ', Like, 'Current Like:', CurLike
+    if (Feedback > 1) write (*,*) r, 'Likelihood: ', Like, 'Current Like:', CurLike, 'Rank', MPIRank
 
     if ((Like /= logZero) .and. (CurLike > Like .or. randexp1() > Like - CurLike)) then
     !Accept
@@ -501,7 +501,7 @@ function WL_Weight(L) result (W)
            output_lines = output_lines +1
            call WriteParams(CurParams, real(mult), CurLike)
           end if
-          if (Feedback > 1) write (*,*) instance, 'Slicing, Current Like:', CurLike
+          if (Feedback > 1) write (*,*) instance, 'Slicing, Current Like:', CurLike, 'Rank', MPIRank
           mult = 1
           if (num_slow /=0) call SliceSampleSlowParam(CurParams, CurLike) 
           if (num_fast /=0) call SliceSampleFastParams(CurParams, CurLike) 
@@ -536,7 +536,7 @@ function WL_Weight(L) result (W)
 
     Like = GetLogLike(Trial) 
 
-    if (Feedback > 1) write (*,*) 'Likelihood: ', Like, 'Current Like:', CurLike
+    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 --git a/source/Makefile b/source/Makefile
index d3219e8..f2333e2 100644
--- a/source/Makefile
+++ b/source/Makefile
@@ -1,119 +1,84 @@
-#You may need to edit the library paths for MKL for Intel
-#Beware of using optmizations that lose accuracy - may give errors when running
+# >>> DESIGNED FOR GMAKE <<<
 
-##Uncomment the next line to include dr7 LRG
-EXTDATA =
-#EXTDATA = LRG
+# Unified Systems makefile for COSMOMC
+# Add FLAGS -DMPI for using MPI
+
+ext=$(shell uname | cut -c1-3)
+
+CC=cc
+
+ifeq ($(ext),IRI)
+F90C= f90
+FFLAGS= -Ofast -mp -n32 -LANG:recursive=ON -lmpi -DMPI
+WMAPFLAGS= $(FFLAGS)
+LAPACKL = -lcomplib.sgimath_mp
+INCLUDE = -I../camb
+CFITSIODIR =
+GSLDIR =
+endif
+
+ifeq ($(ext),Lin)
+F90C=gfortran
+FFLAGS= -O -cpp -fopenmp
+WMAPFLAGS= -O
+LAPACKL = -llapack -lblas
+INCLUDE = -I../camb 
+CFITSIODIR = 
+GSLDIR = /usr
+endif 
+
+ifeq ($(ext),OSF)
+F90C=f90
+FFLAGS= -omp -O -arch host -math_library fast -tune host -fpe1
+WMAPFLAGS= $(FFLAGS)
+LAPACKL  = -lcxml
+INCLUDE = -I../camb 
+CFITSIODIR = 
+GSLDIR =
+endif
 
-#set WMAP empty not to compile with WMAP
-WMAP = /home/aml1005/WMAP7/likelihood_v4
-
-#Only needed for WMAP
-cfitsio = /usr/local/cfitsio/intel10/64/3.040
-
-#GSL only needed for DR7 LRG
-GSLPATH = /home/aml1005/libs/gsl
-
-IFLAG = -I
-INCLUDE=
-
-#Intel MPI (assuming mpif77 set to point to ifort)
-#these settings for ifort 11.1 and higher; may need to add explicit link directory otherwise
-F90C     = mpif90
-FFLAGS = -O2 -ip -W0 -WB -openmp -fpp -DMPI -vec_report0 -mkl=parallel
-LAPACKL = -lmkl_lapack
-
-#HPCF settings. Use Inteal 9 or 10.1+, not 10.0
-#F90C     = mpif90
-#FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp -DMPI -vec_report0
-#LAPACKL = -L/usr/local/Cluster-Apps/intel/mkl/10.2.2.025/lib/em64t -lmkl_lapack -lmkl -lguide -lpthread
-#GSLPATH = /usr/local/Cluster-Apps/gsl/1.9
-#cfitsio = /usr/local/Cluster-Users/cpac/cmb/2.1.0/cfitsio
-#INCLUDE=
-
-#COSMOS: use "module load cosmolib latest"
-#use "runCosmomc" (globally installed) to run, defining required memory usage
-ifeq ($(COSMOHOST),cosmos)
-F90C = ifort
-FFLAGS = -openmp -O3 -w -fpp2 -DMPI
-LAPACKL = -mkl=sequential -lmkl_lapack -lmpi
-cfitsio = $(CFITSIO)
-WMAP = $(COSMOLIB)/WMAP7
-GSLPATH = $(GSL_ROOT)
+ifeq ($(ext),Sun)
+F90C=f90
+FFLAGS= -O4 -xarch=native64 -openmp -ftrap=%none
+WMAPFLAGS= $(FFLAGS)
+LAPACKL = -lsunperf -lfsu
+INCLUDE = -I../camb -M../camb
+CFITSIODIR = 
+GSLDIR =
 endif
 
-#Intel fortran 8, check you have the latest update from the Intel web pages
-#See Makefile_intel for ifc 7.1 or lower (some versions have problems)
-#F90C     = ifort
-#FFLAGS = -O2 -Vaxlib -ip -W0 -WB -openmp -fpp
-#LAPACKL = -L/opt/intel/mkl61/lib/32 -lmkl_lapack -lmkl_ia32 -lguide -lpthread
-
-#G95; make sure LAPACK and MPI libs also compiled with g95
-#F90C     = mpif90
-#FFLAGS = -O2 -cpp -DMPI
-#LAPACKL = /LAPACK/lapack_LINUX.a /LAPACK/blas_LINUX.a
-
-#GFortran: if pre v4.3 add -D__GFORTRAN__
-#F90C     = gfortran
-#FFLAGS = -O2 -ffree-form -x f95-cpp-input
-#LAPACKL = -Wl,-framework -Wl,accelerate
-#may need to delete -Wl,accelerate, and/or add -shared
-
-#SGI, -mp toggles multi-processor. Use -O2 if -Ofast gives problems.
-#Not various versions of the compiler are buggy giving erroneous seg faults with -mp.
-#Version 7.3 is OK, currently version 7.4 is bugged, as are some earlier versions.
-#F90C     = f90
-#LAPACKL = -lcomplib.sgimath
-#FFLAGS  = -Ofast -mp
-
-#Digital/Compaq fortran, -omp toggles multi-processor
-#F90C    = f90
-#FFLAGS  = -omp -O4 -arch host -math_library fast -tune host -fpe1
-#LAPACKL  = -lcxml
-
-#Absoft ProFortran, single processor, set -cpu:[type] for your local system
-#F90C     = f95
-#FFLAGS = -O2 -s -cpu:athlon -lU77 -w -YEXT_NAMES="LCS" -YEXT_SFX="_"
-#LAPACKL =  -llapack -lblas -lg2c
-#IFLAG = -p
-
-#NAGF95, single processor:
-#F90C     = f95
-#FFLAGS = -DNAGF95 -O3
-#LAPACKL = -llapack -lblas -lg2c
-
-#PGF90
-#F90C = pgf90
-#FFLAGS = -O2 -DESCAPEBACKSLASH
-#LAPACKL = -llapack -lblas
-
-
-#Sun, single processor:
-#F90C     = f90
-#FFLAGS = -fast -ftrap=%none
-#LAPACKL = -lsunperf -lfsu
-#LAPACKL = -lsunperf -lfsu -lsocket -lm
-#IFLAG = -M
-
-#Sun MPI
-#F90C = mpf90
-#FFLAGS =  -O4 -openmp -ftrap=%none -dalign -DMPI
-#LAPACKL = -lsunperf -lfsu -lmpi_mt
-#IFLAG = -M
-
-#Sun parallel enterprise:
-#F90C     = f95
-#FFLAGS =  -O4 -xarch=native64 -openmp -ftrap=%none
-#LAPACKL = -lsunperf -lfsu
-#IFLAG = -M
-
-
-#IBM XL Fortran, multi-processor (run "module load lapack" then run "gmake")
-# See also http://cosmocoffee.info/viewtopic.php?t=326
-#F90C     = xlf90_r $(LAPACK)
-#FFLAGS  = -WF,-DIBMXL -qsmp=omp -qsuffix=f=f90:cpp=F90 -O3 -qstrict -qarch=pwr3 -qtune=pwr3
-#INCLUDE = -lessl
-#LAPACKL =
+ifeq ($(ext),AIX)
+F90C    = mpxlf90_r
+FFLAGS  = -O4 -WF,-DIBMXL,-DMPI -qstrict -qsmp=omp -qmaxmem=-1 -qsuffix=f=f90:cpp=F90
+WMAPFLAGS= $(FFLAGS)
+LAPACKL = -lessl
+INCLUDE = -I../camb
+CFITSIODIR = 
+GSLDIR =
+endif
+
+#EXTDATA = LRG
+EXTDATA = LRG
+EXTINCLUDE = -I$(GSLDIR)/include
+EXTOBJS = bsplinepk.o
+
+
+WMAPDIR = ../WMAP
+WMAPINCLUDE = -I$(CFITSIODIR)/include
+WMAPOBJS = read_archive_map.o \
+	read_fits.o \
+	healpix_types.o \
+	br_mod_dist.o \
+	WMAP_7yr_options.o \
+	WMAP_7yr_util.o \
+	WMAP_7yr_gibbs.o \
+	WMAP_7yr_tt_pixlike.o \
+	WMAP_7yr_tt_beam_ptsrc_chisq.o \
+	WMAP_7yr_teeebb_pixlike.o \
+	WMAP_7yr_tetbeebbeb_pixlike.o \
+	WMAP_7yr_likelihood.o
+
+
 
 PROPOSE = propose.o
 CLSFILE = CMB_Cls_simple.o
@@ -121,51 +86,44 @@ CLSFILE = CMB_Cls_simple.o
 #Can use params_H if you prefer more generic parameters
 PARAMETERIZATION = params_CMB.o
 
-F90FLAGS = -DMATRIX_SINGLE $(FFLAGS) $(IFLAG)../camb $(INCLUDE)
-LINKFLAGS = -L../camb -lcamb $(LAPACKL)
+LINKFLAGS =  -L../camb -lcamb  $(LAPACKL)
+
+F90FLAGS = -DMATRIX_SINGLE $(FFLAGS) $(INCLUDE)
 
-DISTFILES = ParamNames.o Matrix_utils.o settings.o IO.o GetDist.o
+DISTFILES = utils.o ParamNames.o Matrix_utils.o settings.o IO.o GetDist.o
 
-OBJFILES=  ParamNames.o Matrix_utils.o settings.o IO.o cmbtypes.o Planck_like.o  \
-	cmbdata.o WeakLen.o bbn.o bao.o lrggettheory.o mpk.o supernovae.o HST.o SDSSLy-a-v3.o \
+
+
+OBJFILES = utils.o ParamNames.o Matrix_utils.o settings.o IO.o cmbtypes.o Planck_like.o \
+	cmbdata.o WeakLen.o bbn.o bao.o lrggettheory.o mpk.o supernovae.o HST.o SDSSLy-a-v3.o\
 	$(CLSFILE) paramdef.o $(PROPOSE) $(PARAMETERIZATION) calclike.o \
 	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
 
-
-ifeq ($(EXTDATA),LRG)
+ifeq ($(EXTDATA),LRG)	
 F90FLAGS += -DDR71RG
-OBJFILES += bsplinepk.o
-GSLINC = -I$(GSLPATH)/include
-LINKFLAGS += -L$(GSLPATH)/lib -lgsl -lgslcblas
-endif
-
-ifneq ($(WMAP),)
-F90FLAGS += $(IFLAG)$(cfitsio)/include $(IFLAG)$(WMAP)
-LINKFLAGS +=  -L$(cfitsio)/lib -L$(WMAP) -lcfitsio
-
-OBJFILES += $(WMAP)/read_archive_map.o \
-	$(WMAP)/read_fits.o \
-	$(WMAP)/healpix_types.o \
-	$(WMAP)/WMAP_7yr_options.o \
-	$(WMAP)/WMAP_7yr_util.o \
-	$(WMAP)/WMAP_7yr_tt_pixlike.o \
-	$(WMAP)/WMAP_7yr_teeebb_pixlike.o \
-	$(WMAP)/WMAP_7yr_likelihood.o \
-	$(WMAP)/WMAP_7yr_gibbs.o \
-	$(WMAP)/WMAP_7yr_tt_beam_ptsrc_chisq.o \
-	$(WMAP)/br_mod_dist.o
+EXTINCLUDE = -I$(GSLDIR)/include
+LINKFLAGS += -lgsl -lgslcblas
+OBJFILES += $(EXTOBJS)
+endif	
+
+ifneq ($(WMAPDIR),)
+F90FLAGS += $(WMAPINCLUDE)
+LINKFLAGS += -lcfitsio
+OBJFILES += $(WMAPOBJS)
 else
 F90FLAGS += -DNOWMAP
 endif
 
-default: cosmomc
 
-all : cosmomc getdist
 
-settings.o: ../camb/libcamb.a
+default: cosmomc.$(ext)
+
+all : cosmomc.$(ext) getdist.$(ext)
+
+settings.o: utils.o
 cmbtypes.o: settings.o
 Planck_like.o: cmbtypes.o
-cmbdata.o: Planck_like.o  $(WMAPOBJS)
+cmbdata.o: Planck_like.o $(WMAPOBJS)
 WeakLen.o: cmbtypes.o
 bbn.o: settings.o
 bao.o: cmbtypes.o
@@ -184,14 +142,20 @@ postprocess.o: calclike.o
 MCMC.o: calclike.o
 driver.o: MCMC.o
 
-export FFLAGS
-export F90C
-
 .f.o:
 	f77 $(F90FLAGS) -c $<
 
 %.o: %.c
-	$(CC) $(GSLINC) -c $*.c
+	$(CC) $(EXTINCLUDE) -c $*.c
+
+%.o: $(WMAPDIR)/%.f90
+	$(F90C) $(WMAPFLAGS) $(WMAPINCLUDE) -c $<
+
+%.o: $(WMAPDIR)/%.F90
+	$(F90C) $(WMAPFLAGS) $(WMAPINCLUDE) -c $<
+
+utils.o: ../camb/utils.F90
+	$(F90C) $(F90FLAGS) -c $<
 
 %.o: %.f90
 	$(F90C) $(F90FLAGS) -c $*.f90
@@ -200,19 +164,12 @@ export F90C
 	$(F90C) $(F90FLAGS) -c $*.F90
 
 
-cosmomc: camb $(OBJFILES)
-	$(F90C) -o ../cosmomc $(OBJFILES) $(LINKFLAGS) $(F90FLAGS)
+cosmomc.$(ext): $(OBJFILES) ../camb/libcamb.a 	
+	$(F90C) -o ../$@ $(OBJFILES) $(LINKFLAGS) $(F90FLAGS)
 
 
-clean: cleancosmomc
-	rm -f ../camb/*.o ../camb/*.obj ../camb/*.mod
-
-cleancosmomc:
+clean:
 	rm -f *.o *.mod *.d *.pc *.obj ../core
 
-
-getdist: camb $(DISTFILES)
-	$(F90C) -o ../getdist $(DISTFILES) $(LINKFLAGS) $(F90FLAGS)
-
-camb:
-	cd ../camb && $(MAKE) --file=Makefile_main libcamb.a
+getdist.$(ext): $(DISTFILES)
+	$(F90C) -o ../$@ $(DISTFILES) $(LINKFLAGS) $(F90FLAGS) 
diff --git a/source/driver.F90 b/source/driver.F90
index d5376fb..3fd1f44 100644
--- a/source/driver.F90
+++ b/source/driver.F90
@@ -57,11 +57,17 @@ program SolveCosmology
 
 #ifdef MPI 
 
-        if (instance /= 0) call DoAbort('With MPI should not have second parameter')
+!wig        if (instance /= 0) call DoAbort('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                                             
-        write (numstr,*) instance
+!wig        write (numstr,*) instance
+ 	instance_shift = MPIrank + instance_shift
+ 	write (numstr,*) instance_shift
+!end wig
         rand_inst = instance
         if (ierror/=MPI_SUCCESS) call DoAbort('MPI fail')
 
@@ -102,6 +108,13 @@ program SolveCosmology
 
         HighAccuracyDefault = Ini_Read_Logical('high_accuracy_default',.false.)
         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
         
         if (Ini_HasKey('highL_unlensed_cl_template')) &
           highL_unlensed_cl_template=  ReadIniFilename(DefIni,'highL_unlensed_cl_template')
diff --git a/source/params_CMB.f90 b/source/params_CMB.f90
index a453def..0634a8a 100644
--- a/source/params_CMB.f90
+++ b/source/params_CMB.f90
@@ -57,24 +57,46 @@
        integer, intent(in) :: in
 
 
-       if (Power_Name == 'power_tilt') then
-
-       P%InitPower%k_0_scalar = pivot_k
-       P%InitPower%k_0_tensor = pivot_k
+!wig
+       if (Power_Name == 'power_sr') then
+         
+          P%InitPower%P_scalar =  cl_norm*CMB%norm(norm_As)
+          P%InitPower%eps1(in) =  CMB%InitPower(1)
+          P%InitPower%eps2(in) =  CMB%InitPower(2)
+          P%InitPower%eps3(in) =  CMB%InitPower(3)
+
+          P%InitPower%logeps1invsigma0(in) = CMB%InitPower(5)
+          P%InitPower%loginvsigma0(in) = P%InitPower%logeps1invsigma0(in) &
+               - log10(P%InitPower%eps1(in))
+          P%InitPower%invsigma0(in) = 10**(P%InitPower%loginvsigma0(in))
+          P%InitPower%Xsigma0(in) = CMB%InitPower(6)
+          P%InitPower%varphi(in)  = CMB%InitPower(7)
+          !CMB%InitPower(7) unused.
+          CMB%InitPower(7) = 1./P%InitPower%invsigma0(in)
+
+       elseif (Power_Name == 'power_dbi') then
+
+          P%InitPower%P_scalar =  cl_norm*CMB%norm(norm_As)
+          P%InitPower%eps1(in) =  10.**(CMB%InitPower(4))
+          P%InitPower%del1(in) = P%InitPower%eps1(in) - CMB%InitPower(1)
+          P%InitPower%eps2(in) =  CMB%InitPower(2) - P%InitPower%del1(in)
+          P%InitPower%eps3(in) =  CMB%InitPower(3)
+
+          P%InitPower%logeps1invsigma0(in) = CMB%InitPower(5)
+          P%InitPower%loginvsigma0(in) = P%InitPower%logeps1invsigma0(in) &
+               - log10(P%InitPower%eps1(in))
+          P%InitPower%invsigma0(in) = 10**(P%InitPower%loginvsigma0(in))
+          P%InitPower%Xsigma0(in) = CMB%InitPower(6)
+          P%InitPower%varphi(in)  = CMB%InitPower(7)
+          !CMB%InitPower(7) unused.
+          CMB%InitPower(7) = 1./P%InitPower%invsigma0(in)
+
+          P%InitPower%gam(in) = 10**(-CMB%InitPower(8))*P%InitPower%eps1(in)
 
-       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 wig
 
     end subroutine SetCAMBInitPower
  
diff --git a/source/settings.f90 b/source/settings.f90
index dd29269..3beb896 100644
--- a/source/settings.f90
+++ b/source/settings.f90
@@ -6,6 +6,13 @@ module settings
   implicit none
 
   real :: AccuracyLevel = 1.
+!wig
+  real :: lAccuracyLevel = 1.
+  real :: lSampleLevel = 1.
+  real :: kSrcSampleLevel = 1.
+  real :: kIntSampleLevel = 1.
+  logical :: CambComputeEveryl = .false.
+!end wig
   !Set to >1 to use CAMB etc on higher accuracy settings. 
   !Does not affect MCMC (except making it all slower)
 
@@ -19,7 +26,10 @@ module settings
 !                    (e.g. beam uncertainty modes, etc, specific to dataset)
 
   integer, parameter :: num_hard =7 
-  integer, parameter :: num_initpower = 3 
+!wig
+!  integer, parameter :: num_initpower = 3 
+  integer, parameter :: num_initpower = 8
+!end wig
   integer, parameter :: num_freq_params = 1
   integer, parameter :: num_norm = 2 + num_freq_params 
   integer, parameter :: num_nuisance_params= 0
@@ -64,6 +74,9 @@ module settings
 
   integer :: num_threads = 0
   integer :: instance = 0
+!wig
+  integer :: instance_shift = 0
+!end wig
   integer :: MPIchains = 1, MPIrank = 0
 
   logical :: Use_LSS = .true.