stellinstall/evspline_8f_source.html

      subroutine evspline(xget,x,nx,ilinx,f,ict,fval,ier)

C

C  use mkspline to set up spline coefficients...

C

C  evaluate a 1d cubic Spline interpolant -- this is C2

C

C  this subroutine calls two subroutines:

C     herm1x  -- find cell containing (xget)

C     fvspline  -- evaluate interpolant function and (optionally) derivatives

C

C  input arguments:

C  ================

C

      integer nx                        ! grid size

      real xget                         ! target of this interpolation

      real x(nx)                        ! ordered x grid

      integer ilinx                     ! ilinx=1 => assume x evenly spaced

C

      real f(0:1,nx)                    ! function data

C

C  f(0,i) = f @ x(i)

C  f(1,i) = d2f/dx2 @ x(i)

C

C      (these are spline coefficients selected for continuous 2-

C      diffentiability, see mkspline.for)

C

      integer ict(3)                    ! code specifying output desired

C

C  ict(1)=1 -- return f  (0, don't)

C  ict(2)=1 -- return df/dx  (0, don't)

C  ict(3)=1 -- return d2f/dx2  (0, don't)

c

c        set ict(1)=3 to get d3f/dx3 (only)

C

C output arguments:

C =================

C

      real fval(*)                      ! output data

      integer ier                       ! error code =0 ==> no error

C

C  fval(1) receives the first output (depends on ict(...) spec)

C  fval(2) receives the second output (depends on ict(...) spec)

C  fval(3) receives the third output (depends on ict(...) spec)

C

C  examples:

C    on input ict = [1,1,1]

C   on output fval= [f,df/dx,d2f/dx2]

C

C    on input ict = [1,0,0]

C   on output fval= [f] ... elements 2 -- 3 never referenced.

C

C    on input ict = [0,0,1]

C   on output fval= [d2f/dx2] ... elements 2 -- 3 never referenced.

C

C  ier -- completion code:  0 means OK

C-------------------

C  local:

C

      integer i(1)                         ! cell indices

C

C  normalized displacement from x(i) grid point.

C    xparam=0 @x(i)  xparam=1 @x(i+1)

C

      real xparam(1)

C

C  cell dimensions and

C  inverse cell dimensions hxi = 1/(x(i+1)-x(i))

C

      real hx(1)

      real hxi(1)

C

C  0 .le. xparam .le. 1

C

C  ** the interface is very similar to herm2ev.for; can use herm2xy **

C---------------------------------------------------------------------

C

      i(1)=0

C

      call herm1x(xget,x,nx,ilinx,i(1),xparam(1),hx(1),hxi(1),ier)

      if(ier.ne.0) return

c

      call fvspline(ict,1,1,fval,i,xparam,hx,hxi,f,nx)

C

      return

      end

C---------------------------------------------------------------------

C  evaluate C1 cubic Hermite function interpolation -- 2d fcn

C   --vectorized-- dmc 10 Feb 1999

C

      subroutine fvspline(ict,ivec,ivecd,

     >   fval,ii,xparam,hx,hxi,fin,nx)

C

C  use mkspline to set up spline coefficients...

C

      integer ict(3)                    ! requested output control

      integer ivec                      ! vector length

      integer ivecd                     ! vector dimension (1st dim of fval)

C

      integer ii(ivec)                  ! target cells (i)

      real xparam(ivec)

                          ! normalized displacements from (i) corners

C

      real hx(ivec)                     ! grid spacing, and

      real hxi(ivec)                    ! inverse grid spacing 1/(x(i+1)-x(i))

C

      real fin(0:1,nx)                  ! interpolant data (cf "evspline")

C

      real fval(ivecd,*)                ! output returned

C

C  for detailed description of fin, ict and fval see subroutine

C  evspline comments.  Note ict is not vectorized; the same output

C  is expected to be returned for all input vector data points.

C

C  note that the index inputs ii and parameter inputs

C     xparam,hx,hxi, are vectorized, and the

C     output array fval has a vector ** 1st dimension ** whose

C     size must be given as a separate argument

C

C  to use this routine in scalar mode, pass in ivec=ivecd=1

C

C---------------

C  Spline evaluation consists of a "mixing" of the interpolant

C  data using the linear functionals xparam, xpi = 1-xparam,

C  xparam**3-xparam, xpi**3-xpi ...

C  and their derivatives as needed.

C

      integer v

C

      real sum

      real sixth

C

      data sixth/0.166666666666666667/

C

C---------------

C   ...in x direction

C

      iadr=0

C

      if(ict(1).le.2) then

C

C  normal call

C

         if(ict(1).eq.1) then

C

C  get desired values:

C

            iadr = iadr + 1

            do v=1,ivec

               i=ii(v)

C

               xp=xparam(v)

               xpi=1.0-xp

               xp2=xp*xp

               xpi2=xpi*xpi

C

               cx=xp*(xp2-1.0)

               cxi=xpi*(xpi2-1.0)

               hx2=hx(v)*hx(v)

C

C  function value:

C

               sum=xpi*fin(0,i) +xp*fin(0,i+1)

               sum=sum+sixth*hx2*(cxi*fin(1,i) +cx*fin(1,i+1))

C

               fval(v,iadr)=sum

            enddo

         endif

C

         if(ict(2).eq.1) then

C

C  df/dx:

C

            iadr=iadr+1

            do v=1,ivec

               i=ii(v)

C

               xp=xparam(v)

               xpi=1.0-xp

               xp2=xp*xp

               xpi2=xpi*xpi

C

               cxd=3.0*xp2-1.0

               cxdi=-3.0*xpi2+1.0


               sum=hxi(v)*(fin(0,i+1)-fin(0,i))

               sum=sum+sixth*hx(v)*(cxdi*fin(1,i) +cxd*fin(1,i+1))

C

               fval(v,iadr)=sum

            enddo

         endif

C

         if(ict(3).eq.1) then

C

C  d2f/dx2:

C

            iadr=iadr+1

            do v=1,ivec

               i=ii(v)

C

               xp=xparam(v)

               xpi=1.0-xp

C

               sum=xpi*fin(1,i) +xp*fin(1,i+1)

               fval(v,iadr)=sum

            enddo

         endif

C

      else

C

C  return fxxx = d3f/dx3

C

         iadr=1

         do v=1,ivec

            i=ii(v)

            fval(v,iadr)=hxi(v)*(fin(1,i+1)-fin(1,i))

         enddo

C

      endif

C

      return

      end