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