diff --git a/README b/README
new file mode 100644
index 0000000..6c059c2
--- /dev/null
+++ b/README
@@ -0,0 +1,67 @@
+**********************************************************************
+Adds fingerprints features and modulations to power law primordial
+power spectra as in arXigv:1106.1635
+**********************************************************************
+
+-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 "feat" 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_tilt.f90
+Implementation of fingerprints'features.
+
+- 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/driver.F90
+2nd instance in mpirun shifts the file names.
+
+- cosmomc/params.ini
+New keys for the above material.
+
+
+
diff --git a/camb/Makefile b/camb/Makefile
index 3a42267..8c94044 100644
--- a/camb/Makefile
+++ b/camb/Makefile
@@ -1,85 +1,136 @@
-#CAMB Makefile
+# >>> 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
 
-#Set FISHER=Y to compile bispectrum fisher matrix code
-FISHER=
+ifeq ($(ext),Lin)
+F90C=gfortran
+FFLAGS= -O -fopenmp
+LFLAGS= 
+endif 
 
-#Edit for your compiler
-#Note there are many old ifort versions, some of which behave oddly
+ifeq ($(ext),OSF)
+F90C=f90
+FFLAGS= -omp -O -arch host -math_library fast -tune host -fpe1
+LFLAGS= 
+endif
 
+ifeq ($(ext),Sun)
+F90C=f90
+FFLAGS= -O4 -xarch=native64 -openmp -ftrap=%none
+LFLAGS= 
+endif
 
-#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 -fast -W0 -WB -fpp2 -vec_report0
-## This is flag is passed to the Fortran compiler allowing it to link C++ if required (not usually):
-F90CRLINK = -cxxlib
-ifneq ($(FISHER),)
-FFLAGS += -mkl
+ifeq ($(ext),AIX)
+F90C=xlf90_r
+FFLAGS= -O4 -qsmp=omp -qmaxmem=-1 -qstrict -qfree=f90 -qsuffix=f=f90:cpp=F90
+LFLAGS= 
 endif
 
-#Gfortran compiler:
-#The options here work in v4.5, delete from RHS in earlier versions (15% slower)
-#if pre v4.3 add -D__GFORTRAN__
-#With v4.6+ try -Ofast -march=native -fopenmp
-#On my machine v4.5 is about 20% slower than ifort
-#F90C     = gfortran
-#FFLAGS =  -O3 -fopenmp -ffast-math -march=native -funroll-loops
-
-
-#Old Intel ifc, add -openmp for multi-processor (some have bugs):
-#F90C     = ifc
-#FFLAGS = -O2 -Vaxlib -ip -W0 -WB -quiet -fpp2
-#some systems can can also add e.g. -tpp7 -xW
-
-#G95 compiler
-#F90C   = g95
-#FFLAGS = -O2
-
-#SGI, -mp toggles multi-processor. Use -O2 if -Ofast gives problems.
-#F90C     = f90
-#FFLAGS  = -Ofast -mp
-
-#Digital/Compaq fortran, -omp toggles multi-processor
-#F90C    = f90
-#FFLAGS  = -omp -O4 -arch host -math_library fast -tune host -fpe1
-
-#Absoft ProFortran, single processor:
-#F90C     = f95
-#FFLAGS = -O2 -cpu:athlon -s -lU77 -w -YEXT_NAMES="LCS" -YEXT_SFX="_"
-
-#NAGF95, single processor:
-#F90C     = f95
-#FFLAGS = -DNAGF95 -O3
-
-#PGF90
-#F90C = pgf90
-#FFLAGS = -O2 -DESCAPEBACKSLASH -Mpreprocess
-
-#Sun V880
-#F90C = mpf90
-#FFLAGS =  -O4 -openmp -ftrap=%none -dalign -DMPI
-
-#Sun parallel enterprise:
-#F90C     = f95
-#FFLAGS =  -O2 -xarch=native64 -openmp -ftrap=%none
-#try removing -openmp if get bus errors. -03, -04 etc are dodgy.
-
-#IBM XL Fortran, multi-processor (run gmake)
-#F90C     = xlf90_r
-#FFLAGS  = -DESCAPEBACKSLASH -DIBMXL -qsmp=omp -qsuffix=f=f90:cpp=F90 -O3 -qstrict -qarch=pwr3 -qtune=pwr3
-
-#Settings for building camb_fits
-#Location of FITSIO and name of library
-FITSDIR       ?= /home/cpac/cpac-tools/lib
-FITSLIB       = cfitsio
-#Location of HEALPIX for building camb_fits
-HEALPIXDIR    ?= /home/cpac/cpac-tools/healpix
-
-ifneq ($(FISHER),)
+
+#Files containing evolution equations initial power spectrum module
+EQUATIONS     = equations
+POWERSPECTRUM = power_tilt
+REIONIZATION = reionization
+RECOMBINATION = recfast
+#RECOMBINATION = cosmorec
+
+
+BISPECTRUM = SeparableBispectrum
+
+DENABLE_FISHER=
+
+#Module doing non-linear scaling
+NONLINEAR     ?= halofit
+
+#Driver program
+DRIVER        ?= inidriver.F90
+#DRIVER        = sigma8.f90
+#DRIVER        = tester.f90
+
+
+CAMBBIN       = camb.$(ext)
+CAMBLIB       = libcamb_$(RECOMBINATION).a
+
+ifneq ($(DENABLE_FISHER),)
 FFLAGS += -DFISHER
