module Numeric.Interpolation.Piece (
   Piece.T,
   linear,
   hermite1,
   ) where

import qualified Numeric.Interpolation.Private.Piece as Piece


{- $setup
>>> import qualified Numeric.Interpolation.Piece as Piece
>>> import qualified Numeric.Interpolation.Private.Piece as PiecePriv
>>> import qualified Test.QuickCheck as QC
>>> import Test.QuickCheck ((==>))
>>>
>>> forAllDistinctPoints ::
>>>    (Show a, QC.Arbitrary a, QC.Testable prop) =>
>>>    ((Rational, a) -> (Rational, a) -> prop) -> QC.Property
>>> forAllDistinctPoints f =
>>>    QC.forAll QC.arbitrary $ \p1@(x1,_) ->
>>>    QC.forAll QC.arbitrary $ \p2@(x2,_) ->
>>>       x1/=x2  ==>  f p1 p2
-}

{- |
prop> forAllDistinctPoints $ \p1 p2 x -> Piece.linear p1 p2 x == Piece.linear p2 p1 x
-}
linear :: (Fractional a) => Piece.T a a a
linear :: T a a a
linear = T a a a
forall a. Fractional a => T a a a
Piece.linear

{- |
Hermite interpolation with one derivative per node.
That is, the interpolating polynomial is cubic.

prop> forAllDistinctPoints $ \p1 p2 x -> Piece.hermite1 p1 p2 x == Piece.hermite1 p2 p1 x
prop> forAllDistinctPoints $ \p1@(x1,y1) p2@(x2,y2) x -> Piece.linear p1 p2 x == let slope = (y2-y1)/(x2-x1) in Piece.hermite1 (x1, (y1,slope)) (x2, (y2,slope)) x
prop> forAllDistinctPoints $ \p1 p2 x -> Piece.hermite1 p1 p2 x == PiecePriv.hermite1 p1 p2 x
-}
hermite1 :: (Fractional a) => Piece.T a a (a, a)
hermite1 :: T a a (a, a)
hermite1 (a
x0,(a
y0,a
dy0)) (a
x1,(a
y1,a
dy1)) a
x =
   let d :: a
d = (a
y1a -> a -> a
forall a. Num a => a -> a -> a
-a
y0) a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
dx10
       dx0 :: a
dx0 = a
xa -> a -> a
forall a. Num a => a -> a -> a
-a
x0
       dx1 :: a
dx1 = a
x1a -> a -> a
forall a. Num a => a -> a -> a
-a
x
       dx10 :: a
dx10 = a
x1a -> a -> a
forall a. Num a => a -> a -> a
-a
x0
   in  (a
y0a -> a -> a
forall a. Num a => a -> a -> a
*a
dx1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
y1a -> a -> a
forall a. Num a => a -> a -> a
*a
dx0 a -> a -> a
forall a. Num a => a -> a -> a
+
        ((a
dy0a -> a -> a
forall a. Num a => a -> a -> a
-a
d) a -> a -> a
forall a. Num a => a -> a -> a
* a
dx1 a -> a -> a
forall a. Num a => a -> a -> a
- (a
dy1a -> a -> a
forall a. Num a => a -> a -> a
-a
d) a -> a -> a
forall a. Num a => a -> a -> a
* a
dx0) a -> a -> a
forall a. Num a => a -> a -> a
* a
dx0 a -> a -> a
forall a. Num a => a -> a -> a
* a
dx1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
dx10)
          a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
dx10