diff --git a/montepython/likelihoods/ska1_lensing/__init__.py b/montepython/likelihoods/ska1_lensing/__init__.py
index ef1cf924..6639c940 100644
--- a/montepython/likelihoods/ska1_lensing/__init__.py
+++ b/montepython/likelihoods/ska1_lensing/__init__.py
@@ -37,14 +37,14 @@ def __init__(self, path, data, command_line):
         if (self.use_halofit):
             self.need_cosmo_arguments(data, {'non linear':'halofit'})
 
-	# Warn if theoretical error and linear cutoff are requested
-	if (self.use_lmax_lincut and self.theoretical_error!=0):
-	    warnings.warn("A lmax cutoff infered from kmax and a theoretical error are requested. This combination is not implemented and in most cases not necessary as the theoretical error should induce a cutoff by itself.")
+        # Warn if theoretical error and linear cutoff are requested
+        if (self.use_lmax_lincut and self.theoretical_error!=0):
+            warnings.warn("A lmax cutoff infered from kmax and a theoretical error are requested. This combination is not implemented and in most cases not necessary as the theoretical error should induce a cutoff by itself.")
 
         # Define array of l values, and initialize them
         # It is a logspace
         # find nlmax in order to reach lmax with logarithmic steps dlnl
-        self.nlmax = np.int(np.log(self.lmax/self.lmin)/self.dlnl)+1
+        self.nlmax = int(np.log(self.lmax/self.lmin)/self.dlnl)+1
         # redefine slightly dlnl so that the last point is always exactly lmax
         self.dlnl = np.log(self.lmax/self.lmin)/(self.nlmax-1)
         self.l = self.lmin*np.exp(self.dlnl*np.arange(self.nlmax))
@@ -69,7 +69,7 @@ def __init__(self, path, data, command_line):
         # first value is already correctly set to 0.
         self.z_bin_edge = np.zeros(self.nbin+1, 'float64')
 
-	total_count = 0.
+        total_count = 0.
         for Bin in xrange(self.nbin-1):
             bin_count = 0.
             z = self.z_bin_edge[Bin]
@@ -79,23 +79,23 @@ def __init__(self, path, data, command_line):
                 bin_count += 0.5*(gd_1+gd_2)*self.dz
                 z += self.dz
             self.z_bin_edge[Bin+1] = z
-	    total_count += bin_count
+            total_count += bin_count
         self.z_bin_edge[self.nbin] = self.zmax
 
         # Fill array of discrete z values
         self.z = np.linspace(0, self.zmax, num=self.nzmax)
 
-	# check normalizations:
+        # check normalizations:
         #print("n_tot = %s" % (n_tot))
-	#int_nz = self.nzmax
-	#int_z = np.linspace(0, self.zmax, num=int_nz)
-	#DDD = np.zeros(int_nz, 'float64')
-	#for nz in xrange(self.nzmax):
-	#    #for nz2 in xrange(int_nz):
-	#    #    DDD[nz2] = self.photo_z_distribution(self.z[nz], int_z[nz2],False)
+        #int_nz = self.nzmax
+        #int_z = np.linspace(0, self.zmax, num=int_nz)
+        #DDD = np.zeros(int_nz, 'float64')
+        #for nz in xrange(self.nzmax):
+        #    #for nz2 in xrange(int_nz):
+        #    #    DDD[nz2] = self.photo_z_distribution(self.z[nz], int_z[nz2],False)
         #    DDD = self.photo_z_distribution(self.z[nz], int_z)
         #    D_tot = scipy.integrate.trapz(DDD,int_z)
-	#    print("%s\t%s\t%s" % (self.z[nz], self.galaxy_distribution(self.z[nz]), D_tot ))
+        #    print("%s\t%s\t%s" % (self.z[nz], self.galaxy_distribution(self.z[nz]), D_tot ))
 
         # Fill distribution for each bin (convolving with photo_z distribution)
         self.eta_z = np.zeros((self.nzmax, self.nbin), 'float64')
@@ -157,7 +157,7 @@ def galaxy_distribution(self, z, array=False):
         If the array flag is set to True, z is then interpretated as an array,
         and not as a single value.
 