+LFLAGS += -llapack -lblas
 EXTCAMBFILES = Matrix_utils.o
 else
 EXTCAMBFILES =
 endif
 
-include ./Makefile_main
+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
+
+
+ifeq ($(RECOMBINATION),cosmorec)
+COSMOREC_PATH =../CosmoRec/
+LFLAGS += -L$(COSMOREC_PATH) -lCosmoRec -lgsl -lgslcblas -lstdc++
+endif
+
+
+ifeq ($(RECOMBINATION),hyrec)
+HYREC_PATH = ../HyRec/
+LFLAGS += -L$(HYREC_PATH) -lhyrec -lstdc++
+endif
+
+default: $(CAMBBIN)
+all: $(CAMBBIN) $(CAMBLIB)
+
+
+subroutines.o: constants.o utils.o
+$(POWERSPECTRUM).o: 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).o: 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
+
+
+$(CAMBBIN): $(CAMBOBJ) $(DRIVER)
+	$(F90C) $(F90FLAGS) $(CAMBOBJ) $(DRIVER) $(LFLAGS) -o $@
+
+$(CAMBLIB): $(CAMBOBJ)
+	ar -r $@ $(CAMBOBJ)
+
+camb_fits: writefits.f90 $(CAMBOBJ) $(DRIVER)
+	$(F90C) $(F90FLAGS) -I$(HEALPIXDIR)/include $(CAMBOBJ) writefits.f90 $(DRIVER) $(HEALPIXLD) -DWRITE_FITS -o $@
+
+
+%.o: %.f90
+	$(F90C) $(F90FLAGS) -c $*.f90
+
+%.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 *.$(ext)
diff --git a/camb/cmbmain.f90 b/camb/cmbmain.f90
index 71468ed..7521ea0 100644
--- a/camb/cmbmain.f90
+++ b/camb/cmbmain.f90
@@ -663,17 +663,24 @@ contains
 
          SourceAccuracyBoost = AccuracyBoost  
          if (CP%WantScalars .and. CP%Reion%Reionization .and. CP%AccuratePolarization) then
-            dlnk0=2._dl/10/SourceAccuracyBoost 
+!feat
+!            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 feat
             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 
+!feat
+!         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 feat
          if (HighAccuracyDefault) dkn2=dkn2/1.2
          if (CP%WantTensors .or. CP%WantVectors) then
               dkn1=dkn1  *0.8_dl
@@ -1096,8 +1103,11 @@ contains
       
 
        qmax_int = min(qmax,max_bessels_etak/CP%tau0)
+!feat
+!       IntSampleBoost=AccuracyBoost  
+       IntSampleBoost=AccuracyBoost*kIntSampleBoost
+!end feat
 
-       IntSampleBoost=AccuracyBoost 
        if (do_bispectrum) then
         IntSampleBoost = IntSampleBoost * 2  
         if (hard_bispectrum) IntSampleBoost = IntSampleBoost * 2  
diff --git a/camb/inidriver.F90 b/camb/inidriver.F90
index b4f893a..5efa6c2 100644
--- a/camb/inidriver.F90
+++ b/camb/inidriver.F90
@@ -281,6 +281,19 @@
           write(*,*) 'Output _params.ini not created as would overwrite input'    
          end if
        end if
+!feat
+           
+       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 feat
 
        call Ini_Close
 
@@ -295,11 +308,28 @@
          write (*,'("Age of universe/GYr  = ",f7.3)') Age  
        end if 
 
