diff --git a/python/cosmotool/__init__.py b/python/cosmotool/__init__.py
index 70fc1f1..bf0df59 100644
--- a/python/cosmotool/__init__.py
+++ b/python/cosmotool/__init__.py
@@ -1 +1,2 @@
 from _cosmotool import *
+from borg import read_borg_vol
diff --git a/python/cosmotool/borg.py b/python/cosmotool/borg.py
new file mode 100644
index 0000000..5ac0f7b
--- /dev/null
+++ b/python/cosmotool/borg.py
@@ -0,0 +1,265 @@
+### 
+### BORG code is from J. Jasche
+###
+import StringIO
+import numpy as np
+from numpy import *
+import os.path
+import array
+import glob
+
+class BorgVolume(object):
+
+    def __init__(self, density, ranges):
+        self._density = density
+        self._ranges = ranges
+
+    property density:
+        def __get__(self):
+            return self._density
+
+    property ranges:
+        def __get__(self):
+            return self._ranges
+
+def build_filelist(fdir):
+    #builds list of all borg density fields which may be distributed over several directories
+    
+    fname_0=glob.glob(fdir[0]+'initial_density_*')
+    fname_1=glob.glob(fdir[0]+'final_density_*')
+       
+    fdir=fdir[1:] #eliminate first element
+       
+    for fd in fdir:
+	fname_0=fname_0+glob.glob(fd+'initial_density_*')
+	fname_1=fname_1+glob.glob(fd+'final_density_*')
+
+	   
+    return fname_0, fname_1
+
+def read_borg_vol(BORGFILE):
+    """ Reading routine for BORG data
+    """
+
+    openfile=open(BORGFILE,'rb')
+
+    period=0
+    N0=0
+    N1=0
+    N2=0
+
+    xmin=0
+    xmax=0
+
+    ymin=0
+    ymax=0
+
+    zmin=0
+    zmax=0
+    nlines=0
+
+    while True:
+      line=openfile.readline()
+      s=line.rstrip('\n')
+      r=s.rsplit(' ')
+    
+      if size(r)==5 : 
+	if r[0] =="define":
+	  if r[1]=="Lattice" : N0=int(r[2])
+	  N1=int(r[3])
+	  N2=int(r[4])
+         
+      
+      if size(r)==11 : 
+	if r[4] =="BoundingBox": xmin=float(r[5])
+	xmax=float(r[6])
+	ymin=float(r[7])
+	ymax=float(r[8])
+	zmin=float(r[9])
+	zmax=float(r[10].rstrip(','))
+  
+      if r[0]=='@1': break
+  
+    ranges=[]
+    ranges.append(xmin)
+    ranges.append(xmax)
+    ranges.append(ymin)
+    ranges.append(ymax)
+    ranges.append(zmin)
+    ranges.append(zmax)
+    
+    #now read data
+    data=np.fromfile(openfile, '>f4')
+    data=data.reshape(N2,N0,N1)
+    return BorgVolume(data,ranges)
+
+def read_spec( fname ):
+    """ Reading routine for ARES spectrum samples
+    """
+    x,y=np.loadtxt( fname ,usecols=(0,1),unpack=True)
+
+    return x , y
+    
+        
+def read_bias_nmean( fname ):
+    """ Reading routine for ARES bias data
+    """
+    x,b0,b1,nmean=np.loadtxt( fname ,usecols=(0,1,2,3),unpack=True)
+
+    return x , b0, b1, nmean
+
+def read_nmean( fname ):
+    """ Reading routine for BORG bias data
+    """
+    x,nmean=np.loadtxt( fname ,usecols=(0,1),unpack=True)
+
+    return x, nmean    
+    
+def get_grid_values(xx,data, ranges):
+    """ return values at grid positions
+    """
+    xmin=ranges[0]
+    xmax=ranges[1]
+    
+    ymin=ranges[2]
+    ymax=ranges[3]
+    
+    zmin=ranges[4]
+    zmax=ranges[5]
+    
+    Lx= xmax-xmin
+    Ly= ymax-ymin
+    Lz= zmax-zmin
+        
+    Nx=shape(data)[0]
+    Ny=shape(data)[1]
+    Nz=shape(data)[2]
+    
+    dx=Lx/float(Nx)
+    dy=Ly/float(Ny)
+    dz=Lz/float(Nz)
+    
+    idx=(xx[:,0]-xmin)/dx
+    idy=(xx[:,1]-ymin)/dz
+    idz=(xx[:,2]-zmin)/dy
+    
+    idx=idx.astype(int)
+    idy=idy.astype(int)
+    idz=idz.astype(int)
+    
+    idflag=np.where( (idx>-1)*(idx<Nx)*(idy>-1)*(idy<Ny)*(idz>-1)*(idz<Nz) )
+        
+    flag=[False]*len(xx)
+    vals=[-999.]*len(xx)
+    
+    flag=np.array(flag)
+    vals=np.array(vals)
+    flag[idflag]=True
+    vals[idflag]=data[idx[idflag],idy[idflag],idz[idflag]]
+        
+    return vals,flag    
+    
+def get_mean_density(fdir, smin, step):
+    """ estimate ensemble mean
+    """
+    print '-'*60
+    print 'Get 3D ensemble mean density field'
+    print '-'*60
+
+    fname0 = fdir + 'initial_density_'+str(0)+'.dat'
+    fname1 = fdir + 'final_density_'+str(0)+'.dat'
+        
+    MEAN0,ranges=read_borg_vol(fname0)
+    MEAN0=MEAN0*0.;
+    VAR0=copy(MEAN0)
+    MEAN1=copy(MEAN0)
+    VAR1=copy(MEAN0)
+    norm=0.
+    
+    idat=smin
+    
+    fname0 = fdir + 'initial_density_'+str(idat)+'.dat'
+    fname1 = fdir + 'final_density_'+str(idat)+'.dat'
+        #and (idat<smin+1000)
+    while((os.path.exists(fname0))):
+	    auxdata0,auxranges0=read_borg_vol(fname0)
+	    auxdata1,auxranges1=read_borg_vol(fname1)
+	    auxx0=auxdata0
+	    auxx1=auxdata1
+	    MEAN0+=auxx0
+	    VAR0+=auxx0**2
+	    MEAN1+=auxx1
+	    VAR1+=auxx1**2
+	    
+	    norm+=1
+	    idat+=step
+	    fname0 = fdir + 'initial_density_'+str(idat)+'.dat'
+	    fname1 = fdir + 'final_density_'+str(idat)+'.dat'
+	    del auxranges0
+	    del auxdata0
+	    del auxranges1
+	    del auxdata1
+	    
+    MEAN0/=norm
+    VAR0/=norm
+    VAR0-=MEAN0**2
+    VAR0=sqrt(fabs(VAR0))
+    
+    MEAN1/=norm
+    VAR1/=norm
+    VAR1-=MEAN1**2
+    VAR1=sqrt(fabs(VAR1))
+    
+    
+    return MEAN0,VAR0,MEAN1,VAR1,ranges
+
+def get_mean_density_fdir(fdir,init,steps):
+    """ estimate ensemble mean
+    """
+    print '-'*60
+    print 'Get 3D ensemble mean density field'
+    print '-'*60
+
+    fname0,fname1=build_filelist(fdir)
+    
+    fname0=fname0[init::steps]
+    fname1=fname1[init::steps]
+    
+    MEAN0,ranges=read_borg_vol(fname0[0])
+    MEAN0=MEAN0*0.;
+    VAR0=copy(MEAN0)
+    MEAN1=copy(MEAN0)
+    VAR1=copy(MEAN0)
+    norm0=0.
+    norm1=0.
+    
+    for fn in fname0:
+	auxdata0,auxranges0=read_borg_vol(fn)
+	MEAN0+=auxdata0
+	VAR0+=auxdata0**2.
+	norm0+=1.
+	del auxranges0
+	del auxdata0
+
+    for fn in fname1:
+	auxdata1,auxranges1=read_borg_vol(fn)
+	MEAN1+=auxdata1
+	VAR1+=auxdata1**2.
+	norm1+=1.
+	del auxranges1
+	del auxdata1
+	
+    MEAN0/=norm0
+    VAR0/=norm0
+    VAR0-=MEAN0**2
+    VAR0=sqrt(fabs(VAR0))
+    
+    MEAN1/=norm1
+    VAR1/=norm1
+    VAR1-=MEAN1**2
+    VAR1=sqrt(fabs(VAR1))
+    
+    
+    return MEAN0,VAR0,MEAN1,VAR1,ranges    
+    
+