      function sigmatn(x1,tau1)
      implicit NONE
      real*4 x1,tau1
      real*4 sigmatn
c--------------------------------------------------------------      
c     Purpose:
c              Find sigmaty at (x1,tau1)
c
c--------------------------------------------------------------      
c     Input :  
c              x1  x-coordinate
c              tau1  tau-coordinate
c-------------------------------------------------------------      

C      include 'Pareg.inc'
      include '../Src/raypar.inc'
      include '../Src/model.inc'

      real*4 ys4(4),yp1s4(4),yp2s4(4),ypps4(4),c(4,4),u,t,sig
      real*4 x1c,tau1c,eps/1.0e-6/

      integer ncol,nlin,num1,num2,num3,num4,i
      
c     correct coordinates relative to sigma model.

      x1c = x1 - xminco
      tau1c = tau1 - tauminco

c     check that point is within the model, if not: make them:

      x1c = amax1(amin1(x1c,xlength*(1.0-eps)),0.0)
      tau1c = amax1(amin1(tau1c,taulength*(1.0-eps)),0.0)

c     find nodes around the point.

      ncol = int(x1c*rdelx) + 1
      nlin = int(tau1c*rdeltau) + 1

c     check that point is within the model, if not: make it:

c      ncol = max0(min0(ncol,numx-1),1)
c      nlin = max0(min0(nlin,numtau-1),1)

      num1 = (ncol-1)*numtau + nlin
      num4 = num1 + 1
      num2 = ncol*numtau + nlin
      num3 = num2 + 1

      ys4(1) = sigmap(num1) 
      ys4(4) = sigmap(num4) 
      ys4(2) = sigmap(num2) 
      ys4(3) = sigmap(num3) 

      yp1s4(1) = dsigdx(num1) 
      yp1s4(4) = dsigdx(num4) 
      yp1s4(2) = dsigdx(num2) 
      yp1s4(3) = dsigdx(num3) 

      yp2s4(1) = dsigdtau(num1) 
      yp2s4(4) = dsigdtau(num4) 
      yp2s4(2) = dsigdtau(num2) 
      yp2s4(3) = dsigdtau(num3) 

      ypps4(1) = d2sigdxtau(num1) 
      ypps4(4) = d2sigdxtau(num4) 
      ypps4(2) = d2sigdxtau(num2) 
      ypps4(3) = d2sigdxtau(num3) 

c     compute coefficients 

      call bcucof(ys4,yp1s4,yp2s4,ypps4,delx,deltau,c)

      t =(x1c - real(ncol-1)*delx)*rdelx
      u =(tau1c - real(nlin-1)*deltau)*rdeltau 

      sig = 0.d0
      do i = 4,1,-1
        sig = t*sig+((c(i,4)*u+c(i,3))*u+c(i,2))*u+c(i,1)
      enddo

      sigmatn = sig
      return
      end

