-----------------------------------------------------------------------------
-- |
-- Module : Diagrams.CubicSpline
-- Copyright : (c) 2011 diagrams-lib team (see LICENSE)
-- License : BSD-style (see LICENSE)
-- Maintainer : diagrams-discuss@googlegroups.com
--
-- A /cubic spline/ is a smooth, connected sequence of cubic curves
-- passing through a given sequence of points. This module implements
-- a straightforward spline generation algorithm based on solving
-- tridiagonal systems of linear equations.
--
-----------------------------------------------------------------------------
module Diagrams.CubicSpline.Internal
(
-- * Solving for spline coefficents
solveCubicSplineDerivatives
, solveCubicSplineDerivativesClosed
, solveCubicSplineCoefficients
) where
import Diagrams.Solve.Tridiagonal
import Data.List
-- | Use the tri-diagonal solver with the appropriate parameters for an open cubic spline.
solveCubicSplineDerivatives :: Fractional a => [a] -> [a]
solveCubicSplineDerivatives (x:xs) = solveTriDiagonal as bs as ds
where
as = replicate (l - 1) 1
bs = 2 : replicate (l - 2) 4 ++ [2]
l = length ds
ds = zipWith f (xs ++ [last xs]) (x:x:xs)
f a b = 3*(a - b)
solveCubicSplineDerivatives _ = error "argument to solveCubicSplineDerivatives must be nonempty"
-- | Use the cyclic-tri-diagonal solver with the appropriate parameters for a closed cubic spline.
solveCubicSplineDerivativesClosed :: Fractional a => [a] -> [a]
solveCubicSplineDerivativesClosed xs = solveCyclicTriDiagonal as bs as ds 1 1
where
as = replicate (l - 1) 1
bs = replicate l 4
l = length xs
xs' = cycle xs
ds = take l $ zipWith f (drop 1 xs') (drop (l - 1) xs')
f a b = 3*(a - b)
-- | Use the cyclic-tri-diagonal solver with the appropriate parameters for a closed cubic spline.
solveCubicSplineCoefficients :: Fractional a => Bool -> [a] -> [[a]]
solveCubicSplineCoefficients closed xs =
[ [x,d,3*(x1-x)-2*d-d1,2*(x-x1)+d+d1]
| (x,x1,d,d1) <- zip4 xs' (tail xs') ds' (tail ds')
]
where
ds | closed = solveCubicSplineDerivativesClosed xs
| otherwise = solveCubicSplineDerivatives xs
close as | closed = as ++ [head as]
| otherwise = as
xs' = close xs
ds' = close ds