diff --git a/pbhmergers.py b/pbhmergers.py
index e9617d7..e017bb9 100644
--- a/pbhmergers.py
+++ b/pbhmergers.py
@@ -63,7 +63,9 @@ def rhocrit(self):
         hubz2 = (self.overden.omega_matter0/aa**3 + self.overden.omega_lambda0) * hubble**2
         #Critical density at redshift in units of kg m^-3
         rhocrit = 3 * hubz2 / (8*math.pi* self.ureg.newtonian_constant_of_gravitation)
-        print "rhocrit = ", rhocrit
+        # print "rhocrit = ", rhocrit
+        # THIS WAS A BUG! THE CORRECT LINE IS:
+        print("rhocrit = ", rhocrit)
         return rhocrit.to_base_units()
 
     def R200(self, mass):
@@ -116,7 +118,9 @@ def cross_section(self, mass):
         gammaint = sigma**(10/7)*scipy.special.gammainc(5/7,vbysig**2)* scipy.special.gamma(5/7)/2
         #We also need to normalise the probability function for v:
         #Integrate[4*Pi*v^2*P[v, sigma, Vvir], {v, 0, Vvir}]
-        probnorm = math.pi**(3/2)*sigma**3*scipy.special.erf(vbysig) - 2*math.pi/3*np.exp(-(vvir**2/sigma**2))*(3*sigma**2*vvir + 2*vvir**3) *0
+        # probnorm = math.pi**(3/2)*sigma**3*scipy.special.erf(vbysig) - 2*math.pi/3*np.exp(-(vvir**2/sigma**2))*(3*sigma**2*vvir + 2*vvir**3) *0
+        # THIS WAS A BUG! THE CORRECT LINE IS:
+        probnorm = math.pi**(3/2)*sigma**3*scipy.special.erf(vbysig) - 2*math.pi/3*np.exp(-(vvir**2/sigma**2))*(3*sigma**2*vvir + 2*vvir**3)
         assert np.all(probnorm.magnitude > 0)
         cross_section = prefac*(gammaint + cutoff)/probnorm
         assert self.ureg.get_dimensionality('[length]**3 [time]**(-1) [mass]**(-2)') == self.ureg.get_dimensionality(cross_section)
@@ -224,7 +228,9 @@ def mergerhalflife(self,mass,threefac=True, bhmass=None):
             threefac = self.threebodyratio(mass)
             threefac = np.max([threefac, np.ones_like(threefac)],axis=0)
             rate *= threefac
-        return 0.5*(mass/bhmass)/rat
+        # return 0.5*(mass/bhmass)/rat
+        # THIS WAS A BUG! THE CORRECT LINE IS:
+        return 0.5*(mass/bhmass)/rate
 
     def bias(self,mass):
         """The formula for halo bias in EPS theory (Mo & White 1996), eq. 13"""
diff --git a/power_specs.py b/power_specs.py
index 1e00346..06bffff 100644
--- a/power_specs.py
+++ b/power_specs.py
@@ -1,9 +1,7 @@
 #!/usr/bin/env python
 # vim: set fileencoding=UTF-8 :
 
-"""
-Various flux derivative stuff
-"""
+"""Various flux derivative stuff."""
 
 import numpy as np
 import math
@@ -14,658 +12,663 @@
 import re
 
 class DataError(Exception):
+    """Do some thing."""
+    
     def __init__(self, value):
+        """Define some value."""
         self.value = value
     def __str__(self):
+        """Define some value."""
         return repr(self.value)
 
-""" A where function to find where a floating point value is equal to another"""
 def wheref(array, value):
-         #Floating point inaccuracy.
-         eps=1e-7
-         return np.where((array > value-eps)*(array < value+eps))
+    """Find where a floating point value is equal to another."""
+    #Floating point inaccuracy.
+    eps=1e-7
+    return np.where((array > value-eps)*(array < value+eps))
 
-"""Just rebins the data"""
 def rebin(data, xaxis,newx):
-        if newx[0] < xaxis[0] or newx[-1]> xaxis[-1]:
-                raise ValueError("A value in newx is beyond the interpolation range")
-        intp=scipy.interpolate.InterpolatedUnivariateSpline(np.log(xaxis),data)
-        newdata=intp(np.log(newx))
-        return newdata
+    """Just rebins the data."""
+    if newx[0] < xaxis[0] or newx[-1]> xaxis[-1]:
+            raise ValueError("A value in newx is beyond the interpolation range")
+    intp=scipy.interpolate.InterpolatedUnivariateSpline(np.log(xaxis),data)
+    newdata=intp(np.log(newx))
+    return newdata
 
-"""Saves the figure, automatically determining file extension"""
 def save_figure(path):
-        bk=matplotlib.backends.backend
-        if path == "":
-                return
-        elif bk == 'TkAgg' or bk == 'Agg' or bk == 'GTKAgg':
-                path = path+".png"
-        elif bk == 'PDF' or bk == 'pdf':
-                path = path+".pdf"
-        elif bk == 'PS' or bk == 'ps':
-                path = path+".ps"
-        return plt.savefig(path)
-
-"""Little function to adjust a table so it has a different central value"""
+    """Save the figure, automatically determining file extension."""
+    bk=matplotlib.backends.backend
+    if path == "":
+            return
+    elif bk == 'TkAgg' or bk == 'Agg' or bk == 'GTKAgg':
+            path = path+".png"
+    elif bk == 'PDF' or bk == 'pdf':
+            path = path+".pdf"
+    elif bk == 'PS' or bk == 'ps':
+            path = path+".ps"
+    return plt.savefig(path)
+
 def corr_table(table, dvecs,table_name):