+!feat
+       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 (*,'("Feature pivot        = ",f12.8)') P%InitPower%krhalf(i)
+          write (*,'("Feature size         = ",f12.8)') P%InitPower%size(i)
+          write (*,'("Feature frequency    = ",f14.8)') P%InitPower%omega(i)
+          write (*,'("Feature phase        = ",f12.8)') P%InitPower%phase(i)
+          write (*,'("Feature power        = ",f12.8)') P%InitPower%p(i)
+       enddo
+!end feat
+
        if (global_error_flag==0) call CAMB_GetResults(P)
        if (global_error_flag/=0) then
         write (*,*) 'Error result '//trim(global_error_message)
         stop
        endif
+
     
         if (P%WantTransfer .and. .not. (P%NonLinear==NonLinear_lens .and. P%DoLensing)) then
          call Transfer_SaveToFiles(MT,TransferFileNames)
diff --git a/camb/modules.f90 b/camb/modules.f90
index 6b0a02f..7330607 100644
--- a/camb/modules.f90
+++ b/camb/modules.f90
@@ -38,8 +38,10 @@
         character(LEN=*), parameter :: version = 'Oct_12'
 
         integer :: FeedbackLevel = 0 !if >0 print out useful information about the model
-
-        logical, parameter :: DebugMsgs=.false. !Set to true to view progress and timing
+!feat
+!        logical, parameter :: DebugMsgs=.false. !Set to true to view progress and timing
+        logical, parameter :: DebugMsgs=.false.
+!end feat
 
         logical, parameter :: DebugEvolution = .false. !Set to true to do all the evolution for all k
 
@@ -189,14 +191,18 @@
       real(dl) :: lSampleBoost=1._dl
           !Increase lSampleBoost to increase sampling in lSamp%l for Cl interpolation
 
-      real(dl) :: AccuracyBoost =1._dl
-
+!feat          
+      real(dl) :: AccuracyBoost =8._dl  
+!end feat
           !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.
+!feat
+      real(sp) :: lAccuracyBoost=2._dl
+!end feat
+
           !Boost number of multipoles integrated in Boltzman heirarchy
 
       integer, parameter :: lmin = 2
@@ -212,6 +218,14 @@
       integer, parameter:: l0max=4000
 
 !     lmax is max possible number of l's evaluated
+
+!feat
+      real(dl) :: kSrcSampleBoost=1._dl
+      real(dl) :: kIntSampleBoost=1._dl
+      logical :: DoClsInterpol=.true.
+      logical :: ComputeEveryl=.true.
+!end feat 
+ 
       integer, parameter :: lmax_arr = l0max
 
       character(LEN=1024) :: highL_unlensed_cl_template = 'HighLExtrapTemplate_lenspotentialCls.dat'
@@ -592,9 +606,28 @@
         integer, intent(IN) :: max_l
         integer lind, lvar, step,top,bot,ls(lmax_arr)
         real(dl) AScale
+    
+!feat        
+        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
 
-        Ascale=scale/lSampleBoost
+           lSet%l0=lind
+           lSet%l = ls(1:lind)
+
+        else
+!end feat
 
+        Ascale=scale/lSampleBoost       
+       
         if (lSampleBoost >=50) then
          !just do all of them
          lind=0
@@ -733,6 +766,10 @@
         lSet%l0=lind
         lSet%l(1:lind) = ls(1:lind)
 
+!feat        
+        endif
+!end feat
+        
       end subroutine initlval
 
       subroutine InterpolateClArr(lSet,iCl, all_Cl, max_ind)
@@ -750,7 +787,10 @@
       if (max_ind > lSet%l0) stop 'Wrong max_ind in InterpolateClArr'
 
       xl = real(lSet%l(1:lSet%l0),dl)
-      call spline(xl,iCL(1),max_ind,cllo,clhi,ddCl(1))
+
+      call spline(xl,iCl(1),max_ind,cllo,clhi,ddCl(1))
+!feat
+      if (DoClsInterpol) then
 
             llo=1
             do il=lmin,lSet%l(max_ind)
@@ -768,6 +808,21 @@
 
             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 feat 
+
       end subroutine InterpolateClArr
 
        subroutine InterpolateClArrTemplated(lSet,iCl, all_Cl, max_ind, template_index)
@@ -2058,9 +2113,16 @@
 
 
              minkh = 1e-4
