1 c
   2 c     http://www.netlib.org/fmm/spline.f
   3 c
   4       subroutine spline(n, x, y, b, c, d)
   5       integer n
   6       double precision x(n), y(n), b(n), c(n), d(n)
   7 c
   8       integer nmi, ib, i
   9       double precision t
  10 c
  11       nmi = n-1
  12       if ( n .lt. 2 ) return
  13       if ( n .lt. 3 ) goto 50
  14 c
  15       d(1= x(2- x(1)
  16       c(2= ( y(2- y(1) ) / d(1)
  17       do 10 i = 2,nmi
  18         d(i) = x(i+1- x(i)
  19         b(i) = 2.0*(d(i-1)+d(i))
  20         c(i+1= (y(i+1)-y(i))/d(i)
  21         c(i) = c(i+1)-c(i)
  22    10 continue
  23 c
  24       b(1= -d(1)
  25       b(n) = -d(n-1)
  26       c(1= 0.0
  27       c(n) = 0.0
  28       if ( n .eq. 3 ) goto 15
  29       c(1= c(3/ (x(4)-x(2)) - c(2)/(x(3)-x(1))
  30       c(n) = c(n-1)/(x(n)-x(n-2)) - c(n-2)/(x(n-1)-x(n-3))
  31       c(1= c(1)*d(1)**2/(x(4)-x(1))
  32       c(n) = -c(n)*d(n-1)**2/(x(n)-x(n-3))
  33 c
  34    15 do 20 i = 2, n
  35       t = d(i-1)/b(i-1)
  36       b(i) = b(i) - t*d(i-1)
  37       c(i) = c(i) - t*c(i-1)
  38    20 continue
  39 c
  40       c(n) = c(n)/b(n)
  41       do 30 ib = 1, nmi
  42         i = n - ib
  43         c(i) = (c(i) - d(i)*c(i+1))/b(i)
  44    30 continue
  45 c
  46       b(n) = (y(n) - y(nmi))/d(nmi) + d(nmi)*(c(nmi) + 2.0*c(n))
  47       do 40 i = 1, nmi
  48         b(i) = (y(i+1- y(i))/d(i) - d(i)*(c(i+1+ 2.0*c(i))
  49         d(i) = (c(i+1)-c(i))/d(i)
  50         c(i) = 3.0*c(i)
  51    40 continue
  52       c(n) = 3.0*c(n)
  53       d(n) = d(n-1)
  54       return
  55 c
  56    50 b(i) = (y(2)-y(1))/(x(2)-x(1))
  57       c(1= 0.0
  58       d(1= 0.0
  59       b(2= b(1)
  60       c(2= 0.0
  61       d(2= 0.0
  62       return
  63       end
  64 c
  65       double precision function seval(n, u, x, y, b, c, d)
  66       integer n
  67       double precision u, x(n), y(n), b(n), c(n), d(n)
  68 c
  69       integer i, j, k
  70       double precision dx
  71       data i/1/
  72       if ( i .ge. n ) i = 1
  73       if ( u .lt. x(i) ) goto 10
  74       if ( u .le. x(i+1) ) goto 30
  75    10 i = 1
  76       j = n + 1
  77    20 k = (i+j)/2
  78       if ( u .lt. x(k) ) j = k
  79       if ( u .ge. x(k) ) i = k
  80       if ( j .gt. i+1 ) goto 20
  81 c
  82    30 dx = u - x(i)
  83       seval = y(i) + dx*(b(i) + dx*(c(i) + dx*d(i)))
  84       return
  85       end
  86 
  87 c same but returns interpolation segment index
  88 c side effect: f=f(u) and fp=f'(u)
  89       integer function seval1(n, u, x, y, b, c, d, f, fp)
  90       integer n
  91       double precision u, x(n), y(n), b(n), c(n), d(n)
  92       double precision f, fp
  93 c
  94       integer i, j, k
  95       double precision dx
  96       data i/1/
  97       if ( i .ge. n ) i = 1
  98       if ( u .lt. x(i) ) goto 10
  99       if ( u .le. x(i+1) ) goto 30
 100    10 i = 1
 101       j = n + 1
 102    20 k = (i+j)/2
 103       if ( u .lt. x(k) ) j = k
 104       if ( u .ge. x(k) ) i = k
 105       if ( j .gt. i+1 ) goto 20
 106 c
 107    30 dx = u - x(i)
 108       f = y(i) + dx*(b(i) + dx*(c(i) + dx*d(i)))
 109       fp = b(i)+dx*(c(i)*2+dx*(d(i)*3))
 110       seval1 = i
 111       return
 112       end