-	Replaced by Tim Sprenger due to switch from euclid to ska
+        Replaced by Tim Sprenger due to switch from euclid to ska
         """
 
         z0 = self.par_zm/self.par_a
@@ -170,58 +170,58 @@ def galaxy_distribution(self, z, array=False):
         return galaxy_dist
 
     def photo_z_distribution(self, z, zph, array=True):
-		"""
-		Photo z distribution
-
-		If the array flag is set to True, zph is then interpretated as an array,
-		and not as a single value.
-
-		Replaced by Tim Sprenger due to switch from euclid to ska
-		"""
-
-		# Note: you must normalize it yourself to one if you want to get nice
-		# plots of the galaxy distribution function in each bin (otherwise, the
-		# spectra will remain correct, but each D_i(x) will loot strangely
-		# normalized when compared to the original D(z)
-		if not array:
-			if (z<=self.z_spec):
-				photo_z_dist =  (1.-self.f_spec)*math.exp(
-					-0.5*((z-zph)/self.sigma_phot/(1.+z))**2)/self.sigma_phot/(
-					1.+z)/math.sqrt(2.*math.pi)
-				if (z==zph):
-					photo_z_dist += self.f_spec/self.zmax*(self.nzmax-1.)
-					if (z==0.):
-						photo_z_dist += self.f_spec/self.zmax*(self.nzmax-1.)
-			elif (z<=self.z_phot):
-				photo_z_dist = math.exp(
-					-0.5*((z-zph)/self.sigma_phot/(1.+z))**2)/self.sigma_phot/(
-					1.+z)/math.sqrt(2.*math.pi)
-			else:
-				photo_z_dist = math.exp(
-					-0.5*((z-zph)/self.sigma_noz/(1.+z))**2)/self.sigma_noz/(
-					1.+z)/math.sqrt(2.*math.pi)
-		else:
-			if (z<=self.z_spec):
-				photo_z_dist =  (1.-self.f_spec)*np.exp(
-					-0.5*((z-zph)/self.sigma_phot/(1.+z))**2)/self.sigma_phot/(
-					1.+z)/math.sqrt(2.*math.pi)
-				for index_z in xrange(len(zph)):
-					if (z==zph[index_z]):
-						photo_z_dist[index_z] += self.f_spec/self.zmax*(self.nzmax-1.)
-						if (z==0.):
-							photo_z_dist[index_z] += self.f_spec/self.zmax*(self.nzmax-1.)
-						break
-			elif (z<=self.z_phot):
-				photo_z_dist = np.exp(
-					-0.5*((z-zph)/self.sigma_phot/(1.+z))**2)/self.sigma_phot/(
-					1.+z)/math.sqrt(2.*math.pi)
-			else:
-				photo_z_dist = np.exp(
-					-0.5*((z-zph)/self.sigma_noz/(1.+z))**2)/self.sigma_noz/(
-					1.+z)/math.sqrt(2.*math.pi)
-
-
-		return photo_z_dist
+                """
+                Photo z distribution
+
+                If the array flag is set to True, zph is then interpretated as an array,
+                and not as a single value.
+
+                Replaced by Tim Sprenger due to switch from euclid to ska
+                """
+
+                # Note: you must normalize it yourself to one if you want to get nice
+                # plots of the galaxy distribution function in each bin (otherwise, the
+                # spectra will remain correct, but each D_i(x) will loot strangely
+                # normalized when compared to the original D(z)
+                if not array:
+                        if (z<=self.z_spec):
+                                photo_z_dist =  (1.-self.f_spec)*math.exp(
+                                        -0.5*((z-zph)/self.sigma_phot/(1.+z))**2)/self.sigma_phot/(
+                                        1.+z)/math.sqrt(2.*math.pi)
+                                if (z==zph):
+                                        photo_z_dist += self.f_spec/self.zmax*(self.nzmax-1.)
+                                        if (z==0.):
+                                                photo_z_dist += self.f_spec/self.zmax*(self.nzmax-1.)
+                        elif (z<=self.z_phot):
+                                photo_z_dist = math.exp(
+                                        -0.5*((z-zph)/self.sigma_phot/(1.+z))**2)/self.sigma_phot/(
+                                        1.+z)/math.sqrt(2.*math.pi)
+                        else:
+                                photo_z_dist = math.exp(
+                                        -0.5*((z-zph)/self.sigma_noz/(1.+z))**2)/self.sigma_noz/(
+                                        1.+z)/math.sqrt(2.*math.pi)
+                else:
+                        if (z<=self.z_spec):
+                                photo_z_dist =  (1.-self.f_spec)*np.exp(
+                                        -0.5*((z-zph)/self.sigma_phot/(1.+z))**2)/self.sigma_phot/(
+                                        1.+z)/math.sqrt(2.*math.pi)
+                                for index_z in xrange(len(zph)):
+                                        if (z==zph[index_z]):
+                                                photo_z_dist[index_z] += self.f_spec/self.zmax*(self.nzmax-1.)
+                                                if (z==0.):
+                                                        photo_z_dist[index_z] += self.f_spec/self.zmax*(self.nzmax-1.)
+                                                break
+                        elif (z<=self.z_phot):
+                                photo_z_dist = np.exp(
+                                        -0.5*((z-zph)/self.sigma_phot/(1.+z))**2)/self.sigma_phot/(
+                                        1.+z)/math.sqrt(2.*math.pi)
+                        else:
+                                photo_z_dist = np.exp(
+                                        -0.5*((z-zph)/self.sigma_noz/(1.+z))**2)/self.sigma_noz/(
+                                        1.+z)/math.sqrt(2.*math.pi)
+
+
+                return photo_z_dist
 
     def loglkl(self, cosmo, data):
 
@@ -251,9 +251,9 @@ def loglkl(self, cosmo, data):
                     fun[1:]+fun[:-1])*(self.r[nr+1:]-self.r[nr:-1]))
                 g[nr, Bin] *= 2.*self.r[nr]*(1.+self.z[nr])
 
-	# compute the maximum l where most contributions are linear
-	# as a function of the lower bin number
-	if self.use_lmax_lincut:
+        # compute the maximum l where most contributions are linear
+        # as a function of the lower bin number
+        if self.use_lmax_lincut:
             lintegrand_lincut_o = np.zeros((self.nzmax, self.nbin, self.nbin), 'float64')
             lintegrand_lincut_u = np.zeros((self.nzmax, self.nbin, self.nbin), 'float64')
             l_lincut = np.zeros((self.nbin, self.nbin), 'float64')
@@ -276,26 +276,26 @@ def loglkl(self, cosmo, data):
                         lintegrand_lincut_u[1:, Bin1, Bin2] +
                         lintegrand_lincut_u[:-1, Bin1, Bin2])*(
                         self.r[1:]-self.r[:-1]))
-	    z_peak = np.zeros((self.nbin, self.nbin), 'float64')
+            z_peak = np.zeros((self.nbin, self.nbin), 'float64')
             for Bin1 in xrange(self.nbin):
                 for Bin2 in xrange(Bin1,self.nbin):
-	    	    z_peak[Bin1, Bin2] = self.zmax
+                    z_peak[Bin1, Bin2] = self.zmax
                     for index_z in xrange(self.nzmax):
-			if (self.r[index_z]>l_lincut[Bin1, Bin2]):
-			    z_peak[Bin1, Bin2] = self.z[index_z]
-			    break
+                        if (self.r[index_z]>l_lincut[Bin1, Bin2]):
+                            z_peak[Bin1, Bin2] = self.z[index_z]
+                            break
                     if self.use_zscaling:
-		        l_lincut[Bin1, Bin2] *= self.kmax_hMpc*cosmo.h()*pow(1.+z_peak[Bin1, Bin2],2./(2.+cosmo.n_s()))
+                        l_lincut[Bin1, Bin2] *= self.kmax_hMpc*cosmo.h()*pow(1.+z_peak[Bin1, Bin2],2./(2.+cosmo.n_s()))
                     else:
-		        l_lincut[Bin1, Bin2] *= self.kmax_hMpc*cosmo.h()
-		l_lincut_mean[Bin1] = np.sum(l_lincut[Bin1, :])/(self.nbin-Bin1)
+                        l_lincut[Bin1, Bin2] *= self.kmax_hMpc*cosmo.h()
+                l_lincut_mean[Bin1] = np.sum(l_lincut[Bin1, :])/(self.nbin-Bin1)
 
-	    #for Bin1 in xrange(self.nbin):
-	        #for Bin2 in xrange(Bin1,self.nbin):
-		    #print("%s\t%s\t%s\t%s" % (Bin1, Bin2, l_lincut[Bin1, Bin2], l_lincut_mean[Bin1]))
+            #for Bin1 in xrange(self.nbin):
+                #for Bin2 in xrange(Bin1,self.nbin):
+                    #print("%s\t%s\t%s\t%s" % (Bin1, Bin2, l_lincut[Bin1, Bin2], l_lincut_mean[Bin1]))
 
-	#for nr in xrange(1, self.nzmax-1):
-	#	print("%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s" % (self.z[nr], g[nr, 0], g[nr, 1], g[nr, 2], g[nr, 3], g[nr, 4], g[nr, 5], g[nr, 6], g[nr, 7], g[nr, 8], g[nr, 9]))
+        #for nr in xrange(1, self.nzmax-1):
+        #        print("%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s" % (self.z[nr], g[nr, 0], g[nr, 1], g[nr, 2], g[nr, 3], g[nr, 4], g[nr, 5], g[nr, 6], g[nr, 7], g[nr, 8], g[nr, 9]))
 
         # Get power spectrum P(k=l/r,z(r)) from cosmological module
         kmin_in_inv_Mpc = self.k_min_h_by_Mpc * cosmo.h()
@@ -318,7 +318,7 @@ def loglkl(self, cosmo, data):
                 else:
                     pk[index_l, index_z] = cosmo.pk(self.l[index_l]/self.r[index_z], self.z[index_z])
 
-		#print("%s\t%s\t%s" %(self.l[index_l], self.z[index_z], pk[index_l, index_z]))
+                #print("%s\t%s\t%s" %(self.l[index_l], self.z[index_z], pk[index_l, index_z]))
 
         # Recover the non_linear scale computed by halofit. If no scale was
         # affected, set the scale to one, and make sure that the nuisance
@@ -329,15 +329,15 @@ def loglkl(self, cosmo, data):
         else:
             k_sigma = cosmo.nonlinear_scale(self.z, self.nzmax)
 
-	if not (cosmo.nonlinear_method == 0):
-	    k_sigma_problem = False
-	    for index_z in xrange(self.nzmax-1):
-	        if (k_sigma[index_z+1]<k_sigma[index_z]) or (k_sigma[index_z+1]>2.5):
-	            k_sigma[index_z+1] = 2.5
-	            k_sigma_problem = True
-	        #print("%s\t%s" % (k_sigma[index_z], self.z[index_z]))
-	    if k_sigma_problem:
-	        warnings.warn("There were unphysical (decreasing in redshift or exploding) values of k_sigma (=cosmo.nonlinear_scale(...)). To proceed they were set to 2.5, the highest scale that seems to be stable.")
+        if not (cosmo.nonlinear_method == 0):
+            k_sigma_problem = False
+            for index_z in xrange(self.nzmax-1):
+                if (k_sigma[index_z+1]<k_sigma[index_z]) or (k_sigma[index_z+1]>2.5):
+                    k_sigma[index_z+1] = 2.5
+                    k_sigma_problem = True
+                #print("%s\t%s" % (k_sigma[index_z], self.z[index_z]))
+            if k_sigma_problem:
+                warnings.warn("There were unphysical (decreasing in redshift or exploding) values of k_sigma (=cosmo.nonlinear_scale(...)). To proceed they were set to 2.5, the highest scale that seems to be stable.")
 
         # Define the alpha function, that will characterize the theoretical
         # uncertainty. Chosen to be 0.001 at low k, raise between 0.1 and 0.2
@@ -504,7 +504,7 @@ def loglkl(self, cosmo, data):
         # TODO parallelize this
         for index, ell in enumerate(ells):
 
-	    if self.use_lmax_lincut:
+            if self.use_lmax_lincut:
                 CutBin = -1
                 for zBin in xrange(self.nbin):
                     if (ell<l_lincut_mean[zBin]):
@@ -514,7 +514,7 @@ def loglkl(self, cosmo, data):
                         break
                 if (CutBin==-1):
                     break
-	    else:
+            else:
                 det_theory = np.linalg.det(Cov_theory[index, :, :])
                 det_observ = np.linalg.det(Cov_observ[index, :, :])
 
@@ -549,14 +549,14 @@ def loglkl(self, cosmo, data):
                             newCov[:, i] = Cov_observ[index,:, i]#MArchi#newCov[:, i] = Cov_observ[:, i]
                             det_theory_plus_error_cross_obs += np.linalg.det(
                                 newCov)
-			try:
-                        	function_vector[k] = (2.*ell+1.)*self.fsky*(det_theory_plus_error_cross_obs/det_theory_plus_error + math.log(det_theory_plus_error/det_observ) - self.nbin ) + dof*vector[k]**2
-			except ValueError:
-				warnings.warn("ska1_lensing: Could not evaluate chi2 including theoretical error with the current parameters. The corresponding chi2 is now set to nan!")
-				break
-				break
-				break
-				chi2 = np.nan
+                        try:
+                                function_vector[k] = (2.*ell+1.)*self.fsky*(det_theory_plus_error_cross_obs/det_theory_plus_error + math.log(det_theory_plus_error/det_observ) - self.nbin ) + dof*vector[k]**2
+                        except ValueError:
+                                warnings.warn("ska1_lensing: Could not evaluate chi2 including theoretical error with the current parameters. The corresponding chi2 is now set to nan!")
+                                break
+                                break
+                                break
+                                chi2 = np.nan
 
 
                     # Computing first
@@ -582,20 +582,20 @@ def loglkl(self, cosmo, data):
 
 
             else:
-		if self.use_lmax_lincut:
+                if self.use_lmax_lincut:
                     det_cross = 0.
                     for i in xrange(0,self.nbin-CutBin):
                         newCov = np.copy(Cov_theory[index, CutBin:, CutBin:])
                         newCov[:, i] = Cov_observ[index, CutBin:, CutBin+i]
                         det_cross += np.linalg.det(newCov)
-		else:
+                else:
                     det_cross = 0.
                     for i in xrange(self.nbin):
                         newCov = np.copy(Cov_theory[index, :, :])
                         newCov[:, i] = Cov_observ[index, :, i]
                         det_cross += np.linalg.det(newCov)
 
-		if self.use_lmax_lincut:
+                if self.use_lmax_lincut:
                     chi2 += (2.*ell+1.)*self.fsky*(det_cross/det_theory + math.log(det_theory/det_observ) - self.nbin+CutBin)
                 else:
                     chi2 += (2.*ell+1.)*self.fsky*(det_cross/det_theory + math.log(det_theory/det_observ) - self.nbin)
diff --git a/montepython/likelihoods/ska1_pk/__init__.py b/montepython/likelihoods/ska1_pk/__init__.py
index c7d5082f..15c204ff 100644
--- a/montepython/likelihoods/ska1_pk/__init__.py
+++ b/montepython/likelihoods/ska1_pk/__init__.py
@@ -4,6 +4,7 @@
 import warnings
 from numpy import newaxis as na
 from math import exp, log, pi, log10
+from montepython import io_mp
 import scipy.integrate
 # Created by Tim Sprenger in 2017 based on euclid_pk
 
@@ -23,17 +24,17 @@ def __init__(self, path, data, command_line):
         self.need_cosmo_arguments(data, {'P_k_max_1/Mpc': 1.5*self.k_cut(self.zmax)})
 
         # Compute non-linear power spectrum if requested
-	if self.use_halofit:
-            	self.need_cosmo_arguments(data, {'non linear':'halofit'})
-		print("Using halofit")
+        if self.use_halofit:
+            self.need_cosmo_arguments(data, {'non linear':'halofit'})
+            print("Using halofit")
 
         # Deduce the dz step from the number of bins and the edge values of z
         self.dz = (self.zmax-self.zmin)/self.nbin
 
-	# Compute new zmin and zmax which are bin centers
-	# Need to be defined as edges if zmin can be close to z=0
-	self.zmin += self.dz/2.
-	self.zmax -= self.dz/2.
+        # Compute new zmin and zmax which are bin centers
+        # Need to be defined as edges if zmin can be close to z=0
+        self.zmin += self.dz/2.
+        self.zmax -= self.dz/2.
 
         # self.z_mean will contain the central values
         self.z_mean = np.linspace(self.zmin, self.zmax, num=self.nbin)
@@ -45,21 +46,21 @@ def __init__(self, path, data, command_line):
 
         # Compute the number of galaxies for each \bar z
         # N_g(\bar z) = int_{\bar z - dz/2}^{\bar z + dz/2} dn/dz dz
-	# dn/dz fit formula from 1412.4700v2: 10^c1*z^c2*e^-c3z
-	# dn/dz = number of galaxies per redshift and deg^2
-	# old approx.: self.N_g = pow(10.,self.c1)*pow(self.z_mean,self.c2)*np.exp(-self.c3*self.z_mean) * self.dz * self.skycov
-	self.N_g = np.zeros(self.nbin)
-	N_tot = 0.0
-	for index_z in xrange(self.nbin):
-		self.N_g[index_z], error = scipy.integrate.quad(self.dndz, self.z_mean[index_z]-self.dz/2., self.z_mean[index_z]+self.dz/2.)
-		assert error/self.N_g[index_z] <= 0.001, ("dndz integration error is bigger than 0.1%")
-		N_tot += self.N_g[index_z]
-
-	# Ntot output
-	#print("\nSKA1: Number of detected galaxies and bias in each redshift bin:")
-	#for index_z in xrange(self.nbin):
-	#	print("z-bin[" + str(self.z_mean[index_z]-self.dz/2.) + "," + str(self.z_mean[index_z]+self.dz/2.) + "]: \tN = %.4g" % (self.N_g[index_z]) + " ,\t b = %.4g" % (b[index_z]))
-	#print("Total number of detected galaxies: N = %.4g\n" % (N_tot))
+        # dn/dz fit formula from 1412.4700v2: 10^c1*z^c2*e^-c3z
+        # dn/dz = number of galaxies per redshift and deg^2
+        # old approx.: self.N_g = pow(10.,self.c1)*pow(self.z_mean,self.c2)*np.exp(-self.c3*self.z_mean) * self.dz * self.skycov
+        self.N_g = np.zeros(self.nbin)
+        N_tot = 0.0
+        for index_z in xrange(self.nbin):
+                self.N_g[index_z], error = scipy.integrate.quad(self.dndz, self.z_mean[index_z]-self.dz/2., self.z_mean[index_z]+self.dz/2.)
+                assert error/self.N_g[index_z] <= 0.001, ("dndz integration error is bigger than 0.1%")
+                N_tot += self.N_g[index_z]
+
+        # Ntot output
+        #print("\nSKA1: Number of detected galaxies and bias in each redshift bin:")
+        #for index_z in xrange(self.nbin):
+        #        print("z-bin[" + str(self.z_mean[index_z]-self.dz/2.) + "," + str(self.z_mean[index_z]+self.dz/2.) + "]: \tN = %.4g" % (self.N_g[index_z]) + " ,\t b = %.4g" % (b[index_z]))
+        #print("Total number of detected galaxies: N = %.4g\n" % (N_tot))
 
         # Define the k values for the integration (from kmin to kmax), at which
         # the spectrum will be computed (and stored for the fiducial model)
@@ -72,7 +73,7 @@ def __init__(self, path, data, command_line):
         # If the file exists, initialize the fiducial values, the spectrum will
         # be read first, with k_size values of k and nbin values of z. Then,
         # H_fid and D_A fid will be read (each with nbin values).
-	# Then V_fid, b_fid and the fiducial errors on real space coordinates follow.
+        # Then V_fid, b_fid and the fiducial errors on real space coordinates follow.
         self.fid_values_exist = False
         self.pk_nl_fid = np.zeros((self.k_size, 2*self.nbin+1), 'float64')
         if self.use_linear_rsd:
@@ -81,8 +82,8 @@ def __init__(self, path, data, command_line):
         self.D_A_fid = np.zeros(2*self.nbin+1, 'float64')
         self.V_fid = np.zeros(self.nbin, 'float64')
         self.b_fid = np.zeros(self.nbin, 'float64')
-	self.sigma_A_fid = np.zeros(self.nbin, 'float64')
-	self.sigma_B_fid = np.zeros(self.nbin, 'float64')
+        self.sigma_A_fid = np.zeros(self.nbin, 'float64')
+        self.sigma_B_fid = np.zeros(self.nbin, 'float64')
 
         fid_file_path = os.path.join(self.data_directory, self.fiducial_file)
         if os.path.exists(fid_file_path):
@@ -113,22 +114,22 @@ def __init__(self, path, data, command_line):
                     self.sigma_A_fid[index_z] = float(line.split()[0])
                     self.sigma_B_fid[index_z] = float(line.split()[1])
                     line = fid_file.readline()
-		self.sigma_NL_fid = float(line)
+                self.sigma_NL_fid = float(line)
 
         # Else the file will be created in the loglkl() function.
         return
 
     def dndz(self, z):
-	# dn/dz/deg^2 fit formula from 1412.4700v2: 10^c1*z^c2*e^-c3z * sky coverage in deg^2
-	dndz = pow(10.,self.c1)*pow(z,self.c2)*np.exp(-self.c3*z)*self.skycov
-	return dndz
+        # dn/dz/deg^2 fit formula from 1412.4700v2: 10^c1*z^c2*e^-c3z * sky coverage in deg^2
+        dndz = pow(10.,self.c1)*pow(z,self.c2)*np.exp(-self.c3*z)*self.skycov
+        return dndz
 
     def k_cut(self, z,h=0.6693,n_s=0.9619):
-	kcut = self.kmax*h
-	# compute kmax according to highest redshift linear cutoff (1509.07562v2)
-	if self.use_zscaling:
-		kcut *= pow(1.+z,2./(2.+n_s))
-	return kcut
+        kcut = self.kmax*h
+        # compute kmax according to highest redshift linear cutoff (1509.07562v2)
+        if self.use_zscaling:
+                kcut *= pow(1.+z,2./(2.+n_s))
+        return kcut
 
     def loglkl(self, cosmo, data):
 
@@ -145,34 +146,34 @@ def loglkl(self, cosmo, data):
 
         # Compute V_survey, for each given redshift bin, which is the volume of
         # a shell times the sky coverage (only fiducial needed):
-	if self.fid_values_exist is False:
-        	V_survey = np.zeros(self.nbin, 'float64')
-        	for index_z in xrange(self.nbin):
-           	 	V_survey[index_z] = 4./3.*pi*self.skycov/41253.*(
-               			r[2*index_z+2]**3-r[2*index_z]**3)
+        if self.fid_values_exist is False:
+                V_survey = np.zeros(self.nbin, 'float64')
+                for index_z in xrange(self.nbin):
+                            V_survey[index_z] = 4./3.*pi*self.skycov/41253.*(
+                                       r[2*index_z+2]**3-r[2*index_z]**3)
 
         # At the center of each bin, compute the bias function,
         # using fit fromula from 1412.4700v2: c4*e^c5z
-	if 'beta_0^SKA1' in self.use_nuisance:
-        	b = self.c4*np.exp(self.c5*self.z_mean*data.mcmc_parameters['beta_1^SKA1']['current']*data.mcmc_parameters['beta_1^SKA1']['scale'])*data.mcmc_parameters['beta_0^SKA1']['current']*data.mcmc_parameters['beta_0^SKA1']['scale']
-	else:
-        	b = self.c4*np.exp(self.c5*self.z_mean)
-
-    	# Compute freq.res. sigma_r = (1+z)^2/H*sigma_nu/nu_21cm, nu in Mhz
-	# Compute ang.res. sigma_perp = 1.22(1+z)^2*D_A*lambda_21cm/Baseline, Baseline in km
-	# combine into exp(-k^2*(mu^2*(sig_r^2-sig_perp^2)+sig_perp^2)) independent of cosmo
-	# used as exp(-k^2*(mu^2*sigma_A+sigma_B)) all fiducial
-	if self.fid_values_exist is False:
-		sigma_A = np.zeros(self.nbin,'float64')
-		sigma_B = np.zeros(self.nbin,'float64')
-		sigma_A = ((1.+self.z_mean[:])**2/H[1::2]*self.delta_nu/np.sqrt(8.*np.log(2.))/self.nu0)**2 -(
-			1./np.sqrt(8.*np.log(2.))*(1+self.z_mean[:])**2 * D_A[1::2]*2.111e-4/self.Baseline)**2
-		sigma_B = (1./np.sqrt(8.*np.log(2.))*(1+self.z_mean[:])**2 * D_A[1::2]*2.111e-4/self.Baseline)**2
-
-	# sigma_NL in Mpc = nonlinear dispersion scale of RSD (1405.1452v2)
-	sigma_NL = 0.0	# fiducial would be 7 but when kept constant that is more constraining than keeping 0
-	if 'sigma_NL' in self.use_nuisance:
-		sigma_NL = data.mcmc_parameters['sigma_NL']['current']*data.mcmc_parameters['sigma_NL']['scale']
+        if 'beta_0^SKA1' in self.use_nuisance:
+                b = self.c4*np.exp(self.c5*self.z_mean*data.mcmc_parameters['beta_1^SKA1']['current']*data.mcmc_parameters['beta_1^SKA1']['scale'])*data.mcmc_parameters['beta_0^SKA1']['current']*data.mcmc_parameters['beta_0^SKA1']['scale']
+        else:
+                b = self.c4*np.exp(self.c5*self.z_mean)
+
+            # Compute freq.res. sigma_r = (1+z)^2/H*sigma_nu/nu_21cm, nu in Mhz
+        # Compute ang.res. sigma_perp = 1.22(1+z)^2*D_A*lambda_21cm/Baseline, Baseline in km
+        # combine into exp(-k^2*(mu^2*(sig_r^2-sig_perp^2)+sig_perp^2)) independent of cosmo
+        # used as exp(-k^2*(mu^2*sigma_A+sigma_B)) all fiducial
+        if self.fid_values_exist is False:
+                sigma_A = np.zeros(self.nbin,'float64')
+                sigma_B = np.zeros(self.nbin,'float64')
+                sigma_A = ((1.+self.z_mean[:])**2/H[1::2]*self.delta_nu/np.sqrt(8.*np.log(2.))/self.nu0)**2 -(
+                        1./np.sqrt(8.*np.log(2.))*(1+self.z_mean[:])**2 * D_A[1::2]*2.111e-4/self.Baseline)**2
+                sigma_B = (1./np.sqrt(8.*np.log(2.))*(1+self.z_mean[:])**2 * D_A[1::2]*2.111e-4/self.Baseline)**2
+
+        # sigma_NL in Mpc = nonlinear dispersion scale of RSD (1405.1452v2)
+        sigma_NL = 0.0        # fiducial would be 7 but when kept constant that is more constraining than keeping 0
+        if 'sigma_NL' in self.use_nuisance:
+                sigma_NL = data.mcmc_parameters['sigma_NL']['current']*data.mcmc_parameters['sigma_NL']['scale']
 
         # If the fiducial model does not exist, recover the power spectrum and
         # store it, then exit.
@@ -203,8 +204,8 @@ def loglkl(self, cosmo, data):
                 for index_z in xrange(self.nbin):
                     fid_file.write('%.8g %.8g\n' % (V_survey[index_z],b[index_z]))
                 for index_z in xrange(self.nbin):
-			fid_file.write('%.8g %.8g\n' % (sigma_A[index_z], sigma_B[index_z]))
-		fid_file.write('%.8g\n' % sigma_NL)
+                        fid_file.write('%.8g %.8g\n' % (sigma_A[index_z], sigma_B[index_z]))
+                fid_file.write('%.8g\n' % sigma_NL)
             print('\n')
             warnings.warn(
                 "Writing fiducial model in %s, for %s likelihood\n" % (
@@ -220,20 +221,20 @@ def loglkl(self, cosmo, data):
         # Compute the beta_fid function, for observed spectrum,
         # beta_fid(k_fid,z) = 1/2b(z) * d log(P_nl_fid(k_fid,z))/d log a
         #                   = -1/2b(z)* (1+z) d log(P_nl_fid(k_fid,z))/dz
-	if self.use_linear_rsd:
+        if self.use_linear_rsd:
             beta_fid = -0.5/self.b_fid*(1+self.z_mean)*np.log(
-            	self.pk_lin_fid[:, 2::2]/self.pk_lin_fid[:, :-2:2])/self.dz
-	else:
+                    self.pk_lin_fid[:, 2::2]/self.pk_lin_fid[:, :-2:2])/self.dz
+        else:
             beta_fid = -0.5/self.b_fid*(1+self.z_mean)*np.log(
-            	self.pk_nl_fid[:, 2::2]/self.pk_nl_fid[:, :-2:2])/self.dz
+                    self.pk_nl_fid[:, 2::2]/self.pk_nl_fid[:, :-2:2])/self.dz
 
         # Compute the tilde P_fid(k_ref,z,mu) = H_fid(z)/D_A_fid(z)**2 ( 1 + beta_fid(k_fid,z)mu^2)^2 P_nl_fid(k_fid,z)exp(-k_fid^2*(mu_fid^2*sigma_A(z)+sigma_B(z)))
         self.tilde_P_fid = np.zeros((self.k_size, self.nbin, self.mu_size),'float64')
         self.tilde_P_fid = self.H_fid[na, 1::2, na]/(
-            	self.D_A_fid[na, 1::2, na])**2*self.b_fid[na,:,na]**2*(
+                    self.D_A_fid[na, 1::2, na])**2*self.b_fid[na,:,na]**2*(
                 1. + beta_fid[:, :, na] * self.mu_fid[na, na, :]**2)**2 * (
-            	self.pk_nl_fid[:, 1::2, na]) * np.exp(-self.k_fid[:,na,na]**2 *
-		(self.mu_fid[na, na, :]**2*(self.sigma_A_fid[na,:,na]+self.sigma_NL_fid**2) + self.sigma_B_fid[na,:,na]))
+                    self.pk_nl_fid[:, 1::2, na]) * np.exp(-self.k_fid[:,na,na]**2 *
+                (self.mu_fid[na, na, :]**2*(self.sigma_A_fid[na,:,na]+self.sigma_NL_fid**2) + self.sigma_B_fid[na,:,na]))
 
         ######################
         # TH PART
@@ -245,11 +246,11 @@ def loglkl(self, cosmo, data):
             for index_z in xrange(2*self.nbin+1):
                 self.k[index_k,index_z,:] = np.sqrt((1.-self.mu_fid[:]**2)*self.D_A_fid[index_z]**2/D_A[index_z]**2 + self.mu_fid[:]**2*H[index_z]**2/self.H_fid[index_z]**2 )*self.k_fid[index_k]
 
-	# Compute values of mu based on fiducial values:
-	# mu^2 = mu_fid^2 / (mu_fid^2 + ((H_fid*D_A_fid)/(H*D_A))^2)*(1 - mu_fid^2))
-	self.mu = np.zeros((self.nbin,self.mu_size),'float64')
-	for index_z in xrange(self.nbin):
-		self.mu[index_z,:] = np.sqrt(self.mu_fid[:]**2/(self.mu_fid[:]**2 + ((self.H_fid[2*index_z+1]*self.D_A_fid[2*index_z+1])/(D_A[2*index_z+1]*H[2*index_z+1]))**2 * (1.-self.mu_fid[:]**2)))
+        # Compute values of mu based on fiducial values:
+        # mu^2 = mu_fid^2 / (mu_fid^2 + ((H_fid*D_A_fid)/(H*D_A))^2)*(1 - mu_fid^2))
+        self.mu = np.zeros((self.nbin,self.mu_size),'float64')
+        for index_z in xrange(self.nbin):
+                self.mu[index_z,:] = np.sqrt(self.mu_fid[:]**2/(self.mu_fid[:]**2 + ((self.H_fid[2*index_z+1]*self.D_A_fid[2*index_z+1])/(D_A[2*index_z+1]*H[2*index_z+1]))**2 * (1.-self.mu_fid[:]**2)))
 
         # Recover the non-linear power spectrum from the cosmological module on all
         # the z_boundaries, to compute afterwards beta. This is pk_nl_th from the
@@ -266,31 +267,31 @@ def loglkl(self, cosmo, data):
         if self.use_linear_rsd:
             pk_lin_th = cosmo.get_pk_cb_lin(self.k,self.z,self.k_size,2*self.nbin+1,self.mu_size)
 
-	if self.UseTheoError :
-        	# Recover the non_linear scale computed by halofit.
-        	#self.k_sigma = np.zeros(2*self.nbin+1, 'float64')
-            	#self.k_sigma = cosmo.nonlinear_scale(self.z,2*self.nbin+1)
-
-        	# Define the theoretical error envelope
-        	self.alpha = np.zeros((self.k_size,self.nbin,self.mu_size),'float64')
-		th_c1 = 0.75056
-		th_c2 = 1.5120
-		th_a1 = 0.014806
-		th_a2 = 0.022047
-       		for index_z in xrange(self.nbin):
-		    k_z = cosmo.h()*pow(1.+self.z_mean[index_z],2./(2.+cosmo.n_s()))
-		    for index_mu in xrange(self.mu_size):
-		        for index_k in xrange(self.k_size):
-		            if self.k[index_k,2*index_z+1,index_mu]/k_z<0.3:
-	 	                self.alpha[index_k,index_z,index_mu] = th_a1*np.exp(th_c1*np.log10(self.k[index_k,2*index_z+1,index_mu]/k_z))
-		            else:
-		                self.alpha[index_k,index_z,index_mu] = th_a2*np.exp(th_c2*np.log10(self.k[index_k,2*index_z+1,index_mu]/k_z))
-
-		# Define fractional theoretical error variance R/P^2
-		self.R_var = np.zeros((self.k_size,self.nbin,self.mu_size),'float64')
-		for index_k in xrange(self.k_size):
-	    	    for index_z in xrange(self.nbin):
-	                self.R_var[index_k,index_z,:] = self.V_fid[index_z]/(2.*np.pi)**2*self.k_CorrLength_hMpc*cosmo.h()/self.z_CorrLength*self.dz*self.k_fid[index_k]**2*self.alpha[index_k,index_z,:]**2
+        if self.UseTheoError :
+                # Recover the non_linear scale computed by halofit.
+                #self.k_sigma = np.zeros(2*self.nbin+1, 'float64')
+                    #self.k_sigma = cosmo.nonlinear_scale(self.z,2*self.nbin+1)
+
+                # Define the theoretical error envelope
+                self.alpha = np.zeros((self.k_size,self.nbin,self.mu_size),'float64')
+                th_c1 = 0.75056
+                th_c2 = 1.5120
+                th_a1 = 0.014806
+                th_a2 = 0.022047
+                for index_z in xrange(self.nbin):
+                    k_z = cosmo.h()*pow(1.+self.z_mean[index_z],2./(2.+cosmo.n_s()))
+                    for index_mu in xrange(self.mu_size):
+                        for index_k in xrange(self.k_size):
+                            if self.k[index_k,2*index_z+1,index_mu]/k_z<0.3:
+                                 self.alpha[index_k,index_z,index_mu] = th_a1*np.exp(th_c1*np.log10(self.k[index_k,2*index_z+1,index_mu]/k_z))
+                            else:
+                                self.alpha[index_k,index_z,index_mu] = th_a2*np.exp(th_c2*np.log10(self.k[index_k,2*index_z+1,index_mu]/k_z))
+
+                # Define fractional theoretical error variance R/P^2
+                self.R_var = np.zeros((self.k_size,self.nbin,self.mu_size),'float64')
+                for index_k in xrange(self.k_size):
+                    for index_z in xrange(self.nbin):
+                        self.R_var[index_k,index_z,:] = self.V_fid[index_z]/(2.*np.pi)**2*self.k_CorrLength_hMpc*cosmo.h()/self.z_CorrLength*self.dz*self.k_fid[index_k]**2*self.alpha[index_k,index_z,:]**2
 
         # Compute the beta function for nl,
         # beta(k,z) = 1/2b(z) * d log(P_nl_th (k,z))/d log a
@@ -315,47 +316,47 @@ def loglkl(self, cosmo, data):
             self.P_shot[index_z] = self.H_fid[2*index_z+1]/(self.D_A_fid[2*index_z+1]**2)*self.V_fid[index_z]/self.N_g[index_z]
 
         # finally compute chi2, for each z_mean
-	if self.use_zscaling==0:
-		# redshift dependent cutoff makes integration more complicated
-        	chi2 = 0.0
-		index_kmax = 0
-		delta_mu = self.mu_fid[1] - self.mu_fid[0] # equally spaced
-		integrand_low = 0.0
-		integrand_hi = 0.0
-
-		for index_z in xrange(self.nbin):
-			# uncomment printers to get contributions from individual redshift bins
-			#printer1 = chi2*delta_mu
-			# uncomment to display max. kmin (used to infer kmin~0.02):
-			#kmin: #print("z=" + str(self.z_mean[index_z]) + " kmin=" + str(34.56/r[2*index_z+1]) + "\tor " + str(6.283/(r[2*index_z+2]-r[2*index_z])))
-			for index_k in xrange(1,self.k_size):
-				if ((self.k_cut(self.z_mean[index_z],cosmo.h(),cosmo.n_s())-self.k_fid[self.k_size-index_k]) > -1.e-6):
-					index_kmax = self.k_size-index_k
-					break
-			integrand_low = self.integrand(0,index_z,0)*.5
-			for index_k in xrange(1,index_kmax+1):
-				integrand_hi = self.integrand(index_k,index_z,0)*.5
-				chi2 += (integrand_hi+integrand_low)*.5*(self.k_fid[index_k]-self.k_fid[index_k-1])
-				integrand_low = integrand_hi
-			chi2 += integrand_low*(self.k_cut(self.z_mean[index_z],cosmo.h(),cosmo.n_s())-self.k_fid[index_kmax])
-			for index_mu in xrange(1,self.mu_size-1):
-				integrand_low = self.integrand(0,index_z,index_mu)
-				for index_k in xrange(1,index_kmax+1):
-					integrand_hi = self.integrand(index_k,index_z,index_mu)
-					chi2 += (integrand_hi+integrand_low)*.5*(self.k_fid[index_k]-self.k_fid[index_k-1])
-					integrand_low = integrand_hi
-				chi2 += integrand_low*(self.k_cut(self.z_mean[index_z],cosmo.h(),cosmo.n_s())-self.k_fid[index_kmax])
-			integrand_low = self.integrand(0,index_z,self.mu_size-1)*.5
-			for index_k in xrange(1,index_kmax+1):
-				integrand_hi = self.integrand(index_k,index_z,self.mu_size-1)*.5
-				chi2 += (integrand_hi+integrand_low)*.5*(self.k_fid[index_k]-self.k_fid[index_k-1])
-				integrand_low = integrand_hi
-			chi2 += integrand_low*(self.k_cut(self.z_mean[index_z],cosmo.h(),cosmo.n_s())-self.k_fid[index_kmax])
-			#printer2 = chi2*delta_mu-printer1
-			#print("%s\t%s" % (self.z_mean[index_z], printer2))
-		chi2 *= delta_mu
-
-	else:
+        if self.use_zscaling==0:
+                # redshift dependent cutoff makes integration more complicated
+                chi2 = 0.0
+                index_kmax = 0
+                delta_mu = self.mu_fid[1] - self.mu_fid[0] # equally spaced
+                integrand_low = 0.0
+                integrand_hi = 0.0
+
+                for index_z in xrange(self.nbin):
+                        # uncomment printers to get contributions from individual redshift bins
+                        #printer1 = chi2*delta_mu
+                        # uncomment to display max. kmin (used to infer kmin~0.02):
+                        #kmin: #print("z=" + str(self.z_mean[index_z]) + " kmin=" + str(34.56/r[2*index_z+1]) + "\tor " + str(6.283/(r[2*index_z+2]-r[2*index_z])))
+                        for index_k in xrange(1,self.k_size):
+                                if ((self.k_cut(self.z_mean[index_z],cosmo.h(),cosmo.n_s())-self.k_fid[self.k_size-index_k]) > -1.e-6):
+                                        index_kmax = self.k_size-index_k
+                                        break
+                        integrand_low = self.integrand(0,index_z,0)*.5
+                        for index_k in xrange(1,index_kmax+1):
+                                integrand_hi = self.integrand(index_k,index_z,0)*.5
+                                chi2 += (integrand_hi+integrand_low)*.5*(self.k_fid[index_k]-self.k_fid[index_k-1])
+                                integrand_low = integrand_hi
+                        chi2 += integrand_low*(self.k_cut(self.z_mean[index_z],cosmo.h(),cosmo.n_s())-self.k_fid[index_kmax])
+                        for index_mu in xrange(1,self.mu_size-1):
+                                integrand_low = self.integrand(0,index_z,index_mu)
+                                for index_k in xrange(1,index_kmax+1):
+                                        integrand_hi = self.integrand(index_k,index_z,index_mu)
+                                        chi2 += (integrand_hi+integrand_low)*.5*(self.k_fid[index_k]-self.k_fid[index_k-1])
+                                        integrand_low = integrand_hi
+                                chi2 += integrand_low*(self.k_cut(self.z_mean[index_z],cosmo.h(),cosmo.n_s())-self.k_fid[index_kmax])
+                        integrand_low = self.integrand(0,index_z,self.mu_size-1)*.5
+                        for index_k in xrange(1,index_kmax+1):
+                                integrand_hi = self.integrand(index_k,index_z,self.mu_size-1)*.5
+                                chi2 += (integrand_hi+integrand_low)*.5*(self.k_fid[index_k]-self.k_fid[index_k-1])
+                                integrand_low = integrand_hi
+                        chi2 += integrand_low*(self.k_cut(self.z_mean[index_z],cosmo.h(),cosmo.n_s())-self.k_fid[index_kmax])
+                        #printer2 = chi2*delta_mu-printer1
+                        #print("%s\t%s" % (self.z_mean[index_z], printer2))
+                chi2 *= delta_mu
+
+        else:
             chi2 = 0.0
             mu_integrand_lo,mu_integrand_hi = 0.0,0.0
             k_integrand  = np.zeros(self.k_size,'float64')