-             dlnkh = 0.02
+!feat
+!             dlnkh = 0.02
+             dlnkh = 1e-3
+!end feat
              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)
+!feat
+             write(*,*)'Transfert_SaveToFiles: dlnkh= points= ',dlnkh,points
+             write(*,*)'You should check with transfer_k_per_logint > 0'
+!end feat
              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 d33bd36..21655f8 100644
--- a/camb/params.ini
+++ b/camb/params.ini
@@ -1,13 +1,13 @@
 #Parameters for CAMB
 
 #output_root is prefixed to output file names
-output_root = test
+output_root = feature
 
 #What to do
 get_scalar_cls = T
 get_vector_cls = F
 get_tensor_cls = F
-get_transfer   = F
+get_transfer   = T
 
 #if do_lensing then scalar_output_file contains additional columns of l^4 C_l^{pp} and l^3 C_l^{pT}
 #where p is the projected potential. Output lensed CMB Cls (without tensors) are in lensed_output_file below.
@@ -84,6 +84,12 @@ tensor_spectral_index(1)  = 0
 initial_ratio(1)          = 1
 #note vector modes use the scalar settings above
 
+#Primordial features
+feature_kscale(1)         = 0.05
+feature_size(1)           = 1.
+feature_power(1)          = 0.5
+feature_pulsation(1)      = 50
+feature_phase(1)          = 0
 
 #Reionization, ignored unless reionization = T, re_redshift measures where x_e=0.5
 reionization         = T
@@ -230,8 +236,8 @@ number_of_threads       = 0
 high_accuracy_default=T
 
 #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
+#Decrease to speed up at cost of worse accuracy. Suggest 3 to 8 for features.
+accuracy_boost          = 2
 
 #Larger to keep more terms in the hierarchy evolution. 
 l_accuracy_boost        = 1
@@ -242,3 +248,14 @@ l_accuracy_boost        = 1
 #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_tilt.f90 b/camb/power_tilt.f90
index c23155f..d1ace2b 100644
--- a/camb/power_tilt.f90
+++ b/camb/power_tilt.f90
@@ -35,7 +35,9 @@
     
       character(LEN=*), parameter :: Power_Name = 'power_tilt'
 
-      integer, parameter :: nnmax= 5 
+!features
+      integer, parameter :: nnmax= 1 
+!end features
        !Maximum possible number of different power spectra to use
 
       Type InitialPowerParams
@@ -50,7 +52,16 @@
        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)
- 
+!features
+       real(dl) :: krhalf(nnmax) !=kr/2
+       real(dl) :: size(nnmax) !=beta x model dep factor
+       real(dl) :: phase(nnmax) != phi + p^2/(p-1) 2 M/H
+       real(dl) :: omega(nnmax) != 2 M/H
+       real(dl) :: p(nnmax)
+
+
+!end features
+
       end Type InitialPowerParams
 
       real(dl) curv  !Curvature contant, set in InitializePowers     
@@ -76,6 +87,13 @@
          AP%k_0_scalar = 0.05
          AP%k_0_tensor = 0.05
          AP%ScalarPowerAmp = 1
+!feat
+         AP%krhalf = 0.07
+         AP%size = 0.
+         AP%phase = 0.
+         AP%omega = 10
+         AP%p = 10
+!end feat
 
        end subroutine SetDefPowerParams
 
@@ -116,14 +134,42 @@
        !< |\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
+        real(dl) ScalarPower,k, lnrat, lnf
+        real(dl) :: angle, amp
+       
         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)   
-     
-!         ScalarPower = ScalarPower * (1 + 0.1*cos( lnrat*30 ) )
- 
+!features
+          lnf = log(k/P%krhalf(in))
+          amp = 0._dl
+
+          if ((P%p(in).eq.1._dl).or.(abs(lnf/P%p(in)).gt.log(huge(1._dl)))) then
+
+             angle = 0._dl
+
+          else
+
+              angle = P%phase(in) + P%omega(in) * (P%p(in)/(P%p(in)-1._dl)) &
+                  *( P%p(in)*exp(lnf/P%p(in))- P%p(in) )
+
+              if ((P%p(in).gt.0._dl).and.(P%p(in).lt.1._dl)) then
+                 if (k.le.P%krhalf(in)) amp = exp(lnf*(-3._dl + 2.5_dl/P%p(in))) * P%size(in)
+              else
+                 if (k.ge.P%krhalf(in)) amp = exp(lnf*(-3._dl + 2.5_dl/P%p(in))) * P%size(in)
+              endif
+
+           endif
+
+          ScalarPower=P%ScalarPowerAmp(in)*exp((P%an(in)-1)*lnrat + P%n_run(in)/2*lnrat**2) &  
+               * (1._dl + amp*sin(angle))
+
+
+          if(ScalarPower.lt.0.) then
+             ScalarPower=0._dl
+          end if
+!end features                                   
+
       end function ScalarPower
 
       