-        new=np.array(table)
-        new[12:,:] = table[12:,:]+2*table[0:12,:]*dvecs
-        pkd="/home/spb41/cosmomc-src/cosmomc/data/lya-interp/"
-        np.savetxt(pkd+table_name,new,("%1.3g","%1.3g","%1.3g","%1.3g","%1.3g","%1.3g","%1.3g","%1.3g","%1.3g","%1.3g","%1.3g"))
-        return new
-
-""" A class to be derived from by flux and matter power_spec classes. Stores various helper methods."""
-class power_spec:
-        # Snapshots
-        Snaps=()
-        #SDSS redshift bins.
-        Zz=np.array([1.0])
-        #SDSS kbins, in s/km units.
-        sdsskbins=np.array([1.0])
-        # Omega_matter
-        om=0.267
-        #Hubble constant
-        H0=0.71
-        #Boxsize in Mpc/h
-        box=60.0
-        #Size of the best-fit box, for testing varying boxsize
-        bfbox=60.0
-        #Some paths
-        knotpos=np.array([1.0])
-        base=""
-        pre=""
-        suf=""
-        ext=""
-        #For plotting
-        ymin=0.8
-        ymax=1.2
-        figprefix="/figure"
-        def __init__(self, Snaps=("snapshot_000", "snapshot_001","snapshot_002","snapshot_003","snapshot_004","snapshot_005","snapshot_006","snapshot_007","snapshot_008","snapshot_009","snapshot_010","snapshot_011"),
-             Zz=np.array([4.2,4.0,3.8,3.6,3.4,3.2,3.0,2.8,2.6,2.4,2.2,2.0]),
+    """Little function to adjust a table so it has a different central value."""
+    new=np.array(table)
+    new[12:,:] = table[12:,:]+2*table[0:12,:]*dvecs
+    pkd="/home/spb41/cosmomc-src/cosmomc/data/lya-interp/"
+    np.savetxt(pkd+table_name,new,("%1.3g","%1.3g","%1.3g","%1.3g","%1.3g","%1.3g","%1.3g","%1.3g","%1.3g","%1.3g","%1.3g"))
+    return new
+
+
+class PowerSpec:
+    """A class to be derived from by flux and matter PowerSpec classes. Stores various helper methods."""
+    
+    # snapshots
+    snaps=()
+    #SDSS redshift bins.
+    zz=np.array([1.0])
+    #SDSS kbins, in s/km units.
+    sdsskbins=np.array([1.0])
+    # Omega_matter
+    om=0.267
+    #hubble constant
+    h0=0.71
+    #Boxsize in Mpc/h
+    box=60.0
+    #Size of the best-fit box, for testing varying boxsize
+    bfbox=60.0
+    #Some paths
+    knotpos=np.array([1.0])
+    base=""
+    pre=""
+    suf=""
+    ext=""
+    #For plotting
+    ymin=0.8
+    ymax=1.2
+    figprefix="/figure"
+    def __init__(self, snaps=("snapshot_000", "snapshot_001","snapshot_002","snapshot_003","snapshot_004","snapshot_005","snapshot_006","snapshot_007","snapshot_008","snapshot_009","snapshot_010","snapshot_011"),
+             zz=np.array([4.2,4.0,3.8,3.6,3.4,3.2,3.0,2.8,2.6,2.4,2.2,2.0]),
              sdsskbins=np.array([0.00141,0.00178,0.00224,0.00282,0.00355,0.00447,0.00562,0.00708,0.00891,0.01122,0.01413,0.01778]),
-             knotpos=np.array([0.07,0.15,0.475, 0.75, 1.19, 1.89,4,25]), om=0.267, H0=0.71,box=60.0,
+             knotpos=np.array([0.07,0.15,0.475, 0.75, 1.19, 1.89,4,25]), om=0.267, h0=0.71,box=60.0,
              base="/home/spb41/Lyman-alpha/MinParametricRecon/runs/",suf="/", ext=".txt"):
-                if len(Snaps) != np.size(Zz):
-                        raise DataError("There are "+str(len(Snaps))+" snapshots, but "+str(np.size(Zz))+"redshifts given.")
-                self.Snaps=Snaps
-                self.Zz=Zz
-                self.sdsskbins=sdsskbins
-                self.om=om
-                self.H0=H0
-                self.box=box
-                self.knotpos=knotpos*H0
-                self.base=base
-                self.suf=suf
-                self.ext=ext
-                return
-
-        """ Get redshift associated with a snapshot """
-        def GetZ(self, snap):
-                ind=np.where(np.array(self.Snaps) == snap)
-                if np.size(ind):
-                        return self.Zz[ind]
-                else:
-                        raise DataError(str(snap)+" does not exist!")
-
-        """ Get snapshot associated with a redshift """
-        def GetSnap(self, redshift):
-                ind=wheref(self.Zz, redshift)
-                if np.size(ind):
-                        return str(np.asarray(self.Snaps)[ind][0])
-                else:
-                        raise DataError("No snapshot at redshift "+str(redshift))
-
-        """Get the k bins at a given redshift in h/Mpc units"""
-        def GetSDSSkbins(self, redshift):
-               return self.sdsskbins*self.Hubble(redshift)/(1.0+redshift)
-#Corr is /sqrt(self.H0)...
-
-        """ Hubble parameter. Hubble(Redshift) """
-        def Hubble(self, zz):
-                return 100*self.H0*math.sqrt(self.om*(1+zz)**3+(1-self.om))
+        """Plot the snapshots."""
+        if len(snaps) != np.size(zz):
+            raise DataError("There are "+str(len(snaps))+" snapshots, but "+str(np.size(zz))+"redshifts given.")
+        self.snaps=snaps
+        self.zz=zz
+        self.sdsskbins=sdsskbins
+        self.om=om
+        self.h0=h0
+        self.box=box
+        self.knotpos=knotpos*h0
+        self.base=base
+        self.suf=suf
+        self.ext=ext
+        return
+
+    def getz(self, snap):
+        """Get redshift associated with a snapshot."""
+        ind=np.where(np.array(self.snaps) == snap)
+        if np.size(ind):
+            return self.zz[ind]
+        else:
+            raise DataError(str(snap)+" does not exist!")
+
+    def getsnap(self, redshift):
+        """Get snapshot associated with a redshift."""
+        ind=wheref(self.zz, redshift)
+        if np.size(ind):
+            return str(np.asarray(self.snaps)[ind][0])
+        else:
+            raise DataError("No snapshot at redshift "+str(redshift))
+
+    def getsdsskbins(self, redshift):
+        """Get the k bins at a given redshift in h/Mpc units."""
+        return self.sdsskbins*self.hubble(redshift)/(1.0+redshift)
+        #Corr is /sqrt(self.h0)...
+
+    def hubble(self, zz):
+        """Hubble parameter."""
+        return 100*self.h0*math.sqrt(self.om*(1+zz)**3+(1-self.om))
         #Conversion factor between s/km and h/Mpc is (1+z)/H(z)
 
-        """ Do correct units conversion to return k and one-d power """
-        def loaddata(self, file, box):
-                #Adjust Fourier convention.
-                flux_power=np.loadtxt(file)
-                scale=self.H0/box
-                k=flux_power[1:,0]*scale*2.0*math.pi
-                PF=flux_power[1:,1]/scale
-                return (k, PF)
-
-        """ Plot comparisons between a bunch of sims on one graph
-        plot_z(Redshift, Sims to use ( eg, A1.14).
-        Note this will clear current figures."""
-        def plot_z(self,Knot,redshift,title="",ylabel="", legend=True):
-                #Load best-fit
-                (simk,BFPk)=self.loadpk(Knot.bstft+self.suf+self.pre+self.GetSnap(redshift)+self.ext,self.bfbox)
-                #Setup figure plot.
-                ind=wheref(self.Zz, redshift)
-                plt.figure(ind[0][0])
-                plt.clf()
-                if title != '':
-                        plt.title(title+" at z="+str(redshift))
-                plt.ylabel(ylabel)
-                plt.xlabel(r"$k\; (\mathrm{Mpc}^{-1})$")
-                line=np.array([])
-                legname=np.array([])
-                for sim in Knot.names:
-                        (k,Pk)=self.loadpk(sim+self.suf+self.pre+self.GetSnap(redshift)+self.ext,self.box)
-                        oi = np.where(simk <= k[-1])
-                        ti = np.where(simk[oi] >= k[0])
-                        relP=rebin(Pk, k, simk[oi][ti])
-                        relP=relP/rebin(BFPk, simk, simk[oi][ti])
-                        line=np.append(line, plt.semilogx(simk[oi][ti]/self.H0,relP,linestyle="-"))
-                        legname=np.append(legname,sim)
-                if legend:
-                        plt.legend(line,legname)
-                        plt.semilogx(self.knotpos,np.ones(len(self.knotpos)),"ro")
-                plt.ylim(self.ymin,self.ymax)
-                plt.xlim(simk[0]*0.8, 10)
-                return
-
-        """ Plot a whole suite of snapshots: plot_all(Knot, outdir) """
-        def plot_all(self, Knot,zzz=np.array([]), out=""):
-                if np.size(zzz) == 0:
-                        zzz=self.Zz    #lolz
-                for z in zzz:
-                        self.plot_z(Knot,z)
-                        if out != "":
-                                save_figure(out+self.figprefix+str(z))
-                return
-
-        """ Plot absolute power spectrum, not relative"""
-        def plot_power(self,path, redshift,colour="black"):
-                (k_g,Pk_g)=self.loadpk(path+self.suf+self.pre+self.GetSnap(redshift)+self.ext,self.box)
-                plt.loglog(k_g,Pk_g, color=colour)
-                plt.xlim(0.01,k_g[-1]*1.1)
-                plt.ylabel("P(k) /(h-3 Mpc3)")
-                plt.xlabel("k /(h MPc-1)")
-                plt.title("Power spectrum at z="+str(redshift))
-                return(k_g, Pk_g)
-
-        """ Plot absolute power for all redshifts """
-        def plot_power_all(self, Knot,zzz=np.array([]), out=""):
-                if np.size(zzz) == 0:
-                        zzz=self.Zz    #lolz
-                for z in zzz:
-                        ind=wheref(self.Zz, z)
-                        plt.figure(ind[0][0])
-                        for sim in Knot.names:
-                                self.plot_power(sim,z)
-                        if out != "":
-                                save_figure(out+self.figprefix+str(z))
-                return
-
-        """ Compare two power spectra directly. Smooths result.
-           plot_compare_two(first P(k), second P(k))"""
-        def plot_compare_two(self, one, onebox, two,twobox,colour=""):
-                (onek,oneP)=self.loadpk(one,onebox)
-                (twok,twoP)=self.loadpk(two,twobox)
-                onei = np.where(onek <= twok[-1])
-                twoi= np.where (onek[onei] >= twok[0])
-                relP=rebin(twoP, twok, onek[onei][twoi])
-                relP=relP/rebin(oneP, onek, onek[onei][twoi])
-                onek=onek[onei][twoi]
-                plt.title("Relative Power spectra "+one+" and "+two)
-                plt.ylabel(r"$P_2(k)/P_1(k)$")
-                plt.xlabel(r"$k\; (h\,\mathrm{Mpc}^{-1})$")
-                if colour == "":
-                        line=plt.semilogx(onek,relP)
-                else:
-                        line=plt.semilogx(onek,relP,color=colour)
-                plt.semilogx(self.knotpos,np.ones(len(self.knotpos)),"ro")
-                ind=np.where(onek < 10)
-                plt.ylim(min(relP[ind])*0.98,max(relP[ind])*1.01)
-                plt.xlim(onek[0]*0.8, 10)
-                return line
-
-        """Get the difference between two simulations on scales probed by
-           the SDSS power spectrum"""
-        def compare_two(self, one, two,redshift):
-                (onek,oneP)=self.loadpk(one,self.bfbox)
-                (twok,twoP)=self.loadpk(two,self.box)
-                onei = np.where(onek <= twok[-1])
-                twoi= np.where (onek[onei] >= twok[0])
-                relP=rebin(twoP, twok, onek[onei][twoi])
-                relP=relP/rebin(oneP, onek, onek[onei][twoi])
-                onek=onek[onei][twoi]
-                sdss=self.GetSDSSkbins(redshift)
-                relP_r=np.ones(np.size(sdss))
-                ind = np.where(sdss > onek[0])
-                relP_r[ind]=rebin(relP,onek,sdss[ind])
-                return relP_r
-
-        """Do the above 12 times to get a correction table"""
-        def compare_two_table(self,onedir, twodir):
-                nk=np.size(self.sdsskbins)
-                nz=np.size(self.Zz)-1
-                table=np.empty([nz,nk])
-                for i in np.arange(0,nz):
-                        sim=self.pre+self.GetSnap(self.Zz[i])+self.ext
-                        table[-1-i,:]=self.compare_two(onedir+sim,twodir+sim,self.Zz[i])
-                return table
-
-        """ Plot a whole redshift range of relative power spectra on the same figure.
-                plot_all(onedir, twodir)
-                Pass onedir and twodir as relative to basedir.
-                ie, for default settings something like
-                best-fit/flux-power/"""
-        def plot_compare_two_sdss(self, onedir, twodir, zzz=np.array([]), out="", title="", ylabel="", ymax=0,ymin=0, colour="",legend=False):
-                if np.size(zzz) == 0:
-                        zzz=self.Zz    #lolz
-                line=np.array([])
-                legname=np.array([])
-                sdss=self.sdsskbins
-                plt.xlabel(r"$k_v\; (s\,\mathrm{km}^{-1})$")
-                for z in zzz:
-                        sim=self.pre+self.GetSnap(z)+self.ext
-                        Pk=self.compare_two(onedir+sim,twodir+sim,z)
-                        if colour == "":
-                                line=np.append(line, plt.semilogx(sdss,Pk))
-                        else:
-                                line=np.append(line, plt.semilogx(sdss,Pk,color=colour))
-                        legname=np.append(legname,"z="+str(z))
-                plt.title(title)
-                if ylabel != "":
-                        plt.ylabel(ylabel)
-                if legend:
-                        plt.legend(line, legname,bbox_to_anchor=(0., 0, 1., .25), loc=3,ncol=3, mode="expand", borderaxespad=0.)
-                if ymax != 0 and ymin !=0:
-                        plt.ylim(ymin,ymax)
-                plt.xlim(sdss[0],sdss[-1])
-                plt.xticks(np.array([sdss[0],3e-3,5e-3,0.01,sdss[-1]]),("0.0014","0.003","0.005","0.01","0.0178"))
-                if out != "":
-                        save_figure(out)
-                return plt.gcf()
-
-        """ Plot a whole redshift range of relative power spectra on the same figure.
-                plot_all(onedir, twodir)
-                Pass onedir and twodir as relative to basedir.
-                ie, for default settings something like
-                best-fit/flux-power/
-                onedir uses bfbox, twodir uses box"""
-        def plot_compare_two_all(self, onedir, twodir, zzz=np.array([]), out="", title="", ylabel="", ymax=0,ymin=0, colour="",legend=False):
-                if np.size(zzz) == 0:
-                        zzz=self.Zz    #lolz
-                line=np.array([])
-                legname=np.array([])
-                for z in zzz:
-                        line=np.append(line, self.plot_compare_two(onedir+self.pre+self.GetSnap(z)+self.ext,self.bfbox,twodir+self.pre+self.GetSnap(z)+self.ext,self.box,colour))
-                        legname=np.append(legname,"z="+str(z))
-                if title == "":
-                        plt.title("Relative Power spectra "+onedir+" and "+twodir)
-                else:
-                        plt.title(title)
-                if ylabel != "":
-                        plt.ylabel(ylabel)
-                if legend:
-                        plt.legend(line, legname,bbox_to_anchor=(0., 0, 1., .25), loc=3,ncol=3, mode="expand", borderaxespad=0.)
-                if ymax != 0 and ymin !=0:
-                        plt.ylim(ymin,ymax)
-                if out != "":
-                        save_figure(out)
-                return plt.gcf()
-
-        """Get a power spectrum in the flat format we use"""
-        """for inputting some cosmomc tables"""
-        def GetFlat(self,dir, si=0.0):
-                Pk_sdss=np.empty([11, 12])
-                #Note this omits the z=2.0 bin
-                #SiIII corr now done on the fly in lya_sdss_viel.f90
-                for i in np.arange(0,np.size(self.Snaps)-1):
-                        scale=self.Hubble(self.Zz[i])/(1.0+self.Zz[i])
-                        (k,Pk)=self.loadpk(dir+self.Snaps[i]+self.ext, self.box)
-                        Fbar=math.exp(-0.0023*(1+self.Zz[i])**3.65)
-                        a=si/(1-Fbar)
-                        #The SiIII correction is kind of oscillatory, so we want
-                        #to average over the whole interval being probed.
-                        sdss=self.sdsskbins
-                        kmids=np.zeros(np.size(sdss)+1)
-                        sicorr=np.empty(np.size(sdss))
-                        for j in np.arange(0,np.size(sdss)-1):
-                                kmids[j+1]=math.exp((math.log(sdss[j+1])+math.log(sdss[j]))/2.)
-                        #Final segment should make no difference
-                        kmids[-1]=2*math.pi/2271+kmids[-2]
-                        #Near zero bad things will happen
-                        sicorr[0]=1+a**2+2*a*math.cos(2271*sdss[0])
-                        for j in np.arange(1,np.size(sdss)):
-                                sicorr[j]=1+a**2+2*a*(math.sin(2271*kmids[j+1])-math.sin(2271*kmids[j]))/(kmids[j+1]-kmids[j])/2271
-                        sdss=self.GetSDSSkbins(self.Zz[i])
-                        Pk_sdss[-i-1,:]=rebin(Pk, k, sdss)*scale*sicorr
-                return Pk_sdss
-
-        """ Calculate the flux derivatives for a single redshift
-        Output: (kbins d2P...kbins dP (flat vector of length 2xkbins))"""
-        def calc_z(self, redshift,s_knot, kbins):
-                #Array to store answers.
-                #Format is: k x (dP, d²P, χ²)
-                kbins=np.array(kbins)
-                nk=np.size(kbins)
-                snap=self.GetSnap(redshift)
-                results=np.zeros(2*nk)
-                if np.size(s_knot.qvals) > 0:
-                        results = np.zeros(4*nk)
-                pdifs=s_knot.pvals-s_knot.p0
-                qdifs=np.array([])
-                #If we have somethign where the parameter is redshift-dependent, eg, gamma.
-                if np.size(s_knot.p0) > 1:
-                        i=wheref(redshift, self.Zz)
-                        pdifs = s_knot.pvals[:,i]-s_knot.p0[i]
-                        if np.size(s_knot.qvals) > 0:
-                                qdifs=s_knot.qvals[:,i] - s_knot.q0[i]
-                npvals=np.size(pdifs)
-                #Load the data
-                (k,PFp0)=self.loadpk(s_knot.bstft+self.suf+snap+self.ext,s_knot.bfbox)
-                #This is to rescale by the mean flux, for generating mean flux tables.
-                ###
-                #tau_eff=0.0023*(1+redshift)**3.65
-                #tmin=0.2*((1+redshift)/4.)**2
-                #tmax=0.5*((1+redshift)/4.)**4
-                #teffs=tmin+s_knot.pvals*(tmax-tmin)/30.
-                #pdifs=teffs/tau_eff-1.
-                ###
-                PowerFluxes=np.zeros((npvals,np.size(k)))
-                for i in np.arange(0,np.size(s_knot.names)):
-                        (k,PowerFluxes[i,:])=self.loadpk(s_knot.names[i]+self.suf+snap+self.ext, s_knot.bfbox)
-                #So now we have an array of data values, which we want to rebin.
-                ind = np.where(kbins >= k[0])
-                difPF_rebin=np.ones((npvals,np.size(kbins)))
-                for i in np.arange(0, npvals):
-                        difPF_rebin[i,ind]=rebin(PowerFluxes[i,:]/PFp0,k,kbins[ind])
-                        #Set the things beyond the range of the interpolator
-                        #equal to the final value.
-                        if ind[0][0] > 0:
-                                difPF_rebin[i,0:ind[0][0]]=difPF_rebin[i,ind[0][0]]
-                #So now we have an array of data values.
-                #Pass each k value to flux_deriv in turn.
-                for k in np.arange(0,np.size(kbins)):
-                        #Format of returned data is:
-                        # y = ax**2 + bx + cz**2 +dz +e xz
-                        # derives = (a,b,c,d,e)
-                        derivs=self.flux_deriv(difPF_rebin[:,k], pdifs,qdifs)
-                        results[k]=derivs[0]
-                        results[nk+k]=derivs[1]
-                        if np.size(derivs) > 2:
-                                results[2*nk+k]=derivs[2]
-                                results[3*nk+k]=derivs[3]
-                return results
-
-        """ Calculate the flux derivatives for all redshifts
-        Input: Sims to load, parameter values, mean parameter value
-        Output: (2*kbins) x (zbins)"""
-        def calc_all(self, s_knot,kbins):
-                flux_derivatives=np.zeros((2*np.size(kbins),np.size(self.Zz)))
-                if np.size(s_knot.qvals) > 1:
-                        flux_derivatives=np.zeros((4*np.size(kbins),np.size(self.Zz)))
-                #Call flux_deriv_const_z for each redshift.
-                for i in np.arange(0,np.size(self.Zz)):
-                        flux_derivatives[:,i]=self.calc_z(self.Zz[i], s_knot,kbins)
-                return flux_derivatives
-
-        """Calculate the flux-derivative for a single redshift and k bin"""
-        def flux_deriv(self, PFdif, pdif, qdif=np.array([])):
-                pdif=np.ravel(pdif)
-                if np.size(pdif) != np.size(PFdif):
-                        raise DataError(str(np.size(pdif))+" parameter values, but "+str(np.size(PFdif))+" P_F values")
-                if np.size(pdif) < 2:
-                        raise DataError(str(np.size(pdif))+" pvals given. Need at least 2.")
-                PFdif=PFdif-1.0
-                if np.size(qdif) > 2:
-                        qdif=np.ravel(qdif)
-                        mat=np.vstack([pdif**2, pdif, qdif**2, qdif] ).T
-                else:
-                        mat=np.vstack([pdif**2, pdif] ).T
-                (derivs, residues,rank, sing)=np.linalg.lstsq(mat, PFdif)
-                return derivs
-
-        """ Get the error on one test simulation at single redshift """
-        def Get_Error_z(self, Sim, bstft,box, derivs, params, redshift,qarams=np.empty([])):
-                #Need to load and rebin the sim.
-                (k, test)=self.loadpk(Sim+self.suf+self.GetSnap(redshift)+self.ext, box)
-                (k,bf)=self.loadpk(bstft+self.suf+self.GetSnap(redshift)+self.ext,box)
-                kbins=self.Getkbins()
-                ind = np.where(kbins >= 1.0*2.0*math.pi*self.H0/box)
-                test2=np.ones(np.size(kbins))
-                test2[ind]=rebin(test/bf,k,kbins[ind])#
-                if ind[0][0] > 0:
-                        test2[0:ind[0][0]]=test2[ind[0][0]]
-                if np.size(qarams) > 0:
-                        guess=derivs.GetPF(params, redshift,qarams)+1.0
-                else:
-                        guess=derivs.GetPF(params, redshift)+1.0
-                return np.array((test2/guess))[0][0]
-
-""" A class written to store the various methods related to calculating of the flux derivatives and plotting of the flux power spectra"""
-class flux_pow(power_spec):
-        figprefix="/flux-figure"
-        kbins=np.array([])
-        def __init__(self, Snaps=("snapshot_000", "snapshot_001","snapshot_002","snapshot_003","snapshot_004","snapshot_005","snapshot_006","snapshot_007","snapshot_008","snapshot_009","snapshot_010","snapshot_011"),
-             Zz=np.array([4.2,4.0,3.8,3.6,3.4,3.2,3.0,2.8,2.6,2.4,2.2,2.0]),
-             sdsskbins=np.array([0.00141,0.00178,0.00224,0.00282,0.00355,0.00447,0.00562,0.00708,0.00891,0.01122,0.01413,0.01778]),
-             knotpos=np.array([0.07,0.15,0.475, 0.75, 1.19, 1.89,4,25]), om=0.266, H0=0.71,box=60.0,kmax=4.0,
-             base="/home/spb41/Lyman-alpha/MinParametricRecon/runs/",bf="best-fit/",suf="flux-power/", ext="_flux_power.txt"):
-                power_spec.__init__(self, Snaps,Zz,sdsskbins,knotpos, om, H0,box,base,suf, ext)
-                (k_bf,Pk_bf)= self.loadpk(bf+suf+"snapshot_000"+self.ext,self.bfbox)
-                ind=np.where(k_bf <= kmax)
-                self.kbins=k_bf[ind]
-
-        def plot_z(self,Sims,redshift,title="Relative Flux Power",ylabel=r"$\mathrm{P}_\mathrm{F}(k,p)\,/\,\mathrm{P}_\mathrm{F}(k,p_0)$", legend=True):
-                power_spec.plot_z(self,Sims,redshift,title,ylabel,legend)
-                if legend:
-                        kbins=self.GetSDSSkbins(redshift)
-                        plt.axvspan(kbins[0], kbins[-1], color="#B0B0B0")
-                plt.ylim(self.ymin,self.ymax)
-                plt.xlim(self.kbins[0]*0.8, 10)
-
-        """Get the kbins to interpolate onto"""
-        def Getkbins(self):
-               return self.kbins
-
-        """Load the SDSS power spectrum"""
-        def MacDonaldPF(self,sdss, zz):
-               psdss=sdss[np.where(sdss[:,0] == zz)][:,1:3]
-               fbar=math.exp(-0.0023*(1+zz)**3.65)
-               #multiply by the hubble parameter to be in 1/(km/s)
-               scale=self.Hubble(zz)/(1.0+zz)
-               PF=psdss[:,1]*fbar**2/scale
-               k=psdss[:,0]*scale
-               return (k, PF)
-
-        """Load a Pk. Different function due to needing to be different for each class"""
-        def loadpk(self, path,box):
-                #Adjust Fourier convention.
-                flux_power=np.loadtxt(self.base+path)
-                scale=self.H0/box
-                k=(flux_power[1:,0]-0.5)*scale*2.0*math.pi
-                PF=flux_power[1:,1]/scale
-                return (k, PF)
-
-
-""" A class to plot matter power spectra """
-class matter_pow(power_spec):
-        ob=0.0
-        #For plotting
-        ymin=0.4
-        ymax=1.6
-        figprefix="/matter-figure"
-        def __init__(self, Snaps=("snapshot_000", "snapshot_001","snapshot_002","snapshot_003","snapshot_004","snapshot_005","snapshot_006","snapshot_007","snapshot_008","snapshot_009","snapshot_010","snapshot_011"),
-             Zz=np.array([4.2,4.0,3.8,3.6,3.4,3.2,3.0,2.8,2.6,2.4,2.2,2.0]),
-             sdsskbins=np.array([0.00141,0.00178,0.00224,0.00282,0.00355,0.00447,0.00562,0.00708,0.00891,0.01122,0.01413,0.01778]),
-             knotpos=np.array([0.07,0.15,0.475, 0.75, 1.19, 1.89,4,25]), om=0.266,ob=0.0449, H0=0.71,box=60.0,
-             base="/home/spb41/Lyman-alpha/MinParametricRecon/runs/",suf="matter-power/", ext=".0", matpre="PK-by-"):
-                power_spec.__init__(self, Snaps,Zz,sdsskbins,knotpos, om, H0,box,base,suf,ext)
-                self.ob=ob
-                self.pre=matpre
-
-        def plot_z(self,Sims,redshift,title="Relative Matter Power",ylabel=r"$\mathrm{P}(k,p)\,/\,\mathrm{P}(k,p_0)$"):
-                power_spec.plot_z(self,Sims,redshift,title,ylabel)
-
-        """Load a Pk. Different function due to needing to be different for each class"""
-        def loadpk(self, path,box):
-                #Load baryon P(k)
-                matter_power=np.loadtxt(self.base+path)
-                scale=self.H0/box
+    def loaddata(self, file, box):
+        """Do correct units conversion to return k and one-d power."""
+        #Adjust Fourier convention.
+        fluxpower=np.loadtxt(file)
+        scale=self.h0/box
+        k=fluxpower[1:,0]*scale*2.0*math.pi
+        pf=fluxpower[1:,1]/scale
+        return (k, pf)
+
+    def plot_z(self,knot,redshift,title="",ylabel="", legend=True):
+        """Plot comparisons between a bunch of sims on one graph."""
+        """plot_z(Redshift, sims to use ( eg, A1.14)."""
+        """Note this will clear current figures."""
+        (simk,bfpk) = self.loadpk(knot.bstft + self.suf + self.pre + self.getsnap(redshift) + self.ext,self.bfbox)
+        #Setup figure plot.
+        ind=wheref(self.zz, redshift)
+        plt.figure(ind[0][0])
+        plt.clf()
+        if title != '':
+            plt.title(title+" at z="+str(redshift))
+        plt.ylabel(ylabel)
+        plt.xlabel(r"$k\; (\mathrm{Mpc}^{-1})$")
+        line=np.array([])
+        legname=np.array([])
+        for sim in knot.names:
+            (k,pk)=self.loadpk(sim+self.suf+self.pre+self.getsnap(redshift)+self.ext,self.box)
+            oi = np.where(simk <= k[-1])
+            ti = np.where(simk[oi] >= k[0])
+            relp=rebin(pk, k, simk[oi][ti])
+            relp=relp/rebin(bfpk, simk, simk[oi][ti])
+            line=np.append(line, plt.semilogx(simk[oi][ti]/self.h0,relp,linestyle="-"))
+            legname=np.append(legname,sim)
+        if legend:
+            plt.legend(line,legname)
+            plt.semilogx(self.knotpos,np.ones(len(self.knotpos)),"ro")
+        plt.ylim(self.ymin,self.ymax)
+        plt.xlim(simk[0]*0.8, 10)
+        return
+
+    def plot_all(self, knot,zzz=np.array([]), out=""):
+        """Plot a whole suite of snapshots."""
+        if np.size(zzz) == 0:
+            zzz=self.zz    #lolz
+        for z in zzz:
+            self.plot_z(knot,z)
+            if out != "":
+                save_figure(out+self.figprefix+str(z))
+        return
+
+    def plot_power(self,path, redshift,colour="black"):
+        """Plot absolute power spectrum, not relative."""
+        (k_g,pk_g) = self.loadpk(path + self.suf + self.pre + self.getsnap(redshift) + self.ext,self.box)
+        plt.loglog(k_g,pk_g, color=colour)
+        plt.xlim(0.01,k_g[-1]*1.1)
+        plt.ylabel("P(k) /(h-3 Mpc3)")
+        plt.xlabel("k /(h MPc-1)")
+        plt.title("Power spectrum at z="+str(redshift))
+        return(k_g, pk_g)
+
+    def plot_power_all(self, knot,zzz=np.array([]), out=""):
+        """Plot absolute power for all redshifts."""
+        if np.size(zzz) == 0:
+            zzz=self.zz    #lolz
+        for z in zzz:
+            ind=wheref(self.zz, z)
+            plt.figure(ind[0][0])
+            for sim in knot.names:
+                self.plot_power(sim,z)
+            if out != "":
+                save_figure(out+self.figprefix+str(z))
+        return
+
+    def plot_compare_two(self, one, onebox, two,twobox,colour=""):
+        """Compare two power spectra directly. Smooths result."""
+        """plot_compare_two(first P(k), second P(k))."""
+        (onek,onep)=self.loadpk(one,onebox)
+        (twok,twop)=self.loadpk(two,twobox)
+        onei = np.where(onek <= twok[-1])
+        twoi= np.where (onek[onei] >= twok[0])
+        relp=rebin(twop, twok, onek[onei][twoi])
+        relp=relp/rebin(onep, onek, onek[onei][twoi])
+        onek=onek[onei][twoi]
+        plt.title("Relative Power spectra "+one+" and "+two)
+        plt.ylabel(r"$P_2(k)/P_1(k)$")
+        plt.xlabel(r"$k\; (h\,\mathrm{Mpc}^{-1})$")
+        if colour == "":
+            line=plt.semilogx(onek,relp)
+        else:
+            line=plt.semilogx(onek,relp,color=colour)
+        plt.semilogx(self.knotpos,np.ones(len(self.knotpos)),"ro")
+        ind=np.where(onek < 10)
+        plt.ylim(min(relp[ind])*0.98,max(relp[ind])*1.01)
+        plt.xlim(onek[0]*0.8, 10)
+        return line
+
+    def compare_two(self, one, two,redshift):
+        """Get the difference between two simulations on scales probed by the SDSS power spectrum."""
+        (onek,onep)=self.loadpk(one,self.bfbox)
+        (twok,twop)=self.loadpk(two,self.box)
+        onei = np.where(onek <= twok[-1])
+        twoi= np.where (onek[onei] >= twok[0])
+        relp=rebin(twop, twok, onek[onei][twoi])
+        relp=relp/rebin(onep, onek, onek[onei][twoi])
+        onek=onek[onei][twoi]
+        sdss=self.getsdsskbins(redshift)
+        relp_r=np.ones(np.size(sdss))
+        ind = np.where(sdss > onek[0])
+        relp_r[ind]=rebin(relp,onek,sdss[ind])
+        return relp_r
+
+    def compare_two_table(self,onedir, twodir):
+        """Do the above 12 times to get a correction table."""
+        nk=np.size(self.sdsskbins)
+        nz=np.size(self.zz)-1
+        table=np.empty([nz,nk])
+        for i in np.arange(0,nz):
+            sim=self.pre+self.getsnap(self.zz[i])+self.ext
+            table[-1-i,:]=self.compare_two(onedir+sim,twodir+sim,self.zz[i])
+        return table
+
+    def plot_compare_two_sdss(self, onedir, twodir, zzz=np.array([]), out="", title="", ylabel="", ymax=0,ymin=0, colour="",legend=False):
+        """Plot a whole redshift range of relative power spectra on the same figure."""
+        """plot_all(onedir, twodir)."""
+        """Pass onedir and twodir as relative to basedir. ie, for default settings something like best-fit/flux-power/."""
+        if np.size(zzz) == 0:
+            zzz=self.zz    #lolz
+        line=np.array([])
+        legname=np.array([])
+        sdss=self.sdsskbins
+        plt.xlabel(r"$k_v\; (s\,\mathrm{km}^{-1})$")
+        for z in zzz:
+            sim=self.pre+self.getsnap(z)+self.ext
+            pk=self.compare_two(onedir+sim,twodir+sim,z)
+            if colour == "":
+                line=np.append(line, plt.semilogx(sdss,pk))
+            else:
+                line=np.append(line, plt.semilogx(sdss,pk,color=colour))
+            legname=np.append(legname,"z="+str(z))
+        plt.title(title)
+        if ylabel != "":
+            plt.ylabel(ylabel)
+        if legend:
+            plt.legend(line, legname,bbox_to_anchor=(0., 0, 1., .25), loc=3,ncol=3, mode="expand", borderaxespad=0.)
+        if ymax != 0 and ymin !=0:
+            plt.ylim(ymin,ymax)
+        plt.xlim(sdss[0],sdss[-1])
+        plt.xticks(np.array([sdss[0],3e-3,5e-3,0.01,sdss[-1]]),("0.0014","0.003","0.005","0.01","0.0178"))
+        if out != "":
+            save_figure(out)
+        return plt.gcf()
+
+    def plot_compare_two_all(self, onedir, twodir, zzz=np.array([]), out="", title="", ylabel="", ymax=0,ymin=0, colour="",legend=False):
+        """Plot a whole redshift range of relative power spectra on the same figure."""
+        """plot_all(onedir, twodir)"""
+        """Pass onedir and twodir as relative to basedir. ie, for default settings something like best-fit/flux-power/ onedir uses bfbox, twodir uses box."""
+        if np.size(zzz) == 0:
+            zzz=self.zz    #lolz
+        line=np.array([])
+        legname=np.array([])
+        for z in zzz:
+            line=np.append(line, self.plot_compare_two(onedir+self.pre+self.getsnap(z)+self.ext,self.bfbox,twodir+self.pre+self.getsnap(z)+self.ext,self.box,colour))
+            legname=np.append(legname,"z="+str(z))
+        if title == "":
+            plt.title("Relative Power spectra "+onedir+" and "+twodir)
+        else:
+            plt.title(title)
+        if ylabel != "":
+            plt.ylabel(ylabel)
+        if legend:
+            plt.legend(line, legname,bbox_to_anchor=(0., 0, 1., .25), loc=3,ncol=3, mode="expand", borderaxespad=0.)
+        if ymax != 0 and ymin !=0:
+            plt.ylim(ymin,ymax)
+        if out != "":
+            save_figure(out)
+        return plt.gcf()
+
+    def getflat(self,dir, si=0.0):
+        """Get a power spectrum in the flat format we use."""
+        """for inputting some cosmomc tables."""
+        pk_sdss=np.empty([11, 12])
+        #Note this omits the z=2.0 bin
+        #SiIII corr now done on the fly in lya_sdss_viel.f90
+        for i in np.arange(0,np.size(self.snaps)-1):
+            scale=self.hubble(self.zz[i])/(1.0+self.zz[i])
+            (k,pk)=self.loadpk(dir+self.snaps[i]+self.ext, self.box)
+            fbar=math.exp(-0.0023*(1+self.zz[i])**3.65)
+            a=si/(1-fbar)
+            #The SiIII correction is kind of oscillatory, so we want
+            #to average over the whole interval being probed.
+            sdss=self.sdsskbins
+            kmids=np.zeros(np.size(sdss)+1)
+            sicorr=np.empty(np.size(sdss))
+            for j in np.arange(0,np.size(sdss)-1):
+                kmids[j+1]=math.exp((math.log(sdss[j+1])+math.log(sdss[j]))/2.)
+            #Final segment should make no difference
+            kmids[-1]=2*math.pi/2271+kmids[-2]
+            #Near zero bad things will happen
+            sicorr[0]=1+a**2+2*a*math.cos(2271*sdss[0])
+            for j in np.arange(1,np.size(sdss)):
+                sicorr[j]=1+a**2+2*a*(math.sin(2271*kmids[j+1])-math.sin(2271*kmids[j]))/(kmids[j+1]-kmids[j])/2271
+            sdss=self.getsdsskbins(self.zz[i])
+            pk_sdss[-i-1,:]=rebin(pk, k, sdss)*scale*sicorr
+        return pk_sdss
+
+    def calc_z(self, redshift,s_knot, kbins):
+        """Calculate the flux derivatives for a single redshift. Output: (kbins d2P...kbins dP (flat vector of length 2xkbins))."""
+        #Array to store answers.
+        #Format is: k x (dP, d²P, χ²)
+        kbins=np.array(kbins)
+        nk=np.size(kbins)
+        snap=self.getsnap(redshift)
+        results=np.zeros(2*nk)
+        if np.size(s_knot.qvals) > 0:
+            results = np.zeros(4*nk)
+        pdifs=s_knot.pvals-s_knot.p0
+        qdifs=np.array([])
+        #If we have somethign where the parameter is redshift-dependent, eg, gamma.
+        if np.size(s_knot.p0) > 1:
+            i=wheref(redshift, self.zz)
+            pdifs = s_knot.pvals[:,i]-s_knot.p0[i]
+            if np.size(s_knot.qvals) > 0:
+                qdifs=s_knot.qvals[:,i] - s_knot.q0[i]
+        npvals=np.size(pdifs)
+        #Load the data
+        (k,pfp0)=self.loadpk(s_knot.bstft+self.suf+snap+self.ext,s_knot.bfbox)
+        #This is to rescale by the mean flux, for generating mean flux tables.
+        ###
+        #tau_eff=0.0023*(1+redshift)**3.65
+        #tmin=0.2*((1+redshift)/4.)**2
+        #tmax=0.5*((1+redshift)/4.)**4
+        #teffs=tmin+s_knot.pvals*(tmax-tmin)/30.
+        #pdifs=teffs/tau_eff-1.
+        ###
+        powerfluxes=np.zeros((npvals,np.size(k)))
+        for i in np.arange(0,np.size(s_knot.names)):
+            (k,powerfluxes[i,:])=self.loadpk(s_knot.names[i]+self.suf+snap+self.ext, s_knot.bfbox)
+        #So now we have an array of data values, which we want to rebin.
+        ind = np.where(kbins >= k[0])
+        difpf_rebin=np.ones((npvals,np.size(kbins)))
+        for i in np.arange(0, npvals):
+            difpf_rebin[i,ind]=rebin(powerfluxes[i,:]/pfp0,k,kbins[ind])
+            #Set the things beyond the range of the interpolator
+            #equal to the final value.
+            if ind[0][0] > 0:
+                difpf_rebin[i,0:ind[0][0]]=difpf_rebin[i,ind[0][0]]
+        #So now we have an array of data values.
+        #Pass each k value to flux_deriv in turn.
+        for k in np.arange(0,np.size(kbins)):
+            #Format of returned data is:
+            # y = ax**2 + bx + cz**2 +dz +e xz
+            # derives = (a,b,c,d,e)
+            derivs=self.flux_deriv(difpf_rebin[:,k], pdifs,qdifs)
+            results[k]=derivs[0]
+            results[nk+k]=derivs[1]
+            if np.size(derivs) > 2:
+                results[2*nk+k]=derivs[2]
+                results[3*nk+k]=derivs[3]
+        return results
+
+    def calc_all(self, s_knot,kbins):
+        """Calculate the flux derivatives for all redshifts. Input: sims to load, parameter values, mean parameter value. Output: (2*kbins) x (zbins)."""
+        flux_derivatives=np.zeros((2*np.size(kbins),np.size(self.zz)))
+        if np.size(s_knot.qvals) > 1:
+            flux_derivatives=np.zeros((4*np.size(kbins),np.size(self.zz)))
+        #Call flux_deriv_const_z for each redshift.
+        for i in np.arange(0,np.size(self.zz)):
+            flux_derivatives[:,i]=self.calc_z(self.zz[i], s_knot,kbins)
+        return flux_derivatives
+
+    def flux_deriv(self, pfdif, pdif, qdif=np.array([])):
+        """Calculate the flux-derivative for a single redshift and k bin."""
+        pdif=np.ravel(pdif)
+        if np.size(pdif) != np.size(pfdif):
+            raise DataError(str(np.size(pdif))+" parameter values, but "+str(np.size(pfdif))+" P_F values")
+        if np.size(pdif) < 2:
+            raise DataError(str(np.size(pdif))+" pvals given. Need at least 2.")
+        pfdif=pfdif-1.0
+        if np.size(qdif) > 2:
+            qdif=np.ravel(qdif)
+            mat=np.vstack([pdif**2, pdif, qdif**2, qdif] ).T
+        else:
+            mat=np.vstack([pdif**2, pdif] ).T
+        (derivs, residues,rank, sing)=np.linalg.lstsq(mat, pfdif)
+        return derivs
+
+    def get_error_z(self, sim, bstft,box, derivs, params, redshift,qarams=np.empty([])):
+        """Get the error on one test simulation at single redshift."""
+        #Need to load and rebin the sim.
+        (k, test)=self.loadpk(sim+self.suf+self.getsnap(redshift)+self.ext, box)
+        (k,bf)=self.loadpk(bstft+self.suf+self.getsnap(redshift)+self.ext,box)
+        kbins=self.getkbins()
+        ind = np.where(kbins >= 1.0*2.0*math.pi*self.h0/box)
+        test2=np.ones(np.size(kbins))
+        test2[ind]=rebin(test/bf,k,kbins[ind])#
+        if ind[0][0] > 0:
+            test2[0:ind[0][0]]=test2[ind[0][0]]
+        if np.size(qarams) > 0:
+            guess=derivs.Getpf(params, redshift,qarams)+1.0
+        else:
+            guess=derivs.Getpf(params, redshift)+1.0
+        return np.array((test2/guess))[0][0]
+
+class FluxPow(PowerSpec):
+    """A class written to store the various methods related to calculating of the flux derivatives and plotting of the flux power spectra."""
+    
+    figprefix="/flux-figure"
+    kbins=np.array([])
+    def __init__(self, snaps=("snapshot_000", "snapshot_001","snapshot_002","snapshot_003","snapshot_004","snapshot_005","snapshot_006","snapshot_007","snapshot_008","snapshot_009","snapshot_010","snapshot_011"),
+         zz=np.array([4.2,4.0,3.8,3.6,3.4,3.2,3.0,2.8,2.6,2.4,2.2,2.0]),
+         sdsskbins=np.array([0.00141,0.00178,0.00224,0.00282,0.00355,0.00447,0.00562,0.00708,0.00891,0.01122,0.01413,0.01778]),
+         knotpos=np.array([0.07,0.15,0.475, 0.75, 1.19, 1.89,4,25]), om=0.266, h0=0.71,box=60.0,kmax=4.0,
+         base="/home/spb41/Lyman-alpha/MinParametricRecon/runs/",bf="best-fit/",suf="flux-power/", ext="_fluxpower.txt"):
+        """Do something with the snapshots."""
+        PowerSpec.__init__(self, snaps,zz,sdsskbins,knotpos, om, h0,box,base,suf, ext)
+        (k_bf,pk_bf)= self.loadpk(bf+suf+"snapshot_000"+self.ext,self.bfbox)
+        ind=np.where(k_bf <= kmax)
+        self.kbins=k_bf[ind]
+
+    def plot_z(self,sims,redshift,title="Relative Flux Power",ylabel=r"$\mathrm{P}_\mathrm{F}(k,p)\,/\,\mathrm{P}_\mathrm{F}(k,p_0)$", legend=True):
+        """Plot comparisons between a bunch of sims on one graph."""
+        PowerSpec.plot_z(self,sims,redshift,title,ylabel,legend)
+        if legend:
+                kbins=self.getsdsskbins(redshift)
+                plt.axvspan(kbins[0], kbins[-1], color="#B0B0B0")
+        plt.ylim(self.ymin,self.ymax)
+        plt.xlim(self.kbins[0]*0.8, 10)
+
+    def getkbins(self):
+        """Get the kbins to interpolate onto."""
+        return self.kbins
+
+    def macdonaldpf(self,sdss, zz):
+        """Load the SDSS power spectrum."""
+        psdss=sdss[np.where(sdss[:,0] == zz)][:,1:3]
+        fbar=math.exp(-0.0023*(1+zz)**3.65)
+        #multiply by the hubble parameter to be in 1/(km/s)
+        scale=self.hubble(zz)/(1.0+zz)
+        pf=psdss[:,1]*fbar**2/scale
+        k=psdss[:,0]*scale
+        return (k, pf)
+
+    def loadpk(self, path,box):
+        """Load a pk. Different function due to needing to be different for each class."""
+        #Adjust Fourier convention.
+        fluxpower=np.loadtxt(self.base+path)
+        scale=self.h0/box
+        k=(fluxpower[1:,0]-0.5)*scale*2.0*math.pi
+        pf=fluxpower[1:,1]/scale
+        return (k, pf)
+
+class MatterPow(PowerSpec):
+    """A class to plot matter power spectra."""
+    
+    ob=0.0
+    #For plotting
+    ymin=0.4
+    ymax=1.6
+    figprefix="/matter-figure"
+    def __init__(self, snaps=("snapshot_000", "snapshot_001","snapshot_002","snapshot_003","snapshot_004","snapshot_005","snapshot_006","snapshot_007","snapshot_008","snapshot_009","snapshot_010","snapshot_011"),
+         zz=np.array([4.2,4.0,3.8,3.6,3.4,3.2,3.0,2.8,2.6,2.4,2.2,2.0]),
+         sdsskbins=np.array([0.00141,0.00178,0.00224,0.00282,0.00355,0.00447,0.00562,0.00708,0.00891,0.01122,0.01413,0.01778]),
+         knotpos=np.array([0.07,0.15,0.475, 0.75, 1.19, 1.89,4,25]), om=0.266,ob=0.0449, h0=0.71,box=60.0,
+         base="/home/spb41/Lyman-alpha/MinParametricRecon/runs/",suf="matter-power/", ext=".0", matpre="pk-by-"):
+        """Do something with the snapshots."""
+        PowerSpec.__init__(self, snaps,zz,sdsskbins,knotpos, om, h0,box,base,suf,ext)
+        self.ob=ob
+        self.pre=matpre
+
+    def plot_z(self,sims,redshift,title="Relative Matter Power",ylabel=r"$\mathrm{P}(k,p)\,/\,\mathrm{P}(k,p_0)$"):
+        """Plot comparisons between a bunch of sims on one graph."""
+        PowerSpec.plot_z(self,sims,redshift,title,ylabel)
+
+    def loadpk(self, path,box):
+        """Load a pk. Different function due to needing to be different for each class."""
+        #Load baryon P(k)
+        matterpower=np.loadtxt(self.base+path)
+        scale=self.h0/box
+        #Adjust Fourier convention.
+        simk=matterpower[1:,0]*scale*2.0*math.pi
+        pkbar=matterpower[1:,1]/scale**3
+        #Load DM P(k)
+        matterpower=np.loadtxt(self.base+re.sub("by","DM",path))
+        pkdm=matterpower[1:,1]/scale**3
+        pk=(pkbar*self.ob+pkdm*(self.om-self.ob))/self.om
+        return (simk,pk)
+
+    def plot_power(self,path, redshift,camb_filename=""):
+        """Plot absolute power spectrum, not relative."""
+        (k_g, pk_g)=PowerSpec.plot_power(self,path,redshift)
+        sigma=2.0
+        pkg=np.loadtxt(self.base+path+self.suf+self.pre+self.getsnap(redshift)+self.ext)
+        samp_err=pkg[1:,2]
+        sqrt_err=np.array(np.sqrt(samp_err))
+        plt.loglog(k_g,pk_g*(1+sigma*(2.0/sqrt_err+1.0/samp_err)),linestyle="-.",color="black")
+        plt.loglog(k_g,pk_g*(1-sigma*(2.0/sqrt_err+1.0/samp_err)),linestyle="-.",color="black")
+        if camb_filename != "":
+                camb=np.loadtxt(camb_filename)
                 #Adjust Fourier convention.
-                simk=matter_power[1:,0]*scale*2.0*math.pi
-                Pkbar=matter_power[1:,1]/scale**3
-                #Load DM P(k)
-                matter_power=np.loadtxt(self.base+re.sub("by","DM",path))
-                PkDM=matter_power[1:,1]/scale**3
-                Pk=(Pkbar*self.ob+PkDM*(self.om-self.ob))/self.om
-                return (simk,Pk)
-
-        """ Plot absolute power spectrum, not relative"""
-        def plot_power(self,path, redshift,camb_filename=""):
-                (k_g, Pk_g)=power_spec.plot_power(self,path,redshift)
-                sigma=2.0
-                pkg=np.loadtxt(self.base+path+self.suf+self.pre+self.GetSnap(redshift)+self.ext)
-                samp_err=pkg[1:,2]
-                sqrt_err=np.array(np.sqrt(samp_err))
-                plt.loglog(k_g,Pk_g*(1+sigma*(2.0/sqrt_err+1.0/samp_err)),linestyle="-.",color="black")
-                plt.loglog(k_g,Pk_g*(1-sigma*(2.0/sqrt_err+1.0/samp_err)),linestyle="-.",color="black")
-                if camb_filename != "":
-                        camb=np.loadtxt(camb_filename)
-                        #Adjust Fourier convention.
-                        k=camb[:,0]*self.H0
-                        #NOW THERE IS NO h in the T anywhere.
-                        Pk=camb[:,1]
-                        plt.loglog(k/self.H0, Pk, linestyle="--")
-                plt.xlim(0.01,k_g[-1]*1.1)
-                return(k_g, Pk_g)
-
-"""The PDF is an instance of the power_spec class. Perhaps poor naming"""
-class flux_pdf(power_spec):
-        def __init__(self, Snaps=("snapshot_006","snapshot_007","snapshot_008","snapshot_009","snapshot_010","snapshot_011"),
-             Zz=np.array([3.0,2.8,2.6,2.4,2.2,2.0]),
-             sdsskbins=np.arange(0,20),
-             knotpos=np.array([]), om=0.266,ob=0.0449, H0=0.71,box=48.0,
-             base="/home/spb41/Lyman-alpha/MinParametricRecon/runs/",suf="flux-pdf/", ext="_flux_pdf.txt",):
-                power_spec.__init__(self, Snaps,Zz,sdsskbins,knotpos, om, H0,box,base,suf,ext)
-
-        def loadpk(self, path, box):
-                flux_pdf = np.loadtxt(self.base+path)
-                return(flux_pdf[:,0], flux_pdf[:,1])
-
-        """ Compare two power spectra directly. Smooths result.
-           plot_compare_two(first P(k), second P(k))"""
-        def plot_compare_two(self, one, onebox, two,twobox,colour=""):
-                (onek,oneP)=self.loadpk(one,onebox)
-                (twok,twoP)=self.loadpk(two,twobox)
-                relP=oneP/twoP
-                plt.title("Relative flux PDF "+one+" and "+two)
-                plt.ylabel(r"$F_2(k)/F_1(k)$")
-                plt.xlabel(r"$Flux$")
-                line=plt.semilogy(onek,relP)
-                return line
-
-        """ Plot absolute power spectrum, not relative"""
-        def plot_power(self,path, redshift):
-                (k,Pdf)=self.loadpk(path+self.suf+self.pre+self.GetSnap(redshift)+self.ext,self.box)
-                plt.semilogy(k,Pdf, color="black", linewidth="1.5")
-                plt.ylabel("P(k) /(h-3 Mpc3)")
-                plt.xlabel("k /(h MPc-1)")
-                plt.title("PDF at z="+str(redshift))
-                return(k, Pdf)
-
-        """ Calculate the flux derivatives for a single redshift
-        Output: (kbins d2P...kbins dP (flat vector of length 2x21))"""
-        def calc_z(self, redshift,s_knot):
-                #Array to store answers.
-                #Format is: k x (dP, d²P, χ²)
-                npvals=np.size(s_knot.pvals)
-                nk=21
-                results=np.zeros(2*nk)
-                pdifs=s_knot.pvals-s_knot.p0
-                #This is to rescale by the mean flux, for generating mean flux tables.
-                ###
-                #tau_eff=0.0023*(1+redshift)**3.65
-                #tmin=0.2*((1+redshift)/4.)**2
-                #tmax=0.5*((1+redshift)/4.)**4
-                #teffs=tmin+s_knot.pvals*(tmax-tmin)/30.
-                #pdifs=teffs/tau_eff-1.
-                ###
-                ured=np.ceil(redshift*5)/5.
-                lred=np.floor(redshift*5)/5.
-                usnap=self.GetSnap(ured)
-                lsnap=self.GetSnap(lred)
-                #Load the data
-                (k,uPFp0)=self.loadpk(s_knot.bstft+self.suf+usnap+self.ext,s_knot.bfbox)
-                uPower=np.zeros((npvals,np.size(k)))
-                for i in np.arange(0,np.size(s_knot.names)):
-                        (k,uPower[i,:])=self.loadpk(s_knot.names[i]+self.suf+usnap+self.ext, s_knot.bfbox)
-                (k,lPFp0)=self.loadpk(s_knot.bstft+self.suf+lsnap+self.ext,s_knot.bfbox)
-                lPower=np.zeros((npvals,np.size(k)))
-                for i in np.arange(0,np.size(s_knot.names)):
-                        (k,lPower[i,:])=self.loadpk(s_knot.names[i]+self.suf+lsnap+self.ext, s_knot.bfbox)
-                PowerFluxes=5*((redshift-lred)*uPower+(ured-redshift)*lPower)
-                PFp0=5*((redshift-lred)*uPFp0+(ured-redshift)*lPFp0)
-                #So now we have an array of data values.
-                #Pass each k value to flux_deriv in turn.
-                for k in np.arange(0,nk):
-                        (dPF, d2PF,chi2)=self.flux_deriv(PowerFluxes[:,k]/PFp0[k], pdifs)
-                        results[k]=d2PF
-                        results[nk+k]=dPF
-                return results
-
-        """Get the kbins to interpolate onto"""
-        def Getkbins(self):
-               return np.arange(0,20,1)+0.5
-        """ Plot comparisons between a bunch of sims on one graph
-        plot_z(Redshift, Sims to use ( eg, A1.14).
-        Note this will clear current figures."""
-        def plot_z(self,Knot,redshift,title="",ylabel=""):
-                #Load best-fit
-                (simk,BFPk)=self.loadpk(Knot.bstft+self.suf+self.pre+self.GetSnap(redshift)+self.ext,self.bfbox)
-                #Setup figure plot.
-                ind=wheref(self.Zz, redshift)
-                plt.figure(ind[0][0])
-                plt.clf()
-                if title != '':
-                        plt.title(title+" at z="+str(redshift),)
-                plt.ylabel(ylabel)
-                plt.xlabel(r"$\mathcal{F}$")
-                line=np.array([])
-                legname=np.array([])
-                for sim in Knot.names:
-                        (k,Pk)=self.loadpk(sim+self.suf+self.pre+self.GetSnap(redshift)+self.ext,self.box)
-                        line=np.append(line, plt.semilogy(simk,Pk/BFPk,linestyle="-", linewidth=1.5))
-                        legname=np.append(legname,sim)
-                plt.legend(line,legname)
-                return
-
-        """Get a power spectrum in the flat format we use"""
-        """for inputting some cosmomc tables"""
-        def GetFlat(self,dir):
-                Pk_sdss=np.empty([11, 12])
-                #For z=2.07 we need to average snap_011 and snap_010
-                z=2.07
-                (k,PF_a)=self.loadpk(dir+self.suf+"snapshot_011"+self.ext, self.box)
-                (k,PF_b)=self.loadpk(dir+self.suf+"snapshot_010"+self.ext, self.box)
-                PF1=(z-2.0)*5*(PF_b-PF_a)+PF_a
-                z=2.52
-                (k,PF_a)=self.loadpk(dir+self.suf+"snapshot_009"+self.ext, self.box)
-                (k,PF_b)=self.loadpk(dir+self.suf+"snapshot_008"+self.ext, self.box)
-                PF2=(z-2.4)*5*(PF_b-PF_a)+PF_a
-                z=2.94
-                (k,PF_a)=self.loadpk(dir+self.suf+"snapshot_007"+self.ext, self.box)
-                (k,PF_b)=self.loadpk(dir+self.suf+"snapshot_006"+self.ext, self.box)
-                PF3=(z-2.8)*5*(PF_b-PF_a)+PF_a
-                PDF = np.array([PF1,PF2,PF3])
-                np.savetxt(sys.stdout,PDF.T("%1.8f","%1.8f","%1.8f"))
-                return (PF1, PF2, PF3)
+                k=camb[:,0]*self.h0
+                #NOW THERE IS NO h in the T anywhere.
+                pk=camb[:,1]
+                plt.loglog(k/self.h0, pk, linestyle="--")
+        plt.xlim(0.01,k_g[-1]*1.1)
+        return(k_g, pk_g)
+
+class FluxPdf(PowerSpec):
+    """The PDF is an instance of the PowerSpec class. Perhaps poor naming."""
+    
+    def __init__(self, snaps=("snapshot_006","snapshot_007","snapshot_008","snapshot_009","snapshot_010","snapshot_011"),
+        zz=np.array([3.0,2.8,2.6,2.4,2.2,2.0]),
+        sdsskbins=np.arange(0,20),
+        knotpos=np.array([]), om=0.266,ob=0.0449, h0=0.71,box=48.0,
+        base="/home/spb41/Lyman-alpha/MinParametricRecon/runs/",suf="flux-pdf/", ext="_fluxpdf.txt",):
+        """Do something with the snapshots."""
+        PowerSpec.__init__(self, snaps,zz,sdsskbins,knotpos, om, h0,box,base,suf,ext)
+
+    def loadpk(self, path, box):
+        """Load a pk. Different function due to needing to be different for each class."""
+        fluxpdf = np.loadtxt(self.base+path)
+        return(fluxpdf[:,0], fluxpdf[:,1])
+
+    def plot_compare_two(self, one, onebox, two,twobox,colour=""):
+        """Compare two power spectra directly. Smooths result."""
+        """plot_compare_two(first P(k), second P(k))."""
+        (onek,onep)=self.loadpk(one,onebox)
+        (twok,twop)=self.loadpk(two,twobox)
+        relp=onep/twop
+        plt.title("Relative flux PDF "+one+" and "+two)
+        plt.ylabel(r"$F_2(k)/F_1(k)$")
+        plt.xlabel(r"$Flux$")
+        line=plt.semilogy(onek,relp)
+        return line
+
+    def plot_power(self,path, redshift):
+        """Plot absolute power spectrum, not relative."""
+        (k,pdf)=self.loadpk(path+self.suf+self.pre+self.getsnap(redshift)+self.ext,self.box)
+        plt.semilogy(k,pdf, color="black", linewidth="1.5")
+        plt.ylabel("P(k) /(h-3 Mpc3)")
+        plt.xlabel("k /(h MPc-1)")
+        plt.title("PDF at z="+str(redshift))
+        return(k, pdf)
+
+    def calc_z(self, redshift,s_knot):
+        """Calculate the flux derivatives for a single redshift. Output: (kbins d2P...kbins dP (flat vector of length 2x21))."""
+        #Array to store answers.
+        #Format is: k x (dP, d²P, χ²)
+        npvals=np.size(s_knot.pvals)
+        nk=21
+        results=np.zeros(2*nk)
+        pdifs=s_knot.pvals-s_knot.p0
+        #This is to rescale by the mean flux, for generating mean flux tables.
+        ###
+        #tau_eff=0.0023*(1+redshift)**3.65
+        #tmin=0.2*((1+redshift)/4.)**2
+        #tmax=0.5*((1+redshift)/4.)**4
+        #teffs=tmin+s_knot.pvals*(tmax-tmin)/30.
+        #pdifs=teffs/tau_eff-1.
+        ###
+        ured=np.ceil(redshift*5)/5.
+        lred=np.floor(redshift*5)/5.
+        usnap=self.getsnap(ured)
+        lsnap=self.getsnap(lred)
+        #Load the data
+        (k,upfp0)=self.loadpk(s_knot.bstft+self.suf+usnap+self.ext,s_knot.bfbox)
+        upower=np.zeros((npvals,np.size(k)))
+        for i in np.arange(0,np.size(s_knot.names)):
+            (k,upower[i,:])=self.loadpk(s_knot.names[i]+self.suf+usnap+self.ext, s_knot.bfbox)
+        (k,lpfp0)=self.loadpk(s_knot.bstft+self.suf+lsnap+self.ext,s_knot.bfbox)
+        lpower=np.zeros((npvals,np.size(k)))
+        for i in np.arange(0,np.size(s_knot.names)):
+            (k,lpower[i,:])=self.loadpk(s_knot.names[i]+self.suf+lsnap+self.ext, s_knot.bfbox)
+        powerfluxes=5*((redshift-lred)*upower+(ured-redshift)*lpower)
+        pfp0=5*((redshift-lred)*upfp0+(ured-redshift)*lpfp0)
+        #So now we have an array of data values.
+        #Pass each k value to flux_deriv in turn.
+        for k in np.arange(0,nk):
+            (dpf, d2pf,chi2)=self.flux_deriv(powerfluxes[:,k]/pfp0[k], pdifs)
+            results[k]=d2pf
+            results[nk+k]=dpf
+        return results
+
+        def getkbins(self):
+            """Get the kbins to interpolate onto."""
+            return np.arange(0,20,1)+0.5
+        
+    def plot_z(self,knot,redshift,title="",ylabel=""):
+        """Plot comparisons between a bunch of sims on one graph."""
+        """plot_z(Redshift, sims to use ( eg, A1.14)."""
+        """Note this will clear current figures."""
+        #Load best-fit
+        (simk,bfpk)=self.loadpk(knot.bstft+self.suf+self.pre+self.getsnap(redshift)+self.ext,self.bfbox)
+        #Setup figure plot.
+        ind=wheref(self.zz, redshift)
+        plt.figure(ind[0][0])
+        plt.clf()
+        if title != '':
+            plt.title(title+" at z="+str(redshift),)
+        plt.ylabel(ylabel)
+        plt.xlabel(r"$\mathcal{F}$")
+        line=np.array([])
+        legname=np.array([])
+        for sim in knot.names:
+            (k,pk)=self.loadpk(sim+self.suf+self.pre+self.getsnap(redshift)+self.ext,self.box)
+            line=np.append(line, plt.semilogy(simk,pk/bfpk,linestyle="-", linewidth=1.5))
+            legname=np.append(legname,sim)
+        plt.legend(line,legname)
+        return
+
+    def getflat(self,dir):
+        """Get a power spectrum in the flat format we use."""
+        """for inputting some cosmomc tables."""
+        # pk_sdss=np.empty([11, 12])
+        #For z=2.07 we need to average snap_011 and snap_010
+        z=2.07
+        (k,pf_a)=self.loadpk(dir+self.suf+"snapshot_011"+self.ext, self.box)
+        (k,pf_b)=self.loadpk(dir+self.suf+"snapshot_010"+self.ext, self.box)
+        pf1=(z-2.0)*5*(pf_b-pf_a)+pf_a
+        z=2.52
+        (k,pf_a)=self.loadpk(dir+self.suf+"snapshot_009"+self.ext, self.box)
+        (k,pf_b)=self.loadpk(dir+self.suf+"snapshot_008"+self.ext, self.box)
+        pf2=(z-2.4)*5*(pf_b-pf_a)+pf_a
+        z=2.94
+        (k,pf_a)=self.loadpk(dir+self.suf+"snapshot_007"+self.ext, self.box)
+        (k,pf_b)=self.loadpk(dir+self.suf+"snapshot_006"+self.ext, self.box)
+        pf3=(z-2.8)*5*(pf_b-pf_a)+pf_a
+        pdf = np.array([pf1,pf2,pf3])
+        np.savetxt(sys.stdout,pdf.T("%1.8f","%1.8f","%1.8f"))
+        return (pf1, pf2, pf3)
 
 if __name__=='__main__':
-        flux=flux_pow()
-        matter=matter_pow()
-        fpdf=flux_pdf()
-        A_knot=knot(("A0.54/","A0.74/","A0.84/","A1.04/", "A1.14/","A1.34/"), (0.54,0.74,0.84,1.04,1.14,1.34),0.94,"best-fit/", 60)
-        AA_knot=knot(("AA0.54/","AA0.74/","AA1.14/","AA1.34/"), (0.54,0.74,1.14,1.34),0.94,"boxcorr400/", 120)
-        B_knot=knot(("B0.33/","B0.53/","B0.73/", "B0.83/", "B1.03/","B1.13/", "B1.33/"), (0.33,0.53,0.73,0.83,1.03, 1.13,1.33),0.93,"best-fit/", 60)
-        C_knot=knot(("C0.11/", "C0.31/","C0.51/","C0.71/","C1.11/","C1.31/","C1.51/"),(0.11, 0.31,0.51,0.71, 1.11,1.31,1.51),0.91,"bf2/", 60)
-        D_knot=knot(("D0.50/","D0.70/","D1.10/","D1.20/", "D1.30/", "D1.50/", "D1.70/"),(0.50, 0.70,1.10,1.20,1.30, 1.50, 1.70),0.90,"bfD/", 48)
-        interp=flux_interp(flux, (AA_knot, B_knot, C_knot, D_knot))
+        # flux=FluxPow()
+        # matter=MatterPow()
+        # fpdf=fluxpdf()
+        # A_knot=knot(("A0.54/","A0.74/","A0.84/","A1.04/", "A1.14/","A1.34/"), (0.54,0.74,0.84,1.04,1.14,1.34),0.94,"best-fit/", 60)
+        # AA_knot=knot(("AA0.54/","AA0.74/","AA1.14/","AA1.34/"), (0.54,0.74,1.14,1.34),0.94,"boxcorr400/", 120)
+        # B_knot=knot(("B0.33/","B0.53/","B0.73/", "B0.83/", "B1.03/","B1.13/", "B1.33/"), (0.33,0.53,0.73,0.83,1.03, 1.13,1.33),0.93,"best-fit/", 60)
+        # C_knot=knot(("C0.11/", "C0.31/","C0.51/","C0.71/","C1.11/","C1.31/","C1.51/"),(0.11, 0.31,0.51,0.71, 1.11,1.31,1.51),0.91,"bf2/", 60)
+        # D_knot=knot(("D0.50/","D0.70/","D1.10/","D1.20/", "D1.30/", "D1.50/", "D1.70/"),(0.50, 0.70,1.10,1.20,1.30, 1.50, 1.70),0.90,"bfD/", 48)
+        # interp=flux_interp(flux, (AA_knot, B_knot, C_knot, D_knot))
 #Thermal parameters for models
 
 #Models where gamma is varied
@@ -692,6 +695,6 @@ def GetFlat(self,dir):
         bf2ag=np.array([1.0184905 ,1.0332849 ,1.0404092 ,1.0538235 ,1.0598905 ,1.0695964 ,1.0771414 ,1.0827464 ,1.0904784 ,1.0964669 ,1.1009428,1.1068357])
         bf2at=np.array([26914.856, 27497.555, 28407.030, 28357.785, 28919.376, 28610.038, 28478.749, 28484.856, 27841.160, 27335.328, 26933.453,26099.643])/1e3
 
-        G_knot=knot(("bf2z/","bf2b/","bf2c/","bf2d/", "bf2T15/","bf2T35/","bf2T45/","bf2a/"),(bf2zg,bf2bg,bf2cg,bf2dg,T15g,T35g,T45g,bf2ag),bf2g,"bf2/",60, (bf2zt,bf2bt,bf2ct,bf2dt,T15t,T35t,T45t,bf2at),bf2t)
-        g_int=flux_interp(flux,(G_knot))
+        # G_knot=knot(("bf2z/","bf2b/","bf2c/","bf2d/", "bf2T15/","bf2T35/","bf2T45/","bf2a/"),(bf2zg,bf2bg,bf2cg,bf2dg,T15g,T35g,T45g,bf2ag),bf2g,"bf2/",60, (bf2zt,bf2bt,bf2ct,bf2dt,T15t,T35t,T45t,bf2at),bf2t)
+        # g_int=flux_interp(flux,(G_knot))
 
diff --git a/schwarzchild.ipynb b/schwarzchild.ipynb
new file mode 100644
index 0000000..fbaf01c
--- /dev/null
+++ b/schwarzchild.ipynb
@@ -0,0 +1,165 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "da9f47ba-2fb4-4ac5-9581-4f6f05b5d2e4",
+   "metadata": {},
+   "source": [
+    "Problem 1.3.a."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "b7ccde90-d65f-408f-82f7-ba12bcf9c895",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Schwarzschild radius of the Sun: 2953.9779243533935 meter\n"
+     ]
+    }
+   ],
+   "source": [
+    "import pint\n",
+    "ureg  = pint.UnitRegistry()\n",
+    "\n",
+    "G     = 6.67430e-11 * ureg.meter**3 / (ureg.kilogram * ureg.second**2)\n",
+    "M_sun = 1.989e30 * ureg.kilogram\n",
+    "c     = 2.998e8 * ureg.meter / ureg.second\n",
+    "\n",
+    "# Schwarzschild radius\n",
+    "R_s   = (2 * G * M_sun) / c**2\n",
+    "\n",
+    "print(f\"Schwarzschild radius of the Sun: {R_s.to(ureg.meter)}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a9afa71f-4643-4c4c-9256-245e7a4072a8",
+   "metadata": {},
+   "source": [
+    "Problem 1.3.b."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "e093b079-b8d9-48ca-8341-01a663c3fd3a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Using nested for loops\n",
+    "def matrix_multiply_loops(A, B):\n",
+    "    result = [[0 for _ in range(len(B[0]))] for _ in range(len(A))]\n",
+    "    for i in range(len(A)):\n",
+    "        for j in range(len(B[0])):\n",
+    "            for k in range(len(B)):\n",
+    "                result[i][j] += A[i][k] * B[k][j]\n",
+    "    return result\n",
+    "\n",
+    "# Using a list comprehension\n",
+    "def matrix_multiply_list_comp(A, B):\n",
+    "    return [[sum(A[i][k] * B[k][j] for k in range(len(B))) for j in range(len(B[0]))] for i in range(len(A))]\n",
+    "\n",
+    "# Using the built-in numpy matrix multiplication routines\n",
+    "import numpy as np\n",
+    "def matrix_multiply_numpy(A, B):\n",
+    "    return np.dot(A, B)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b0b38b84-d86c-42ed-b306-84d170baa952",
+   "metadata": {},
+   "source": [
+    "Problem 1.3.c."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "id": "2a2329c0-aa74-4814-8aa1-7b7737f3c73e",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Method: nested for loops\n",
+      "Result:  True\n",
+      "Time taken: 0.42284202575683594 seconds\n",
+      "Method: a list comprehension\n",
+      "Result:  True\n",
+      "Time taken with list comprehension: 0.3502190113067627 seconds\n",
+      "Method: the built-in numpy matrix multiplication routines\n",
+      "Result:  True\n",
+      "Time taken with NumPy: 0.0010690689086914062 seconds\n"
+     ]
+    }
+   ],
+   "source": [
+    "A     = np.random.rand(200, 200) # Create a large random matrix\n",
+    "A_inv = np.linalg.inv(A)         # Compute its inverse\n",
+    "I     = np.eye(200)              # Identity matrix for comparison\n",
+    "\n",
+    "import time\n",
+    "\n",
+    "# Test with using nested for loops\n",
+    "print('Method: nested for loops')\n",
+    "start        = time.time()\n",
+    "result_loops = matrix_multiply_loops(A.tolist(), A_inv.tolist())\n",
+    "end          = time.time()\n",
+    "print('Result: ', np.allclose(result_loops, I))\n",
+    "print(f'Time taken: {end - start} seconds')\n",
+    "\n",
+    "# Test with using a list comprehension\n",
+    "print('Method: a list comprehension')\n",
+    "start            = time.time()\n",
+    "result_list_comp = matrix_multiply_list_comp(A.tolist(), A_inv.tolist())\n",
+    "end              = time.time()\n",
+    "print('Result: ', np.allclose(result_list_comp, I))\n",
+    "print(f\"Time taken with list comprehension: {end - start} seconds\")\n",
+    "\n",
+    "# Test with using the built-in numpy matrix multiplication routines\n",
+    "print('Method: the built-in numpy matrix multiplication routines')\n",
+    "start        = time.time()\n",
+    "result_numpy = matrix_multiply_numpy(A, A_inv)\n",
+    "end          = time.time()\n",
+    "print('Result: ', np.allclose(result_numpy, I))\n",
+    "print(f\"Time taken with NumPy: {end - start} seconds\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c85da8d3-450f-49a8-b4fa-0412bc91e0c9",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}