diff --git a/geostatspy/geostats.py b/geostatspy/geostats.py
index 51103ed..033feb1 100644
--- a/geostatspy/geostats.py
+++ b/geostatspy/geostats.py
@@ -38,7 +38,7 @@ def backtr(df,vcol,vr,vrg,zmin,zmax,ltail,ltpar,utail,utpar):
     :param utpar: upper tail extrapolation parameter
     :return: TODO
     """    
-    
+    #comment
     EPSLON=1.0e-20
     nd = len(df); nt = len(vr) # number of data to transform and number of data in table
     backtr = np.zeros(nd)
@@ -2196,6 +2196,222 @@ def nscore(
 
     return ns, vr, wt_ns
 
+# *** all parameters have been copied over form kb2d with respective z parameter added if necessary
+# we can probably delete a couple of these if we don't do block kriging
+def kt3d (
+    df,
+    xcol,
+    ycol,
+    zcol,
+    vcol,
+    tmin,
+    tmax,
+    nx,
+    xmn,
+    xsiz,
+    ny,
+    ymn,
+    ysiz,
+    nz,
+    zmn,
+    zsiz,
+    nxdis,
+    nydis,
+    nzdis,
+    ndmin,
+    ndmax,
+    radius,
+    ktype,
+    skmean,
+    vario,
+):
+    UNEST = -999.
+    EPSLON = 1.0e-10
+    PMX = 9999.0    
+    MAXSAM = ndmax + 1
+    MAXDIS = nxdis * nydis * nzdis
+    MAXKD = MAXSAM + 1
+    MAXKRG = MAXKD * MAXKD
+
+# Allocate the needed memory:  
+    
+    # used for block kriging displacements
+    xdb = np.zeros(MAXDIS)
+    ydb = np.zeros(MAXDIS)
+    zdb = np.zeros(MAXDIS)
+    # temp arrays to hold nearest neighbor search results
+    xa = np.zeros(MAXSAM)
+    ya = np.zeros(MAXSAM)
+    za = np.zeros(MAXSAM)
+    vra = np.zeros(MAXSAM)
+    # result grid meshes for estimated values and estimation variance
+    kmap = np.zeros((nx,ny,nz))
+    vmap = np.zeros((nx,ny,nz))
+
+    # trim values outside tmin and tmax
+    df_extract = df.loc[(df[vcol] >= tmin) & (df[vcol] <= tmax)]   
+    nd = len(df_extract)
+    ndmax = min(ndmax,nd)
+    x = df_extract[xcol].values
+    y = df_extract[ycol].values
+    z = df_extract[zcol].values
+    vr = df_extract[vcol].values
+
+# set up tree for nearest neighbor search
+    dp = list((z[i], y[i], x[i]) for i in range(0,nd))
+    data_locs = np.column_stack((z, y, x))
+    tree = sp.cKDTree(data_locs, leafsize=16, compact_nodes=True, copy_data=False, balanced_tree=True) #why this error?
+
+# Summary statistics for the data after trimming
+    avg = vr.mean()
+    stdev = vr.std()
+    ss = stdev**2.0 # variance
+    vrmin = vr.min()
+    vrmax = vr.max()
+
+# load the variogram
+    nst = vario['nst'] # num structures
+    cc = np.zeros(nst) # variance contribution
+    aa = np.zeros(nst) # major range
+    it = np.zeros(nst) # type vario structure
+    aa = np.zeros(nst) # major range
+    ang_azi = np.zeros(nst) # azimuth
+    ang_dip = np.zeros(nst) # dip
+    anis = np.zeros(nst) # anistropy ratio w/ minor range
+    anis_v = np.zeros(nst) # anistropy ratio w/ vertical range
+    
+    # only works w/ nst == 1 or nst == 2
+    c0 = vario['nug'] 
+    it[0] = vario['it1'] 
+    cc[0] = vario['cc1'] 
+    ang_azi[0] = vario['azi1']
+    ang_dip[0] = vario['dip1']
+    aa[0] = vario['hmax']
+    anis[0] = vario['hmed1']/vario['hmax1']
+    anis_v[0] = vario['hmin1']/vario['hmax1']
+    if nst == 2:
+        cc[1] = vario['cc2']
+        it[1] = vario['it2']
+        ang_azi[1] = vario['azi2']
+        ang_dip[1] = vario['dip1']
+        aa[1] = vario['hmax']
+        anis[1] = vario['hmed2']/vario['hmax2']
+        anis_v[1] = vario['hmin2']/vario['hmax2']
+
+    rotmat, maxcov = setup_rotmat_3D(c0, nst, it, cc, ang_azi, ang_dip, PMX)
+    cbb = maxcov
+    unbias = maxcov
+
+    # limited to simple kriging for now
+    ktype = 0
+
+    # TODO Set up discretization points per block (block kriging)
+
+    # main loop over all points:
+    nk = 0
+    ak = 0.0
+    vk = 0.0
+    for iz in range(0,nz):
+        zloc = zmn + (iz-0)*zsiz 
+        for iy in range(0,ny):
+            yloc = ymn + (iy-0)*ysiz
+            for ix in range (0,nz):     #triple nested for loop that loops through points and adds to matrices
+                xloc = xmn + (ix-0)*xsiz
+                current_node = (zloc, yloc,xloc) # xloc, yloc, zloc centroid of current cell
+
+                # Find the nearest samples within each octant:
+                na = -1   # accounting for 0 as first index
+                dist, nums = tree.query(current_node,ndmax) # use kd tree for fast nearest data search
+                # remove any data outside search radius
+                na = len(dist)
+                nums = nums[dist<radius]
+                dist = dist[dist<radius] 
+                na = len(dist)        
+
+                # Is there enough samples?
+                if na + 1 < ndmin:   # accounting for min index of 0
+                    est  = UNEST
+                    estv = UNEST
+                    # not enough samples, mark as "unestimated"
+                    print('UNEST at ' + str(ix) + ',' + str(iy) + ','+str(iz)) 
+                else:
+                    # Put coordinates and values of neighborhood samples into xa,ya,za,vra:
+                    for ia in range(0,na):
+                        jj = int(nums[ia])
+                        xa[ia]  = x[jj] #why this error?
+                        ya[ia]  = y[jj]
+                        za[ia]  = z[jj]
+                        vra[ia] = vr[jj]
+                    
+                    # Handle the situation of only one sample:
+                    if na == 0:  # accounting for min index of 0 - one sample case na = 0
+                        cb1 = cova3(xa[0],ya[0],za[0],xa[0],ya[0],za[0],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov)
+
+                        xx  = xa[0] - xloc
+                        yy  = ya[0] - yloc
+                        zz  = za[0] - zloc
+
+                    # Establish Right Hand Side Covariance:
+                        cb = cova3(xx,yy,zz,xdb[0],ydb[0],zdb[i],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov)
+                        
+                        if ktype == 0:
+                            s[0] = cb/cbb
+                            est  = s[0]*vra[0] + (1.0-s[0])*skmean
+                            estv = cbb - s[0] * cb
+                        # else:
+                        #     est  = vra[0]
+                        #     estv = cbb - 2.0*cb + cb1
+                    else:
+                    # Solve the Kriging System with more than one sample:
+                        neq = na + ktype # accounting for first index of 0
+                        # print('NEQ' + str(neq))
+                        nn  = (neq + 1)*neq/2   
+
+                        # Set up kriging matrices:
+                        a = np.zeroes((na,na))
+                        r = np.zeroes(na)
+                        rr = np.zeroes(na)
+
+                        # Establish Left Hand Side Covariance Matrix:
+
+                        for j in range(0,na):
+                            for i in range(0,na): 
+                                # add covariance into left-hand square matrix
+                                a[j][i] = cova3(xa[i],ya[i],za[i],xa[j],ya[j],za[j],nst,c0,PMX,cc,aa,it,ang_azi,anis,rotmat,maxcov) 
+                            # add covariance into right-hand column matrix
+                            r[j] = cova3(xloc,yloc,zloc,xa[j],ya[j],za[j],nst,c0,PMX,cc,aa,it,ang_azi,anis,rotmat,maxcov) 
+
+                        # solve the kriging system
+                        s = ksol_numpy(neq,a,r)
+                        ising = 0
+
+                    # Write a warning if the matrix is singular:
+                        if ising != 0:
+                            print('WARNING KT3D: singular matrix')
+                            est  = UNEST
+                            estv = UNEST
+                        else:
+                    # Compute the estimate and the kriging variance:
+                            est  = 0.0
+                            estv = cbb
+                            sumw = 0.0
+                            # if ktype == 1: 
+                            #     estv = estv - (s[na])*unbias
+                            for i in range(0,na):                          
+                                sumw = sumw + s[i]
+                                est  = est  + s[i]*vra[i]
+                                estv = estv - s[i]*rr[i]
+                            if ktype == 0: 
+                                est = est + (1.0-sumw)*skmean
+                # add estimation variance and estimated value to maps
+                kmap[nz-iz-1, ny-iy-1,ix] = est
+                vmap[nz-iz-1, ny-iy-1,ix] = estv
+                if est > UNEST:
+                    nk = nk + 1
+                    ak = ak + est
+                    vk = vk + est*est
+
+    return kmap, vmap   
 
 def kb2d(
     df,
@@ -2258,8 +2474,11 @@ def kb2d(
     
 # load the variogram
     nst = vario['nst']
-    cc = np.zeros(nst); aa = np.zeros(nst); it = np.zeros(nst)
-    ang = np.zeros(nst); anis = np.zeros(nst)
+    cc = np.zeros(nst);
+    aa = np.zeros(nst); 
+    it = np.zeros(nst) 
+    ang = np.zeros(nst); 
+    anis = np.zeros(nst) 
     
     c0 = vario['nug']; 
     cc[0] = vario['cc1']; it[0] = vario['it1']; ang[0] = vario['azi1']; 
@@ -2267,6 +2486,7 @@ def kb2d(
     if nst == 2:
         cc[1] = vario['cc2']; it[1] = vario['it2']; ang[1] = vario['azi2']; 
         aa[1] = vario['hmaj2']; anis[1] = vario['hmin2']/vario['hmaj2'];
+
     
 # Allocate the needed memory:   
     xdb = np.zeros(MAXDIS)
@@ -2299,7 +2519,7 @@ def kb2d(
 # Summary statistics for the data after trimming
     avg = vr.mean()
     stdev = vr.std()
-    ss = stdev**2.0
+    ss = stdev**2.0 # variance
     vrmin = vr.min()
     vrmax = vr.max()
 
@@ -2324,6 +2544,7 @@ def kb2d(
             xdb[i] = xloc
             ydb[i] = yloc
 
+
 # Initialize accumulators:
     cbb  = 0.0
     rad2 = radius*radius
@@ -2353,7 +2574,7 @@ def kb2d(
         yloc = ymn + (iy-0)*ysiz  
         for ix in range(0,nx):
             xloc = xmn + (ix-0)*xsiz
-            current_node = (yloc,xloc)
+            current_node = (yloc,xloc) # xloc, yloc centroid of cell
         
 # Find the nearest samples within each octant: First initialize
 # the counter arrays:
@@ -2421,7 +2642,7 @@ def kb2d(
                         for i in range(0,na):  # was j - want full matrix                    
                             iin = iin + 1
                             a[iin] = cova2(xa[i],ya[i],xa[j],ya[j],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov) 
-                        if ktype == 1:
+                        if ktype == 1: 
                             iin = iin + 1
                             a[iin] = unbias
                         xx = xa[j] - xloc
@@ -2429,8 +2650,9 @@ def kb2d(
 
 # Establish Right Hand Side Covariance:
                         if ndb <= 1:
+                            # ndb <= 1 -> use single value
                             cb = cova2(xx,yy,xdb[0],ydb[0],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov)
-                        else:
+                        else: 
                             cb  = 0.0
                             for j1 in range(0,ndb):    
                                 cb = cb + cova2(xx,yy,xdb[j1],ydb[j1],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov)
@@ -4470,9 +4692,6 @@ def cova3(x1, y1, z1, x2, y2, z2, nst, c0, pmx, cc, aa, it, anis, anis_v, rotmat
             cova3_ = cova3_ + cov1
     return cova3_	
 	
-
-
-	
 	
 def gamv_3D(
     df,
diff --git a/geostatspy/testing.ipynb b/geostatspy/testing.ipynb
new file mode 100644
index 0000000..9844bc2
--- /dev/null
+++ b/geostatspy/testing.ipynb
@@ -0,0 +1,389 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 102,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "testing the kt3d function using some data from Geostats Guy\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(\"testing the kt3d function using some data from Geostats Guy\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 103,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import geostatspy.GSLIB as GSLIB    \n",
+    "# GSLIB utilies, visualization and wrapper\n",
+    "import geostatspy.geostats as geostats\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import matplotlib.pyplot as plt\n",
+    "import scipy.spatial as sp  # for fast nearest neighbor search"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 104,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "\n",
+    "\n",
+    "# *** all parameters have been copied over form kb2d with respective z parameter added if necessary\n",
+    "# we can probably delete a couple of these if we don't do block kriging\n",
+    "def kt3d (\n",
+    "    df,\n",
+    "    xcol,\n",
+    "    ycol,\n",
+    "    zcol,\n",
+    "    vcol,\n",
+    "    tmin,\n",
+    "    tmax,\n",
+    "    nx,\n",
+    "    xmn,\n",
+    "    xsiz,\n",
+    "    ny,\n",
+    "    ymn,\n",
+    "    ysiz,\n",
+    "    nz,\n",
+    "    zmn,\n",
+    "    zsiz,\n",
+    "    nxdis,\n",
+    "    nydis,\n",
+    "    nzdis,\n",
+    "    ndmin,\n",
+    "    ndmax,\n",
+    "    radius,\n",
+    "    ktype,\n",
+    "    skmean,\n",
+    "    vario,\n",
+    "):\n",
+    "    UNEST = -999.\n",
+    "    EPSLON = 1.0e-10\n",
+    "    # VERSION = 2.907\n",
+    "    # first = True\n",
+    "    PMX = 9999.0    \n",
+    "    MAXSAM = ndmax + 1\n",
+    "    MAXDIS = nxdis * nydis * nzdis\n",
+    "    MAXKD = MAXSAM + 1\n",
+    "    MAXKRG = MAXKD * MAXKD\n",
+    "\n",
+    "# Allocate the needed memory:   \n",
+    "    xdb = np.zeros(MAXDIS)\n",
+    "    ydb = np.zeros(MAXDIS)\n",
+    "    zdb = np.zeros(MAXDIS)\n",
+    "    xa = np.zeros(MAXSAM)\n",
+    "    ya = np.zeros(MAXSAM)\n",
+    "    za = np.zeros(MAXSAM)\n",
+    "    vra = np.zeros(MAXSAM)\n",
+    "    dist = np.zeros(MAXSAM)\n",
+    "    nums = np.zeros(MAXSAM)\n",
+    "    r = np.zeros(MAXKD)\n",
+    "    rr = np.zeros(MAXKD)\n",
+    "    s = np.zeros(MAXKD)\n",
+    "    a = np.zeros(MAXKRG)\n",
+    "    kmap = np.zeros((nx,ny,nz))\n",
+    "    vmap = np.zeros((nx,ny,nz))\n",
+    "\n",
+    "    # trim values outside tmin and tmax\n",
+    "    df_extract = df.loc[(df[vcol] >= tmin) & (df[vcol] <= tmax)]  \n",
+    "    assert df[vcol].max() <= tmax and df[vcol].min() >= tmin\n",
+    "    nd = len(df_extract)\n",
+    "    ndmax = min(ndmax,nd)\n",
+    "    x = df_extract[xcol].values\n",
+    "    y = df_extract[ycol].values\n",
+    "    z = df_extract[zcol].values\n",
+    "    vr = df_extract[vcol].values\n",
+    "\n",
+    "\n",
+    "    # Allocate the needed memory:   \n",
+    "    xdb = np.zeros(MAXDIS)\n",
+    "    ydb = np.zeros(MAXDIS)\n",
+    "    zdb = np.zeros(MAXDIS)\n",
+    "    xa = np.zeros(MAXSAM)\n",
+    "    ya = np.zeros(MAXSAM)\n",
+    "    za = np.zeros(MAXSAM)\n",
+    "    vra = np.zeros(MAXSAM)\n",
+    "    dist = np.zeros(MAXSAM)\n",
+    "    nums = np.zeros(MAXSAM)\n",
+    "    r = np.zeros(MAXKD)\n",
+    "    rr = np.zeros(MAXKD)\n",
+    "    s = np.zeros(MAXKD)\n",
+    "    a = np.zeros(MAXKRG)\n",
+    "    kmap = np.zeros((nx,ny,nz))\n",
+    "    vmap = np.zeros((nx,ny,nz)) \n",
+    "\n",
+    "# set up tree for nearest neighbor search\n",
+    "    dp = list((z[i], y[i], x[i]) for i in range(0,nd))\n",
+    "    data_locs = np.column_stack((z, y, x))\n",
+    "    tree = sp.cKDTree(data_locs, leafsize=16, compact_nodes=True, copy_data=False, balanced_tree=True) #why this error?\n",
+    "\n",
+    "# Summary statistics for the data after trimming\n",
+    "    avg = vr.mean()\n",
+    "    stdev = vr.std()\n",
+    "    ss = stdev**2.0 # variance\n",
+    "    vrmin = vr.min()\n",
+    "    vrmax = vr.max()\n",
+    "\n",
+    "# load the variogram\n",
+    "    nst = vario['nst'] # num structures\n",
+    "    cc = np.zeros(nst) # variance contribution\n",
+    "    aa = np.zeros(nst) # major range\n",
+    "    it = np.zeros(nst) # type vario structure\n",
+    "    aa = np.zeros(nst) # major range\n",
+    "    ang_azi = np.zeros(nst) # azimuth\n",
+    "    ang_dip = np.zeros(nst) # dip\n",
+    "    anis = np.zeros(nst) # anistropy ratio w/ minor range\n",
+    "    anis_v = np.zeros(nst) # anistropy ratio w/ vertical range\n",
+    "    \n",
+    "    # only works w/ nst == 1 or nst == 2\n",
+    "    print(\"vario\")\n",
+    "    print(vario)\n",
+    "    c0 = vario['nug'] \n",
+    "    it[0] = vario['it1'] \n",
+    "    cc[0] = vario['cc1'] \n",
+    "    ang_azi[0] = vario['azi1']\n",
+    "    ang_dip[0] = vario['dip1']\n",
+    "    aa[0] = vario['hmax']\n",
+    "    anis[0] = vario['hmed1']/vario['hmax1']\n",
+    "    anis_v[0] = vario['hmin1']/vario['hmax1']\n",
+    "    if nst == 2:\n",
+    "        cc[1] = vario['cc2']\n",
+    "        it[1] = vario['it2']\n",
+    "        ang_azi[1] = vario['azi2']\n",
+    "        ang_dip[1] = vario['dip1']\n",
+    "        aa[1] = vario['hmax']\n",
+    "        anis[1] = vario['hmed2']/vario['hmax2']\n",
+    "        anis_v[1] = vario['hmin2']/vario['hmax2']\n",
+    "\n",
+    "    #rotmat, maxcov = setup_rotmat_3D(c0, nst, it, cc, ang_azi, ang_dip, PMX)\n",
+    "    #cbb = maxcov\n",
+    "\n",
+    "    # TODO Set up discretization points per block (bock kriging)\n",
+    "\n",
+    "    # need to set up numPy cube grid\n",
+    "    cubeGrid = np.ndarray(shape=(nx,ny,nz), dtype=float, order='C') \n",
+    "    print(cubeGrid)\n",
+    "    # main loop over all points:\n",
+    "    nk = 0\n",
+    "    ak = 0.0\n",
+    "    vk = 0.0\n",
+    "    for iz in range(0,nz):\n",
+    "        zloc = zmn + (iz-0)*zsiz \n",
+    "        for iy in range(0,ny):\n",
+    "            yloc = ymn + (iy-0)*ysiz\n",
+    "            for ix in range (0,nz):     #triple nested for loop that loops through points and adds to matrices\n",
+    "                xloc = xmn + (ix-0)*xsiz\n",
+    "                current_node = (zloc, yloc,xloc) # xloc, yloc, zloc centroid of current cell\n",
+    "\n",
+    "\n",
+    "    # Find the nearest samples within each octant: First initialize\n",
+    "# the counter arrays:\n",
+    "                na = -1   # accounting for 0 as first index\n",
+    "                dist.fill(1.0e+20)\n",
+    "                nums.fill(-1)\n",
+    "                dist, nums = tree.query(current_node,ndmax) # use kd tree for fast nearest data search\n",
+    "                # remove any data outside search radius\n",
+    "                na = len(dist)\n",
+    "                nums = nums[dist<radius]\n",
+    "                dist = dist[dist<radius] \n",
+    "                na = len(dist)        \n",
+    "\n",
+    "                # Is there enough samples?\n",
+    "                if na + 1 < ndmin:   # accounting for min index of 0\n",
+    "                    est  = UNEST\n",
+    "                    estv = UNEST\n",
+    "                    print('UNEST at ' + str(ix) + ',' + str(iy) + ','+str(iz)) \n",
+    "                else:\n",
+    "\n",
+    "                    # Put coordinates and values of neighborhood samples into xa,ya,za,vra:\n",
+    "                    for ia in range(0,na):\n",
+    "                        jj = int(nums[ia])\n",
+    "                        xa[ia]  = x[jj] #why this error?\n",
+    "                        ya[ia]  = y[jj]\n",
+    "                        za[ia]  = z[jj]\n",
+    "                        vra[ia] = vr[jj] #what is vra??\n",
+    "                    \n",
+    "                    # Handle the situation of only one sample:\n",
+    "                    if na == 0:  # accounting for min index of 0 - one sample case na = 0\n",
+    "                        cb1 = geostats.cova3(xa[0],ya[0],za[0],xa[0],ya[0],za[0],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov)\n",
+    "\n",
+    "                        xx  = xa[0] - xloc\n",
+    "                        yy  = ya[0] - yloc\n",
+    "                        zz  = za[0] - zloc\n",
+    "\n",
+    "                # Is there enough samples?\n",
+    "                if na + 1 < ndmin:   # accounting for min index of 0\n",
+    "                    est  = UNEST\n",
+    "                    estv = UNEST\n",
+    "                    print('UNEST at ' + str(ix) + ',' + str(iy) + ','+str(iz)) \n",
+    "                else:\n",
+    "                    # Put coordinates and values of neighborhood samples into xa,ya,za,vra:\n",
+    "                    for ia in range(0,na):\n",
+    "                        jj = int(nums[ia])\n",
+    "                        xa[ia]  = x[jj] #why this error?\n",
+    "                        ya[ia]  = y[jj]\n",
+    "                        za[ia]  = z[jj]\n",
+    "                        vra[ia] = vr[jj] #what is vra??\n",
+    "                    \n",
+    "                    # Handle the situation of only one sample:\n",
+    "                    if na == 0:  # accounting for min index of 0 - one sample case na = 0\n",
+    "                        cb1 = geostats.cova3(xa[0],ya[0],za[0],xa[0],ya[0],za[0],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov)\n",
+    "\n",
+    "                        xx  = xa[0] - xloc\n",
+    "                        yy  = ya[0] - yloc\n",
+    "                        zz  = za[0] - zloc\n",
+    "\n",
+    "                    # Establish Right Hand Side Covariance:\n",
+    "                        if ndb <= 1:\n",
+    "                            cb = geostats.cova3(xx,yy,zz,xdb[0],ydb[0],zdb[i],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov)\n",
+    "                        else:    \n",
+    "                            cb  = 0.0\n",
+    "                            for i in range(0,ndb):                  \n",
+    "                                cb = cb + geostats.cova3(xx,yy,zz,xdb[i],ydb[i],zdb[i],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov)\n",
+    "                                dx = xx - xdb(i)\n",
+    "                                dy = yy - ydb(i)\n",
+    "                                dz = zz - zdb(i)\n",
+    "                                if (dx*dx+dy*dy+dz*dz) < EPSLON:\n",
+    "                                    cb = cb - c0\n",
+    "                                cb = cb / real(ndb)\n",
+    "                        if ktype == 0:\n",
+    "                            s[0] = cb/cbb\n",
+    "                            est  = s[0]*vra[0] + (1.0-s[0])*skmean\n",
+    "                            estv = cbb - s[0] * cb\n",
+    "                        else:\n",
+    "                            est  = vra[0]\n",
+    "                            estv = cbb - 2.0*cb + cb1\n",
+    "                    else:\n",
+    "                    # Solve the Kriging System with more than one sample:\n",
+    "                        neq = na + ktype # accounting for first index of 0\n",
+    "                        # print('NEQ' + str(neq))\n",
+    "                        nn  = (neq + 1)*neq/2   \n",
+    "\n",
+    "                        # Set up kriging matrices:\n",
+    "                        a = np.zeroes((na,na))\n",
+    "                        r = np.zeroes(na)\n",
+    "\n",
+    "                        # Establish Left Hand Side Covariance Matrix:\n",
+    "\n",
+    "                        for j in range(0,na):\n",
+    "                            for i in range(0,na): #loops through and adds covariances into matrix\n",
+    "                                a[j][i] = geostats.cova3(xa[i],ya[i],za[i],xa[j],ya[j],za[j],nst,c0,PMX,cc,aa,it,ang_azi,anis,rotmat,maxcov) \n",
+    "                            r[j] = geostats.cova3(xloc,yloc,zloc,xa[j],ya[j],za[j],nst,c0,PMX,cc,aa,it,ang_azi,anis,rotmat,maxcov) \n",
+    "\n",
+    "# Write a warning if the matrix is singular:\n",
+    "                        if ising != 0:\n",
+    "                            print('WARNING KT3D: singular matrix')\n",
+    "                            print('              for block' + str(ix) + ',' + str(iy) + ',' + str(iz))\n",
+    "                            est  = UNEST\n",
+    "                            estv = UNEST\n",
+    "                        else:\n",
+    "\n",
+    "# Compute the estimate and the kriging variance:\n",
+    "                            est  = 0.0\n",
+    "                            estv = cbb\n",
+    "                            sumw = 0.0\n",
+    "                            if ktype == 1: \n",
+    "                                estv = estv - (s[na])*unbias\n",
+    "                            for i in range(0,na):                          \n",
+    "                                sumw = sumw + s[i]\n",
+    "                                est  = est  + s[i]*vra[i]\n",
+    "                                estv = estv - s[i]*rr[i]\n",
+    "                            if ktype == 0: \n",
+    "                                est = est + (1.0-sumw)*skmean\n",
+    "                kmap[nz-iz-1, ny-iy-1,ix] = est\n",
+    "                vmap[nz-iz-1, ny-iy-1,ix] = estv\n",
+    "                if est > UNEST:\n",
+    "                    nk = nk + 1\n",
+    "                    ak = ak + est\n",
+    "                    vk = vk + est*est\n",
+    "# TODO Set up kriging matrices:\n",
+    "\n",
+    "    return kmap, vmap   \n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 105,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "skmean_por = 0.10; skmean_perm = 65.0      # simple kriging mean (used if simple kriging is selected below)\n",
+    "ktype = 1                                  # kriging type, 0 - simple, 1 - ordinary\n",
+    "radius = 100                               # search radius for neighbouring data\n",
+    "nxdis = 1; nydis = 1 ; nzdis = 1         # number of grid discretizations for block kriging (not tested)\n",
+    "ndmin = 0; ndmax = 10                      # minimum and maximum data for an estimate\n",
+    "tmin = 0.0                                 # minimum property value\n",
+    "xsiz = 10; ysiz = 10 ; zsiz = 10                   # cell size\n",
+    "nx = 100; ny = 100; nz = 100                      # number of cells\n",
+    "xmn = 5; ymn = 5 ; zmn = 5                       # grid origin, location center of lower left cell\n",
+    "por_vario = GSLIB.make_variogram(nug=0.0,nst=1,it1=1,cc1=1.0,azi1=0,hmaj1=100,hmin1=100) # porosity variogram\n",
+    "por_vario[\"dip1\"] = 0\n",
+    "por_vario[\"hmax\"] = 0\n",
+    "por_vario[\"hmed1\"] = 0\n",
+    "por_vario[\"hmax1\"] = 1\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 106,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df = pd.read_csv(\"completions_production.csv\")\n",
+    "por_kmap, por_vmap = kt3d ( df, 'X', 'Y', 'Z', 'Porosity',  tmin, tmax, nx,\n",
+    "    xmn,  xsiz,  ny,  ymn,   ysiz,  nz,  zmn,  zsiz,  nxdis, nydis,  nzdis, ndmin, ndmax, radius, ktype, skmean_por, por_vario)\n",
+    "print(por_kmap)\n",
+    "print(por_vmap)\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "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.7.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}