@@ -168,6 +214,23 @@
          num=num+1
          Keys(num) = 's_pivot'
          Vals(num) = P%k_0_scalar
+!features
+         num=num+1
+         Keys(num) = 'f_pivot'
+         Vals(num) = P%krhalf(in)
+         num=num+1
+         Keys(num) = 'size'
+         Vals(num) = P%size(in)
+         num=num+1
+         Keys(num) = 'omega'
+         Vals(num) = P%omega(in)
+         num=num+1
+         Keys(num) = 'power'
+         Vals(num) = P%p(in)
+         num=num+1
+         Keys(num) = 'phase'
+         Vals(num) = P%phase(in)
+!end features
          end if
          if (Tens) then
          num=num+1
@@ -202,7 +265,13 @@
 
               InitPower%an(i) = Ini_Read_Double_Array_File(Ini,'scalar_spectral_index', i)
               InitPower%n_run(i) = Ini_Read_Double_Array_File(Ini,'scalar_nrun',i,0._dl) 
-              
+!features              
+              InitPower%krhalf(i) = Ini_Read_Double_Array_File(Ini,'feature_kscale',i,0.02_dl) 
+              InitPower%size(i) = Ini_Read_Double_Array_File(Ini,'feature_size',i, 1._dl)
+              InitPower%omega(i) = Ini_Read_Double_Array_File(Ini,'feature_pulsation', i)
+              InitPower%phase(i) = Ini_Read_Double_Array_File(Ini,'feature_phase',i,0._dl) 
+              InitPower%p(i) = Ini_Read_Double_Array_File(Ini,'feature_power', i)
+!end features
               if (WantTensors) then
                  InitPower%ant(i) = Ini_Read_Double_Array_File(Ini,'tensor_spectral_index',i)
                  InitPower%rat(i) = Ini_Read_Double_Array_File(Ini,'initial_ratio',i)
diff --git a/params.ini b/params.ini
index 6baac2d..092b656 100644
--- a/params.ini
+++ b/params.ini
@@ -1,7 +1,7 @@
 #Sample parameters for cosmomc in default parameterization
 
 #Root name for files produced
-file_root = chains/test
+file_root = chains/feat
 
 #action = 0:  MCMC, action=1: postprocess .data file, action=2: find best fit point only
 action = 0
@@ -72,11 +72,11 @@ 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
@@ -141,7 +141,7 @@ 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 = F
+checkpoint = T
 
 #whether to stop on CAMB error, or continue ignoring point
 stop_on_error=  T
@@ -167,7 +167,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 = 2
+l_accuracy_level = 1
+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
@@ -197,8 +202,14 @@ param[ns] = 0.95 0.5 1.5 0.02 0.01
 param[nt] = 0 0 0 0 0
 param[nrun] = 0 0 0 0 0
 
-#log[10^10 A_s]
-param[logA] = 3 2.7 4 0.01 0.01
-param[r] = 0 0 0 0 0
+param[logkf] = -1.3 -1.3 -1.3 0 0
+param[size] = 0 0 0 0 0
+param[pow] = 0 0 0 0 0
+param[freq] = 0 0 0 0 0
+param[phase] = 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..895d4fa 100644
--- a/params_CMB.paramnames
+++ b/params_CMB.paramnames
@@ -8,7 +8,12 @@ w             w			#constant equation of state parameter for scalar field dark en
 ns            n_s		#beware that pivot scale can change in .ini file
 nt            n_t
 nrun          n_{run}
-logA          log[10^{10} A_s]
+logkf	      log(k_r/2)	#feature scale
+size	      A_{\omega}	#feature size
+pow	      p	      		#power dependency
+freq	      \omega		#oscillation frequency
+phase	      \psi		#phase shift
+lnA           ln[10^{10} A_s]
 r             r			#ratio of tensor to scalar primordial amplitudes at pivot scale
 asz           A_{SZ}		#SZ template amplitude, as in WMAP
 omegal*       \Omega_\Lambda
