* add missing routine (sorry!)

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEOS/EOS_Omni/EOS_Omni@10 8e189c6b-2ab8-4400-aa02-70a9cfce18b9

2 files changed, 134 insertions, 1 deletions

diff --git a/src/nuc_eos/linterp.f b/src/nuc_eos/linterp.f
new file mode 100644
index 0000000..282d090
--- /dev/null
+++ b/src/nuc_eos/linterp.f
@@ -0,0 +1,132 @@
+      SUBROUTINE intp3d ( x, y, z, f, kt, ft, nx, ny, nz, xt, yt, zt,
+     .                   d1, d2, d3 )     
+c
+      implicit none
+c                                                          
+c---------------------------------------------------------------------
+c
+c     purpose: interpolation of a function of three variables in an
+c              equidistant(!!!) table.
+c
+c     method:  8-point Lagrange linear interpolation formula          
+c
+c     x        input vector of first  variable
+c     y        input vector of second variable
+c     z        input vector of third  variable
+c
+c     f        output vector of interpolated function values
+c
+c     kt       vector length of input and output vectors
+c
+c     ft       3d array of tabulated function values
+c     nx       x-dimension of table
+c     ny       y-dimension of table
+c     nz       z-dimension of table
+c     xt       vector of x-coordinates of table
+c     yt       vector of y-coordinates of table
+c     zt       vector of z-coordinates of table
+c
+c     d1       centered derivative of ft with respect to x
+c     d2       centered derivative of ft with respect to y
+c     d3       centered derivative of ft with respect to z
+c     Note that d? only make sense when intp3d is called with kt=1
+c---------------------------------------------------------------------
+c
+c 
+
+c
+      integer kt,nx,ny,nz,ktx
+      double precision x(kt),y(kt),z(kt),f(kt)
+      double precision xt(nx),yt(ny),zt(nz)
+      double precision ft(nx,ny,nz)
+      double precision d1,d2,d3
+c
+c
+      PARAMETER   (ktx = 400)
+      double precision  fh(ktx,8), delx(ktx), dely(ktx), delz(ktx),
+     &            a1(ktx), a2(ktx), a3(ktx), a4(ktx),
+     &            a5(ktx), a6(ktx), a7(ktx), a8(ktx)
+
+      double precision dx,dy,dz,dxi,dyi,dzi,dxyi,dxzi,dyzi,dxyzi
+      integer n,ix,iy,iz
+
+      IF (kt .GT. ktx)  STOP'***KTX**'
+c
+c
+c------  determine spacing parameters of (equidistant!!!) table
+c
+      dx    = (xt(nx) - xt(1)) / FLOAT(nx-1)
+      dy    = (yt(ny) - yt(1)) / FLOAT(ny-1)
+      dz    = (zt(nz) - zt(1)) / FLOAT(nz-1)
+c
+      dxi   = 1. / dx
+      dyi   = 1. / dy
+      dzi   = 1. / dz
+c
+      dxyi  = dxi * dyi
+      dxzi  = dxi * dzi
+      dyzi  = dyi * dzi
+c
+      dxyzi = dxi * dyi * dzi
+c
+c
+c------- loop over all points to be interpolated
+c
+      DO  n = 1, kt                                            
+c
+c------- determine location in (equidistant!!!) table 
+c                                                                  
+         ix = 2 + INT( (x(n) - xt(1) - 1.e-10) * dxi )
+         iy = 2 + INT( (y(n) - yt(1) - 1.e-10) * dyi )
+         iz = 2 + INT( (z(n) - zt(1) - 1.e-10) * dzi )
+c                                                     
+         ix = MAX( 2, MIN( ix, nx ) )
+         iy = MAX( 2, MIN( iy, ny ) )
+         iz = MAX( 2, MIN( iz, nz ) )
+c
+c         write(*,*) iy-1,iy,iy+1
+c
+c------- set-up auxiliary arrays for Lagrange interpolation
+c                                                                 
+         delx(n) = xt(ix) - x(n)
+         dely(n) = yt(iy) - y(n)
+         delz(n) = zt(iz) - z(n)
+c                                                                   
+         fh(n,1) = ft(ix  , iy  , iz  )                             
+         fh(n,2) = ft(ix-1, iy  , iz  )                             
+         fh(n,3) = ft(ix  , iy-1, iz  )                             
+         fh(n,4) = ft(ix  , iy  , iz-1)                             
+         fh(n,5) = ft(ix-1, iy-1, iz  )                             
+         fh(n,6) = ft(ix-1, iy  , iz-1)                             
+         fh(n,7) = ft(ix  , iy-1, iz-1)                             
+         fh(n,8) = ft(ix-1, iy-1, iz-1)                             
+c              
+c------ set up coefficients of the interpolation polynomial and 
+c       evaluate function values 
+c                                                    
+         a1(n) = fh(n,1)                             
+         a2(n) = dxi   * ( fh(n,2) - fh(n,1) )       
+         a3(n) = dyi   * ( fh(n,3) - fh(n,1) )       
+         a4(n) = dzi   * ( fh(n,4) - fh(n,1) )       
+         a5(n) = dxyi  * ( fh(n,5) - fh(n,2) - fh(n,3) + fh(n,1) )
+         a6(n) = dxzi  * ( fh(n,6) - fh(n,2) - fh(n,4) + fh(n,1) )
+         a7(n) = dyzi  * ( fh(n,7) - fh(n,3) - fh(n,4) + fh(n,1) )
+         a8(n) = dxyzi * ( fh(n,8) - fh(n,1) + fh(n,2) + fh(n,3) +
+     &                     fh(n,4) - fh(n,5) - fh(n,6) - fh(n,7) )
+c
+         d1 = -a2(n)
+         d2 = -a3(n)
+         d3 = -a4(n)
+         f(n)  = a1(n) +  a2(n) * delx(n)                         
+     &                 +  a3(n) * dely(n)                         
+     &                 +  a4(n) * delz(n)                         
+     &                 +  a5(n) * delx(n) * dely(n)               
+     &                 +  a6(n) * delx(n) * delz(n)               
+     &                 +  a7(n) * dely(n) * delz(n)               
+     &                 +  a8(n) * delx(n) * dely(n) * delz(n)     
+c
+      ENDDO
+c                                                                    
+      RETURN                                                         
+      END                                                            
+
diff --git a/src/nuc_eos/make.code.defn b/src/nuc_eos/make.code.defn
index 9be5ee3..59c8118 100644
--- a/src/nuc_eos/make.code.defn
+++ b/src/nuc_eos/make.code.defn
@@ -1,5 +1,6 @@
 SRCS = eosmodule.F90 nuc_eos.F90 bisection.F90 \
-       findtemp.F90 linterp_many.F90 readtable.F90
+       findtemp.F90 linterp_many.F90 readtable.F90 \
+       linterp.f
 
 SUBDIRS =