diff --git a/src/gray.f b/src/gray.f
index f2c08e2..fa226e7 100644
--- a/src/gray.f
+++ b/src/gray.f
@@ -505,6 +505,7 @@ c
      .'N Npl ki alpha tau Pt Jcd dIds nh iohkw index_rt'
       write(8,*) ' #istep j k xt yt zt rt psin'
       write(9,*) ' #istep j k xt yt zt rt psin'
+      write(82,*) ' #istep j k xt yt zwspl zwparab'
       write(17,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1'
       write(33,*) ' #i j k sst x y R z psi tauv Npl alpha index_rt'
       write(12,*) ' #i sst w1 w2 phiw rci1 rci2 phir errw errr '//
@@ -593,6 +594,8 @@ c
       common/mode/sox
       common/angles/alpha0,beta0
       common/scal/iscal
+      
+      common/waist/w0csi,w0eta
 c
       open(2,file='gray.data',status= 'unknown')
 c
@@ -5669,6 +5672,22 @@ c     gg=F(u)/u with F(u) as in Cohen paper
       complex*16 ui,aac,bbc,ccc,ddc,aak,bbk,cck
       parameter(ui=(0.0d0,1.0d0))
       dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
+      parameter(nrmax=(jmx-1)*kmx+1)
+      dimension xtiv(nrmax),ytiv(nrmax),zwjv(nrmax),w(nrmax)
+c      parameter(nxmax=2*jmx-1)
+      parameter(nxmax=2*kmx)
+      parameter(kspl=3,nxest=nxmax+4)
+      parameter(km=kspl+1,nu=nxest-km)
+      parameter(ne=nxest,kb1=kspl*nu+km,kb2=kb1+nu-kspl)
+      parameter(lwrk1=nu*nu*(2+kb1+kb2)+
+     .                2*(2*nu+km*(nrmax+ne)+ne-2*kspl)+kb2+1)
+      parameter(lwrk2=nu*nu*(kb2+1)+kb2)
+      parameter(kwrk=nrmax+(nxest-kspl-km)*(nxest-kspl-km))
+      dimension xgridv(nxmax),ygridv(nxmax)
+      dimension zwint(nxmax*nxmax),ccexp(nxmax*nxmax)
+      dimension tx(nxest),ty(nxest),wrk1(lwrk1),wrk2(lwrk2),iwrk(kwrk)
+      parameter(lwrkbsp=8*nxmax,liwrkbsp=2*nxmax)
+      dimension wrkbsp(lwrkbsp),iwrkbsp(liwrkbsp)
 c
       common/nray/nrayr,nrayth
       common/wrk/ywrk,ypwrk
@@ -5678,6 +5697,7 @@ c
       common/ss/st
       common/dnpar/dnpara
       common/atjki/tauv,alphav
+      common/waist/w0csi,w0eta
 c
       x4m=0.0d0
       x3ym=0.0d0
@@ -5692,6 +5712,9 @@ c
       y2zwm=0.0d0
       z2wm=0.0d0
       z2rm=0.0d0
+c initialize grid dimension for spline interpolation
+      xmaxgrid=2*max(w0csi,w0eta)
+      iray=0
 c
       do j=1,nrayr     
         kktx=nrayth
@@ -5699,7 +5722,6 @@ c
         do k=1,kktx
           zwj=(dble(j-1)*rwmax/dble(nrayr-1))**2+
      .    0.5*(tauv(j,k,istep)-tauv(1,1,istep))
-  
 c
           dx=ywrk(1,j,k)-ywrk(1,1,1)
           dy=ywrk(2,j,k)-ywrk(2,1,1)
@@ -5723,7 +5745,14 @@ c
           xti= dx*csps1-dy*snps1
           yti=(dx*snps1+dy*csps1)*csth1-dz*snth1
           zti=(dx*snps1+dy*csps1)*snth1+dz*csth1
-          rti=sqrt(xti**2+yti**2)      
+          rti=sqrt(xti**2+yti**2)     
+c store x,y,z values for spline interpolation
+          iray=iray+1
+          xtiv(iray)=xti
+          ytiv(iray)=yti 
+          zwjv(iray)=zwj
+c initialize grid dimension for spline interpolation
+          xmaxgrid=max(xmaxgrid,rti)
 c
           dirxt= (dirx*csps1-diry*snps1)/dir
           diryt=((dirx*snps1+diry*csps1)*csth1-dirz*snth1)/dir
@@ -5880,6 +5909,58 @@ c
 c     in common dnpara=dnpar/2 to match  westerhof Delta function
 c     Delta -> exp[-(...)^2/(2*DeltaQ)]
 c
+
+
+c     Spline fit
+      if(nrayr.gt.1) then
+
+      npoints=iray
+      xmin=-xmaxgrid
+      xmax=xmaxgrid
+      ymin=-xmaxgrid
+      ymax=xmaxgrid
+      nxgrid=2*nrayr-1
+      dx=(xmax-xmin)/(nxgrid-1)
+      do i=1,nxgrid
+        xgridv(i)=xmin+(i-1)*dx
+        ygridv(i)=xgridv(i)
+      end do
+
+c     call interp spline
+      iopt=0
+      sspl=1.0d-3
+      eps=1.0d-7
+      do iray=1,npoints
+        w(iray) = 1.0d0
+      end do
+      call surfit(iopt,npoints,xtiv,ytiv,zwjv,w,xmin,xmax,ymin,ymax,
+     .            kspl,kspl,sspl,nxest,nxest,nxest,eps,
+     .            nkntx,tx,nknty,ty,ccexp,resid,wrk1,lwrk1,wrk2,lwrk2,
+     .            iwrk,kwrk,ierr)
+      if (ierr.gt.0) then
+        print*, 'surfit:',istep,nkntx,nknty,ierr,resid
+        do j=1,nxgrid
+          do i=1,nxgrid
+            zwint(nxgrid*(i-1)+j)=0.0d0
+          end do
+        end do
+      else
+c     call eval spline
+        call bispev(tx,nkntx,ty,nknty,ccexp,kspl,kspl,xgridv,nxgrid,
+     .     ygridv,nxgrid,zwint,wrkbsp,lwrkbsp,iwrkbsp,liwrkbsp,ierr)
+        if (ierr.ne.0) print*, 'bispev:',istep,ierr
+      end if
+      do j=1,nxgrid
+        do i=1,nxgrid
+          ij=nxgrid*(i-1)+j
+          write(82,111) istep,i,j,xgridv(i),ygridv(j),zwint(ij),
+     .      aaw*xgridv(i)**2+ccw*xgridv(i)*ygridv(j)+bbw*ygridv(j)**2
+        end do
+        write(82,*) ''
+      end do
+      write(82,*) ''
+      end if
+
       write(12,99) istep,st,
      . wcsi,weta,phiwdeg,rcicsi,rcieta,phirdeg,errw,errr,
      . dk1,dk2,dkpar,phik*180.0d0/pi,dnpar