@@ -17,4 +22,4 @@ 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
+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 db9d905..e3659ed 100644
--- a/source/CMB_Cls_simple.f90
+++ b/source/CMB_Cls_simple.f90
@@ -404,6 +404,28 @@ contains
    integer zix
    real redshifts(matter_power_lnzsteps)
 
+!feat
+   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 feat
+
         Threadnum =num_threads
         w_lam = -1
         call CAMB_SetDefParams(P)
@@ -436,9 +458,11 @@ contains
             P%Transfer%high_precision=.true.
             P%Transfer%kmax=P%Transfer%kmax + 0.2
            end if
-          AccuracyBoost = AccuracyLevel
-          lAccuracyBoost = AccuracyLevel
-          lSampleBoost = AccuracyLevel
+!feat
+!          AccuracyBoost = AccuracyLevel
+!          lAccuracyBoost = AccuracyLevel
+!          lSampleBoost = AccuracyLevel
+!end feat
           P%AccurateReionization = .true.
         end if
                 
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 1e571c0..93f4837 100644
--- a/source/Makefile
+++ b/source/Makefile
@@ -1,122 +1,72 @@
-#use "make RECOMBINATION=cosmorec" to build with CosmoRec rather than RECFAST default
-
-#You may need to edit the library paths for MKL for Intel
-
-##Uncomment the next line to include dr7 LRG
-EXTDATA =
-#EXTDATA = LRG
-
-RECOMBINATION ?=recfast
-
-#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
-#Can add -xHost if your cluster is uniform, or specify specific processor optimizations -x...
-#If getdist gives segfaults remove openmp when compiling getdist
-F90C     = mpif90
-FFLAGS = -O3 -W0 -WB -openmp -fpp -DMPI -vec_report0 -mkl=parallel
-LAPACKL = -lmpi
-
-#GFortran: defaults for v4.5; if pre v4.3 add -D__GFORTRAN__
-#F90C     = gfortran
-#in earlier versions use FFLAGS = -O2 -ffree-form -x f95-cpp-input -D__GFORTRAN__
-#FFLAGS =  -O3 -fopenmp -ffree-form -x f95-cpp-input -ffast-math -march=native -funroll-loops
-#LAPACKL = -Wl,-framework -Wl,accelerate
-#commented above is (I think) for Mac; this is standard linux (sudo apt-get install liblapack-dev)
-#LAPACKL = -lblas -llapack
-
-#HPCF settings. Use Intel 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 -fast -w -fpp2 -DMPI
-LAPACKL = -mkl=sequential -lmkl_lapack -lmpi
-cfitsio = $(CFITSIO)
-WMAP = $(COSMOLIB)/WMAP7
-GSLPATH = $(GSL_ROOT)
+# >>> DESIGNED FOR GMAKE <<<
+
+ext=$(shell uname | cut -c1-3)
+
+CC=cc
+
+ifeq ($(ext),IRI)
+F90C= f90
+FFLAGS= -Ofast -mp -n32 -LANG:recursive=ON -lmpi -DMPI
+LAPACKL = -lcomplib.sgimath_mp
+INCLUDE = -I../camb
+LFLAGS=
+endif
+
+ifeq ($(ext),Lin)
+F90C=gfortran
+FFLAGS= -O -cpp -fopenmp
+LAPACKL = -llapack -lblas
+INCLUDE = -I../camb
+LFLAGS=
+endif 
+
+ifeq ($(ext),OSF)
+F90C=f90
+FFLAGS= -omp -O -arch host -math_library fast -tune host -fpe1
+LAPACKL  = -lcxml
+INCLUDE = -I../camb 
+LFLAGS=
+endif
+
+
+ifeq ($(ext),Sun)
+F90C=f90
+FFLAGS= -O4 -xarch=native64 -openmp -ftrap=%none
+LAPACKL = -lsunperf -lfsu
+INCLUDE = -I../camb -M../camb
+LFLAGS=
+endif
+
+ifeq ($(ext),AIX)
+F90C    = mpxlf90_r
+FFLAGS  = -O4 -WF,-DIBMXL,-DMPI -qstrict -qsmp=omp -qmaxmem=-1 -qsuffix=f=f90:cpp=F90
+LAPACKL = -lessl
+INCLUDE = -I../camb
+LFLAGS=
 endif
 
-#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
-
-#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 -Mpreprocess
-#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 =
+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
+
+
+RECOMBINATION =recfast
 
 PROPOSE = propose.o
 CLSFILE = CMB_Cls_simple.o
