program cspline

! uses numerical recipes routines spline and splint.

implicit real*8 (a - h, o - z)
dimension x (1000), y (1000), deriv2 (1000)
dimension xtemp (1000), ytemp (1000)
character*40 input1, input2, output

write (*, '(5x, a)') 'input spectrum file.'
read (*, 10) input1
write (*, '(5x, a)') 'input spline position file.'
read (*, 10) input2
write (*,'(5x,a)') 'output file.'
read (*, 10) output
10    format (a)
open (unit = 1, file = input1, status = 'old')
open (unit = 2, file = input2, status = 'old')
open (unit = 3, file = output, status='unknown')

i = 1
20 read (1, *, end = 30) x (i), y (i)
  i = i + 1
  go to 20
30 npoints = i - 1

if (x (npoints) .lt. x (1)) then
  do i = 1, npoints
    xtemp (i) = x (i)
    ytemp (i) = y (i)
  end do
  do i = 1, npoints
    x (i) = xtemp (npoints + 1 - i)
    y (i) = ytemp (npoints + 1 - i)
  end do
end if

call spline (x, y, npoints, deriv2)

40 read (2, *, end = 50) xcalc
  call splint (x, y, deriv2, npoints, xcalc, ycalc)
  write (3, '(f11.5, 1pe13.5)') xcalc, ycalc
  go to 40
50 continue

close (unit = 1)
close (unit = 2)
close (unit = 3)
stop
end
!
subroutine spline (x, y, n, y2)

! modified to always use "natural" boundary conditions

implicit real*8 (a - h, o - z)
parameter (nmax = 1000)
dimension x (n), y (n), y2 (n), u (nmax)
y2 (1) = 0.
u (1) = 0.
do i = 2, n - 1
  sig = (x (i) - x (i - 1)) / (x (i + 1) - x (i - 1))
  p = sig * y2 (i - 1) + 2.
  y2 (i) = (sig - 1.) / p
  u (i) = (6.* ((y (i + 1) - y (i)) / (x (i + 1) - x (i)) - (y (i) - &
  y (i - 1)) / (x (i) - x (i - 1))) / (x (i + 1) - x (i - 1)) - sig * &
  u (i - 1)) / p
end do
qn = 0.
un = 0.
y2 (n) = (un - qn * u (n - 1)) / (qn * y2 (n - 1) + 1.)
do k = n - 1, 1, -1
  y2 (k) = y2 (k) * y2 (k + 1) + u (k)
end do
return
end
!
subroutine splint (xa, ya, y2a, n, x, y)

implicit real*8 (a - h, o - z)
dimension xa (n), ya (n), y2a (n)
klo = 1
khi = n
1 if (khi - klo .gt. 1) then
  k = (khi + klo) / 2
  if (xa (k) .gt. x) then
    khi = k
  else
    klo = k
  end if
  go to 1
end if
h = xa (khi) - xa (klo)
if (h .eq. 0.) pause 'bad xa input.'
a = (xa (khi) - x) / h
b = (x - xa (klo)) / h
y = a * ya (klo) + b * ya (khi) + ((a**3 - a) * y2a (klo) + (b**3 - b) * &
y2a (khi)) * (h**2) / 6.

return
end