diff -r -N -c -b -B cosmomc/camb/inidriver.F90 cosmomc_addon/camb/inidriver.F90
*** cosmomc/camb/inidriver.F90	2005-03-31 00:37:14.000000000 +0100
--- cosmomc_addon/camb/inidriver.F90	2006-04-04 10:38:28.000000000 +0100
***************
*** 171,193 ****
               stop 'Unknown recombination'
             end if
  
             P%InitPower%nn = Ini_Read_Int('initial_power_num')
             if (P%InitPower%nn>nnmax) stop 'Too many initial power spectra - increase nnmax in InitialPower'
             P%InitPower%rat(:) = 1
             do i=1, P%InitPower%nn
                write (numstr,*) i 
                numstr=adjustl(numstr)
                P%InitPower%an(i) = &
                     Ini_Read_Double('scalar_spectral_index('//trim(numstr)//')')
  
                P%InitPower%n_run(i) = &
                     Ini_Read_Double('scalar_nrun('//trim(numstr)//')',0._dl)
      
                if (P%WantTensors) then
                   P%InitPower%ant(i) = Ini_Read_Double('tensor_spectral_index('//trim(numstr)//')')
                   P%InitPower%rat(i) = Ini_Read_Double('initial_ratio('//trim(numstr)//')')
                end if              
! 
                P%InitPower%ScalarPowerAmp(i) = Ini_Read_Double('scalar_amp('//trim(numstr)//')',1.d0) 
                !Always need this as may want to set tensor amplitude even if scalars not computed
             end do
--- 171,214 ----
               stop 'Unknown recombination'
             end if
  
+ !addon
+             P%InitPower%WantExternalScalar = Ini_Read_Logical('use_external_scalar',.false.)
+             P%InitPower%WantExternalTensor = Ini_Read_Logical('use_external_tensor',.false.)
+ !end addon
+ 
             P%InitPower%nn = Ini_Read_Int('initial_power_num')
             if (P%InitPower%nn>nnmax) stop 'Too many initial power spectra - increase nnmax in InitialPower'
             P%InitPower%rat(:) = 1
             do i=1, P%InitPower%nn
                write (numstr,*) i 
                numstr=adjustl(numstr)
+ !addon
+               if (P%WantScalars .or. P%WantTransfer) then
+                  if (P%InitPower%WantExternalScalar) then
+                     P%InitPower%ScalarPowerFileName(i) & 
+                          = Ini_Read_String('scalar_input_file('//trim(numstr)//')')
+                  else
+ 
                P%InitPower%an(i) = &
                     Ini_Read_Double('scalar_spectral_index('//trim(numstr)//')')
  
                P%InitPower%n_run(i) = &
                     Ini_Read_Double('scalar_nrun('//trim(numstr)//')',0._dl)
      
+                  endif
+               endif
+ 
                if (P%WantTensors) then
+ 
+                  if (P%InitPower%WantExternalTensor) then
+                     P%InitPower%TensorPowerFileName(i) & 
+                          = Ini_Read_String('tensor_input_file('//trim(numstr)//')')
+                  else
                   P%InitPower%ant(i) = Ini_Read_Double('tensor_spectral_index('//trim(numstr)//')')
                   P%InitPower%rat(i) = Ini_Read_Double('initial_ratio('//trim(numstr)//')')
+                  endif
                end if              
! !end addon
                P%InitPower%ScalarPowerAmp(i) = Ini_Read_Double('scalar_amp('//trim(numstr)//')',1.d0) 
                !Always need this as may want to set tensor amplitude even if scalars not computed
             end do
***************
*** 247,252 ****
--- 268,282 ----
         AccuracyBoost  = Ini_Read_Double('accuracy_boost',1.d0)
         lAccuracyBoost = Ini_Read_Double('l_accuracy_boost',1.d0)
         lSampleBoost   = Ini_Read_Double('l_sample_boost',1.d0)
+ !addon     
+        DoClsInterpol  = Ini_Read_Logical('do_cls_interpol',.false.)
+        ComputeEveryl  = Ini_Read_Logical('compute_every_l',.false.)
+   
+        if (FeedbackLevel > 0) then
+           write (*,*)'do_cls_interpol is',DoClsInterpol
+           write (*,*)'compute_every_l is',ComputeEveryl          
+        endif
+ !end addon
  
         if (outroot /= '') then
           call Ini_SaveReadValues(trim(outroot) //'params.ini',1)
***************
*** 290,296 ****
  #endif
          end if
  
!              
    
          end program driver
  
--- 320,329 ----
  #endif
          end if
  
! !addon
!         if (P%InitPower%WantExternalScalar.or.P%InitPower%WantExternalTensor) &
!              call FreeExternalData(P%InitPower)
! !end addon
    
          end program driver
  
diff -r -N -c -b -B cosmomc/camb/Makefile cosmomc_addon/camb/Makefile
*** cosmomc/camb/Makefile	2006-03-14 11:07:24.000000000 +0000
--- cosmomc_addon/camb/Makefile	2006-04-04 10:38:38.000000000 +0100
***************
*** 1,57 ****
! #CAMB Makefile
! #Edit for your compiler
! #Note there are many ifc versions, some of which behave oddly
! 
! #Intel 8 , -openmp toggles mutli-processor:
! F90C     = ifort
! FFLAGS = -openmp -ip -O2 -Vaxlib -W0 -WB -fpp2
! 
! # Intel 9 on IA-64 (eg. COSMOS)
! # (do "module load icomp90" before compiling)
! #F90C = ifort
! #FFLAGS = -openmp -fpp2 -w -O3 -ip -mP2OPT_hlo_prefetch=F
! 
! #Intel ifc, add -openmp for multi-processor (some have bugs):
! #F90C     = ifc
! #FFLAGS = -O2 -Vaxlib -ip -W0 -WB -quiet -fpp2
! #some systems can can also add e.g. -tpp7 -xW
! 
! #SGI, -mp toggles multi-processor. Use -O2 if -Ofast gives problems.
! #F90C     = f90
! #FFLAGS  = -Ofast -mp
! 
! #Digital/Compaq fortran, -omp toggles multi-processor
! #F90C    = f90
! #FFLAGS  = -omp -O4 -arch host -math_library fast -tune host -fpe1
! 
! #Absoft ProFortran, single processor:
! #F90C     = f95
! #FFLAGS = -O2 -cpu:athlon -s -lU77 -w -YEXT_NAMES="LCS" -YEXT_SFX="_"
! 
! #NAGF95, single processor:
! #F90C     = f95
! #FFLAGS = -DNAGF95 -O3
! 
! #G95 compiler
! #F90C   = g95
! #FFLAGS = -O2
! 
! #Sun, single processor:
! #F90C     = f90
! #FFLAGS = -O5
! 
! #Sun V880
! #F90C = mpf90
! #FFLAGS =  -O4 -openmp -ftrap=%none -dalign -DMPI
! 
! #Sun parallel enterprise:
! #F90C     = f95
! #FFLAGS =  -O2 -xarch=native64 -openmp -ftrap=%none
! #try removing -openmp if get bus errors. -03, -04 etc are dodgy.
! 
! #IBM XL Fortran, multi-processor (run gmake)
! #F90C     = xlf90_r
! #FFLAGS  = -qsmp=omp -qsuffix=f=f90:cpp=F90 -O3 -qstrict -qarch=pwr3 -qtune=pwr3
  
  
  #Files containing evolution equations initial power spectrum module
--- 1,33 ----
! # >>> DESIGNED FOR GMAKE <<<
! 
! #CAMB System unified Makefile
! 
! ext=$(shell uname | cut -c1-3)
! 
! ifeq ($(ext),IRI)
! F90C=f90
! FFLAGS = -O3 -mp -64
! endif
! 
! ifeq ($(ext),Lin)
! F90C=ifort
! FFLAGS= -O -openmp  -fpp
! endif 
! 
! ifeq ($(ext),OSF)
! F90C=f90
! FFLAGS= -omp -O -arch host -math_library fast -tune host -fpe1
! endif
! 
! ifeq ($(ext),Sun)
! F90C=f90
! FFLAGS= -O4 -xarch=native64 -openmp -ftrap=%none
! endif
! 
! ifeq ($(ext),AIX)
! F90C=xlf90_r
! FFLAGS= -O -q64 -qsmp=omp -qthreaded  -qmaxmem=-1 -qfree=f90 -qsuffix=f=f90:cpp=F90
! endif
  
  
  #Files containing evolution equations initial power spectrum module
***************
*** 83,99 ****
  CAMBOBJ       = subroutines.o inifile.o $(POWERSPECTRUM).o recfast.o modules.o \
  	bessels.o $(EQUATIONS).o $(NONLINEAR).o lensing.o cmbmain.o camb.o
  
! default: camb
  
! all: camb $(CAMBLIB)
  
! camb: $(CAMBOBJ) $(DRIVER)
  	$(F90C) $(F90FLAGS) $(CAMBOBJ) $(DRIVER) -o $@
  
  $(CAMBLIB): $(CAMBOBJ)
  	ar -r $@ $?
  
! camb_fits: writefits.f90 $(CAMBOBJ) $(DRIVER)
  	$(F90C) $(F90FLAGS) -I$(HEALPIXDIR)/include $(CAMBOBJ) writefits.f90 $(DRIVER) $(HEALPIXLD) -DWRITE_FITS -o $@
  
  %.o: %.f90
--- 59,75 ----
  CAMBOBJ       = subroutines.o inifile.o $(POWERSPECTRUM).o recfast.o modules.o \
  	bessels.o $(EQUATIONS).o $(NONLINEAR).o lensing.o cmbmain.o camb.o
  
! default: camb.$(ext)
  
! all: camb.$(ext) $(CAMBLIB)
  
! camb.$(ext): $(CAMBOBJ) $(DRIVER)
  	$(F90C) $(F90FLAGS) $(CAMBOBJ) $(DRIVER) -o $@
  
  $(CAMBLIB): $(CAMBOBJ)
  	ar -r $@ $?
  
! camb_fits.$(ext): writefits.f90 $(CAMBOBJ) $(DRIVER)
  	$(F90C) $(F90FLAGS) -I$(HEALPIXDIR)/include $(CAMBOBJ) writefits.f90 $(DRIVER) $(HEALPIXLD) -DWRITE_FITS -o $@
  
  %.o: %.f90
diff -r -N -c -b -B cosmomc/camb/modules.f90 cosmomc_addon/camb/modules.f90
*** cosmomc/camb/modules.f90	2006-03-16 11:07:20.000000000 +0000
--- cosmomc_addon/camb/modules.f90	2006-04-04 10:38:50.000000000 +0100
***************
*** 186,192 ****
        integer, parameter:: l0max=4000
  
  !     lmax is max possible number of l's evaluated
!       integer, parameter :: lmax_arr = 100+l0max/7
   
          contains
        
--- 186,197 ----
        integer, parameter:: l0max=4000
  
  !     lmax is max possible number of l's evaluated
! !addon
! !      integer, parameter :: lmax_arr = 100+l0max/7
!       integer, parameter :: lmax_arr = 100 + l0max
!       logical :: DoClsInterpol=.true.
!       logical :: ComputeEveryl=.true.
! !end addon
   
          contains
        
***************
*** 539,544 ****
--- 544,567 ----
          integer lind, lvar, step,top,bot,ls(lmax_arr)
          real(dl) AScale
      
+ !addon    
+         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 addon
          Ascale=scale/lSampleBoost       
        
          lind=0
***************
*** 641,646 ****
--- 664,673 ----
          lSet%l0=lind
          lSet%l(1:lind) = ls(1:lind)
          
+ !addon
+         endif
+ !end addon
+         
        end subroutine initlval
  
        subroutine InterpolateClArr(lSet,iCl, all_Cl, max_ind)
***************
*** 660,665 ****
--- 687,694 ----
        xl = real(lSet%l(1:lSet%l0),dl)
        call spline(xl,iCl(1),max_ind,cllo,clhi,ddCl(1))
       
+ !addon
+       if (DoClsInterpol) then
              llo=1
              do il=lmin,lSet%l(max_ind)
                 xi=il
***************
*** 676,681 ****
--- 705,727 ----
                
              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 addon
+ 
+ 
        end subroutine InterpolateClArr
  
   
diff -r -N -c -b -B cosmomc/camb/params.ini cosmomc_addon/camb/params.ini
*** cosmomc/camb/params.ini	2005-03-30 16:05:28.000000000 +0100
--- cosmomc_addon/camb/params.ini	2006-04-04 10:38:57.000000000 +0100
***************
*** 64,69 ****
--- 64,76 ----
  
  #Initial power spectrum, amplitude, spectral index and running
  initial_power_num         = 1
+ 
+ #read the power spectra form ascii files in the form (Mpc^-1): k P(k)
+ use_external_scalar       = F
+ use_external_tensor       = F
+ scalar_input_file(1)      = scalar_power.dat
+ tensor_input_file(1)      = tensor_power.dat
+ 
  scalar_amp(1)             = 2.3e-9
  scalar_spectral_index(1)  = 1
  scalar_nrun(1)            = 0
***************
*** 100,105 ****
--- 107,118 ----
  #Matter power spectrum output against k/h in units of h^3 Mpc^{-3}
  transfer_matterpower(1) = matterpower.dat
  
+ #Use spline interpolation for the Cls
+ do_cls_interpol = T
+ 
+ #Remove undersampling of the Cls: if true lSampleBoost has no effect
+ compute_every_l = F
+ 
  
  #Output files not produced if blank. make camb_fits to use use the FITS setting.
  scalar_output_file = scalCls.dat
diff -r -N -c -b -B cosmomc/camb/power_tilt.f90 cosmomc_addon/camb/power_tilt.f90
*** cosmomc/camb/power_tilt.f90	2005-03-31 00:37:14.000000000 +0100
--- cosmomc_addon/camb/power_tilt.f90	2006-02-20 20:35:53.000000000 +0000
***************
*** 51,60 ****
         real(dl) k_0_scalar, k_0_tensor
         real(dl) ScalarPowerAmp(nnmax)
   
        end Type InitialPowerParams
  
!       real(dl) curv  !Curvature contant, set in InitializePowers     
  
        
        Type(InitialPowerParams) :: P
    
--- 51,80 ----
         real(dl) k_0_scalar, k_0_tensor
         real(dl) ScalarPowerAmp(nnmax)
   
+ !addon 
+        logical :: WantExternalScalar
+        logical :: WantExternalTensor
+        character(80) :: ScalarPowerFileName(nnmax),TensorPowerFileName(nnmax)
+ !end addon      
+ 
        end Type InitialPowerParams
  
! !addon
!       type ExternalData
! !only pointers may be defered shape in derived type
!          real(dl) :: ScalarMomentumMin, ScalarMomentumMax
!          real(dl), dimension(:), pointer :: PtrScalarPower, PtrScalarLnMomentum
!          real(dl), dimension(:), pointer :: PtrScalarPower2
!          real(dl) :: TensorMomentumMin, TensorMomentumMax
!          real(dl), dimension(:), pointer :: PtrTensorPower, PtrTensorLnMomentum
!          real(dl), dimension(:), pointer :: PtrTensorPower2
!       end type ExternalData
  
+       type(ExternalData), dimension(:), pointer  :: PtrExternalData => null()
+ !end addon      
+ 
+ 
+       real(dl) curv  !Curvature contant, set in InitializePowers     
        
        Type(InitialPowerParams) :: P
    
***************
*** 63,68 ****
--- 83,91 ----
        public  InitializePowers, ScalarPower, TensorPower
        public nnmax,InitialPowerParams,Power_Descript, Power_Name, SetDefPowerParams
  
+ !addon
+       public FreeExternalData
+ !end addon
      contains
         
    
***************
*** 77,82 ****
--- 100,110 ----
           AP%k_0_scalar = 0.05
           AP%k_0_tensor = 0.05
           AP%ScalarPowerAmp = 1
+ !addon
+ !         nullify(PtrExternalData)
+          AP%WantExternalScalar = .false.
+          AP%WantExternalTensor = .false.
+ !end addon        
  
         end subroutine SetDefPowerParams
  
***************
*** 87,92 ****
--- 115,128 ----
  
          real(dl) acurv
  
+ !addon
+         integer :: i,ioflag,j
+         integer :: nPowerUnit=111
+         integer :: nData
+         real(dl), save :: sbuffer1,sbuffer2
+         real(dl), parameter :: spl_large=1.e40_dl     
+ !end addon
+ 
          P = AParamSet
          if (P%nn > nnmax) then
             write (*,*) 'To use ',P%nn,'power spectra you need to increase'
***************
*** 97,105 ****
--- 133,295 ----
  
  !Write implementation specific code here...
  
+ !addon
+ 
+ !create and point to P%nn ExternalData objects (pointers and reals)
+         if ((P%WantExternalScalar).or.(P%WantExternalTensor)) then
+            if (.not.associated(PtrExternalData)) then
+               allocate(PtrExternalData(P%nn)) 
+            else
+               write(*,*)'InitializePowers: PtrExternalData already Initialized!'
+               return
+            endif
+         endif
+          
+         if (P%WantExternalScalar) then
+                         
+            do i=1,P%nn
+                
+               write(*,*)'using primordial scalar power spectrum: ',P%ScalarPowerFileName(i)
+                
+               open(unit=nPowerUnit,file=P%ScalarPowerFileName(i),action='read' &
+                    ,form='formatted',status='old')             
+ !count the number of records
+               nData=0
+               do                                    
+                  read(nPowerUnit,*,iostat=ioflag)                  
+                  if (ioflag.lt.0) exit
+                  nData=nData + 1                                    
+               enddo
+ !create and point to
+               allocate(PtrExternalData(i)%PtrScalarPower(nData) & 
+                    ,PtrExternalData(i)%PtrScalarLnMomentum(nData))
+               
+               rewind(nPowerUnit)
+               nData=1
+               do                                   
+                  read(nPowerUnit,*,iostat=ioflag) sbuffer1,sbuffer2
+                  if (ioflag.lt.0) exit
+                  PtrExternalData(i)%PtrScalarLnMomentum(nData) = log(sbuffer1)
+                  PtrExternalData(i)%PtrScalarPower(nData) = sbuffer2      
+                  nData=nData+1
+               enddo
+               nData=nData - 1 
+               close(nPowerUnit,iostat=ioflag)
+             
+ !keep trace of the accessible range              
+               PtrExternalData(i)%ScalarMomentumMin = exp(PtrExternalData(i)%PtrScalarLnMomentum(1))
+               PtrExternalData(i)%ScalarMomentumMax = exp(PtrExternalData(i)%PtrScalarLnMomentum(nData))
+               
+ 
+ !2nd order derivatives allocation
+               allocate(PtrExternalData(i)%PtrScalarPower2(nData))
+ !splining in log k
+               call spline(PtrExternalData(i)%PtrScalarLnMomentum(:) &
+                    ,PtrExternalData(i)%PtrScalarPower(:),nData, spl_large &
+                    , spl_large,PtrExternalData(i)%PtrScalarPower2(:))
+               
+            enddo
+ 
+         endif
+ 
+ 
+ !the same for tensors
+ 
+         if (P%WantExternalTensor) then
+             
+            do i=1,P%nn
+               
+               write(*,*)'using primordial tensor power spectrum: ',P%TensorPowerFileName(i)
+ 
+               open(unit=nPowerUnit,file=P%TensorPowerFileName(i) &
+                    ,action='read',form='formatted',status='old')     
+ !count the number of records
+               nData=0
+               do                                    
+                  read(nPowerUnit,*,iostat=ioflag)                  
+                  if (ioflag.lt.0) exit
+                  nData=nData + 1                                    
+               enddo
+ !store the data
+               allocate(PtrExternalData(i)%PtrTensorPower(nData) & 
+                    ,PtrExternalData(i)%PtrTensorLnMomentum(nData))
+                
+               rewind(nPowerUnit)
+               nData=1
+               do                                   
+                  read(nPowerUnit,*,iostat=ioflag) sbuffer1,sbuffer2
+                  if (ioflag.lt.0) exit
+                  PtrExternalData(i)%PtrTensorLnMomentum(nData) = log(sbuffer1)
+                  PtrExternalData(i)%PtrTensorPower(nData) = sbuffer2      
+                  nData=nData+1
+               enddo
+               nData=nData - 1 
+               close(nPowerUnit)
+ 
+ !keep trace of the accessible range              
+               PtrExternalData(i)%TensorMomentumMin = exp(PtrExternalData(i)%PtrTensorLnMomentum(1))
+               PtrExternalData(i)%TensorMomentumMax = exp(PtrExternalData(i)%PtrTensorLnMomentum(nData))
+               
+ 
+ !2nd order derivatives
+               allocate(PtrExternalData(i)%PtrTensorPower2(nData))
+ !splining in log k
+               call spline(PtrExternalData(i)%PtrTensorLnMomentum(:) &
+                    ,PtrExternalData(i)%PtrTensorPower(:),nData, spl_large &
+                    , spl_large,PtrExternalData(i)%PtrTensorPower2(:))
+               
+            enddo
+ 
+             
+         endif
+          
         end subroutine InitializePowers
         
  
+ 
+        subroutine FreeExternalData(AP)
+          type (InitialPowerParams), intent(in) :: AP
+          integer :: i
+ 
+          if (AP%nn.ne.size(PtrExternalData,1)) then
+             write(*,*)'InitialPower.FreeExternalData : wrong initial power params'
+             stop
+          endif
+ 
+          if (AP%WantExternalScalar) then
+             do i=1, AP%nn               
+ !free heap arrays and disassociate ExternalData pointers            
+                deallocate(PtrExternalData(i)%PtrScalarPower)
+                deallocate(PtrExternalData(i)%PtrScalarPower2)
+                deallocate(PtrExternalData(i)%PtrScalarLnMomentum)
+ !should still point to ExternalData objects
+                if (.not.associated(PtrExternalData)) stop 'InitialPower.FreeExternalData : ooops'
+             enddo
+          endif
+          if (AP%WantExternalTensor) then
+             do i=1, AP%nn               
+                deallocate(PtrExternalData(i)%PtrTensorPower)
+                deallocate(PtrExternalData(i)%PtrTensorPower2)
+                deallocate(PtrExternalData(i)%PtrTensorLnMomentum)
+                if (.not.associated(PtrExternalData)) stop 'InitialPower.FreeExternalData : ooops'
+             enddo
+          endif
+ !Free remaining objects and disassociate pointer to ExternalData
+          if (associated(PtrExternalData)) then
+             deallocate(PtrExternalData)
+             PtrExternalData => null()
+          else
+             if (AP%WantExternalScalar.or.AP%WantExternalTensor) then
+ !should never happen               
+                stop 'InitialPower.FreeExternalData : ooops'
+             endif
+          endif
+ 
+        end subroutine FreeExternalData
+ 
+ !end addon
+ 
+ 
        function ScalarPower(k,in)
  
         !"in" gives the index of the power to return for this k
***************
*** 120,128 ****
--- 310,340 ----
          real(dl) ScalarPower,k, lnrat
          integer in
  
+ !addon
+         integer :: nData       
+ 
+         if (P%WantExternalScalar) then
+            nData=size(PtrExternalData(in)%PtrScalarLnMomentum,1)
+ 
+            if ((k.ge.PtrExternalData(in)%ScalarMomentumMin) &
+                 .and.(k.le.PtrExternalData(in)%ScalarMomentumMax)) then
+               call splinterpol(PtrExternalData(in)%PtrScalarLnMomentum &
+                    ,PtrExternalData(in)%PtrScalarPower &
+                    ,PtrExternalData(in)%PtrScalarPower2,nData,log(k),ScalarPower)
+ !confusing              ScalarPower = P%ScalarPowerAmp(in)*ScalarPower
+            else
+               write(*,*)'InitPower.ScalarPower() : out of range: k= ',k
+               ScalarPower=0_dl              
+            endif
+ 
+         else
+ 
            lnrat = log(k/P%k_0_scalar)
            ScalarPower=P%ScalarPowerAmp(in)*exp((P%an(in)-1)*lnrat + P%n_run(in)/2*lnrat**2)   
   
+         endif
+ !end addon
+ 
        end function ScalarPower
  
        
***************
*** 143,149 ****
--- 355,384 ----
     
          integer in
  
+ !addon
+         integer :: nData
+         
+ 
+         if (P%WantExternalTensor) then
+            nData=size(PtrExternalData(in)%PtrTensorLnMomentum,1)
+            if ((k.ge.PtrExternalData(in)%TensorMomentumMin) &
+                 .and.(k.le.PtrExternalData(in)%TensorMomentumMax)) then
+               call splinterpol(PtrExternalData(in)%PtrTensorLnMomentum &
+                    ,PtrExternalData(in)%PtrTensorPower &
+                    ,PtrExternalData(in)%PtrTensorPower2,nData,log(k),TensorPower)
+ 
+ !confusing              TensorPower = P%rat(in)*P%ScalarPowerAmp(in)*TensorPower
+            else
+               write(*,*)'InitPower.TensorPower() : out of range: k= ',k
+               TensorPower=0_dl               
+            endif
+         else
+ 
          TensorPower=P%rat(in)*P%ScalarPowerAmp(in)*exp(P%ant(in)*log(k/P%k_0_tensor))
+ 
+         endif
+ !end addon
+         
          if (curv < 0) TensorPower=TensorPower*tanh(PiByTwo*sqrt(-k**2/curv-3)) 
  
        end function TensorPower
diff -r -N -c -b -B cosmomc/camb/subroutines.f90 cosmomc_addon/camb/subroutines.f90
*** cosmomc/camb/subroutines.f90	2005-03-31 00:37:14.000000000 +0100
--- cosmomc_addon/camb/subroutines.f90	2006-04-04 10:39:04.000000000 +0100
***************
*** 962,964 ****
--- 962,997 ----
  !  end abort action
  !
        end subroutine dverk
+ 
+ 
+ !addon
+ !the numrec splint renamed
+ 
+ SUBROUTINE splinterpol(xa,ya,y2a,n,x,y)
+   use Precision
+   INTEGER n
+   REAL(dl) x,y,xa(n),y2a(n),ya(n)
+   INTEGER k,khi,klo
+   REAL(dl) a,b,h
+   klo=1
+   khi=n
+ 1 if (khi-klo.gt.1) then
+      k=(khi+klo)/2
+      if(xa(k).gt.x)then
+         khi=k
+      else
+         klo=k
+      endif
+      goto 1
+   endif
+   h=xa(khi)-xa(klo)
+   if (h.eq.0.) stop 'bad xa input in splint'
+   a=(xa(khi)-x)/h
+   b=(x-xa(klo))/h
+   y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))&
+        *(h**2)/6.
+   return
+ END SUBROUTINE splinterpol
+ 
+ 
+ !end addon
diff -r -N -c -b -B cosmomc/distparams.ini cosmomc_addon/distparams.ini
*** cosmomc/distparams.ini	2006-03-17 23:39:38.000000000 +0000
--- cosmomc_addon/distparams.ini	2006-04-04 10:34:35.000000000 +0100
***************
*** 32,37 ****
--- 32,38 ----
  
  plot_meanlikes = T
  shade_meanlikes = T
+ td_meanlikes = F
  
  # if non-zero, output _thin file, thinned by thin_factor
  thin_factor = 0
diff -r -N -c -b -B cosmomc/params.ini cosmomc_addon/params.ini
*** cosmomc/params.ini	2006-04-03 14:00:41.000000000 +0100
--- cosmomc_addon/params.ini	2006-04-04 10:34:35.000000000 +0100
***************
*** 117,123 ****
  #If true, generate checkpoint files and terminated runs can be restarted using exactly the same command
  #and chains continued from where they stopped
  #With checkpoint=T note you must delete all chains/file_root.* files if you want new chains with an old file_root
! checkpoint = F
  
  
  #CAMB parameters
--- 117,123 ----
  #If true, generate checkpoint files and terminated runs can be restarted using exactly the same command
  #and chains continued from where they stopped
  #With checkpoint=T note you must delete all chains/file_root.* files if you want new chains with an old file_root
! checkpoint = T
  
  
  #CAMB parameters
***************
*** 133,138 ****
--- 133,141 ----
  #increase accuracy_level to run CAMB on higher accuracy
  #(default is about 0.3%, accuracy_level=2 around 0.1% at high l)
  accuracy_level = 1
+ l_accuracy_level = 1
+ l_sample_level = 1
+ camb_compute_every_l = F
  
  #If action = 1
  redo_likelihoods = T
diff -r -N -c -b -B cosmomc/params.ini.orig cosmomc_addon/params.ini.orig
*** cosmomc/params.ini.orig	1970-01-01 01:00:00.000000000 +0100
--- cosmomc_addon/params.ini.orig	2006-04-04 10:34:22.000000000 +0100
***************
*** 0 ****
--- 1,181 ----
+ #Sample parameters for cosmomc in default parameterization
+ 
+ #Root name for files produced
+ file_root = chains/test
+ 
+ #action = 0:  MCMC, action=1: postprocess .data file, action=2: find best fit point only
+ action = 0
+ 
+ #Maximum number of chain steps
+ samples = 200000
+ 
+ #Feedback level ( 2=lots,1=chatty,0=none)
+ feedback = 1
+ 
+ #Temperature at which to Monte-Carlo
+ temperature = 1
+ 
+ #filenames for CMB datasets
+ #Note you may need to change lmax in cmbtypes.f90 to use small scales (e.g. lmax=2100)
+ cmb_numdatasets = 1
+ cmb_dataset1 = WMAP
+ 
+ cmb_dataset2 = data/ACBAR_lge800_.dataset
+ cmb_dataset3 = data/CBI2.0_lge600.dataset
+ cmb_dataset4 = data/B03_NA_21July05.newdat
+ 
+ #filenames for matter power spectrum datasets, incl twodf
+ mpk_numdatasets = 1
+ mpk_dataset1 = twodf
+ #mpk_dataset1 = data/sdss.dataset
+ 
+ #if true, use HALOFIT for non-linear corrections (astro-ph/0207664).
+ #note lyman-alpha (lya) code assumes linear spectrum
+ nonlinear_pk = F
+ 
+ use_CMB = T
+ use_HST = F
+ use_mpk = F
+ use_clusters = F
+ use_BBN = F
+ use_Age_Tophat_Prior = T
+ use_SN = F
+ use_lya = F
+ use_min_zre = 0
+ 
+ #Force computation of sigma_8 even if use_mpk = F
+ get_sigma8 = F
+ 
+ #1: Simple Metropolis, 2: slice sampling, 3: slice sampling fast parameters, 4: directional gridding
+ sampling_method = 1
+ 
+ #if sampling_method =4, iterations per gridded direction
+ directional_grid_steps = 20
+ 
+ #use fast-slow parameter distinctions to speed up 
+ #(note for basic models WMAP3 code is only ~3x as fast as CAMB)
+ use_fast_slow = F
+ 
+ #Can use covariance matrix for proposal density, otherwise use settings below
+ #Covariance matrix can be produced using "getdist" program.
+ propose_matrix = params_CMB.covmat
+ 
+ #If propose_matrix is blank (first run), can try to use numerical Hessian to 
+ #estimate a good propose matrix. As a byproduct you also get an approx best fit point
+ estimate_propose_matrix = F
+ 
+ #Tolerance on log likelihood to use when estimating best fit point
+ delta_loglike = 2
+ 
+ #Scale of proposal relative to covariance; 2.4 is recommended by astro-ph/0405462 for Gaussians
+ #If propose_matrix is much broader than the new distribution, make proportionately smaller
+ #Generally make smaller if your acceptance rate is too low
+ propose_scale = 2.4
+ 
+ #Increase to oversample fast parameters more, e.g. if space is odd shape
+ oversample_fast = 1
+ 
+ #if non-zero number of steps between sample info dumped to file file_root.data
+ indep_sample = 0
+ 
+ #number of samples to disgard at start; usually set to zero and remove later
+ burn_in = 0
+ 
+ #If zero set automatically
+ num_threads = 0
+ 
+ #MPI mode multi-chain options (recommended)
+ #MPI_Converge_Stop is a (variance of chain means)/(mean of variances) parameter that can be used to stop the chains
+ #Set to a negative number not to use this feature. Does not guarantee good accuracy of confidence limits.
+ MPI_Converge_Stop = 0.03
+ 
+ #Do initial period of slice sampling; may be good idea if 
+ #cov matrix or widths are likely to be very poor estimates
+ MPI_StartSliceSampling  = F
+ 
+ #Can optionally also check for convergence of confidence limits (after MPI_Converge_Stop reached)
+ #Can be good idea as small value of MPI_Converge_Stop does not (necessarily) imply good exploration of tails
+ MPI_Check_Limit_Converge = F
+ 
+ #if MPI_Check_Limit_Converge = T, give tail fraction to check (checks both tails):
+ MPI_Limit_Converge = 0.025
+ #permitted quantile chain variance in units of the standard deviation (small values v slow):
+ MPI_Limit_Converge_Err = 0.2
+ #which parameter's tails to check. If zero, check all parameters:
+ MPI_Limit_Param = 0
+ 
+ #if MPI_LearnPropose = T, the proposal density is continally updated from the covariance of samples so far (since burn in)
+ MPI_LearnPropose = T
+ #can set a value of converge at which to stop updating covariance (so that it becomes rigorously Markovian)
+ #e.g. MPI_R_StopProposeUpdate = 0.4 will stop updating when (variance of chain means)/(mean of variances) < 0.4
+ MPI_R_StopProposeUpdate = 0
+ 
+ 
+ #if blank this is set from system clock
+ rand_seed = 
+ 
+ #If true, generate checkpoint files and terminated runs can be restarted using exactly the same command
+ #and chains continued from where they stopped
+ #With checkpoint=T note you must delete all chains/file_root.* files if you want new chains with an old file_root
+ checkpoint = F
+ 
+ 
+ #CAMB parameters
+ #If we are including tensors
+ compute_tensors = F
+ #Initial power spectrum amplitude point (Mpc^{-1})
+ pivot_k = 0.05
+ #If using tensors, enforce n_T = -A_T/(8A_s)
+ inflation_consistency = F
+ 
+ #Whether the CMB should be lensed (slows a lot unless also computing matter power)
+ CMB_lensing = F
+ #increase accuracy_level to run CAMB on higher accuracy
+ #(default is about 0.3%, accuracy_level=2 around 0.1% at high l)
+ accuracy_level = 1
+ 
+ #If action = 1
+ redo_likelihoods = T
+ redo_theory = F
+ redo_cls = F
+ redo_pk = F
+ redo_skip = 0
+ redo_outroot = 
+ redo_thin = 1
+ redo_add = F
+ redo_from_text = F
+ #If large difference in log likelihoods may need to offset to give sensible weights
+ #for exp(difference in likelihoods)
+ redo_likeoffset =  0
+ 
+ #parameter start center, min, max, start width, st. dev. estimate
+ #ombh2
+ param1 = 0.0223 0.005 0.1 0.001 0.001
+ #omdmh2
+ param2 = 0.105 0.01 0.99 0.01 0.01
+ #theta
+ param3 = 1.04 0.5 10 0.002 0.002
+ #tau
+ param4 = 0.09 0.01 0.8 0.03 0.03
+ 
+ #omk
+ param5 = 0 0 0 0 0
+ #nufrac
+ param6 = 0 0 0 0 0
+ #w
+ param7 = -1 -1 -1 0 0
+ 
+ #n_s
+ param8 = 0.95 0.5 1.5 0.03 0.01
+ #n_t
+ param9 = 0 0 0 0 0
+ #n_run
+ param10 = 0 0 0 0 0
+ 
+ #log[10^10 A_s]
+ param11 = 3 2.7 4 0.01 0.01
+ #amp_ratio
+ param12 = 0 0 0 0 0
+ #unused fast parameter
+ param13 = 0 0 0 0 0
+ 
diff -r -N -c -b -B cosmomc/source/CMB_Cls_simple.f90 cosmomc_addon/source/CMB_Cls_simple.f90
*** cosmomc/source/CMB_Cls_simple.f90	2006-03-17 16:18:11.000000000 +0000
--- cosmomc_addon/source/CMB_Cls_simple.f90	2006-04-04 10:34:35.000000000 +0100
***************
*** 266,271 ****
--- 266,289 ----
     type(CAMBParams)  P 
     integer zix
  
+ !addon
+    AccuracyBoost = AccuracyLevel
+    lAccuracyBoost = lAccuracyLevel
+    lSampleBoost = lSampleLevel
+    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(*,*)
+ !end addon
+ 
          Threadnum =num_threads
          w_lam = -1
          call CAMB_SetDefParams(P)
***************
*** 294,302 ****
              P%Transfer%high_precision=.true.
              P%Transfer%k_per_logint=3*nint(AccuracyLevel)
            end if
!           AccuracyBoost = AccuracyLevel
!           lAccuracyBoost = AccuracyLevel
!           lSampleBoost = AccuracyLevel
            P%AccurateReionization = .true.
          end if
                  
--- 312,323 ----
              P%Transfer%high_precision=.true.
              P%Transfer%k_per_logint=3*nint(AccuracyLevel)
            end if
! !addon
! !moved up for screen output
! !          AccuracyBoost = AccuracyLevel
! !          lAccuracyBoost = AccuracyLevel
! !          lSampleBoost = AccuracyLevel
! !end addon
            P%AccurateReionization = .true.
          end if
                  
diff -r -N -c -b -B cosmomc/source/driver.F90 cosmomc_addon/source/driver.F90
*** cosmomc/source/driver.F90	2006-03-17 21:11:02.000000000 +0000
--- cosmomc_addon/source/driver.F90	2006-04-04 10:34:35.000000000 +0100
***************
*** 50,61 ****
  
  #ifdef MPI 
  
!         if (instance /= 0) call DoStop('With MPI should not have second parameter')
          call mpi_comm_rank(mpi_comm_world,MPIrank,ierror)
  
  
          instance = MPIrank +1 !start at 1 for chains                                             
!         write (numstr,*) instance
          rand_inst = instance
          if (ierror/=MPI_SUCCESS) call DoStop('MPI fail')
  
--- 50,69 ----
  
  #ifdef MPI 
  
! !addon
! !        if (instance /= 0) call DoStop('With MPI should not have second parameter')
!         instance_shift=instance
!         write(*,*)'MPI file names shifted by: ',instance_shift
! !end addon
          call mpi_comm_rank(mpi_comm_world,MPIrank,ierror)
  
  
          instance = MPIrank +1 !start at 1 for chains                                             
! !addon
! !        write (numstr,*) instance
!         instance_shift = MPIrank + instance_shift
!         write (numstr,*) instance_shift
! !end addon
          rand_inst = instance
          if (ierror/=MPI_SUCCESS) call DoStop('MPI fail')
  
***************
*** 77,82 ****
--- 85,95 ----
  
          propose_scale = Ini_Read_Real('propose_scale',2.4)
          AccuracyLevel = Ini_Read_Real('accuracy_level',1.)
+ !addon
+         lAccuracyLevel = Ini_Read_Real('l_accuracy_level',1.)
+         lSampleLevel = Ini_Read_Real('l_sample_level',1.)
+         CambComputeEveryl = Ini_Read_Logical('camb_compute_every_l',.false.)
+ !end addon
  
          checkpoint = Ini_Read_Logical('checkpoint',.false.)
          if (checkpoint) flush_write = .true.
diff -r -N -c -b -B cosmomc/source/GetDist.f90 cosmomc_addon/source/GetDist.f90
*** cosmomc/source/GetDist.f90	2006-04-02 23:49:21.000000000 +0100
--- cosmomc_addon/source/GetDist.f90	2006-04-04 10:34:35.000000000 +0100
***************
*** 67,72 ****
--- 67,75 ----
          logical force_twotail 
          logical smoothing, plot_meanlikes, plot_NDcontours
          logical shade_meanlikes, make_single_samples
+ !addon
+         logical td_meanlikes
+ !end addon
          integer single_thin,cust2DPlots(max_cols**2)
          integer ix_min(max_cols),ix_max(max_cols)
          real center(max_cols),width(max_cols)   
***************
*** 1146,1160 ****
   !             write (aunit,fmt) (/(I, I=ix_min(j), ix_max(j), 1)/)*width(j)+center(j)
                if (DoShade) then
                if (shade_meanlikes) then
                  write (aunit,'(a)') 'load (fullfile(''plot_data'',''' // trim(plotfile) //'_likes''));'
                  write (aunit,*) 'contourf(x1,x2,'//trim(plotfile)//'_likes,64);'
                else     
                  write (aunit,*) 'contourf(x1,x2,pts,64);'
                end if
! 
                write (aunit,*) 'set(gca,''climmode'',''manual''); shading flat; hold on;'
!               end if
!               
             
                if (num_contours /= 0) then
  
--- 1149,1176 ----
   !             write (aunit,fmt) (/(I, I=ix_min(j), ix_max(j), 1)/)*width(j)+center(j)
                if (DoShade) then
                if (shade_meanlikes) then
+ !addon
                  write (aunit,'(a)') 'load (fullfile(''plot_data'',''' // trim(plotfile) //'_likes''));'
+                  if (td_meanlikes) then
+                     write (aunit,*) 'surfl(x1,x2,'//trim(plotfile)//'_likes);'
+                  else
                  write (aunit,*) 'contourf(x1,x2,'//trim(plotfile)//'_likes,64);'
+                  endif
+ !end addon
                else     
                  write (aunit,*) 'contourf(x1,x2,pts,64);'
                end if
! !addon
! !              write (aunit,*) 'set(gca,''climmode'',''manual''); shading flat; hold on;'
!               if (td_meanlikes) then
!                  write (aunit,*) 'set(gca,''climmode'',''manual'');shading interp; hold on;'
!                  write (aunit,*) 'grid off;box on;camproj(''perspective'');'
!                  write (aunit,*) 'set(gca,''climmode'',''manual'');alpha(0.9);axis tight'
!               else
                write (aunit,*) 'set(gca,''climmode'',''manual''); shading flat; hold on;'
!               endif
! !end addon
!               endif
             
                if (num_contours /= 0) then
  
***************
*** 1449,1454 ****
--- 1465,1473 ----
            write (*,*) 'WARNING: Mapping params - .covmat file is new params.'
  
      shade_meanlikes = Ini_Read_Logical('shade_meanlikes',.false.)   
+ !addon
+     td_meanlikes = Ini_Read_Logical('td_meanlikes',.false.)  
+ !end addon
  
      matlab_col = Ini_Read_String('matlab_colscheme')
      smoothing = Ini_Read_Logical('smoothing',.true.)
***************
*** 1955,1961 ****
--- 1974,1987 ----
                write (51,*)  'set(gcf, ''PaperUnits'',''inches'');'
       
                if (plot_row < plot_col .and. plot_row*plot_col>9) then
+ !addon
+ !            write (51,*) 'set(gcf, ''PaperPosition'',[ 0 0 8 10]);';
+          if (td_meanlikes) then
+             write (51,*) 'set(gcf, ''PaperPosition'',[ 0 0 7 9]);';
+          else
                  write (51,*) 'set(gcf, ''PaperPosition'',[ 0 0 8 10]);';
+          endif
+ !end addon           
                else
                  write (51,*) 'set(gcf, ''PaperPosition'',[ 0 0 8 8]);';
                end if
diff -r -N -c -b -B cosmomc/source/GetDist.f90.orig cosmomc_addon/source/GetDist.f90.orig
*** cosmomc/source/GetDist.f90.orig	1970-01-01 01:00:00.000000000 +0100
--- cosmomc_addon/source/GetDist.f90.orig	2006-04-04 10:34:22.000000000 +0100
***************
*** 0 ****
--- 1,2197 ----
+ ! Program to process .txt chain files produced by CosmoMC (cosmologist.info/cosmomc)
+ ! Calculates statistics, thins, produces 1D and 2D marginalized plots and scatter plots
+ ! and writes out useful information.
+ ! Option to adjust priors and map parameters - write code below
+ 
+ !Output files are (depending on options)
+ 
+ ! * file_root.margestats - the mean, limits and stddev from the 1D marginalized distributions
+ ! * file_root.likestats  - the best fit, and limits from the N-D confidence region
+ ! * file_root.m          - MatLab .m file to produce 1D marginalized plots
+ ! * file_root.sm         - SuperMongo .sm file to produce 1D marginalized plots
+ ! * file_root_2D.m       - MatLab .m file to produce 2D marginalized plots
+ ! * file_root_3D.m       - MatLab .m file to produce 2D sample plots colored by a third parameter
+ ! * file_root_tri.m      - MatLab .m file to produce triangle plots (1D on diagonal, 2D off-diagonal)
+ ! * file_root_thin.txt   - thinned merged versions of input files
+ ! * file_root_single.txt - File of weight-1 single samples generated by probabilistic importance sampling
+ ! * file_root.covmat     - covariance matrix for parameters
+ ! * file_root.corr       - covariance matrix for parameters normalized to have unit variance
+ ! * file_root.PCA        - human-readable file giving details of PCA, constrained parameters, etc
+ ! * file_root.converge   - Various statistics to help assess chain convergence/sampling error
+ 
+ ! You can easily edit the .m and .sm files produced to produce custom layouts of plots, etc.
+ ! Data for plots are exported to the plot_data folder
+ 
+ !This version April 2006
+ !March 03: fixed bug computing the limits in the .margestats file
+ !May 03: Added support for triangle plots 
+ !Dec 03: Support for non-chain samples, auto-correlation convergence, auto_label,
+ !         split tests to quanify sampling error on upper and lower quantiles
+ !Jan 05: MatLab 7 options, various minor changes
+ !Jul 05: Added limits info to .margestats output
+ !Apr 06: Fixed some version confusions
+ module MCSamples
+  use settings
+  implicit none
+ 
+         integer, parameter :: gp = KIND(1.d0)
+ 
+         real(gp), dimension(:,:), target, allocatable :: coldata  !(col_index, row_index)
+         integer, dimension(:), allocatable :: thin_ix  
+         integer thin_rows
+ 
+         integer, parameter :: max_rows = 500000
+         integer, parameter :: max_cols = 500
+         integer, parameter :: max_chains = 100
+         integer, parameter :: max_split_tests = 4 
+ 
+         integer chain_indices(max_chains), num_chains_used
+ 
+         integer nrows, ncols
+         real numsamp
+         integer ND_cont1, ND_cont2
+         integer covmat_dimension
+         integer colix(max_cols), num_vars  !Parameters with non-blank labels
+         real mean(max_cols), sddev(max_cols)
+         real, dimension(:,:), allocatable :: corrmatrix
+         character(LEN=120) rootname
+         character(LEN=120) labels(max_cols)
+         logical has_limits(max_cols), has_limits_top(max_cols),has_limits_bot(max_cols)
+         integer num_contours 
+         real matlab_version
+         real contours(2), contour_levels(2)
+         real mean_mult, max_mult
+         integer :: indep_thin = 0 
+         integer chain_numbers(max_chains)
+         logical isused(max_cols)
+         logical force_twotail 
+         logical smoothing, plot_meanlikes, plot_NDcontours
+         logical shade_meanlikes, make_single_samples
+         integer single_thin,cust2DPlots(max_cols**2)
+         integer ix_min(max_cols),ix_max(max_cols)
+         real center(max_cols),width(max_cols)   
+         real, dimension(:,:), allocatable :: TheBins, bins2D, bin2Dlikes, bin2Dmax
+         real meanlike, maxlike, maxmult
+         logical BW,do_shading
+         character(LEN=80) ComparePlots(20)
+         integer Num_ComparePlots
+         logical ::   prob_label = .false.      
+             
+ contains
+ 
+   subroutine AdjustPriors
+    !Can adjust the multiplicity of each sample in coldata(1, rownum) for new priors
+    !Be careful as this code is parameterisation dependent
+    integer i
+    real tau, chisq
+ 
+    stop 'You need to write the AdjustPriors subroutine in GetDist.F90 first!'
+ 
+    write (*,*) 'Adjusting priors'
+    do i=0, nrows-1
+ !E.g. tau prior
+       tau = coldata(4+2,i)
+       chisq = (tau - 0.1)**2/0.03**2
+       coldata(1,i) = coldata(1,i)*exp(-chisq/2) 
+       coldata(2,i) = coldata(2,i) + chisq/2
+      
+    end do
+    
+ 
+   end subroutine AdjustPriors
+ 
+   subroutine MapParameters(invars)
+     real invars(1:ncols)
+ ! map parameters in invars: eg. invars(3)=invars(3)*invars(4) 
+ 
+ !    invars(2+13)=invars(17+2)*exp(-invars(2+4))
+    
+     stop 'Need to write MapParameters routine first'  
+   
+   end subroutine MapParameters
+ 
+   subroutine CoolChain(cool)
+     real, intent(in) :: cool
+     integer i
+     real maxL, newL
+     
+     write (*,*) 'Cooling chains by ', cool
+     MaxL = minval(coldata(2,0:nrows-1))
+     do i=0, nrows-1
+      newL = coldata(2,i)*cool
+      coldata(1,i) = coldata(1,i)*exp(-(newL - coldata(2,i)) - MaxL*(1-cool) )
+      coldata(2,i) = newL
+    end do
+ 
+   end subroutine CoolChain
+ 
+   subroutine DeleteZeros
+     integer i,ii
+ 
+     ii=0
+     do i=0, nrows-1
+      if (coldata(1,i)/=0) then
+        coldata(:,ii) = coldata(:,i)
+        ii=ii+1
+      end if
+      end do
+      if (ii==0) stop 'Prior has removed all models!'
+      if (ii /= nrows) write (*,*) 'Prior removed ',nrows-ii,'models'
+      nrows = ii
+           
+   end subroutine DeleteZeros
+ 
+ 
+ 
+   subroutine SortColData(bycol)
+  !Sort coldata in order of likelihood
+      Type(double_pointer), dimension(:), allocatable :: rowptrs
+      integer, intent(in) :: bycol
+      real(gp), dimension(:,:), pointer :: tmp 
+      integer i
+ 
+      allocate(tmp(ncols,0:nrows-1))
+   
+      allocate(rowptrs(0:nrows-1))
+      do i = 0, nrows -1
+       rowptrs(i)%p => coldata(1:,i)
+      end do
+ 
+     call QuickSortArr(rowptrs(0), 1, nrows, bycol)
+ 
+ !Make new table using sorted pointers   
+     do i = 0, nrows-1
+     tmp(1:ncols,i) = rowptrs(i)%p(1:ncols)
+     end do
+ 
+     deallocate(rowptrs)
+     do i=0, nrows-1
+      !Write this out explicitly to avoid stack overflow problems
+      coldata(1:ncols,i) = tmp(1:ncols,i)
+     end do
+     deallocate(tmp)
+ 
+   end  subroutine SortColData
+ 
+ 
+   subroutine MakeSingleSamples(single_thin)
+ 
+   !Make file of weight-1 samples by choosing samples with probability given by their weight
+ 
+     use Random
+     integer i, single_thin
+     character(LEN=20) :: fmt
+     
+     real maxmult
+     Feedback = 0          
+     call initRandom()
+     fmt = trim(numcat('(',ncols))//'E16.7)'
+     open(unit=50,file=trim(rootname)//'_single.txt',form='formatted',status='replace')
+     maxmult = maxval(coldata(1,0:nrows-1))
+     do i= 0, nrows -1
+      if (ranmar() <= coldata(1,i)/maxmult/single_thin) write (50,fmt) 1.0, coldata(2:ncols,i)          
+     end do         
+     close(50)
+ 
+   end subroutine MakeSingleSamples
+ 
+   subroutine WriteThinData(fname,cool)
+     character(LEN=*), intent(in) :: fname
+     real,intent(in) :: cool
+  
+     character(LEN=20) :: fmt
+     integer i
+     real MaxL, NewL      
+ 
+     if (cool /= 1) write (*,*) 'Cooled thinned output with temp: ', cool
+   
+     MaxL = minval(coldata(2,0:nrows-1))
+ 
+     fmt = trim(numcat('(',ncols))//'E16.7)'
+     open(unit=50,file=fname,form='formatted',status='replace')
+ 
+     do i=0, thin_rows -1   
+       if (cool/=1) then
+         newL = coldata(2,thin_ix(i))*cool
+         write (50,fmt) exp(-(newL - coldata(2,thin_ix(i))) - MaxL*(1-cool) ), newL, &
+                coldata(3:ncols,thin_ix(i))
+       else 
+         write (50,fmt) 1., coldata(2:ncols,thin_ix(i)) 
+       end if
+     end do
+     write (*,*) 'Wrote ',thin_rows, 'thinned samples'
+     close(50)
+ 
+   end subroutine WriteThinData
+ 
+ 
+   subroutine ThinData(fac, ix1,ix2)
+   !Make thinned samples
+  
+     integer, intent(in) :: fac
+     integer, intent(in), optional :: ix1,ix2
+     integer i, tot, nout, nend, mult
+           
+     if (allocated(thin_ix)) deallocate(thin_ix)
+     allocate(thin_ix(0:nint(numsamp)/fac))
+               
+     tot= 0
+     nout=0
+     i = 0
+     if (present(ix1)) i=ix1
+     nend = nrows
+     if (present(ix2)) nend = ix2+1
+  
+     mult = coldata(1,i)
+     do while (i< nend)
+         
+     if (abs(nint(coldata(1,i)) - coldata(1,i)) > 1e-4) &
+         stop 'non-integer weights in ThinData'
+ 
+     if (mult + tot < fac) then
+        tot = tot + mult
+        i=i+1
+        if (i< nend) mult = nint(coldata(1,i))
+     else
+         thin_ix(nout) = i
+         nout= nout+1
+        if (mult == fac - tot) then  
+           i=i+1
+           if (i< nend) mult = nint(coldata(1,i))
+        else
+          mult = mult - (fac -tot)
+        end if
+        tot = 0
+     end if
+ 
+     end do
+        
+     thin_rows = nout
+ 
+     close(50)
+ 
+   end subroutine ThinData
+ 
+   
+       subroutine GetCovMatrix
+         integer i, j
+         real mean(max_cols)
+         real scale
+         real, dimension(:,:), allocatable :: covmatrix
+                 
+         allocate(corrmatrix(ncols-2,ncols-2))
+     
+         corrmatrix = 0
+ 
+         do i=1, ncols-2
+          if (isused(i+2)) then
+                mean(i) = sum(coldata(1,0:nrows-1)*coldata(i+2,0:nrows-1))/numsamp
+          end if
+         end do
+ 
+          do i = 1, ncols-2
+           if (isused(i+2)) then
+              do j = i, ncols-2
+              if (isused(j+2)) then
+              corrmatrix(i,j) = sum(coldata(1,0:nrows-1)*((coldata(2+i,0:nrows-1)-mean(i))* &
+                                            (coldata(2+j,0:nrows-1)-mean(j))))/numsamp
+              corrmatrix(j,i) = corrmatrix(i,j)
+              end if
+              end do
+            end if
+          end do 
+ 
+              if (covmat_dimension /= 0) then
+               if (covmat_dimension > ncols -2) stop 'covmat_dimension larger than number of parameters'
+               write (*,*) 'Writing covariance matrix for ',covmat_dimension,' parameters'
+               allocate(covmatrix(covmat_dimension,covmat_dimension))
+               covmatrix=corrmatrix(1:covmat_dimension,1:covmat_dimension)
+ 
+               call WriteSqMatrix(trim(rootname) //'.covmat',covmatrix, covmat_dimension)
+               deallocate(covmatrix)
+              end if
+ 
+               do i=1, ncols-2
+                if (corrmatrix(i,i) > 0) then
+                scale = sqrt(corrmatrix(i,i))
+                corrmatrix(i,:) = corrmatrix(i,:)/scale
+                corrmatrix(:,i) = corrmatrix(:,i)/scale
+                end if
+               end do
+              call WriteSqMatrix(trim(rootname) //'.corr',corrmatrix, ncols -2)
+              
+           end subroutine GetCovMatrix
+ 
+ 
+           function MostCorrelated2D(i1,i2, direc)
+             integer, intent(in) :: i1,i2, direc
+             integer MostCorrelated2D
+             !Find which parameter is most correllated with the degeneracy in ix1, ix2
+             integer pars(2)
+             real mat2D(2,2), evals(2), u(2,2)
+             real corrs(2,ncols-2)
+           
+             if (direc /= 0 .and. direc /= -1) stop 'Invalid 3D color parameter'
+ 
+             pars(1)=i1
+             pars(2)=i2
+             mat2D = corrmatrix(pars,pars)
+             u = mat2D
+ 
+             call Diagonalize(u, evals, 2)
+             corrs = matmul(transpose(u), corrmatrix(pars,:))
+             corrs(:,pars) = 0
+           
+             MostCorrelated2D = MaxIndex(abs(corrs(2+direc,colix(1:num_vars)-2)), num_vars)
+             MostCorrelated2D = colix(MostCorrelated2D) -2
+          
+           end function MostCorrelated2D
+ 
+           subroutine GetFractionIndices(fraction_indices,n)
+              integer, intent(in) :: n
+              integer num, fraction_indices(*), i
+              real tot, aim
+ 
+              tot = 0
+              aim=numsamp/n
+              num = 1
+              fraction_indices(1) = 0
+              fraction_indices(n+1) = nrows
+              do i=0, nrows-1
+               tot = tot + coldata(1,i)
+               if (tot > aim) then
+                num=num+1
+                fraction_indices(num) = i
+                if (num==n) exit              
+                aim=aim+numsamp/n
+               end if   
+              end do
+ 
+           end subroutine GetFractionIndices
+ 
+            function ConfidVal(ix,limfrac,upper,ix1,ix2)
+              integer, intent(IN) :: ix
+              real, intent(IN) :: limfrac
+              logical, intent(IN) :: upper
+              integer, intent(IN), optional :: ix1,ix2
+              real ConfidVal, samps
+              real try, lasttry, try_t, try_b
+              integer l,t 
+           !find upper and lower bounds
+             
+              l=0
+              t=nrows-1
+              if (present(ix1)) l=ix1
+              if (present(ix2)) t = ix2
+              try_b = minval(coldata(ix,l:t))
+              try_t = maxval(coldata(ix,l:t))
+ 
+              samps = sum(coldata(1,l:t))
+ 
+              lasttry = -1
+              if (upper) then
+              do 
+                try = sum(coldata(1,l:t),mask = coldata(ix,l:t) > (try_b + try_t)/2) 
+                if (try > samps*limfrac) then
+                   try_b = (try_b+try_t)/2
+                else
+                   try_t = (try_b+try_t)/2
+                end if
+                if (try == lasttry) exit
+                lasttry = try
+              end do
+              else
+                do 
+                try = sum(coldata(1,l:t),mask = coldata(ix,l:t) < (try_b + try_t)/2) 
+                if (try > samps*limfrac) then
+                   try_t = (try_b+try_t)/2
+                else
+                   try_b = (try_b+try_t)/2
+                end if
+                if (try == lasttry) exit
+                lasttry = try
+                 end do
+              end if
+ 
+              ConfidVal = try_t
+             end function ConfidVal
+ 
+           subroutine PCA(pars,n,param_map,normparam_num)
+             !Perform principle component analysis
+             !In other words, get eigenvectors and eigenvalues for normalized variables
+             !with optional (log) mapping
+             integer, intent(in) :: n, pars(n), normparam_num
+             character(LEN=*), intent(IN) ::  param_map
+             integer i, j, normparam, locarr(1:1)
+             character(LEN =60) fmt, tmpS
+ 
+             real PCmean(n), sd(n), newmean(n), newsd(n)
+             real evals(n)
+             real corrmatrix(n,n), u(n,n)
+             real, dimension(:,:), allocatable :: PCdata
+             character(LEN=100) PClabs(n), div, expo
+             logical doexp
+ 
+            
+              write (*,*) 'Doing PCA for ',n,' parameters'
+              call CreateTxtFile(trim(rootname) //'.PCA', 40)
+              write (40,*) 'PCA for parameters: '
+ 
+              if (normparam_num /=0) then
+                 normparam = IndexOf(normparam_num,pars,n)
+                 if (normparam ==0) stop 'Invalid PCA normalization parameter'
+              else
+                 normparam = 0
+              end if
+ 
+              allocate(PCdata(n,0:nrows-1))
+              
+              PCdata(:,:) = coldata(2+pars,0:nrows-1)
+              fmt = trim(numcat('',n))//'f8.4)'
+ 
+              doexp =.false.
+              do i = 1, n
+                 if (param_map(i:i)=='L') then 
+                    doexp = .true.
+                    PCdata(i,:) = log(PCdata(i,:))
+                    PClabs(i) = 'ln(' //trim(labels(pars(i)+2))//')'
+                 elseif (param_map(i:i)=='M') then 
+                    doexp = .true.
+                    PCdata(i,:) = log(-PCdata(i,:))
+                    PClabs(i) = 'ln(-' //trim(labels(pars(i)+2))//')'
+                 else 
+                    PClabs(i) = trim(labels(pars(i)+2))
+                 end if
+                 write (40,*) pars(i),':',trim(PClabs(i))
+                 PCmean(i) = sum(coldata(1,0:nrows-1)*PCdata(i,:))/numsamp
+                 PCdata(i,:) = PCdata(i,:) - PCmean(i)
+                 sd(i) = sqrt(sum(coldata(1,0:nrows-1)*PCdata(i,:)**2)/numsamp)
+                 if (sd(i)/=0) PCdata(i,:) = PCdata(i,:)/sd(i)
+                 corrmatrix(i,i) = 1
+              end do
+ 
+              write (40,*) ''
+              write (40,*) 'Correlation matrix for reduced parameters'
+              do i = 1, n
+               do j = i, n
+               corrmatrix(i,j) = sum(coldata(1,0:nrows-1)*PCdata(i,:)*PCdata(j,:))/numsamp
+                   corrmatrix(j,i) = corrmatrix(i,j)
+           end do
+                 write (40,'(1I4,'': '','//trim(fmt)) pars(i), corrmatrix(i,:)
+          end do 
+             
+              u = corrmatrix
+              call Diagonalize(u, evals, n)
+ 
+              write (40,*) ''
+              write (40,*) 'e-values of correlation matrix'
+              do i = 1, n
+                 write (40,'(''PC'',1I2,'': '','//trim(fmt)) i, evals(i)
+              end do
+ 
+              write (40,*) ''
+              write (40,*) 'e-vectors'
+              
+              do i = 1, n
+                 write (40,'(1I3,'': '','//trim(fmt)) pars(i), u(i,:)
+              end do
+ 
+           
+               if (normparam /= 0) then
+                 !Set so parameter normparam has exponent 1
+                 do i=1, n
+                    u(:,i) = u(:,i)/u(normparam,i)*sd(normparam)
+                 end do
+               else
+                 !Normalize so main component has exponent 1
+                 do i=1, n
+                    locarr(1:1) = maxloc(abs(u(:,i)))
+                    u(:,i) = u(:,i)/u(locarr(1),i)*sd(locarr(1))
+                 end do
+               end if
+                      
+              do i = 0, nrows -1
+                PCdata(:,i) = matmul(transpose(u), PCdata(:,i))
+                if (doexp) PCdata(:,i) = exp(PCdata(:,i))
+              end do
+ 
+              write (40,*) ''
+              write (40,*) 'Principle components'
+             
+              do i = 1, n
+                 tmpS = trim(numcat('PC',i))//' (e-value: '//trim(RealToStr(evals(i)))//')'
+                 !ifc gives recursive IO error if you use a write within a write, 
+                 !even if separate or internal files
+                 write (40,*) trim(tmpS)
+                 do j=1,n
+                    if (param_map(j:j)=='L' .or. param_map(j:j)=='M' ) then
+                    expo = RealToStr(1/sd(j)*u(j,i) )
+                    if (param_map(j:j)=='M') then
+                       div = RealToStr( -exp(PCmean(j)))
+                     else
+                       div = RealToStr( exp(PCmean(j)))
+                     end if
+                     write(40,*) '['//trim(RealToStr(u(j,i)))//']   ('//trim(labels(pars(j)+2)) &
+                           //'/'//trim(div)//')^{'//trim(expo)//'}'
+                    else
+                       expo = RealToStr(sd(j)/u(j,i))
+                       if (doexp) then
+                        write(40,*) '['//trim(RealToStr(u(j,i)))//']    exp(('// &
+                         trim(labels(pars(j)+2))//'-'//trim(RealToStr(PCmean(j)))// ')/'// trim(expo)//')'
+                       else
+                         write(40,*) '['//trim(RealToStr(u(j,i)))//']   ('// &
+                        trim(labels(pars(j)+2))//'-'//trim(RealToStr(PCmean(j)))// ')/'// trim(expo)
+                       end if
+                    end if
+ 
+                 end do
+                 newmean(i) = sum(coldata(1,0:nrows-1)*PCdata(i,:))/numsamp
+                 newsd(i) = sqrt(sum(coldata(1,0:nrows-1)*(PCdata(i,:)-newmean(i))**2)/numsamp) 
+                 write (40,*) '          = '//trim(RealToStr(newmean(i)))// ' +- '//trim(RealToStr(newsd(i)))
+                 write (40,'('' ND limits: '',4f9.3)') minval(PCdata(i,0:ND_cont1)),maxval(PCdata(i,0:ND_cont1)), &
+                                              minval(PCdata(i,0:ND_cont2)),maxval(PCdata(i,0:ND_cont2))
+                 write (40,*) ''
+                 
+              end do
+ 
+ !Find out how correlated these components are with other parameters
+               write (40,*) 'Correlations of principle components'
+    
+               write (40,trim(numcat('(''    '',',n))//'I8)') (I, I=1, n)
+ 
+               do i=1, n
+                PCdata(i,:) = (PCdata(i,:) - newmean(i)) /newsd(i)
+               end do
+ 
+                do j=1,n
+                   write (40,'(''PC'',1I2)', advance = 'no') j
+                   do i=1,n
+                    write (40,'(1f8.3)', advance ='no') sum(coldata(1,0:nrows-1)*PCdata(i,:)* &
+                        PCdata(j,:))/numsamp
+                   end do
+                   write (40,*) ''
+                end do
+ 
+                do j=1, num_vars
+                    write (40,'(1I4)', advance = 'no') colix(j)-2
+                    do i=1,n
+                     write (40,'(1f8.3)', advance ='no') sum(coldata(1,0:nrows-1)*PCdata(i,:)* &
+                       (coldata(colix(j),0:nrows-1)-mean(j))/sddev(j))/numsamp
+                    end do
+                    write (40,*) '  ('//trim(labels(colix(j)))//')'
+                end do
+ 
+ 
+              close(40)
+              deallocate(PCdata)
+ 
+           end subroutine PCA
+ 
+          subroutine GetUsedCols
+           integer j
+                  
+            do j = 3, ncols
+             isused(j) = any(coldata(j,0:nrows-1)/=coldata(j,0)) 
+            end do  
+          end subroutine GetUsedCols
+ 
+          subroutine DoConvergeTests
+           real chain_means(max_chains,max_cols),chain_samp(max_chains)
+           real between_chain_var(max_cols), in_chain_var(max_cols)
+           integer frac(max_split_tests+1), split_n, chain_start(max_chains)
+           integer i,j,k,jj,kk,ix, maxoff
+           real split_tests(max_split_tests)
+           real mean(max_cols), fullmean(max_cols),fullvar(max_cols)
+           real, parameter :: cutfrac = 0.5
+           real usedsamps,evals(max_cols),sc,R, maxsamp
+           real, dimension(:,:), allocatable :: cov, meanscov
+           integer usedvars(max_cols), num, thin_fac(max_chains), markov_thin(max_chains)
+           real(gp) u, g2, fitted, focus, alpha, beta, probsum, tmp1
+           real confid
+           integer i1,i2,i3, nburn(max_chains), hardest, endb,hardestend
+           integer tran(2,2,2), tran2(2,2), off    
+           integer, dimension(:), allocatable :: binchain
+           real, parameter :: epsilon = 0.001
+           double precision, dimension(:,:), allocatable :: corrs
+           character(LEN=10) :: typestr
+           integer autocorr_thin
+ 
+ ! Get statistics for individual chains, and do split tests on the samples
+ 
+           call CreateTxtFile(trim(rootname) //'.converge', 40)
+ 
+          if (num_chains_used > 1) write (*,*) 'Number of chains used =  ',num_chains_used
+   
+           chain_indices(num_chains_used+1) = nrows
+           do i=1, num_chains_used
+                chain_start(i) =  chain_indices(i) + nint((chain_indices(i+1)-chain_indices(i))*cutfrac) 
+                chain_samp(i) = sum(coldata(1,chain_start(i):chain_indices(i+1)-1)) 
+           end do
+           usedsamps = sum(chain_samp(1:num_chains_used))
+           maxsamp = maxval(chain_samp(1:num_chains_used))
+           num = 0  
+           do j = 3, ncols
+            if (isused(j)) then
+             mean(j)=0
+             fullmean(j) =  sum(coldata(1,0:nrows-1)*coldata(j,0:nrows-1))/numsamp
+            
+             if (num_chains_used> 1) then
+              num=num+1
+              usedvars(num)=j
+              do i=1, num_chains_used
+                  mean(j) = mean(j) + sum(coldata(1,chain_start(i):chain_indices(i+1)-1)*&
+                            coldata(j,chain_start(i):chain_indices(i+1)-1))
+              end do
+              mean(j) = mean(j)/usedsamps
+             end if
+            end if
+           end do
+ 
+           do j = 3, ncols
+              if (isused(j)) &
+                   fullvar(j)=  sum(coldata(1,0:nrows-1)*(coldata(j,0:nrows-1)-fullmean(j))**2)/numsamp
+           end do
+ 
+ 
+           if (num_chains_used > 1) then
+           
+           write (40,*) ''
+           write(40,*)  'Variance test convergence stats using last half chains'
+           write (40,*) 'param var(chain mean)/mean(chain var)'
+           write (40,*) ''
+ 
+           do j = 3, ncols
+            between_chain_var(j) = 0
+            in_chain_var(j) = 0
+   
+            if (isused(j)) then
+ 
+            if (num_chains_used > 1) then
+ !Get stats for individual chains - the variance of the means over the mean of the variances
+            do i=1, num_chains_used
+              chain_means(i,j) =   sum(coldata(1,chain_start(i):chain_indices(i+1)-1)* &
+                    coldata(j,chain_start(i):chain_indices(i+1)-1))/chain_samp(i) 
+ 
+              between_chain_var(j) = between_chain_var(j) + &
+                ! chain_samp(i)/maxsamp* & !Weight for different length chains
+                   (chain_means(i,j) - mean(j))**2
+ 
+              in_chain_var(j) = in_chain_var(j) +  & !chain_samp(i)/maxsamp *&
+                sum(coldata(1,chain_start(i):chain_indices(i+1)-1)* &
+                (coldata(j,chain_start(i):chain_indices(i+1)-1)-chain_means(i,j))**2)
+    
+            end do
+            between_chain_var(j) = between_chain_var(j)/(num_chains_used-1) !(usedsamps/maxsamp -1)
+            in_chain_var(j) = in_chain_var(j)/usedsamps
+            write (40,'(1I3,f9.4,"  '//trim(labels(j))//'")') j-2, & 
+                 between_chain_var(j) /in_chain_var(j)
+            end if
+ 
+            end if
+            end do          
+     
+           end if
+      
+           if (num_chains_used > 1 .and. covmat_dimension>0) then
+            !Assess convergence in the var(mean)/mean(var) in the worst eigenvalue
+            !c.f. Brooks and Gelman 1997          
+        
+           do while (usedvars(num) > covmat_dimension+2) 
+              num= num - 1
+           end do
+        
+           allocate(meanscov(num,num))
+           allocate(cov(num,num))
+      
+           do jj=1,num
+            j=usedvars(jj)
+            do kk=jj, num
+             k = usedvars(kk)
+             meanscov(jj,kk) =0
+             cov(jj,kk)= 0
+              do i= 1, num_chains_used
+ 
+               cov(jj,kk) = cov(jj,kk) + &
+               !sqrt(chain_samp(jj)*chain_samp(kk))/maxsamp * &
+                  sum(coldata(1,chain_start(i):chain_indices(i+1)-1)* &
+                     (coldata(j,chain_start(i):chain_indices(i+1)-1)-chain_means(i,j))* &
+                     (coldata(k,chain_start(i):chain_indices(i+1)-1)-chain_means(i,k)))
+ 
+               meanscov(jj,kk) = meanscov(jj,kk)+ & !sqrt(chain_samp(jj)*chain_samp(kk))/maxsamp * &
+                          (chain_means(i,j)-mean(j))*(chain_means(i,k)-mean(k))
+              end do 
+              meanscov(kk,jj) = meanscov(jj,kk)
+              cov(kk,jj) = cov(jj,kk)              
+                   
+             end do
+           end do
+           meanscov = meanscov/(num_chains_used-1) !(usedsamps/maxsamp -1)
+           cov = cov / usedsamps
+ 
+           do jj=1,num
+                 sc = sqrt(cov(jj,jj))
+                 cov(jj,:) = cov(jj,:) / sc
+                 cov(:,jj) = cov(:,jj) / sc
+                 meanscov(jj,:) = meanscov(jj,:) /sc
+                 meanscov(:,jj) = meanscov(:,jj) /sc                 
+           end do
+  
+          call Diagonalize(meanscov, evals, num)
+          cov = matmul(matmul(transpose(meanscov),cov),meanscov)
+          write (40,*) ''
+          write (40,*) 'var(mean)/mean(var) for eigenvalues of covariance of means'
+          R = 0
+          do jj=1,num
+            write (40,'(1I3,f9.4)') jj,evals(jj)/cov(jj,jj)       
+            R = max(R,evals(jj)/cov(jj,jj))
+          end do
+ !R is essentially the Gelman and Rubin statistic
+          write (*,'(" var(mean)/mean(var), 1/2 chains, worst e-value: R-1 = ",f9.4)') R      
+         
+          deallocate(cov,meanscov)                 
+          end if
+      
+ 
+ 
+ !Do tests for robustness under using splits of the samples
+ !Return the rms ([change in upper/lower quantile]/[standard deviation]) 
+ !when data split into 2, 3,.. sets
+           write (40,*) ''
+           write (40,*)  'Split tests: rms_n([delta(upper/lower quantile)]/sd) n={2,3,4}:'
+           write(40,*) 'i.e. mean sample splitting change in the quantiles in units of the st. dev.'
+           write (40,*) '' 
+           do j = 3, ncols
+            if (isused(j)) then
+            do endb =0,1
+            do split_n = 2,max_split_tests
+             call GetFractionIndices(frac,split_n)
+              split_tests(split_n) = 0
+              confid =  ConfidVal(j,(1-contours(num_contours))/2,endb==0,0,nrows-1)
+             do i=1,split_n
+ 
+               split_tests(split_n) = split_tests(split_n) + &
+                 (ConfidVal(j,(1-contours(num_contours))/2,endb==0,frac(i),frac(i+1)-1)-confid)**2
+   
+             end do !i
+              split_tests(split_n) = sqrt(split_tests(split_n)/split_n/fullvar(j))
+            end do 
+             if (endb==0) then
+               typestr = 'upper'
+             else
+               typestr = 'lower'
+             end if       
+             write (40,'(1I3,'//trim(IntToStr(max_split_tests-1)) // 'f9.4,"  ' &
+                 //trim(labels(j))//' '//trim(typestr)//'")') j-2, split_tests(2:max_split_tests)
+ 
+            end do !endb
+            end if
+           end do
+ 
+ ! Now do Raftery and Lewis method
+ ! See http://www.stat.washington.edu/tech.reports/raftery-lewis2.ps
+ ! Raw non-importance sampled chains only
+          if (all(abs(coldata(1,0:nrows-1) - nint(coldata(1,0:nrows-1)))<1e-4)) then
+            nburn = 0
+ 
+            do ix=1, num_chains_used
+            thin_fac(ix) = nint(maxval(coldata(1,chain_indices(ix):chain_indices(ix+1)-1)))
+ 
+            do j = 3, covmat_dimension+2
+             if (isused(j) .and. (force_twotail .or. .not. has_limits(j))) then          
+               do endb =0,1
+              !Get binary chain depending on whether above or below confidence value
+               u = ConfidVal(j,(1-contours(num_contours))/2,endb==0,&
+                            chain_indices(ix),chain_indices(ix+1)-1)
+               do !thin_fac
+               call ThinData(thin_fac(ix),chain_indices(ix),chain_indices(ix+1)-1)
+               if (thin_rows < 2) exit  
+               allocate(binchain(0:thin_rows-1))
+               where (coldata(j,thin_ix(0:thin_rows-1)) >= u)
+                 binchain = 1
+               elsewhere
+                 binchain = 2
+               endwhere 
+               tran = 0
+               !Estimate transitions probabilities for 2nd order process
+               do i = 2, thin_rows-1
+                tran(binchain(i-2),binchain(i-1),binchain(i)) = &
+                 tran(binchain(i-2),binchain(i-1),binchain(i)) +1 
+               end do                 
+               deallocate(binchain)
+    
+               !Test whether 2nd order is better than Markov using BIC statistic
+               g2 = 0
+               do i1=1,2
+                do i2=1,2
+                 do i3=1,2
+                   if (tran(i1,i2,i3)/=0) then
+                    fitted = dble(tran(i1,i2,1) + tran(i1,i2,2)) * &
+                      (tran(1,i2,i3) + tran(2,i2,i3))  / dble( tran(1,i2,1) + &
+                      tran(1,i2,2) + tran(2,i2,1) + tran(2,i2,2) )
+                    focus = dble( tran(i1,i2,i3) )
+                    g2 = g2 + log( focus / fitted ) * focus
+ 
+                   end if 
+                 end do !i1
+                end do !i2
+               end do !i3
+               g2 = g2 * 2 
+               if (g2 - log( dble(thin_rows-2) ) * 2 < 0) exit
+               thin_fac(ix) = thin_fac(ix) + 1
+               end do !thin_fac
+        
+             !Get Markov transition probabilities for binary processes
+             if (sum(tran(:,1,2))==0 .or. sum(tran(:,2,1))==0) then
+               thin_fac(ix) = 0
+               goto 203
+             end if
+             alpha = sum(tran(:,1,2))/dble(sum(tran(:,1,1))+sum(tran(:,1,2)))
+             beta =  sum(tran(:,2,1))/dble(sum(tran(:,2,1))+sum(tran(:,2,2)))
+             probsum = alpha + beta
+             tmp1 = log(probsum * epsilon / max(alpha,beta))/ log( dabs(1.0d0 - probsum) )
+             if (int( tmp1 + 1 ) * thin_fac(ix) > nburn(ix)) then
+               nburn(ix) = int( tmp1 + 1 ) * thin_fac(ix)
+               hardest = j
+               hardestend = endb
+             end if
+             end do
+             end if
+             end do !j
+        
+              markov_thin(ix) = thin_fac(ix)
+ 
+             !Get thin factor to have independent samples rather than Markov
+               hardest = max(hardest,1)
+               u = ConfidVal(hardest,(1-contours(num_contours))/2,hardestend==0)
+               thin_fac(ix) = thin_fac(ix) + 1             
+               do !thin_fac
+               call ThinData(thin_fac(ix),chain_indices(ix),chain_indices(ix+1)-1)
+               if (thin_rows < 2) exit  
+               allocate(binchain(0:thin_rows-1))
+               where (coldata(hardest,thin_ix(0:thin_rows-1)) > u)
+                 binchain = 1
+               elsewhere
+                 binchain = 2
+               endwhere 
+               tran2 = 0
+               !Estimate transitions probabilities for 2nd order process
+               do i = 1, thin_rows-1
+                tran2(binchain(i-1),binchain(i)) = &
+                 tran2(binchain(i-1),binchain(i)) +1 
+               end do                 
+               deallocate(binchain)
+    
+               !Test whether independence is better than Markov using BIC statistic
+               g2 = 0
+               do i1=1,2
+                do i2=1,2
+                   if (tran2(i1,i2)/=0) then
+                     fitted = dble( (tran2(i1,1) + tran2(i1,2)) * (tran2(1,i2) + &
+                         tran2(2,i2)) ) / dble(thin_rows -1)
+                     focus = dble( tran2(i1,i2) )
+                     if (fitted <= 0 .or. focus <= 0) then
+                      write (*,*) 'Raftery and Lewis estimator had problems'
+                      return
+                     end if
+                     g2 = g2 + dlog( focus / fitted ) * focus
+                     end if 
+                 end do !i1
+               end do !i2
+               g2 = g2 * 2
+ 
+               if (g2 - log( dble(thin_rows-1) ) < 0) exit
+               thin_fac(ix) = thin_fac(ix) + 1
+               end do !thin_fac
+               
+ 203           if (thin_rows < 2) thin_fac(ix) = 0
+  
+              end do !chains
+  
+              write (40,*) ''
+              write (40,*) 'Raftery&Lewis statistics'
+              write (40,*) ''
+              write (40,*) 'chain  markov_thin  indep_thin    nburn'
+              do ix = 1, num_chains_used
+               if (thin_fac(ix)==0) then
+                  write (40,'(1I4,"      Not enough samples")') chain_numbers(ix)
+               else
+                  write (40,'(1I4,3I12)') chain_numbers(ix), markov_thin(ix), &
+                                          thin_fac(ix), nburn(ix)
+               end if
+              end do
+                                     
+             if (any(thin_fac(1:num_chains_used)==0)) then
+               write (*,*) 'RL: Not enough samples to estimate convergence stats'
+             else        
+               call writeS ('RL: Thin for Markov:         '//&
+                  Trim(IntToStr(maxval(markov_thin(1:num_chains_used)))))
+               indep_thin = maxval(thin_fac(1:num_chains_used))
+               call writeS('RL: Thin for indep samples:  '// &
+                trim(IntToStr(indep_thin)))            
+               call WriteS('RL: Estimated burn in steps: '//&
+                   trim(IntToStr(maxval(nburn(1:num_chains_used))))//' ('//&
+                   trim(IntToStr(nint(maxval(nburn(1:num_chains_used))/mean_mult)))//' rows)')
+ 
+             end if
+ 
+ !!Get correlation lengths
+            write (40,*) ''
+            write (40,*) 'Parameter auto-correlations as function of step separation'
+            write (40,*) ''
+            if (indep_thin ==0) then
+              autocorr_thin = 20
+             elseif (indep_thin <= 30) then
+              autocorr_thin =  5
+             else
+              autocorr_thin = 5* (indep_thin/30)
+            end if  
+ 
+             call ThinData(autocorr_thin,0,nrows-1)
+             maxoff = min(15,thin_rows/(autocorr_thin*num_chains_used))
+             allocate(Corrs(ncols,maxoff))
+             corrs = 0
+             do off =1,maxoff                 
+                 do i=off, thin_rows-1
+                  do j = 3, ncols
+                   if (isused(j)) &
+                    corrs(j,off) = corrs(j,off) + (coldata(j,thin_ix(i))-fullmean(j))* &
+                          (coldata(j,thin_ix(i-off)) - fullmean(j))
+                  end do
+                 end do
+                 do j = 3, ncols
+                   if (isused(j)) &
+                      corrs(j,off) = corrs(j,off)/(thin_rows-off)/fullvar(j)
+                 end do
+   
+              end do
+           
+              write (40,'("   ",'//trim(IntToStr(maxoff)) // 'I8)')  &
+                   (/(I, I=autocorr_thin, maxoff*autocorr_thin,autocorr_thin)/)     
+             
+              do j = 3, ncols
+               if (isused(j)) then                        
+                write (40,'(1I3,'//trim(IntToStr(maxoff)) // 'f8.3,"  ' &
+                 //trim(labels(j))//'")') j-2, corrs(j,:)
+               end if
+              end do
+         
+             deallocate(Corrs)
+ 
+ 
+ 
+          end if
+      
+          close(40)
+ 
+          end subroutine DoConvergeTests
+          
+ 
+ 
+         subroutine Get2DPlotData(j,j2)
+          integer, intent(in) :: j,j2
+          integer i,ix1,ix2,wx,wy
+          real binweight, dist, norm, maxbin
+          real try_b, try_t,try_sum, try_last
+          character(LEN=120) :: plotfile, filename,numstr
+       
+               bins2D = 0
+               bin2Dlikes = 0
+               bin2Dmax = 1e30        
+ 
+               do i = 0, nrows-1
+                 ix1=nint((coldata(colix(j),i)-center(j))/width(j))
+                 ix2=nint((coldata(colix(j2),i)-center(j2))/width(j2))
+        
+                 if (smoothing) then
+                 do wx = max(ix_min(j),ix1-2),min(ix_max(j),ix1+2)
+                    dist = (coldata(colix(j),i) - (center(j) + wx*width(j)))**2/width(j)**2
+                  do wy = max(ix_min(j2),ix2-2),min(ix_max(j2),ix2+2)
+                    binweight = coldata(1,i)*exp( - (&
+                       dist +(coldata(colix(j2),i) - (center(j2) + wy*width(j2)))**2/width(j2)**2 )/2)
+                    bins2D(wx,wy) = bins2D(wx,wy) + binweight
+                    bin2Dlikes(wx,wy) = bin2Dlikes(wx,wy) + binweight*exp(meanlike - coldata(2,i))
+                  end do
+                 end do
+                 else
+                   bins2D(ix1,ix2) = bins2D(ix1,ix2) + coldata(1,i) 
+                   bin2Dlikes(ix1,ix2) = bin2Dlikes(ix1,ix2) + coldata(1,i)*exp(meanlike-coldata(2,i))
+                 end if
+               bin2Dmax(ix1,ix2) = min(bin2Dmax(ix1,ix2),real(coldata(2,i)))
+               end do
+ 
+              
+               do ix1=ix_min(j),ix_max(j)
+                do ix2 =ix_min(j2),ix_max(j2) 
+                   if (bins2D(ix1,ix2) >0) bin2Dlikes(ix1,ix2) = bin2Dlikes(ix1,ix2)/bins2D(ix1,ix2)
+                   if (plot_NDcontours) then
+                     !Map maximum likelihood in each bin into a significance value from the full N-D distribution
+                     if (bin2DMax(ix1,ix2) == 1e30) then
+                      bin2DMax(ix1,ix2) = 0
+                       else
+                       bin2DMax(ix1,ix2) = sum(coldata(1,0:nrows-1), mask = coldata(2,0:nrows-1) > bin2DMax(ix1,ix2))/numsamp
+                     end if
+                  end if
+                 end do 
+               end do
+      
+ 
+               if (has_limits(colix(j)) .or. has_limits(colix(j2))) then
+ 
+                !Fix up underweighting near edges. Note this makes edge pixels noisier.
+                do ix1 = ix_min(j), ix_max(j)
+                  do ix2 = ix_min(j2),ix_max(j2)
+                     if (ix1 ==ix_min(j) .and. has_limits_bot(colix(j))) bins2D(ix1,ix2) = bins2D(ix1,ix2)*2
+                     if (ix1 ==ix_max(j) .and. has_limits_top(colix(j))) bins2D(ix1,ix2) = bins2D(ix1,ix2)*2
+                     if (ix2 ==ix_min(j2).and. has_limits_bot(colix(j2))) bins2D(ix1,ix2) = bins2D(ix1,ix2)*2
+                     if (ix2 ==ix_max(j2).and. has_limits_top(colix(j2))) bins2D(ix1,ix2) = bins2D(ix1,ix2)*2
+ 
+                     if (ix1 ==ix_min(j)+1.and. has_limits_bot(colix(j))) bins2D(ix1,ix2) = bins2D(ix1,ix2)/0.84
+                     if (ix1 ==ix_max(j)-1.and. has_limits_top(colix(j))) bins2D(ix1,ix2) = bins2D(ix1,ix2)/0.84
+                     if (ix2 ==ix_min(j2)+1.and. has_limits_bot(colix(j2))) bins2D(ix1,ix2) = bins2D(ix1,ix2)/0.84
+                     if (ix2 ==ix_max(j2)-1.and. has_limits_top(colix(j2))) bins2D(ix1,ix2) = bins2D(ix1,ix2)/0.84
+                  end do
+                end do
+ 
+              end if
+ 
+ ! Get contour containing contours(:) of the probability
+ 
+              allocate(TheBins(ix_max(j)-ix_min(j)+1,ix_max(j2)-ix_min(j2)+1))
+              TheBins = bins2D(ix_min(j):ix_max(j),ix_min(j2):ix_max(j2))
+                 
+              norm = sum(TheBins)
+              
+              do ix1 = 1, num_contours
+              try_t = maxval(TheBins)
+              try_b = 0
+              try_last = -1
+              do
+                try_sum = sum(TheBins,mask = TheBins < (try_b + try_t)/2) 
+                if (try_sum > (1-contours(ix1))*norm) then
+                   try_t = (try_b+try_t)/2
+                else
+                   try_b = (try_b+try_t)/2
+                end if
+                if (try_sum == try_last) exit
+                try_last = try_sum
+          
+              end do
+              contour_levels(ix1) = (try_b+try_t)/2
+ 
+              end do
+              deallocate(TheBins)
+ 
+ 
+               plotfile = numcat(trim(numcat(trim(rootname)//'_2D_',colix(j)-2))//'_',colix(j2)-2)
+               filename = 'plot_data/'//trim(plotfile)
+               open(unit=49,file=filename,form='formatted',status='replace')
+               do ix1 = ix_min(j), ix_max(j)
+                   write (49,trim(numcat('(',ix_max(j2)-ix_min(j2)+1))//'E16.7)') bins2D(ix1,ix_min(j2):ix_max(j2))
+               end do
+               close(49)
+ 
+             open(unit=49,file=trim(filename)//'_cont',form='formatted',status='replace')
+             write(numstr,*) contour_levels(1:num_contours)
+             if (num_contours==1) numstr = trim(numstr)//' '//trim(numstr)
+             write (49,*) trim(numstr)
+             close(49)
+ 
+           if (shade_meanlikes) then
+                   open(unit=49,file=trim(filename)//'_likes',form='formatted',status='replace')
+                   maxbin = maxval(bin2Dlikes(ix_min(j):ix_max(j),ix_min(j2):ix_max(j2)))
+                   do ix1 = ix_min(j), ix_max(j)
+                     write (49,trim(numcat('(',ix_max(j2)-ix_min(j2)+1))//'E16.7)') &
+                            bin2Dlikes(ix1,ix_min(j2):ix_max(j2))/maxbin
+                   end do
+                  close(49)
+           end if
+ 
+           if (plot_NDcontours) then
+                   open(unit=49,file=trim(filename)//'_confid',form='formatted',status='replace')
+                   do ix1 = ix_min(j), ix_max(j)
+                     write (49,trim(numcat('(',ix_max(j2)-ix_min(j2)+1))//'E16.7)') &
+                            bin2Dmax(ix1,ix_min(j2):ix_max(j2))
+ 
+                   end do
+                  close(49)
+           end if
+      
+         end subroutine Get2DPlotData
+ 
+         function PlotContMATLAB(aunit,aroot,j,j2, DoShade)
+          character(LEN=*), intent(in) :: aroot
+          integer, intent(in) :: j, j2,aunit
+          logical, intent(in) :: DoShade
+          character(LEN=120) fname,numstr,plotfile
+          logical PlotContMATLAB
+ 
+              PlotContMATLAB= .false. 
+              plotfile = numcat(trim(numcat(trim(aroot)//'_2D_',colix(j)-2))//'_',colix(j2)-2)
+              if (.not. FileExists('plot_data/'//trim(plotfile))) return
+              PlotContMATLAB= .true.
+               write (aunit,'(a)') 'pts=load (fullfile(''plot_data'',''' // trim(plotfile) //'''));'
+      
+               fname =trim(trim(aroot)//trim(numcat('_p',colix(j2)-2)))
+               write (aunit,'(a)') 'tmp = load (fullfile(''plot_data'',''' // trim(fname)//  '.dat''));'
+               write (aunit,*) 'x1 = tmp(:,1);'
+ 
+               fname =trim(trim(aroot)//trim(numcat('_p',colix(j)-2)))
+               write (aunit,'(a)') 'tmp = load (fullfile(''plot_data'',''' // trim(fname)//  '.dat''));'
+               write (aunit,*) 'x2 = tmp(:,1);'
+     
+ !              fmt = trim(numcat('(''x1 = ['',',ix_max(j2)-ix_min(j2)+1))//'E16.7,''];'')'
+ !              write (aunit,fmt) (/(I, I=ix_min(j2), ix_max(j2), 1)/)*width(j2)+center(j2)
+ 
+  !             fmt = trim(numcat('(''x2 = ['',',ix_max(j)-ix_min(j)+1))//'E16.7,''];'')'
+  !             write (aunit,fmt) (/(I, I=ix_min(j), ix_max(j), 1)/)*width(j)+center(j)
+               if (DoShade) then
+               if (shade_meanlikes) then
+                 write (aunit,'(a)') 'load (fullfile(''plot_data'',''' // trim(plotfile) //'_likes''));'
+                 write (aunit,*) 'contourf(x1,x2,'//trim(plotfile)//'_likes,64);'
+               else     
+                 write (aunit,*) 'contourf(x1,x2,pts,64);'
+               end if
+ 
+               write (aunit,*) 'set(gca,''climmode'',''manual''); shading flat; hold on;'
+               end if
+               
+            
+               if (num_contours /= 0) then
+ 
+                write (aunit,'(a)') 'cnt = load (fullfile(''plot_data'',''' // trim(plotfile) //'_cont''));'
+   
+                if (plot_NDcontours) then
+                 write (aunit,'(a)') 'load (fullfile(''plot_data'',''' // trim(plotfile) //'_confid''));'
+                 write(numstr,*) 1-contours(1:num_contours)  
+                 if (num_contours==1) numstr = trim(numstr)//' '//trim(numstr)
+                write (aunit,'(a)') 'contour(x1,x2,'//trim(plotfile)//'_confid,['//trim(numstr)//'],''r'');'
+                
+                end if
+               end if
+ 
+                       
+         end function PlotContMATLAB
+ 
+ 
+         subroutine Write2DPlotMATLAB(aunit,j,j2, DoLabelx,DoLabely)
+           integer, intent(in) :: aunit,j,j2
+           logical, intent(in) :: DoLabelx,DoLabely
+           character(LEN=120) fmt
+           integer i
+           
+ 
+               if (PlotContMATLAB(aunit,rootname,j,j2,do_shading) .and. &
+                   num_contours /= 0) then
+ 
+                write (aunit,'(a)') '[C h] = contour(x1,x2,pts,cnt,lineM{1});'
+                write (aunit,'(a)') 'set(h,''LineWidth'',lw1);'
+                write (aunit,*) 'hold on; axis manual; '
+        
+                do i = 1, Num_ComparePlots
+                 fmt = trim(numcat('lineM{',i+1)) // '}'
+                 if (PlotContMATLAB(aunit,ComparePlots(i),j,j2,.false.)) &
+                  write (aunit,'(a)') '[C h] = contour(x1,x2,pts,cnt,'//trim(fmt)//');'
+                  write (aunit,'(a)') 'set(h,''LineWidth'',lw2);'
+              
+                end do
+ 
+               end if
+ 
+               write (aunit,*) 'hold off; set(gca,''Layer'',''top'',''FontSize'',axes_fontsize);'
+               fmt = ''',''FontSize'',lab_fontsize);'
+               if (DoLabelx) write (aunit,*) 'xlabel('''//trim(labels(colix(j2)))//trim(fmt)
+               if (DoLabely) write (aunit,*) 'ylabel(''',trim(labels(colix(j)))//trim(fmt)
+  
+          end subroutine  Write2DPlotMATLAB
+ 
+ 
+          subroutine WriteMatLabInit(unit,sm)
+           integer, intent(in) :: unit
+           logical, intent(in) :: sm
+           if (sm) then
+            write(unit,*) 'lab_fontsize = 9; axes_fontsize = 9;'
+           else
+            write(unit,*) 'lab_fontsize = 12; axes_fontsize = 12;'
+           end if
+           if (BW) then
+            if (plot_meanlikes) then
+             write(unit,*) 'lineM = {''-k'',''-r'',''-b'',''-m'',''-g'',''-y''};';
+             write(unit,*) 'lineL = {'':k'','':r'','':b'','':m'','':g'','':y''};';
+             write(unit,*) 'lw1=3;lw2=1;' !Line Widths     
+            else
+             write(unit,*) 'lineM = {''-k'',''--r'',''-.b'','':m'',''--g'',''-.y''};';
+             write(unit,*) 'lw1=1;lw2=1;' !Line Widths     
+            end if
+            else
+             write(unit,*) 'lineM = {''-k'',''-r'',''-b'',''-m'',''-g'',''-y''};';
+             write(unit,*) 'lineL = {'':k'','':r'','':b'','':m'','':g'','':y''};';
+            write(unit,*) 'lw1=1;lw2=1;' !Line Widths     
+           end if
+ 
+          end subroutine WriteMatLabInit
+ 
+      subroutine Write1DplotMatLab(aunit,j)
+       integer, intent(in) :: aunit, j
+       character(LEN=120) fname
+        character(LEN=60) fmt
+        integer ix1
+ 
+         fname =trim(trim(rootname)//trim(numcat('_p',colix(j)-2)))
+         write (aunit,*) 'load (fullfile(''plot_data'',''' // trim(fname)//  '.dat''));'
+         write (aunit,*) 'pts='//trim(fname)//';'
+         write (aunit,*) 'plot(pts(:,1),pts(:,2),lineM{1},''LineWidth'',lw1);'
+         write (aunit,*) 'axis([-Inf,Inf,0,1.1]);axis manual;'
+         write (aunit,*) 'set(gca,''FontSize'',axes_fontsize); hold on;'
+ 
+         if (plot_meanlikes) then
+          write (aunit,*) 'load (fullfile(''plot_data'',''' // trim(fname)//  '.likes''));'
+          write (aunit,*) 'pts='//trim(fname)//';'
+          write (aunit,*) 'plot(pts(:,1),pts(:,2),lineL{1},''LineWidth'',lw1);'
+         end if
+ 
+ 
+        do ix1 = 1, Num_ComparePlots
+          fname =  trim(trim(ComparePlots(ix1))//trim(numcat('_p',colix(j)-2)))    
+          if (FileExists('plot_data/' // trim(fname)//'.dat')) then
+ 
+           write (aunit,*) 'pts = load (fullfile(''plot_data'',''' // trim(fname)// '.dat''));'
+           fmt = trim(numcat('lineM{',ix1+1)) // '}'        
+           write (aunit,*) 'plot(pts(:,1),pts(:,2),'//trim(fmt)//',''LineWidth'',lw2);'
+     
+          if (plot_meanlikes) then
+ 
+             write (aunit,*) 'pts=load (fullfile(''plot_data'',''' //trim(fname)//'.likes''));'
+            fmt = trim(numcat('lineL{',ix1+1)) // '}'        
+             write (aunit,*) 'plot(pts(:,1),pts(:,2),'//trim(fmt)//',''LineWidth'',lw2);'
+ 
+          end if
+ 
+          end if
+        end do
+ 
+      end subroutine Write1DplotMatLab
+  
+ end module MCSamples
+ 
+ subroutine CheckMatlabAxes(afile)
+  integer, intent(in) :: afile
+ 
+   write (afile,*) 'ls =get(gca,''XTick'');sz=size(ls,2);'
+   write (afile,*) 'if(sz > 4)'
+   write (afile,*) ' set(gca,''XTick'',ls(:,1:2:sz));'
+   write (afile,*) 'end;'
+ 
+ end subroutine CheckMatlabAxes
+ 
+ 
+ 
+ program GetDist
+ 
+         use IniFile
+         use MCSamples
+         implicit none
+ 
+         real  bincounts(-1000:1000),binlikes(-1000:1000),binmaxlikes(-1000:1000)
+     
+         character(LEN=80) InputFile, numstr
+         character(LEN=300) fname
+         character(LEN=10000) InLine  ! ,LastL,OutLine
+         real invars(max_cols)
+         integer ix, ix1,i,ix2,ix3, nbins, wx
+         real ignorerows
+    
+         logical bad, adjust_priors
+ 
+         real  amax, amin
+         character(LEN=120) filename, infile, in_root,out_root,sm_file
+         character(LEN=120) matlab_col, InS1,InS2,fmt, tmpstr
+         integer plot_row, plot_col, chain_num, first_chain,chain_ix, plot_num
+        
+         integer x,y,j,j2, num_2D_plots, num_cust2D_plots
+         integer num_3D_plots
+         character(LEN=80) contours_str, plot_3D(20), bin_limits
+         integer plot_x, plot_y, thin_factor, outliers
+         real limmin(max_cols),limmax(max_cols), maxbin
+         integer plot_2D_param, j2min
+         real try_b, try_t, binweight
+         real cont_lines(max_cols,2,2), dist, distweight, limfrac 
+         integer bestfit_ix
+         integer chain_exclude(max_chains), num_exclude
+         logical map_params
+         logical :: triangle_plot = .false.
+         real counts
+         real cool
+         integer PCA_num, PCA_normparam
+         integer PCA_params(max_cols)
+         integer plotparams_num, plotparams(max_cols)
+         character(LEN=max_cols) PCA_func
+         logical Done2D(max_cols,max_cols)
+         logical samples_are_chains, no_plots
+ 
+         real thin_cool
+         logical no_tests, auto_label
+ 
+ 
+         InputFile = GetParam(1)
+         if (InputFile == '') stop 'No parameter input file'
+ 
+         call Ini_Open(InputFile, 1, bad, .false.)
+         if (bad) stop 'Error opening parameter file'
+         Ini_fail_on_not_found = .true.
+ 
+         in_root = Ini_Read_String('file_root')
+         rootname = ExtractFileName(in_root)
+         chain_num = Ini_Read_Int('chain_num')
+ 
+         ncols = Ini_Read_Int('columnnum',0)
+ 
+     
+         allocate(coldata(ncols,0:max_rows))
+ 
+         nbins = Ini_Read_Int('num_bins')
+ 
+         ignorerows = Ini_Read_Double('ignore_rows',0.d0)
+ 
+         adjust_priors = Ini_Read_Logical('adjust_priors',.false.)
+     
+         Ini_fail_on_not_found = .false.
+ 
+         matlab_version = Ini_Read_Real('matlab_version',6.)
+ 
+         labels(1) = 'mult'
+         labels(2) = 'likelihood'
+         auto_label = Ini_Read_Logical('auto_label',.false.)
+         do ix=3, ncols
+            write (numstr,*) ix-2
+            if (auto_label) then
+              labels(ix) = trim(adjustl(numstr))
+            else
+              labels(ix) = Ini_Read_String('lab'//trim(adjustl(numstr)))
+            end if 
+         end do
+ 
+       
+   
+         prob_label = Ini_Read_Logical('prob_label',.false.)
+         triangle_plot = Ini_Read_Logical('triangle_plot',.false.)
+         
+         samples_are_chains = Ini_Read_Logical('samples_are_chains',.true.)
+ 
+         no_tests = Ini_Read_Logical('no_tests',.false.)
+         no_plots = Ini_Read_Logical('no_plots',.false.)
+ 
+         thin_factor = Ini_Read_Int('thin_factor',0)
+         thin_cool = Ini_read_Real('thin_cool',1.)
+ 
+         first_chain = Ini_Read_Int('first_chain',1)
+ 
+         make_single_samples = Ini_Read_logical('make_single_samples',.false.)
+         single_thin = Ini_Read_Int('single_thin',1)
+         cool = Ini_Read_Real('cool',1.)
+ 
+         Num_ComparePlots = Ini_Read_Int('compare_num',0)
+         do ix = 1, Num_ComparePlots
+            ComparePlots(ix) = ExtractFileName(Ini_Read_String(numcat('compare',ix)))
+         end do
+ 
+         has_limits_top = .false.
+         has_limits_bot = .false.
+         
+         bin_limits = Ini_Read_String('all_limits')
+         do ix=3, ncols
+            write (numstr,*) ix-2
+            if (bin_limits /= '') then
+             InLine = bin_limits
+            else
+             InLine = Ini_Read_String('limits'//trim(adjustl(numstr)))
+            end if
+            if (InLine /= '') then
+               read (InLine,*) InS1, InS2 
+               if (trim(adjustl(InS1)) /= 'N') then
+                has_limits_bot(ix) = .true.
+                read(InS1,*) limmin(ix)
+               end if
+               if (trim(adjustl(InS2)) /= 'N') then
+                 has_limits_top(ix) = .true.
+                 read(InS2,*) limmax(ix)
+               end if
+            end if
+         end do
+ 
+         has_limits = has_limits_top .or. has_limits_bot
+ 
+         plot_2D_param = Ini_Read_Int('plot_2D_param',0)
+         if (plot_2D_param /= 0) then
+         plot_2D_param = plot_2D_param + 2
+         num_cust2D_plots = 0
+     else
+          !Use custom array of specific plots
+      num_cust2D_plots = Ini_Read_Int('plot_2D_num',0)
+      do ix = 1, num_cust2D_plots
+         InLine = Ini_Read_String(numcat('plot',ix))
+         read (InLine, *) ix1,ix2
+         cust2DPLots(ix) = ix1+2 + (ix2+2)*100 
+      end do
+     end if
+ 
+     InS1 =Ini_Read_String('exclude_chain')
+     num_exclude = 0
+     chain_exclude = 0
+     read (InS1, *, end =20) chain_exclude
+    
+ 20  do while (chain_exclude(num_exclude+1)/=0) 
+      num_exclude = num_exclude + 1
+     end do 
+     
+ 
+     map_params = Ini_Read_Logical('map_params',.false.)
+     if (map_params)  &
+           write (*,*) 'WARNING: Mapping params - .covmat file is new params.'
+ 
+     shade_meanlikes = Ini_Read_Logical('shade_meanlikes',.false.)   
+ 
+     matlab_col = Ini_Read_String('matlab_colscheme')
+     smoothing = Ini_Read_Logical('smoothing',.true.)
+ 
+     out_root = Ini_Read_String('out_root')
+     if (out_root /= '') then
+      rootname = out_root
+      write (*,*) 'producing files with root '//trim(out_root)
+     end if
+     sm_file = trim(rootname) // '.sm'
+ 
+ 
+     contours(1) = Ini_Read_Real('contour1',0.)
+     contours(2) = Ini_Read_Real('contour2',0.)
+     contours_str = trim(Ini_Read_String('contour1')) // ' ; '// &
+              trim(Ini_Read_String('contour2'))
+     num_contours = count(contours/= 0)
+     if (contours(1)==0) contours(1) = contours(2)
+ 
+     force_twotail = Ini_Read_Logical('force_twotail',.false.)  
+     if (force_twotail) write (*,*) 'Computing two tail limits'
+ 
+     covmat_dimension = Ini_Read_Int('cov_matrix_dimension',0)
+     if (covmat_dimension == -1) covmat_dimension = ncols-2
+     
+ 
+     plot_meanlikes = Ini_Read_Logical('plot_meanlikes',.false.)
+     plot_NDcontours = Ini_Read_Logical('plot_NDcontours',.false.)
+ 
+         PCA_num = Ini_Read_Int('PCA_num',0)
+         if (PCA_num /= 0) then
+            if (PCA_num <2) stop 'Can only do PCA for 2 or more parameters'
+            InLine = Ini_Read_String('PCA_params')
+            PCA_func =Ini_Read_String('PCA_func')
+             !Characters representing functional mapping
+            if (PCA_func == '') PCA_func(1:PCA_num) = 'N' !no mapping
+            if (InLine == 'ALL' .or. InLine == 'all') then
+              PCA_params(1:PCA_num) = (/ (i, i=1,PCA_num)/)
+            else
+              read (InLine,*) PCA_params(1:PCA_num)
+            end if
+            PCA_NormParam = Ini_Read_Int('PCA_normparam',0)
+         end if
+ 
+         plotparams_num = Ini_Read_Int('plotparams_num',0)
+         if (plotparams_num /= 0) then
+           InLine = Ini_Read_String('plotparams')
+           read(InLine,*) plotparams(1:plotparams_num)
+         end if
+ 
+         num_3D_plots = Ini_Read_Int('num_3D_plots',0)
+         do ix =1, num_3D_plots
+            plot_3D(ix) = Ini_Read_String(numcat('3D_plot',ix))
+         end do
+ 
+         BW = Ini_Read_Logical('B&W',.false.)
+         do_shading = Ini_Read_Logical('do_shading',.true.)
+         call Ini_Close
+ 
+         allocate(bins2D(-1000:1000,-1000:1000),bin2Dlikes(-1000:1000,-1000:1000),bin2Dmax(-1000:1000,-1000:1000))
+  
+ 
+ !Read in the chains     
+         nrows = 0
+         num_chains_used =0
+ 
+         do chain_ix = first_chain, first_chain-1 + max(1,chain_num)
+ 
+          do ix = 1, num_exclude
+            if (chain_exclude(ix) == chain_ix) goto 30
+          end do  
+       
+          if (chain_num == 0) then
+             infile = trim(in_root) // '.txt'
+          else
+             write (numstr,*) chain_ix
+             infile = trim(in_root) //'_'//trim(adjustl(numstr))// '.txt'
+          end if
+ 
+          write (*,*) 'reading ' // trim(infile)
+          open(unit=50,file=infile,form='formatted',status='old', err=28)
+ 
+           if (ignorerows >=1) then
+            do ix = 1, nint(ignorerows)
+              read (50,'(a)',end=1) InLine
+            end do
+           end if
+ 
+          num_chains_used = num_chains_used + 1
+          if (num_chains_used > max_chains) stop 'Increase max_chains in GetDist'
+          chain_indices(num_chains_used) = nrows
+          chain_numbers(num_chains_used) = chain_ix
+           do
+            read (50,'(a)',end=1) InLine
+            if (samples_are_chains) then
+             read(InLine,*,end=1, err=2) invars(1:ncols) 
+            else
+             read(InLine,*,end=1, err=2) invars(3:ncols) 
+             invars(1) = 1
+             invars(2) = 1
+            end if
+ 
+            if (map_params) call MapParameters(invars)
+            coldata(1:ncols, nrows) = invars(1:ncols)
+            
+            nrows = nrows + 1
+            if (nrows > max_rows) stop 'need to increase max_rows'
+ 
+           end do    
+  2       write (*,*) 'error reading line ', nrows + int(ignorerows), ' - skipping rest of file'
+      
+  1       close(50)
+          if (ignorerows<1 .and. ignorerows/=0) then
+            i = chain_indices(num_chains_used)
+            j = nint((nrows-i-1)*ignorerows)
+            do ix = i,nrows-j-1
+             coldata(:,ix) = coldata(:,ix+j)
+            end do
+            nrows = nrows - j
+          end if
+          goto 30
+  28      write (*,'(" chain ",1I4," missing")') chain_ix
+  30      end do
+ 
+           if (nrows == 0) stop 'No un-ignored rows! (check number of chains/burn in)'
+ 
+ 
+           if (cool /= 1) call CoolChain(cool)
+ !Adjust weights if requested
+           if (adjust_priors) then
+              call AdjustPriors
+           end if
+ 
+           call GetUsedCols !See which parameters are fixed
+ 
+           mean_mult = sum(coldata(1,0:nrows-1))/nrows
+   
+           max_mult = (mean_mult*nrows)/min(nrows/2,500)
+           outliers = count(coldata(1,0:nrows-1) > max_mult)
+           if (outliers /=0) write (*,*) 'outlier fraction ', real(outliers)/nrows
+ 
+           maxmult = maxval(coldata(1,0:nrows-1))
+           numsamp = sum(coldata(1,0:nrows-1))
+          
+          if (.not. no_tests) call DoConvergeTests
+          if (adjust_priors) call DeleteZeros  
+ 
+          write (*,*) 'mean input multiplicity = ',mean_mult
+ 
+ !Output thinned data if requested
+ !Must do this with unsorted output  
+           if (thin_factor /= 0) then 
+              call ThinData(thin_factor)
+              call WriteThinData(trim(rootname)//'_thin.txt',thin_cool)
+           end if
+ 
+ !Produce file of weight-1 samples if requested
+ 
+        if (num_3D_plots/=0 .and. .not. make_single_samples ) then
+             make_single_samples = .true.
+             single_thin = max(1,nint(numsamp/maxmult)/2000)
+         end if
+  
+       if (make_single_samples) call MakeSingleSamples(single_thin)
+ 
+ !Sort data in order of likelihood of points
+ 
+           call SortColData(2)
+       
+           numsamp = sum(coldata(1,0:nrows-1))
+ 
+ 
+ !Get ND confidence region (index into sorted coldata)
+           counts = 0
+       ND_cont1=-1; ND_cont2=-1
+       do j=0, nrows -1
+         counts = counts + coldata(1,j)
+         if (counts > numsamp*contours(1) .and. ND_cont1==-1) then
+          ND_cont1 = j
+              if (ND_cont2 /= -1) exit
+         end if
+         if (counts > numsamp*contours(2) .and. ND_cont2==-1) then
+          ND_cont2 = j
+         end if
+       end do          
+ 
+ !Only use variables whose label's are not empty (and in list of plotparams if plotparams_num /= 0)
+           num_vars = 0
+           
+           do ix = 3,ncols
+ 
+            if (labels(ix) /= '' .and. isused(ix)) then
+                 if (plotparams_num == 0 .or. any(plotparams(1:plotparams_num)==ix-2)) then
+                  num_vars = num_vars + 1
+                  colix(num_vars) = ix
+                  mean(num_vars) = sum(coldata(1,0:nrows-1)*coldata(ix,0:nrows-1))/numsamp
+                  sddev(num_vars)  = sqrt(sum(coldata(1,0:nrows-1)*(coldata(ix,0:nrows-1) &
+                            -mean(num_vars))**2)/numsamp)
+                 end if
+              end if
+           end do
+ 
+           triangle_plot = triangle_plot .and. (num_vars > 1)
+ 
+           write (*,*) 'using ',nrows,' rows, processing ',num_vars,' parameters'
+           if (indep_thin/=0) then
+            write (*,*) 'Approx indep samples: ', nint(numsamp/indep_thin)
+           else
+            write (*,*) 'effective number of samples (assuming indep): ', nint(numsamp/maxmult)
+           end if
+ !Get covariance matrix and correlation matrix
+ 
+           call GetCovMatrix
+ 
+           if (PCA_num>0) call PCA(PCA_params,PCA_num,PCA_func, PCA_NormParam)
+ 
+ !Find best fit, and mean likelihood     
+       bestfit_ix = 0 !Since we have sorted the lines
+       maxlike = coldata(2,bestfit_ix)
+       write (*,*) 'Best fit -Ln(like) = ', maxlike
+ 
+       if (coldata(2,nrows-1) - maxlike < 30) then
+        meanlike = log(sum(exp((coldata(2,0:nrows-1) -maxlike))*coldata(1,0:nrows-1)) & 
+                   / numsamp) + maxlike
+        write (*,*) 'Ln(mean 1/like) = ', meanlike
+       end if
+ 
+       meanlike = sum(coldata(2,0:nrows-1)*coldata(1,0:nrows-1)) / numsamp 
+       write (*,*) 'mean(-Ln(like)) = ', meanlike
+       
+       meanlike = -log(sum(exp(-(coldata(2,0:nrows-1) -maxlike))*coldata(1,0:nrows-1)) & 
+                   / numsamp) + maxlike
+       write (*,*) '-Ln(mean like)  = ', meanlike
+ 
+ 
+ !Output files for 1D plots
+ 
+           plot_col =  nint(sqrt(num_vars/1.4))
+           plot_row = (num_vars +plot_col-1)/plot_col 
+    
+           !So we don't forget to run sm, delete current .ps file
+       call DeleteFile(trim(rootname)//'.ps')
+     
+     
+           open(unit=51,file=trim(rootname) // '.m',form='formatted',status='replace')
+          !MatLab file for 1D plots
+           call WriteMatLabInit(51,num_vars>3)
+ 
+           open(unit=50,file=sm_file,form='formatted',status='replace')
+        !SuperMongo .sm file for 1D plots
+           write (50,*) 'DEFINE TeX_strings 1'
+           write (50,*) 'device postfile '//trim(rootname)//'.ps'
+           write (50,*) 'DEFINE default_font rm'
+           if (num_vars > 6) then
+            write (50,*) 'expand 0.8'
+            write (50, *)'define x_gutter 0.6'
+            write (50,*) 'define y_gutter 0.8'
+           else
+            write (50,*) 'expand 1.1'
+           endif
+ 
+           write (50,*) 'location 2500 31000 3500 31000'
+           if (BW) then
+             write (50,*) 'lweight 4'
+           else
+              write (50,*) 'lweight 3'
+           end if
+ 
+          if (triangle_plot) then
+           open(unit=52,file=trim(rootname) // '_tri.m',form='formatted',status='replace')
+           call WriteMatLabInit(52,num_vars>4)
+          end if
+ 
+           cont_lines = 0
+ 
+ !Do 1D bins
+           do j = 1,num_vars
+              ix = colix(j)
+              bincounts = 0
+              binlikes = 0
+              binmaxlikes = 0
+ 
+             
+              do ix1 = 1, num_contours
+               limfrac = (1-contours(ix1))
+               if (.not. has_limits(ix) .or. force_twotail) then
+                 limfrac = limfrac/2 !two tail
+               end if
+               cont_lines(j,2,ix1) = ConfidVal(ix,limfrac,.true.)
+               cont_lines(j,1,ix1) = ConfidVal(ix,limfrac,.false.)
+              end do !contour lines
+ 
+ 
+              if (has_limits_bot(ix)) then
+               amin = limmin(ix)
+              else
+               amin = ConfidVal(ix,0.001,.false.)
+              end if
+ 
+              if (has_limits_top(ix)) then
+               amax = limmax(ix)
+              else
+               amax = ConfidVal(ix,0.001,.true.)
+              end if
+ 
+ 
+              width(j) = (amax-amin)/(nbins+1)
+              if (width(j)==0) cycle
+              if (amin <= 0) then 
+                center(j) = amax 
+              else
+                center(j) = amin 
+              end if
+ 
+              ix_min(j) = nint((amin - center(j))/width(j))
+              ix_max(j) = nint((amax - center(j))/width(j))
+ 
+              if (smoothing .and. .not. has_limits_bot(ix) .and. amin > 0) ix_min(j) = ix_min(j)-1
+              if (smoothing .and. .not. has_limits_top(ix)) ix_max(j) = ix_max(j)+1
+ 
+              do i = 0, nrows-1
+                 ix2=nint((coldata(ix,i)-center(j))/width(j))
+                 if (smoothing .and. ix_min(j) /= ix_max(j)) then
+                  do wx = max(ix_min(j),ix2-2),min(ix_max(j),ix2+2)
+                    dist = coldata(colix(j),i) - (center(j) + wx*width(j))
+                    distweight = exp(- dist**2/width(j)**2/2)
+                    binweight = coldata(1,i)*distweight
+                    bincounts(wx) = bincounts(wx)  + binweight
+                    binlikes(wx) = binlikes(wx) + binweight*exp(meanlike-coldata(2,i))
+                    binmaxlikes(wx) = max(binmaxlikes(wx), &
+                                 real(exp(maxlike-coldata(2,i))*distweight))
+                  end do
+                else
+                   bincounts(ix2) = bincounts(ix2) + coldata(1,i)
+                   binlikes(ix2) = binlikes(ix2) + coldata(1,i)*exp(meanlike-coldata(2,i))
+                   binmaxlikes(ix2) = max(binmaxlikes(ix2), real(exp(maxlike-coldata(2,i))))
+                end if
+              end do
+ 
+ 
+          do i=ix_min(j),ix_max(j)
+           if (bincounts(i) >0) binlikes(i) = binlikes(i)/bincounts(i)
+          end do
+ 
+           if (smoothing .and. ix_min(j) /= ix_max(j)) then
+             
+              !account for underweighting near edges
+               if (has_limits_bot(ix)) then
+                 bincounts(ix_min(j))   = bincounts(ix_min(j))*2
+                 bincounts(ix_min(j)+1) = bincounts(ix_min(j)+1)/0.84
+               end if
+               if (has_limits_top(ix)) then
+                 bincounts(ix_max(j))   = bincounts(ix_max(j))*2
+                 bincounts(ix_max(j)-1) = bincounts(ix_max(j)-1)/0.84
+               end if
+  
+           end if
+ 
+ 
+           if (.not. no_plots) then
+ 
+           fname =trim(trim(rootname)//trim(numcat('_p',colix(j)-2)))
+           filename = 'plot_data/'// trim(fname)
+           open(unit=49,file=trim(filename)//'.dat',form='formatted',status='replace')
+           maxbin = maxval(bincounts(ix_min(j):ix_max(j)))
+           if (maxbin==0) then
+               fname =  IntToStr(colix(j)-2)
+               write (*,*) 'no samples in bin, param:'//trim(fname)
+               stop
+            end if
+           do i = ix_min(j), ix_max(j)
+              write (49,'(3E16.7)') center(j) + i*width(j), bincounts(i)/maxbin, binmaxlikes(i)
+           end do
+           if (ix_min(j) == ix_max(j)) write (49,'(1E16.7,'' 0'')') center(j) + ix_min*width(j)
+           close(49)
+ 
+           !Detect ends which don't go to zero 
+           if (.not. force_twotail) then
+            if (bincounts(ix_max(j)) > maxbin*.15) cont_lines(j,2,:) = amax
+            if (bincounts(ix_min(j)) > maxbin*.15) cont_lines(j,1,:) = amin
+           end if
+                    
+           if (plot_meanlikes) then
+  
+            open(unit=49,file=trim(filename)//'.likes',form='formatted',status='replace')       
+            maxbin = maxval(binlikes(ix_min(j):ix_max(j)))
+            do i = ix_min(j), ix_max(j)
+                  write (49,'(2E16.7)') center(j) + i*width(j), binlikes(i)/maxbin
+                end do
+            close(49)
+           end if
+ 
+           plot_x = mod(j+plot_col-1,plot_col)+1
+           plot_y = plot_row - (j+plot_col-1)/plot_col + 1 
+ 
+            write (51,*) 'subplot(',plot_row,',',plot_col,',',j,')'
+            call Write1DplotMatLab(51,j);
+            if (matlab_version < 7) then
+             if (plot_col*plot_row> 4) call CheckMatlabAxes(51)
+              write (51,*) ' p = get(gca,''Position'');'
+             if (plot_col*plot_row > 8) then            
+              write (51,*) 'set(gca,''Position'',[p(1) (p(2)+p(4)/4) p(3) p(4)*4/5]);'
+             else
+              write (51,*) 'set(gca,''Position'',[p(1) (p(2)+p(4)/8) p(3) p(4)*8/9]);'
+             end if
+            end if
+            if (prob_label)  write (51,*) 'ylabel(''Probability'')'
+            write(51,*)  'xlabel('''//   trim(labels(ix))//''',''FontSize'',lab_fontsize);'
+            write (51,*) 'set(gca,''ytick'',[]);hold off;'
+ 
+            if (triangle_plot) then
+             write (52,*) 'subplot(',num_vars,',',num_vars,',',(j-1)*num_vars+j,')'
+             call Write1DplotMatLab(52,j);
+             if (num_vars > 2) call CheckMatlabAxes(52)
+             if (prob_label)  write (52,*) 'ylabel(''Probability'')'
+             if (j==num_vars) then
+              write(52,*)  'xlabel('''//   trim(labels(ix))//''',''FontSize'',lab_fontsize);'
+             end if
+             write (52,*) 'set(gca,''ytick'',[]);hold off;'
+            end if
+               
+             write (50,'(''window'',4I4)') plot_col,plot_row, plot_x, plot_y
+             write (50,*) 'data "' // trim(filename) //'.dat"'
+             write (50,*) 'read { x 1 y 2 mx 3}'
+             if (ix_min(j) == ix_max(j)) then
+               write (50,'(''limits '',2E16.7,'' y'')') center(j) - abs(center(j))/2-0.1, &
+                           center(j) + abs(center(j))/2+0.1
+             else
+              if (has_limits(ix)) then
+               write (50,'(''limits '',2E16.7,'' 0 1.1'')') amin, amax
+              else
+               write (50,*) 'limits x 0 1.1'
+              end if
+              end if
+           
+            write (50,*) 'box 1 3 3 3'
+            
+            write (50,*) 'connect x y'
+ !      write (50,*) 'connect x mx '
+        !uncomment to plot maximum likelihood values for each bin
+ 
+          if (plot_meanlikes) then
+          
+           write (50,*) 'data "' // trim(filename)//'.likes"'
+           write (50,*) 'read { x 1 y 2 }'
+           write (50,*) 'ltype 1'
+           write (50,*) 'connect x y'
+          end if
+          
+            do ix1 = 1, Num_ComparePlots
+              fname =  trim(trim(ComparePlots(ix1))//trim(numcat('_p',colix(j)-2)))    
+              if (FileExists('plot_data/' // trim(fname)//'.dat')) then
+ 
+               write (50,*) 'data "plot_data/' // trim(fname) //'.dat"'
+               write (50,*) 'read { x 1 y 2 mx 3}'
+               write (50,*) 'ltype 0'
+               if (BW) then
+                write (50,*) 'lweight 2'
+                if (ix1>1) write (50,*) 'ltype 2'
+               else
+                write (50,*) 'ctype ',ix1+2
+               end if
+              write (50,*) 'connect x y'
+              if (plot_meanlikes) then
+                 write (50,*) 'data "plot_data/'//trim(fname)//'.likes"'
+                 write (50,*) 'read { x 1 y 2}'
+                 write (50,*) 'ltype 1'
+                 write (50,*) 'connect x y'
+              end if
+              end if
+            end do
+            if (Num_ComparePlots/=0) write (50,*) 'ctype black'
+            if (Num_ComparePlots/=0 .and. BW) write (50,*) 'lweight 4'
+ 
+ 
+   !         write (50,*) 'ltype 1'
+   !         do ix1 = 1, num_contours
+   !         do i2 =1,2
+   !         write (numstr,*) cont_lines(j,i2,ix1)
+   !         write (50,*) 'set x = {'//trim(numstr)//' '//trim(numstr)//'}'
+   !         write (50,*) 'set y = { 0 1.1}'
+   !!         write (50,*) 'connect x y'
+   !uncomment this if you want dooted lines on 1D plots for condifence regions
+   !Note 1D confidence regions are computed from full samples using different definition to
+   !2D confidence regions
+ 
+   !         end do
+   !         end do
+ 
+            
+            write (50,*) 'ltype 0'
+            write (50,*) 'xlabel '// trim(labels(ix))
+            write (50,*)
+            write (50,*) '#Next...'
+ 
+             end if !no plots
+           end do
+ 
+    
+           if (.not. no_plots) then
+               write (50,*) 'quit'
+               close(50)
+             
+               write (51,*)  'set(gcf, ''PaperUnits'',''inches'');'
+      
+               if (plot_row < plot_col .and. plot_row*plot_col>9) then
+                 write (51,*) 'set(gcf, ''PaperPosition'',[ 0 0 8 10]);';
+               else
+                 write (51,*) 'set(gcf, ''PaperPosition'',[ 0 0 8 8]);';
+               end if
+               write (51,*) 'print -dpsc2 '//trim(rootname)//'.ps;'
+               close(51)
+           end if
+                  
+ 
+ !do 2D bins and produce matlab file
+       
+           if (plot_2D_param == 0 .and. num_cust2D_plots==0 .and. .not. no_plots) then
+              !In this case output the most correlated variable combinations
+           
+              write (*,*) 'doing 2D plots for most correlated variables'
+              try_t = 1e5
+ 
+              x=0;y=0;  
+             
+              num_cust2D_plots = 12
+               do j = 1, num_cust2D_plots
+                try_b = -1e5
+                do ix1 =1 ,num_vars
+                 do ix2 = ix1+1,num_vars 
+                    if (abs(corrmatrix(colix(ix1)-2,colix(ix2)-2)) < try_t .and. &
+                         abs(corrmatrix(colix(ix1)-2,colix(ix2)-2)) > try_b) then
+                       try_b = abs(corrmatrix(colix(ix1)-2,colix(ix2)-2))
+                       x = ix1; y = ix2;
+                     end if
+                 end do
+                end do
+                if (try_b == -1e5) then
+                   num_cust2D_plots = j-1
+                   exit
+                end if
+                try_t = try_b
+ 
+              cust2Dplots(j) = colix(x) + colix(y)*100
+              end do
+ 
+           end if
+ 
+ 
+       if (num_cust2D_plots == 0) then
+            num_2D_plots = 0
+         
+            do j= 1, num_vars
+              if (ix_min(j) /= ix_max(j)) then
+              do j2 = j+1, num_vars
+               if (ix_min(j2) /= ix_max(j2)) then
+                 if (plot_2D_param==0 .or. plot_2D_param == colix(j) .or. plot_2D_param == colix(j2)) &
+                       num_2D_plots = num_2D_plots + 1
+               end if
+             end do
+             end if
+            end do
+       else
+        num_2D_plots = num_cust2D_plots
+       end if
+  
+          done2D= .false.
+          if (num_2D_plots > 0 .and. .not. no_plots) then
+ 
+           write (*,*) 'Producing ',num_2D_plots,' 2D plots'
+           filename = trim(rootname)//'_2D.m'
+           open(unit=50,file=filename,form='formatted',status='replace')
+           call WriteMatLabInit(50,num_2D_plots >=7)
+ 
+           plot_col = nint(sqrt(num_2D_plots/1.4))
+           plot_row = (num_2D_plots +plot_col-1)/plot_col 
+           plot_num = 0
+     
+          
+            do j= 1, num_vars
+            if (ix_min(j) /= ix_max(j)) then
+             
+             if (plot_2D_param/=0 .or. num_cust2D_plots /= 0) then
+                if (colix(j) == plot_2D_param) cycle
+                j2min = 1
+             else
+                j2min = j+1
+             end if
+ 
+            do j2 = j2min, num_vars
+               if (ix_min(j2) /= ix_max(j2)) then
+ 
+               if (plot_2D_param/=0 .and. colix(j2) /= plot_2D_param) cycle
+               if (num_cust2D_plots /= 0 .and.  &
+                 count(cust2Dplots(1:num_cust2D_plots) == colix(j)*100 + colix(j2))==0) cycle
+       
+               plot_num = plot_num + 1
+               done2D(j,j2) = .true.
+               call Get2DPlotData(j,j2)
+           
+               write (50,'(''subplot('',1I5,'','',1I5,'','',1I5,'');'')') plot_row,plot_col,plot_num
+               call Write2DPlotMATLAB(50,j,j2,.true.,.true.)
+               if (plot_row*plot_col > 4) call CheckMatlabAxes(50)
+ 
+               end if
+            end do
+            end if
+           end do
+ 
+               
+           write (50,*)  'set(gcf, ''PaperUnits'',''inches'');'
+               if (matlab_col/='') write (50,*) trim(matlab_col)
+           if (num_2D_plots < 5) then
+                 write (50,*) 'set(gcf, ''PaperPosition'',[ 0 0 6 6]);';
+           else
+             write (50,*) 'set(gcf, ''PaperPosition'',[ 0 0 8 10]);';
+           end if
+ 
+           write (50,*) 'print -dpsc2 '//trim(rootname)//'_2D.ps;'
+           close(50)
+ 
+           end if
+ 
+          if (triangle_plot .and. .not. no_plots) then
+           !Add the off-diagonal 2D plots
+           do j = 1, num_vars
+             do j2 = j+1, num_vars
+               if (.not. Done2D(j2,j)) call Get2DPlotData(j2,j)
+               write (52,'(''subplot('',1I5,'','',1I5,'','',1I5,'');'')') num_vars, num_vars, &
+                    (j2-1)*num_vars + j                   
+               call Write2DPlotMATLAB(52,j2,j,j2==num_vars, j==1)
+               if (num_vars > 2) call CheckMatlabAxes(52)
+ 
+            end do
+           end do
+           write (52,*)  'set(gcf, ''PaperUnits'',''inches'');'
+           write (52,*) 'set(gcf, ''PaperPosition'',[ 0 0 8 8]);';
+           write (52,*) 'print -dpsc2 '//trim(rootname)//'_tri.ps;'
+           close(52)
+   
+          end if
+          
+ 
+ !Do 3D plots (i.e. 2D scatter plots with coloured points)
+ 
+           if (num_3D_plots /=0 .and. .not. no_plots) then
+               write (*,*) 'producing ',num_3D_plots, '2D colored scatter plots'
+               filename = trim(rootname)//'_3D.m'
+               open(unit=50,file=filename,form='formatted',status='replace')
+               write (50,*) 'colormap(''jet'');'
+               if (num_3D_plots ==1 ) then
+                 write(50,*) 'lab_fontsize = 12; axes_fontsize = 12;'
+               else
+                 write(50,*) 'lab_fontsize = 9; axes_fontsize = 9;'
+               end if
+               if (mod(num_3D_plots,2)==0 .and. num_3D_plots < 11) then
+                plot_col = num_3D_plots/2
+               else
+                plot_col =  nint(sqrt(1.*num_3D_plots))
+               end if
+               plot_row = (num_3D_plots +plot_col-1)/plot_col 
+ 
+            
+               write(50,*) 'pts = load('''//trim(rootname)//'_single.txt'');'
+              
+               do j=1, num_3D_plots
+                  read(plot_3D(j),*) ix1, ix2, ix3  !x, y, color
+                  if (ix3<1) ix3 = MostCorrelated2D(ix1,ix2,ix3)
+                  
+                  write (50,'(''subplot('',1I5,'','',1I5,'','',1I5,'');'')') plot_row,plot_col,j   
+                  write (50,*) '%Do params ',ix1,ix2,ix3
+                  write (50,*) 'scatter(pts(:,',2+ix1,'),pts(:,',2+ix2,'),3,pts(:,',2+ix3,'));'
+                  fmt = ''',''FontSize'',lab_fontsize);'
+                  write (50,*) 'xlabel('''//trim(labels(ix1+2))//trim(fmt)
+                  write (50,*) 'ylabel('''//trim(labels(ix2+2))//trim(fmt)
+                  write (50,*) 'set(gca,''FontSize'',axes_fontsize); pa = gca;'
+                  write (50,*) 'hbar = colorbar(''horiz'');axes(hbar);'
+                  write (50,*) 'xlabel('''//trim(labels(ix3+2))//trim(fmt)
+                  write (50,*) 'set(gca,''FontSize'',axes_fontsize);'
+                  if (num_3D_plots > 2 .and. matlab_version < 7) then
+                    write (50,*) ' p = get(pa,''Position'');'
+                    write (50,*) 'set(pa,''Position'',[p(1) (p(2)+p(4)/8) p(3) p(4)]);'
+                  end if
+                end do
+            
+               write (50,*)  'set(gcf, ''PaperUnits'',''inches'');'
+      
+               tmpstr = RealToStr((plot_col*9.)/plot_row)
+               write (50,*) 'set(gcf, ''PaperPosition'',[ 0 0 '//  &
+                             trim(tmpstr) //' 8]);'
+ 
+               write (50,*) 'print -dpsc2 '//trim(rootname)//'_3D.ps;'
+               close(50)
+           end if
+     
+    
+ 
+ !write out stats
+ !Marginalized
+ 
+          open(unit=50,file=trim(rootname)//'.margestats',form='formatted',status='replace')
+           write(50,'(a)') 'param  mean          sddev         lower1        upper1        lower2        upper2'
+           open(unit=51,file='plot_data/'//trim(rootname)//'_params',form='formatted',status='replace')
+           do j=1, num_vars
+              write(50,'(1I5,6E14.6,"   '//trim(labels(colix(j)))//'")') colix(j)-2, mean(j), sddev(j), &
+                   cont_lines(j,1:2,1),cont_lines(j,1:2,2)
+              write (51,*) colix(j)-2, trim(labels(colix(j)))
+           end do
+           close(51)
+           write (50,*) '' 
+           write (50,'(a)') 'Limits are: ' // trim(contours_str)
+           if (.not. force_twotail) then
+               do j=1, num_vars
+                if (has_limits(colix(j))) then
+                  write(50,'(1I5," one tail;  '//trim(labels(colix(j)))//'")') colix(j)-2
+                 else
+                  write(50,'(1I5," two tail;  '//trim(labels(colix(j)))//'")') colix(j)-2
+                 end if
+               end do
+           else
+               write (50,*) 'All limits are two tail'
+           end if
+           close(50)
+ 
+ 
+ !Limits from global likelihood
+           open(unit=50,file=trim(rootname)//'.likestats',form='formatted',status='replace')
+           write (50,*) 'Best fit sample -log(Like) = ',coldata(2,bestfit_ix)
+           write (50,*) ''
+           write(50,'(a)') 'param  bestfit       lower1        upper1        lower2        upper2'
+           
+           do j=1, num_vars
+      
+           write(50,'(1I5,5E14.6,"   '//trim(labels(colix(j)))//'")') colix(j)-2, coldata(colix(j),bestfit_ix),&
+           minval(coldata(colix(j),0:ND_cont1)), &
+           maxval(coldata(colix(j),0:ND_cont1)), &
+           minval(coldata(colix(j),0:ND_cont2)), &
+           maxval(coldata(colix(j),0:ND_cont2))
+ 
+           end do
+           close(50)
+ 
+ 
+ 
+ end program GetDist
+ 
diff -r -N -c -b -B cosmomc/source/Makefile cosmomc_addon/source/Makefile
*** cosmomc/source/Makefile	2006-04-03 13:37:37.000000000 +0100
--- cosmomc_addon/source/Makefile	2006-04-04 10:35:03.000000000 +0100
***************
*** 1,160 ****
! #Settings here are for Intel ifort on CITA cluster; uncomment others below
! #You may need to edit the library paths for MKL for Intel
! #Beware of using optmizations which lose accuracy - may give errors when running
! #With Visual Fortan open cosmomc.dsw workspace file instead.
! 
! #Edit for CFITSIO, HEALPIX directories
! cfitsio = /home/pascal/include/cfitsio
! healpix = /home/pascal/include/Healpix_2.00
! WMAP3 = /home/antlewis/WMAP3
! 
! #Intel MPI (assuming mpif77 set to point to ifc)
! #change -lmkl_ia32 to -lmkl_p3 for MKL versions earlier than 6 (6 needed for ifc 8+) 
! F90C     = mpif77
! FFLAGS = -O2 -Vaxlib -ip -W0 -WB -openmp -fpp -DMPI
! LAPACKL = -L/opt/intel/mkl61/lib/32 -lmkl_lapack -lmkl_ia32 -lguide -lpthread
! INCLUDE = -I../camb 
  
! # Use with Intel 9 on the Altix (eg. COSMOS)
! # (do "module load icomp90" before compiling, or add to .bashrc)
! #F90C = ifort
! #FFLAGS = -O3 -w -fpp2 -DMPI 
! #LAPACKL = -lscs -lguide -lmpi -ldl
! #INCLUDE = -I../camb
! 
! #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
! #INCLUDE = -I../camb 
! 
! #Intel fortran 8 on COSMOS
! #F90C     = ifort
! #FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp2 -lmpi -lscs  -ldl -DMPI
! #LAPACKL = -L/opt/intel/mkl80/lib/64 -lmkl_lapack -lmkl_ipf -limf
! #INCLUDE = -I../camb
! 
! #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
! #INCLUDE = -I../camb 
! #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
! #INCLUDE = -I../camb 
! 
! #Intel Itanium efc, -openmp toggles multi-processor:
! #F90C     = efc
! #FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp
! #LAPACKL = -L/opt/intel/mkl/lib/64 -lmkl_lapack -lmkl_ipf -lguide -lpthread -static
! #INCLUDE = -I../camb 
! 
! #Intel Itanium efc MPI:
! # mpirun -np 4 ./cosmomc params.ini
! #F90C     = efc
! #FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp -lmpi -lscs  -ldl -DMPI
! #LAPACKL = -L/opt/intel/mkl/lib/64 -lmkl_lapack -lmkl_ipf -limf
! #INCLUDE = -I../camb
! 
! #G95; make sure LAPACK and MPI libs also compiled with g95
! #F90C     = mpif90
! #FFLAGS = -O2 -DMPI
! #LAPACKL = /LAPACK/lapack_LINUX.a /LAPACK/blas_LINUX.a
! #INCLUDE = -I../camb
! 
! #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 
! #INCLUDE = -I../camb -p../camb
! 
! #NAGF95, single processor:
! #F90C     = f95
! #FFLAGS = -DNAGF95 -O3
! #LAPACKL = -llapack -lblas -lg2c
! #INCLUDE = -I../camb 
! 
! #Sun, single processor:
! #F90C     = f90
! #FFLAGS = -O5
! #LAPACKL = -lsunperf -lfsu
! #INCLUDE = -I../camb -M../camb
! 
! #Sun MPI
! #F90C = mpf90
! #FFLAGS =  -O4 -openmp -ftrap=%none -dalign -DMPI
! #LAPACKL = -lsunperf -lfsu -lmpi_mt
! #INCLUDE = -I../camb -M../camb
! 
! #Sun parallel enterprise:
! #F90C     = f95
! #FFLAGS =  -O4 -xarch=native64 -openmp -ftrap=%none
! #LAPACKL = -lsunperf -lfsu
! #INCLUDE = -I../camb -M../camb
! 
! #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 = -I../camb -lessl
! #LAPACKL =
  
  PROPOSE = propose.o
  CLSFILE = CMB_Cls_simple.o
  
  #Can use params_H if you prefer more generic parameters
  PARAMETERIZATION = params_CMB.o
  
! WMAP3_inc = -I$(cfitsio)/include -I$(healpix)/include -I$(WMAP3)
! WMAP3_lib = -L$(cfitsio)/lib -L$(healpix)/lib -L$(WMAP3) \
! 	-lhealpix -lcfitsio
  
! CLSLIB  = $(WMAP3_lib) -L../camb -lcamb $(LAPACKL) 
  
! F90FLAGS = $(FFLAGS) $(WMAP3_inc) $(INCLUDE)
  
  DISTFILES = utils.o settings.o GetDist.o
  
! WMAPOBJS = $(WMAP3)/read_archive_map.o \
! 	$(WMAP3)/read_fits.o \
! 	$(WMAP3)/WMAP_3yr_options.o \
! 	$(WMAP3)/WMAP_3yr_util.o \
! 	$(WMAP3)/WMAP_3yr_tt_pixlike.o \
! 	$(WMAP3)/WMAP_3yr_tt_beam_and_ptsrc_corr.o \
! 	$(WMAP3)/WMAP_3yr_teeebb_pixlike.o \
! 	$(WMAP3)/WMAP_3yr_likelihood.o
  
  OBJFILES= $(WMAPOBJS) utils.o settings.o cmbtypes.o cmbdata.o WeakLen.o twodf.o mpk.o \
! 	supernovae.o SDSSLy-a-v3.o\
  	$(CLSFILE) paramdef.o $(PROPOSE) $(PARAMETERIZATION) calclike.o \
  	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
  
! default: cosmomc
  
! all : cosmomc getdist
  
  .f.o:
  	f77 $(F90FLAGS) -c $<
  
  %.o: %.f90
! 	$(F90C) $(F90FLAGS) -c $*.f90
  
  %.o: %.F90
! 	$(F90C) $(F90FLAGS) -c $*.F90
! 
  
- cosmomc: $(OBJFILES) 	
- 	$(F90C) -o ../cosmomc $(OBJFILES) $(CLSLIB) $(F90FLAGS)
  
  
  clean:
  	rm -f *.o *.mod *.d *.pc *.obj ../core
  
! getdist: $(DISTFILES)
! 	$(F90C) -o ../getdist $(DISTFILES) $(CLSLIB) $(F90FLAGS) 
--- 1,103 ----
! # >>> DESIGNED FOR GMAKE <<<
! 
! # Unified Systems makefile for COSMOMC
! # Add FLAGS -DMPI for using MPI
  
! ext=$(shell uname | cut -c1-3)
! 
! ifeq ($(ext),IRI)
! F90C= f90
! FFLAGS= -Ofast -mp -n32 -LANG:recursive=ON -lmpi -DMPI
! WMAPFLAGS= $(FFLAGS)
  #LAPACKL = -lcomplib.sgimath
! LAPACKL = -lcomplib.sgimath_mp
! INCLUDE = -I../camb
! CFITSIODIR = /people/ringeval/usr
! HEALPIXDIR = /people/ringeval/usr
! endif
! 
! ifeq ($(ext),Lin)
! F90C=ifort
! FFLAGS= -O -openmp -fpp
! WMAPFLAGS= -O -save -fpp
! LAPACKL = -lmkl -lmkl_lapack
! INCLUDE = -I../camb 
! CFITSIODIR = /opt/intel/cfitsio-3.006
! HEALPIXDIR = /opt/intel/healpix-2.0.1
! endif 
! 
! ifeq ($(ext),OSF)
! F90C=f90
! FFLAGS= -omp -O -arch host -math_library fast -tune host -fpe1
! WMAPFLAGS= $(FFLAGS)
! LAPACKL  = -lcxml
! INCLUDE = -I../camb 
! endif
  
! ifeq ($(ext),Sun)
! F90C=f90
! FFLAGS= -O4 -xarch=native64 -openmp -ftrap=%none
! WMAPFLAGS= $(FFLAGS)
! LAPACKL = -lsunperf -lfsu
! INCLUDE = -I../camb -M../camb
! endif
! 
! ifeq ($(ext),AIX)
! F90C    = mpxlf90_r
! FFLAGS  = -O4 -WF,-DIBMXL,-DMPI -qstrict -qsmp=omp -qmaxmem=-1 -qsuffix=f=f90:cpp=F90
! WMAPFLAGS= $(FFLAGS)
! INCLUDE = -I../camb -lessl
! LAPACKL =
! endif
  
+ WMAPDIR = ../WMAP
  PROPOSE = propose.o
  CLSFILE = CMB_Cls_simple.o
  
  #Can use params_H if you prefer more generic parameters
  PARAMETERIZATION = params_CMB.o
  
! WMAPINCLUDE =  -I$(CFITSIODIR)/include -I$(HEALPIXDIR)/include
! # -I$(WMAPDIR)
  
! CLSLIB  = -L$(HEALPIXDIR)/lib -L$(CFITSIODIR)/lib \
! 	-L../camb -lcamb -lhealpix -lcfitsio $(LAPACKL)
  
! F90FLAGS = $(FFLAGS) $(INCLUDE)
  
  DISTFILES = utils.o settings.o GetDist.o
  
! WMAPOBJS = read_archive_map.o read_fits.o WMAP_3yr_options.o WMAP_3yr_util.o \
! 	WMAP_3yr_tt_pixlike.o WMAP_3yr_tt_beam_and_ptsrc_corr.o \
! 	WMAP_3yr_teeebb_pixlike.o WMAP_3yr_likelihood.o
  
  OBJFILES= $(WMAPOBJS) utils.o settings.o cmbtypes.o cmbdata.o WeakLen.o twodf.o mpk.o \
! 	supernovae.o lya.o\
  	$(CLSFILE) paramdef.o $(PROPOSE) $(PARAMETERIZATION) calclike.o \
  	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
  
! default: cosmomc.$(ext)
  
! all : cosmomc.$(ext) getdist.$(ext)
  
  .f.o:
  	f77 $(F90FLAGS) -c $<
  
+ %.o: $(WMAPDIR)/%.f90
+ 	$(F90C) $(WMAPFLAGS) $(WMAPINCLUDE) -c $<
+ 
  %.o: %.f90
! 	$(F90C) $(F90FLAGS) -c $<
  
  %.o: %.F90
! 	$(F90C) $(F90FLAGS) -c $<
  
  
+ cosmomc.$(ext): $(OBJFILES) 	
+ 	$(F90C) -o ../$@ $(OBJFILES) $(CLSLIB) $(F90FLAGS)
  
  clean:
  	rm -f *.o *.mod *.d *.pc *.obj ../core
  
! getdist.$(ext): $(DISTFILES)
! 	$(F90C) -o ../$@ $(DISTFILES) $(CLSLIB) $(F90FLAGS) 
diff -r -N -c -b -B cosmomc/source/Makefile.orig cosmomc_addon/source/Makefile.orig
*** cosmomc/source/Makefile.orig	1970-01-01 01:00:00.000000000 +0100
--- cosmomc_addon/source/Makefile.orig	2006-04-04 10:34:22.000000000 +0100
***************
*** 0 ****
--- 1,160 ----
+ #Settings here are for Intel ifort on CITA cluster; uncomment others below
+ #You may need to edit the library paths for MKL for Intel
+ #Beware of using optmizations which lose accuracy - may give errors when running
+ #With Visual Fortan open cosmomc.dsw workspace file instead.
+ 
+ #Edit for CFITSIO, HEALPIX directories
+ cfitsio = /home/pascal/include/cfitsio
+ healpix = /home/pascal/include/Healpix_2.00
+ WMAP3 = /home/antlewis/WMAP3
+ 
+ #Intel MPI (assuming mpif77 set to point to ifc)
+ #change -lmkl_ia32 to -lmkl_p3 for MKL versions earlier than 6 (6 needed for ifc 8+) 
+ F90C     = mpif77
+ FFLAGS = -O2 -Vaxlib -ip -W0 -WB -openmp -fpp -DMPI
+ LAPACKL = -L/opt/intel/mkl61/lib/32 -lmkl_lapack -lmkl_ia32 -lguide -lpthread
+ INCLUDE = -I../camb 
+ 
+ # Use with Intel 9 on the Altix (eg. COSMOS)
+ # (do "module load icomp90" before compiling, or add to .bashrc)
+ #F90C = ifort
+ #FFLAGS = -O3 -w -fpp2 -DMPI 
+ #LAPACKL = -lscs -lguide -lmpi -ldl
+ #INCLUDE = -I../camb
+ 
+ #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
+ #INCLUDE = -I../camb 
+ 
+ #Intel fortran 8 on COSMOS
+ #F90C     = ifort
+ #FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp2 -lmpi -lscs  -ldl -DMPI
+ #LAPACKL = -L/opt/intel/mkl80/lib/64 -lmkl_lapack -lmkl_ipf -limf
+ #INCLUDE = -I../camb
+ 
+ #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
+ #INCLUDE = -I../camb 
+ #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
+ #INCLUDE = -I../camb 
+ 
+ #Intel Itanium efc, -openmp toggles multi-processor:
+ #F90C     = efc
+ #FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp
+ #LAPACKL = -L/opt/intel/mkl/lib/64 -lmkl_lapack -lmkl_ipf -lguide -lpthread -static
+ #INCLUDE = -I../camb 
+ 
+ #Intel Itanium efc MPI:
+ # mpirun -np 4 ./cosmomc params.ini
+ #F90C     = efc
+ #FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp -lmpi -lscs  -ldl -DMPI
+ #LAPACKL = -L/opt/intel/mkl/lib/64 -lmkl_lapack -lmkl_ipf -limf
+ #INCLUDE = -I../camb
+ 
+ #G95; make sure LAPACK and MPI libs also compiled with g95
+ #F90C     = mpif90
+ #FFLAGS = -O2 -DMPI
+ #LAPACKL = /LAPACK/lapack_LINUX.a /LAPACK/blas_LINUX.a
+ #INCLUDE = -I../camb
+ 
+ #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 
+ #INCLUDE = -I../camb -p../camb
+ 
+ #NAGF95, single processor:
+ #F90C     = f95
+ #FFLAGS = -DNAGF95 -O3
+ #LAPACKL = -llapack -lblas -lg2c
+ #INCLUDE = -I../camb 
+ 
+ #Sun, single processor:
+ #F90C     = f90
+ #FFLAGS = -O5
+ #LAPACKL = -lsunperf -lfsu
+ #INCLUDE = -I../camb -M../camb
+ 
+ #Sun MPI
+ #F90C = mpf90
+ #FFLAGS =  -O4 -openmp -ftrap=%none -dalign -DMPI
+ #LAPACKL = -lsunperf -lfsu -lmpi_mt
+ #INCLUDE = -I../camb -M../camb
+ 
+ #Sun parallel enterprise:
+ #F90C     = f95
+ #FFLAGS =  -O4 -xarch=native64 -openmp -ftrap=%none
+ #LAPACKL = -lsunperf -lfsu
+ #INCLUDE = -I../camb -M../camb
+ 
+ #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 = -I../camb -lessl
+ #LAPACKL =
+ 
+ PROPOSE = propose.o
+ CLSFILE = CMB_Cls_simple.o
+ 
+ #Can use params_H if you prefer more generic parameters
+ PARAMETERIZATION = params_CMB.o
+ 
+ WMAP3_inc = -I$(cfitsio)/include -I$(healpix)/include -I$(WMAP3)
+ WMAP3_lib = -L$(cfitsio)/lib -L$(healpix)/lib -L$(WMAP3) \
+ 	-lhealpix -lcfitsio
+ 
+ CLSLIB  = $(WMAP3_lib) -L../camb -lcamb $(LAPACKL) 
+ 
+ F90FLAGS = $(FFLAGS) $(WMAP3_inc) $(INCLUDE)
+ 
+ DISTFILES = utils.o settings.o GetDist.o
+ 
+ WMAPOBJS = $(WMAP3)/read_archive_map.o \
+ 	$(WMAP3)/read_fits.o \
+ 	$(WMAP3)/WMAP_3yr_options.o \
+ 	$(WMAP3)/WMAP_3yr_util.o \
+ 	$(WMAP3)/WMAP_3yr_tt_pixlike.o \
+ 	$(WMAP3)/WMAP_3yr_tt_beam_and_ptsrc_corr.o \
+ 	$(WMAP3)/WMAP_3yr_teeebb_pixlike.o \
+ 	$(WMAP3)/WMAP_3yr_likelihood.o
+ 
+ OBJFILES= $(WMAPOBJS) utils.o settings.o cmbtypes.o cmbdata.o WeakLen.o twodf.o mpk.o \
+ 	supernovae.o SDSSLy-a-v3.o\
+ 	$(CLSFILE) paramdef.o $(PROPOSE) $(PARAMETERIZATION) calclike.o \
+ 	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
+ 
+ default: cosmomc
+ 
+ all : cosmomc getdist
+ 
+ .f.o:
+ 	f77 $(F90FLAGS) -c $<
+ 
+ %.o: %.f90
+ 	$(F90C) $(F90FLAGS) -c $*.f90
+ 
+ %.o: %.F90
+ 	$(F90C) $(F90FLAGS) -c $*.F90
+ 
+ 
+ cosmomc: $(OBJFILES) 	
+ 	$(F90C) -o ../cosmomc $(OBJFILES) $(CLSLIB) $(F90FLAGS)
+ 
+ 
+ clean:
+ 	rm -f *.o *.mod *.d *.pc *.obj ../core
+ 
+ getdist: $(DISTFILES)
+ 	$(F90C) -o ../getdist $(DISTFILES) $(CLSLIB) $(F90FLAGS) 
diff -r -N -c -b -B cosmomc/source/Makefile.rej cosmomc_addon/source/Makefile.rej
*** cosmomc/source/Makefile.rej	1970-01-01 01:00:00.000000000 +0100
--- cosmomc_addon/source/Makefile.rej	2006-04-04 10:34:35.000000000 +0100
***************
*** 0 ****
--- 1,265 ----
+ ***************
+ *** 1,160 ****
+ ! #Settings here are for Intel ifort on CITA cluster; uncomment others below
+ ! #You may need to edit the library paths for MKL for Intel
+ ! #Beware of using optmizations which lose accuracy - may give errors when running
+ ! #With Visual Fortan open cosmomc.dsw workspace file instead.
+ ! 
+ ! #Edit for CFITSIO, HEALPIX directories
+ ! cfitsio = /home/pascal/include/cfitsio
+ ! healpix = /home/pascal/include/Healpix_2.00
+ ! WMAP3 = /home/antlewis/WMAP3
+ ! 
+ ! #Intel MPI (assuming mpif77 set to point to ifc)
+ ! #change -lmkl_ia32 to -lmkl_p3 for MKL versions earlier than 6 (6 needed for ifc 8+) 
+ ! F90C     = mpif77
+ ! FFLAGS = -O2 -Vaxlib -ip -W0 -WB -openmp -fpp -DMPI
+ ! LAPACKL = -L/opt/intel/mkl61/lib/32 -lmkl_lapack -lmkl_ia32 -lguide -lpthread
+ ! INCLUDE = -I../camb 
+   
+ ! # Use with Intel 9 on the Altix (eg. COSMOS)
+ ! # (do "module load icomp90" before compiling, or add to .bashrc)
+ ! #F90C = ifort
+ ! #FFLAGS = -O3 -w -fpp2 -DMPI 
+ ! #LAPACKL = -lscs -lguide -lmpi -ldl
+ ! #INCLUDE = -I../camb
+ ! 
+ ! #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
+ ! #INCLUDE = -I../camb 
+ ! 
+ ! #Intel fortran 8 on COSMOS
+ ! #F90C     = ifort
+ ! #FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp2 -lmpi -lscs  -ldl -DMPI
+ ! #LAPACKL = -L/opt/intel/mkl80/lib/64 -lmkl_lapack -lmkl_ipf -limf
+ ! #INCLUDE = -I../camb
+ ! 
+ ! #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
+ ! #INCLUDE = -I../camb 
+ ! #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
+ ! #INCLUDE = -I../camb 
+ ! 
+ ! #Intel Itanium efc, -openmp toggles multi-processor:
+ ! #F90C     = efc
+ ! #FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp
+ ! #LAPACKL = -L/opt/intel/mkl/lib/64 -lmkl_lapack -lmkl_ipf -lguide -lpthread -static
+ ! #INCLUDE = -I../camb 
+ ! 
+ ! #Intel Itanium efc MPI:
+ ! # mpirun -np 4 ./cosmomc params.ini
+ ! #F90C     = efc
+ ! #FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp -lmpi -lscs  -ldl -DMPI
+ ! #LAPACKL = -L/opt/intel/mkl/lib/64 -lmkl_lapack -lmkl_ipf -limf
+ ! #INCLUDE = -I../camb
+ ! 
+ ! #G95; make sure LAPACK and MPI libs also compiled with g95
+ ! #F90C     = mpif90
+ ! #FFLAGS = -O2 -DMPI
+ ! #LAPACKL = /LAPACK/lapack_LINUX.a /LAPACK/blas_LINUX.a
+ ! #INCLUDE = -I../camb
+ ! 
+ ! #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 
+ ! #INCLUDE = -I../camb -p../camb
+ ! 
+ ! #NAGF95, single processor:
+ ! #F90C     = f95
+ ! #FFLAGS = -DNAGF95 -O3
+ ! #LAPACKL = -llapack -lblas -lg2c
+ ! #INCLUDE = -I../camb 
+ ! 
+ ! #Sun, single processor:
+ ! #F90C     = f90
+ ! #FFLAGS = -O5
+ ! #LAPACKL = -lsunperf -lfsu
+ ! #INCLUDE = -I../camb -M../camb
+ ! 
+ ! #Sun MPI
+ ! #F90C = mpf90
+ ! #FFLAGS =  -O4 -openmp -ftrap=%none -dalign -DMPI
+ ! #LAPACKL = -lsunperf -lfsu -lmpi_mt
+ ! #INCLUDE = -I../camb -M../camb
+ ! 
+ ! #Sun parallel enterprise:
+ ! #F90C     = f95
+ ! #FFLAGS =  -O4 -xarch=native64 -openmp -ftrap=%none
+ ! #LAPACKL = -lsunperf -lfsu
+ ! #INCLUDE = -I../camb -M../camb
+ ! 
+ ! #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 = -I../camb -lessl
+ ! #LAPACKL =
+   
+   PROPOSE = propose.o
+   CLSFILE = CMB_Cls_simple.o
+   
+   #Can use params_H if you prefer more generic parameters
+   PARAMETERIZATION = params_CMB.o
+   
+ ! WMAP3_inc = -I$(cfitsio)/include -I$(healpix)/include -I$(WMAP3)
+ ! WMAP3_lib = -L$(cfitsio)/lib -L$(healpix)/lib -L$(WMAP3) \
+ ! 	-lhealpix -lcfitsio
+   
+ ! CLSLIB  = $(WMAP3_lib) -L../camb -lcamb $(LAPACKL) 
+   
+ ! F90FLAGS = $(FFLAGS) $(WMAP3_inc) $(INCLUDE)
+   
+   DISTFILES = utils.o settings.o GetDist.o
+   
+ ! WMAPOBJS = $(WMAP3)/read_archive_map.o \
+ ! 	$(WMAP3)/read_fits.o \
+ ! 	$(WMAP3)/WMAP_3yr_options.o \
+ ! 	$(WMAP3)/WMAP_3yr_util.o \
+ ! 	$(WMAP3)/WMAP_3yr_tt_pixlike.o \
+ ! 	$(WMAP3)/WMAP_3yr_tt_beam_and_ptsrc_corr.o \
+ ! 	$(WMAP3)/WMAP_3yr_teeebb_pixlike.o \
+ ! 	$(WMAP3)/WMAP_3yr_likelihood.o
+   
+   OBJFILES= $(WMAPOBJS) utils.o settings.o cmbtypes.o cmbdata.o WeakLen.o twodf.o mpk.o \
+   	supernovae.o lya.o\
+   	$(CLSFILE) paramdef.o $(PROPOSE) $(PARAMETERIZATION) calclike.o \
+   	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
+   
+ ! default: cosmomc
+   
+ ! all : cosmomc getdist
+   
+   .f.o:
+   	f77 $(F90FLAGS) -c $<
+   
+   %.o: %.f90
+ ! 	$(F90C) $(F90FLAGS) -c $*.f90
+   
+   %.o: %.F90
+ ! 	$(F90C) $(F90FLAGS) -c $*.F90
+ ! 
+   
+ - cosmomc: $(OBJFILES) 	
+ - 	$(F90C) -o ../cosmomc $(OBJFILES) $(CLSLIB) $(F90FLAGS)
+   
+   
+   clean:
+   	rm -f *.o *.mod *.d *.pc *.obj ../core
+   
+ ! getdist: $(DISTFILES)
+ ! 	$(F90C) -o ../getdist $(DISTFILES) $(CLSLIB) $(F90FLAGS) 
+ --- 1,102 ----
+ ! # >>> DESIGNED FOR GMAKE <<<
+ ! 
+ ! # Unified Systems makefile for COSMOMC
+ ! # Add FLAGS -DMPI for using MPI
+   
+ ! ext=$(shell uname | cut -c1-3)
+ ! 
+ ! ifeq ($(ext),IRI)
+ ! F90C= f90
+ ! FFLAGS= -Ofast -mp -n32 -LANG:recursive=ON -lmpi -DMPI
+ ! WMAPFLAGS= $(FFLAGS)
+   #LAPACKL = -lcomplib.sgimath
+ ! LAPACKL = -lcomplib.sgimath_mp
+ ! INCLUDE = -I../camb
+ ! CFITSIODIR = /people/ringeval/usr
+ ! HEALPIXDIR = /people/ringeval/usr
+ ! endif
+ ! 
+ ! ifeq ($(ext),Lin)
+ ! F90C=ifort
+ ! FFLAGS= -O -openmp -fpp
+ ! WMAPFLAGS= -O -save -fpp
+ ! LAPACKL = -lmkl -lmkl_lapack
+ ! INCLUDE = -I../camb 
+ ! CFITSIODIR = /opt/intel/cfitsio-3.006
+ ! HEALPIXDIR = /opt/intel/healpix-2.0.1
+ ! endif 
+ ! 
+ ! ifeq ($(ext),OSF)
+ ! F90C=f90
+ ! FFLAGS= -omp -O -arch host -math_library fast -tune host -fpe1
+ ! WMAPFLAGS= $(FFLAGS)
+ ! LAPACKL  = -lcxml
+ ! INCLUDE = -I../camb 
+ ! endif
+   
+ ! ifeq ($(ext),Sun)
+ ! F90C=f90
+ ! FFLAGS= -O4 -xarch=native64 -openmp -ftrap=%none
+ ! WMAPFLAGS= $(FFLAGS)
+ ! LAPACKL = -lsunperf -lfsu
+ ! INCLUDE = -I../camb -M../camb
+ ! endif
+ ! 
+ ! ifeq ($(ext),AIX)
+ ! F90C    = mpxlf90_r
+ ! FFLAGS  = -O4 -WF,-DIBMXL,-DMPI -qstrict -qsmp=omp -qmaxmem=-1 -qsuffix=f=f90:cpp=F90
+ ! WMAPFLAGS= $(FFLAGS)
+ ! INCLUDE = -I../camb -lessl
+ ! LAPACKL =
+ ! endif
+   
+ + WMAPDIR = ../WMAP
+   PROPOSE = propose.o
+   CLSFILE = CMB_Cls_simple.o
+   
+   #Can use params_H if you prefer more generic parameters
+   PARAMETERIZATION = params_CMB.o
+   
+ ! WMAPINCLUDE =  -I$(CFITSIODIR)/include -I$(HEALPIXDIR)/include -I$(WMAPDIR)
+   
+ ! CLSLIB  = -L$(HEALPIXDIR)/lib -L$(CFITSIODIR)/lib \
+ ! 	-L../camb -lcamb -lhealpix -lcfitsio $(LAPACKL)
+   
+ ! F90FLAGS = $(FFLAGS) $(INCLUDE)
+   
+   DISTFILES = utils.o settings.o GetDist.o
+   
+ ! WMAPOBJS = read_archive_map.o read_fits.o WMAP_3yr_options.o WMAP_3yr_util.o \
+ ! 	WMAP_3yr_tt_pixlike.o WMAP_3yr_tt_beam_and_ptsrc_corr.o \
+ ! 	WMAP_3yr_teeebb_pixlike.o WMAP_3yr_likelihood.o
+   
+   OBJFILES= $(WMAPOBJS) utils.o settings.o cmbtypes.o cmbdata.o WeakLen.o twodf.o mpk.o \
+   	supernovae.o lya.o\
+   	$(CLSFILE) paramdef.o $(PROPOSE) $(PARAMETERIZATION) calclike.o \
+   	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
+   
+ ! default: cosmomc.$(ext)
+   
+ ! all : cosmomc.$(ext) getdist.$(ext)
+   
+   .f.o:
+   	f77 $(F90FLAGS) -c $<
+   
+ + %.o: $(WMAPDIR)/%.f90
+ + 	$(F90C) $(WMAPFLAGS) $(WMAPINCLUDE) -c $<
+ + 
+   %.o: %.f90
+ ! 	$(F90C) $(F90FLAGS) -c $<
+   
+   %.o: %.F90
+ ! 	$(F90C) $(F90FLAGS) -c $<
+   
+   
+ + cosmomc.$(ext): $(OBJFILES) 	
+ + 	$(F90C) -o ../$@ $(OBJFILES) $(CLSLIB) $(F90FLAGS)
+   
+   clean:
+   	rm -f *.o *.mod *.d *.pc *.obj ../core
+   
+ ! getdist.$(ext): $(DISTFILES)
+ ! 	$(F90C) -o ../$@ $(DISTFILES) $(CLSLIB) $(F90FLAGS) 
diff -r -N -c -b -B cosmomc/source/MCMC.f90 cosmomc_addon/source/MCMC.f90
*** cosmomc/source/MCMC.f90	2006-03-17 16:18:11.000000000 +0000
--- cosmomc_addon/source/MCMC.f90	2006-04-04 10:34:35.000000000 +0100
***************
*** 313,320 ****
      end if
  
      Like = GetLogLike(Trial) 
!     if (Feedback > 1) write (*,*) 'Likelihood: ', Like, 'Current Like:', CurLike
! 
      accpt = (Like /= logZero) .and. (CurLike > Like .or. randexp1() > Like - CurLike) 
               !Include the min() so that compilers not doing optimal compilation don't complain
   
--- 313,322 ----
      end if
  
      Like = GetLogLike(Trial) 
! !addon    if (Feedback > 1) write (*,*) 'Likelihood: ', Like, 'Current Like:', CurLike
!         if (Feedback > 1) write (*,*) 'Likelihood: ', Like, 'Current Like:'&
!              , CurLike, 'Rank', MPIRank
! !end addon
      accpt = (Like /= logZero) .and. (CurLike > Like .or. randexp1() > Like - CurLike) 
               !Include the min() so that compilers not doing optimal compilation don't complain
   
diff -r -N -c -b -B cosmomc/source/settings.f90 cosmomc_addon/source/settings.f90
*** cosmomc/source/settings.f90	2006-03-17 16:18:11.000000000 +0000
--- cosmomc_addon/source/settings.f90	2006-04-04 10:34:35.000000000 +0100
***************
*** 5,10 ****
--- 5,15 ----
    implicit none
  
    real :: AccuracyLevel = 1.
+ !addon
+   real :: lAccuracyLevel = 1.
+   real :: lSampleLevel = 1.
+   logical :: CambComputeEveryl = .false.
+ !end addon
    !Set to >1 to use CAMB etc on higher accuracy settings. 
    !Does not affect MCMC (except making it all slower)
  
***************
*** 43,48 ****
--- 48,56 ----
  
    integer :: num_threads = 0
    integer :: instance = 0
+ !addon
+   integer :: instance_shift = 0
+ !end addon
    integer :: MPIchains = 1, MPIrank = 0
  
    logical :: Use_LSS = .true.