@@ -124,10 +74,11 @@ 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_$(RECOMBINATION) $(LAPACKL) $(F90CRLINK) 
+F90FLAGS = -DMATRIX_SINGLE $(FFLAGS) -I../camb $(INCLUDE)
+LFLAGS = -L../camb -lcamb_$(RECOMBINATION) $(LAPACKL)
+
+DISTFILES = utils.o ParamNames.o Matrix_utils.o settings.o IO.o GetDist.o
 
-DISTFILES = ParamNames.o Matrix_utils.o settings.o IO.o GetDist.o
 
 LIKEFILES= calclike.o
 
@@ -137,57 +88,45 @@ OBJFILES = ParamNames.o Matrix_utils.o settings.o IO.o cmbtypes.o Planck_like.o
 	EstCovmat.o PowellConstrainedMinimize.o minimize.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
+LFLAGS += -lgsl -lgslcblas
+OBJFILES += $(EXTOBJS)
+endif	
 
-F90CRLINK =
 
 ifeq ($(RECOMBINATION),cosmorec)
-## This is flag is passed to the Fortran compiler allowing it to link C++ (uncomment the right one).
-# GCC (gfortran/g++)
-COSMOREC_PATH ?= ../CosmoRec/
-F90CRLINK = -lstdc++ -L$(COSMOREC_PATH) -lCosmoRec -L$(GSLPATH)/lib -lgsl -lgslcblas
-# Intel Compilers (ifort/icpc)
-#F90CRLINK = -cxxlib -L$(COSMOREC_PATH) -lCosmoRec -L$(GSLPATH)/lib -lgsl -lgslcblas
+COSMOREC_PATH = ../CosmoRec/
 FFLAGS +=  -DCOSMOREC
+LFLAGS = -lstdc++ -L$(COSMOREC_PATH) -lCosmoRec -lgsl -lgslcblas
 endif
 
+
 ifeq ($(RECOMBINATION),hyrec)
-HYREC_PATH ?= ../HyRec/
-F90CRLINK += -L$(HYREC_PATH) -lhyrec
+HYREC_PATH = ../HyRec/
+LFLAGS += -L$(HYREC_PATH) -lhyrec
 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
+
+ifneq ($(WMAPDIR),)
+F90FLAGS += $(WMAPINCLUDE)
+OBJFILES += $(WMAPOBJS)
+LFLAGS += -lgsl -lgslcblas -lcfitsio
 else
 F90FLAGS += -DNOWMAP
 endif
 
-default: cosmomc
 
-all : cosmomc getdist
+
+
+default: cosmomc.$(ext)
+
+all : cosmomc.$(ext) getdist.$(ext)
 
 settings.o: ../camb/libcamb_$(RECOMBINATION).a
 cmbtypes.o: settings.o
 Planck_like.o: cmbtypes.o
-cmbdata.o: Planck_like.o  $(WMAPOBJS)
+cmbdata.o: Planck_like.o
 WeakLen.o: cmbtypes.o
 bbn.o: settings.o
 mpk.o: cmbtypes.o lrggettheory.o
@@ -205,14 +144,20 @@ MCMC.o: calclike.o
 driver.o: EstCovmat.o MCMC.o minimize.o
 minimize.o: PowellConstrainedMinimize.o calclike.o
 
-export FFLAGS
-export F90C
-
 .f.o:
 	f77 $(F90FLAGS) -c $<
 
 %.o: %.c
-	$(CC) $(GSLINC) -c $*.c
+	$(CC) $(EXTINCLUDE) -c $*.c
+
+%.o: $(WMAPDIR)/%.f90
+	$(F90C) $(FFLAGS) $(WMAPINCLUDE) -c $<
+
+%.o: $(WMAPDIR)/%.F90
+	$(F90C) $(FFLAGS) $(WMAPINCLUDE) -c $<
+
+utils.o: ../camb/utils.F90
+	$(F90C) $(F90FLAGS) -c $<
 
 %.o: %.f90
 	$(F90C) $(F90FLAGS) -c $*.f90
@@ -221,19 +166,12 @@ export F90C
 	$(F90C) $(F90FLAGS) -c $*.F90
 
 
