c     Copyright (C) 2008-2014 Vincent Eymet
c
c     KDISTRIBUTION is free software; you can redistribute it and/or modify
c     it under the terms of the GNU General Public License as published by
c     the Free Software Foundation; either version 3, or (at your option)
c     any later version.

c     KDISTRIBUTION is distributed in the hope that it will be useful,
c     but WITHOUT ANY WARRANTY; without even the implied warranty of
c     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
c     GNU General Public License for more details.

c     You should have received a copy of the GNU General Public License
c     along with KDISTRIBUTION; if not, see <http://www.gnu.org/licenses/>
c
      double precision function S(z1,z2,x1,x2,y1,y2,x)
      implicit none
      include 'max.inc'
c     
c     Purpose: to compute the spline function between data points x1 and x2
c     
c     Inputs:
c       + z1: first coefficient for the interval
c       + z2: second coefficient for the interval
c       + x1: value of x at the beginning of the interval
c       + x2: value of x at the end of the interval
c       + y1: value of y at the beginning of the interval
c       + y2: value of y at the end of the interval
c       + x: value of x for which S has to be computed
c
      double precision z1,z2,x1,x2,y1,y2,x
c     label
      integer strlen
      character*(Nchar_mx) label
      label='subroutine S'

c     Debug
      if (x2.lt.x1) then
         write(*,*) 'Error from function ',label(1:strlen(label)),':'
         write(*,*) 'x1=',x1,' x2=',x2
         stop
      endif
      if (x.lt.x1) then
         write(*,*) 'Error from function ',label(1:strlen(label)),':'
         write(*,*) 'x1=',x1,' x=',x
         stop
      endif
      if (x2.lt.x) then
         write(*,*) 'Error from function ',label(1:strlen(label)),':'
         write(*,*) 'x=',x,' x2=',x2
         stop
      endif
c     Debug
      
      S=(z1*(x2-x)**3.0D+0+z2*(x-x1)**3.0D+0)/(6.0D+0*(x2-x1))
     &     +(y2/(x2-x1)-(x2-x1)/6.0D+0*z2)*(x-x1)
     &     +(y1/(x2-x1)-(x2-x1)/6.0D+0*z1)*(x2-x)
      
      return
      end