Berechnung eines interpolierenden kubischen Splines

module mysubs

contains

  subroutine init(n, t, y, x, u, v, w, nd, hd, c)
    use kinds
    use fileutil
    implicit none
    real(kind=REAL8), dimension(:), pointer :: t, y
    real(kind=REAL8), dimension(:), pointer :: u, v, w, nd, hd, c
    integer, intent(out)          :: n
    real(kind=REAL8), intent(out) :: x

    integer, parameter  :: fileid = 1
    character(len=256)  :: fileName
    integer             :: status, i

    write (*,*) 'Berechnet ein interpoliertes kubisches natuerliches Spline'
    write (*,*)
    write (*,*) 'Welche Datei enthaelt die Punktepaare?&
                & (In splines.dat liegt ein Beispiel.)'
    read (*,*) fileName

    open(unit=fileId, file=fileName, action='read', iostat=status)
    if (status /= 0) then
       write (*,*) 'Kann Datei ', fileName, 'nicht öffnen.'
       stop
    end if

    n = anzahlReals2(fileId) - 1
    allocate(t(0:n), y(0:n))
    allocate(u(0:n), v(0:n), w(0:n), nd(0:n), hd(0:n), c(0:n))

    read(fileid, *, iostat=status) (t(i), y(i), i=0, n)
    if (status /= 0) then
       write (*,*) 'Konnte Array nicht lesen.'
       stop
    end if
    close(unit=fileid)

    write (*,*) 'Insgesamt ', n + 1, ' Stuetzstellen gelesen.'
    write (*,*)
    write (*, '(A)', advance='no') &
         & 'An welcher Stelle soll das Spline berechnet werden? '
    read (*,*) x
    write (*,*)

    return
  end subroutine init


  subroutine Randbedingungen(c, u, v, w, y, t)
    use kinds
    implicit none
    real(kind=REAL8), dimension(:), pointer :: c, u, v, w, y, t

    real(kind=REAL8)  :: h1, hn, d1, dn, b1, b2
    real(kind=REAL8)  :: alpha, beta, det
    real(kind=REAL8)  :: a11, a12, a21, a22
    integer           :: i, n

    n = size(t) - 1
    h1 = t(1) - t(0)
    hn = t(n) - t(n - 1)
    d1 = (y(1) - y(0)) / h1
    dn = (y(n) - y(n - 1)) / hn
    b1 = 3 * d1 - u(1)
    b2 = 3 * dn - u(n - 1)
    a11 = 2 * v(1)
    a12 = w(1)
    a21 = v(n - 1)
    a22 = 2 + w(n - 1)
    det = a11 * a22 - a12 * a21
    alpha = (b1 * a22 - a12 * b2) / det
    beta = (a11 * b2 - b1 * a21) / det
    c(0) = alpha

    do i = 1, n - 1
       c(i) = u(i) + alpha * v(i) + beta * w(i)
    end do
    c(n) = beta

    return
  end subroutine Randbedingungen
  

  function Spline(x, t, y, c)
    use kinds
    implicit none
    real(kind=REAL8), dimension(:), pointer :: t, y, c
    real(kind=REAL8), intent(in) :: x
    real(kind=REAL8)             :: Spline

    real(kind=REAL8)  :: h, delta, a0, a1, a2, a3
    integer           :: k, n
 
    n = size(t) - 1
    k = n + 1
    do
       k = k - 1
       if (x > t(k)) exit
    end do
    k = k + 1
    h = t(k) - t(k - 1)
    delta = (y(k) - y(k - 1)) / h
    a0 = y(k)
    a1 = c(k)
    a2 = (c(k) - delta) / h
    a3 = (c(k) - 2 * delta + c(k - 1)) / (h*h)

    Spline = ((a3 * (x - t(k - 1)) + a2) * (x - t(k)) + a1) * (x- t(k)) + a0
    return
  end function Spline


  subroutine result(x, t, y, c)
    use kinds
    implicit none
    real(kind=REAL8), dimension(:), pointer :: t, y, c
    real(kind=REAL8), intent(in) :: x

    write (*,'(A, F12.5)') 'Das Ergenis hat den Wert ', Spline(x, t, y, c)
    return
  end subroutine result

end module mysubs


program Splines
  ! Berechnung eines interpolierten kubischen natuerlichen Splines}

  use kinds
  use mysubs
  implicit none

  real(kind=REAL8), dimension(:), pointer :: t, y
  real(kind=REAL8), dimension(:), pointer :: u, v, w, nd, hd, c
  real(kind=REAL8)  :: x, h0, h1, a, b, q
  integer           :: n, i

  call init(n, t, y, x, u, v, w, nd, hd, c)

  h0 = t(1) - t(0)
  a = (y(1) - y(0)) / (h0*h0)
  v(1) = -1 / h0
  do i = 1, n - 1
     h1 = t(i + 1) - t(i)
     b = (y(i + 1) - y(i)) / (h1*h1)
     u(i) = 3 * (a + b)
     hd(i) = 2 * (1 / h0 + 1 / h1)
     nd(i) = 1 / h1
     a = b
     h0 = h1
  end do
  w(n - 1) = -1 / h1

  ! Vorwaertselimination
  do i = 1, n - 2
     q = -nd(i) / hd(i)
     hd(i + 1) = hd(i + 1) + nd(i) * q
     u(i + 1) = u(i + 1) + q * u(i)
     v(i + 1) = q * v(i)
  end do

  ! Rueckwaertselimination
  do i = n - 1, 2, -1
     q = hd(i)
     u(i) = u(i) / q
     v(i) = v(i) / q
     w(i) = w(i) / q
     q = nd(i - 1)
     u(i - 1) = u(i - 1) - q * u(i)
     v(i - 1) = v(i - 1) - q * v(i)
     w(i - 1) = -q * w(i)
  end do

  q = hd(1)
  u(1) = u(1) / q
  v(1) = v(1) / q
  w(1) = w(1) / q
  call Randbedingungen(c, u, v, w, y, t)
  call result(x, t, y, c)

  stop
end program Splines