-cosmomc: camb $(OBJFILES)
-	$(F90C) -o ../cosmomc $(OBJFILES) $(LINKFLAGS) $(F90FLAGS)
-
-
-clean: cleancosmomc
-	rm -f ../camb/*.o ../camb/*.obj ../camb/*.mod
-
-cleancosmomc:
-	rm -f *.o *.mod *.d *.pc *.obj ../core
+cosmomc.$(ext): $(OBJFILES) ../camb/libcamb_$(RECOMBINATION).a 	
+	$(F90C) -o ../$@ $(OBJFILES) $(F90FLAGS) $(LFLAGS) 
 
 
-getdist: camb $(DISTFILES)
-	$(F90C) -o ../getdist $(DISTFILES) $(LINKFLAGS) $(F90FLAGS)
+getdist.$(ext): $(DISTFILES)
+	$(F90C) -o ../$@ $(DISTFILES) $(F90FLAGS) $(LFLAGS) 
 
-camb:
-	cd ../camb && $(MAKE) --file=Makefile_main libcamb_$(RECOMBINATION).a RECOMBINATION=$(RECOMBINATION)
+clean:
+	rm -f *.o *.mod *.d *.pc *.obj ../*.$(ext) ../core
\ No newline at end of file
diff --git a/source/driver.F90 b/source/driver.F90
index 2e39ef6..486666e 100644
--- a/source/driver.F90
+++ b/source/driver.F90
@@ -61,11 +61,17 @@ program SolveCosmology
 
 #ifdef MPI 
 
-        if (instance /= 0) call DoAbort('With MPI should not have second parameter')
+!feat        if (instance /= 0) call DoAbort('With MPI should not have second parameter')
+        instance_shift=instance
+        write(*,*)'MPI file names shifted by: ',instance_shift
+!end feat
         call mpi_comm_rank(mpi_comm_world,MPIrank,ierror)
-
+        
         instance = MPIrank +1 !start at 1 for chains                                             
-        write (numstr,*) instance
+!feat        write (numstr,*) instance
+ 	instance_shift = MPIrank + instance_shift
+ 	write (numstr,*) instance_shift
+!end feat
         rand_inst = instance
         if (ierror/=MPI_SUCCESS) call DoAbort('MPI fail')
 
@@ -106,6 +112,13 @@ program SolveCosmology
 
         HighAccuracyDefault = Ini_Read_Logical('high_accuracy_default',.false.)
         AccuracyLevel = Ini_Read_Real('accuracy_level',1.)
+!feat
+        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 feat
         
         if (Ini_HasKey('highL_unlensed_cl_template')) then
           highL_unlensed_cl_template=  ReadIniFilename(DefIni,'highL_unlensed_cl_template') 
diff --git a/source/params_CMB.f90 b/source/params_CMB.f90
index a453def..d90d3bd 100644
--- a/source/params_CMB.f90
+++ b/source/params_CMB.f90
@@ -68,6 +68,13 @@
        P%InitPower%an(in) = CMB%InitPower(1)
        P%InitPower%ant(in) = CMB%InitPower(2)
        P%InitPower%n_run(in) = CMB%InitPower(3)
+!feat
+       P%InitPower%krhalf(in) = 10._dl**CMB%InitPower(4)
+       P%InitPower%size(in) = CMB%InitPower(5)
+       P%InitPower%p(in) = CMB%InitPower(6)
+       P%InitPower%omega(in) = CMB%InitPower(7)
+       P%InitPower%phase(in) = CMB%InitPower(8)       
+!end feat       
        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)
diff --git a/source/settings.f90 b/source/settings.f90
index 90d84d0..deab0c4 100644
--- a/source/settings.f90
+++ b/source/settings.f90
@@ -6,6 +6,13 @@ module settings
   implicit none
 
   real :: AccuracyLevel = 1.
+!feat
+  real :: lAccuracyLevel = 1.
+  real :: lSampleLevel = 1.
+  real :: kSrcSampleLevel = 1.
+  real :: kIntSampleLevel = 1.
+  logical :: CambComputeEveryl = .false.
+!end feat
   !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 
+!feat
+!  integer, parameter :: num_initpower = 3 
+  integer, parameter :: num_initpower = 8
+!end feat
   integer, parameter :: num_freq_params = 1
   integer, parameter :: num_norm = 2 + num_freq_params 
   integer, parameter :: num_nuisance_params= 0
@@ -66,6 +76,9 @@ module settings
 
   integer :: num_threads = 0
   integer :: instance = 0
+!feat
+  integer :: instance_shift = 0
+!end feat
   integer :: MPIchains = 1, MPIrank = 0
 
   logical :: Use_LSS = .true.