-----------------------------------------------------------------------------
-- |
-- 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 :: forall a. Fractional a => [a] -> [a]
solveCubicSplineDerivatives (a
x:[a]
xs) = forall a. Fractional a => [a] -> [a] -> [a] -> [a] -> [a]
solveTriDiagonal [a]
as [a]
bs [a]
as [a]
ds
  where
    as :: [a]
as = forall a. Int -> a -> [a]
replicate (Int
l forall a. Num a => a -> a -> a
- Int
1) a
1
    bs :: [a]
bs = a
2 forall a. a -> [a] -> [a]
: forall a. Int -> a -> [a]
replicate (Int
l forall a. Num a => a -> a -> a
- Int
2) a
4 forall a. [a] -> [a] -> [a]
++ [a
2]
    l :: Int
l  = forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
ds
    ds :: [a]
ds = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
f ([a]
xs forall a. [a] -> [a] -> [a]
++ [forall a. [a] -> a
last [a]
xs]) (a
xforall a. a -> [a] -> [a]
:a
xforall a. a -> [a] -> [a]
:[a]
xs)
    f :: a -> a -> a
f a
a a
b = a
3forall a. Num a => a -> a -> a
*(a
a forall a. Num a => a -> a -> a
- a
b)

solveCubicSplineDerivatives [a]
_ = forall a. HasCallStack => [Char] -> a
error [Char]
"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 :: forall a. Fractional a => [a] -> [a]
solveCubicSplineDerivativesClosed [a]
xs = forall a. Fractional a => [a] -> [a] -> [a] -> [a] -> a -> a -> [a]
solveCyclicTriDiagonal [a]
as [a]
bs [a]
as [a]
ds a
1 a
1
  where
    as :: [a]
as = forall a. Int -> a -> [a]
replicate (Int
l forall a. Num a => a -> a -> a
- Int
1) a
1
    bs :: [a]
bs = forall a. Int -> a -> [a]
replicate Int
l a
4
    l :: Int
l  = forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
xs
    xs' :: [a]
xs' = forall a. [a] -> [a]
cycle [a]
xs
    ds :: [a]
ds = forall a. Int -> [a] -> [a]
take Int
l forall a b. (a -> b) -> a -> b
$ forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
f (forall a. Int -> [a] -> [a]
drop Int
1 [a]
xs') (forall a. Int -> [a] -> [a]
drop (Int
l forall a. Num a => a -> a -> a
- Int
1) [a]
xs')
    f :: a -> a -> a
f a
a a
b = a
3forall a. Num a => a -> a -> a
*(a
a forall a. Num a => a -> a -> a
- a
b)

-- | Use the cyclic-tri-diagonal solver with the appropriate parameters for a closed cubic spline.
solveCubicSplineCoefficients :: Fractional a => Bool -> [a] -> [[a]]
solveCubicSplineCoefficients :: forall a. Fractional a => Bool -> [a] -> [[a]]
solveCubicSplineCoefficients Bool
closed [a]
xs =
    [ [a
x,a
d,a
3forall a. Num a => a -> a -> a
*(a
x1forall a. Num a => a -> a -> a
-a
x)forall a. Num a => a -> a -> a
-a
2forall a. Num a => a -> a -> a
*a
dforall a. Num a => a -> a -> a
-a
d1,a
2forall a. Num a => a -> a -> a
*(a
xforall a. Num a => a -> a -> a
-a
x1)forall a. Num a => a -> a -> a
+a
dforall a. Num a => a -> a -> a
+a
d1]
    | (a
x,a
x1,a
d,a
d1) <- forall a b c d. [a] -> [b] -> [c] -> [d] -> [(a, b, c, d)]
zip4 [a]
xs' (forall a. [a] -> [a]
tail [a]
xs') [a]
ds' (forall a. [a] -> [a]
tail [a]
ds')
    ]
  where
    ds :: [a]
ds | Bool
closed    = forall a. Fractional a => [a] -> [a]
solveCubicSplineDerivativesClosed [a]
xs
       | Bool
otherwise = forall a. Fractional a => [a] -> [a]
solveCubicSplineDerivatives [a]
xs
    close :: [a] -> [a]
close [a]
as | Bool
closed    = [a]
as forall a. [a] -> [a] -> [a]
++ [forall a. [a] -> a
head [a]
as]
             | Bool
otherwise = [a]
as
    xs' :: [a]
xs' = forall a. [a] -> [a]
close [a]
xs
    ds' :: [a]
ds' = forall a. [a] -> [a]
close [a]
ds