This is an alternative site for discovering Elm packages. You may be looking for the official Elm package site instead.

# Interpolate.Cubic

This library interpolates cubic splines for one-dimensional sets of data.

It computes "natural splines", which means the second derivative at the endpoints is set to zero.

For more information on the mathematics used, check out this paper.

# Creating splines

withRange : Float -> Float -> List Float -> Spline

Compute a spline, given the minimum and maximum values of `x` and a list of data for `f(x)`. The data should be evenly spaced and in order of increasing `x`.

For example, if we have the data

``````f(2) = 1
f(3) = 5.2
f(4) = 3.2
f(5) = 0.8
``````

Then we would generate a spline by calling

``````fSpline = withRange 2 6 [ 1, 5.2, 3.2, 0.8 ]
``````

If there is only one data point, then the spline will be a horizontal line passing through that point. If the data is empty, the spline will be zero everywhere.

withDelta : Float -> Float -> List Float -> Spline

Same as `withRange`, except instead of passing the endpoint as the second argument, you pass the x-distance between data points.

``````fSpline = withDelta 2 1 [ 1, 5.2, 3.2, 0.8 ]
-- equivalent to withRange 2 6 [ 1, 5.2, 3.2, 0.8 ]
``````
type Spline = Spline { coefficients : Array Coefficients , start : Float , dx : Float }

# Interpolating

valueAt : Float -> Spline -> Float

Given `x1` and a spline `f(x)`, return `f(x1)`

slopeAt : Float -> Spline -> Float

Return the first derivative of the curve at the given point

concavityAt : Float -> Spline -> Float

Return the second derivative of the curve at the given point

curveAt : Float -> Spline -> LocalCurve

Return a record containing the value, slope, and concavity of the curve at the given point

type alias LocalCurve = { value : Float , slope : Float , concavity : Float }
``````module Interpolate.Cubic
exposing
( Spline
, LocalCurve
, withRange
, withDelta
, valueAt
, slopeAt
, concavityAt
, curveAt
)

{-| This library interpolates cubic splines for one-dimensional sets of data.

It computes "natural splines", which means the second derivative at the endpoints
is set to zero.

For more information on the mathematics used, check out
[this paper.](http://web.archive.org/web/20090408054627/http://online.redwoods.cc.ca.us/instruct/darnold/laproj/Fall98/SkyMeg/Proj.PDF)

# Creating splines
@docs withRange, withDelta, Spline

# Interpolating
@docs valueAt, slopeAt, concavityAt, curveAt, LocalCurve
-}

import Array exposing (Array)

{-| Given `x1` and a spline `f(x)`, return `f(x1)`
-}
valueAt : Float -> Spline -> Float
valueAt =
evaluate cubic

{-| Return the first derivative of the curve at the given point
-}
slopeAt : Float -> Spline -> Float
slopeAt =
evaluate firstDerivative

{-| Return the second derivative of the curve at the given point
-}
concavityAt : Float -> Spline -> Float
concavityAt =
evaluate secondDerivative

{-| Return a record containing the value, slope, and concavity of the curve
at the given point
-}
curveAt : Float -> Spline -> LocalCurve
curveAt =
let
curve x coeff =
{ value = cubic x coeff
, slope = firstDerivative x coeff
, concavity = secondDerivative x coeff
}
in
evaluate curve

cubic : Float -> Coefficients -> Float
cubic x { a, b, c, d } =
a * x ^ 3 + b * x ^ 2 + c * x + d

firstDerivative : Float -> Coefficients -> Float
firstDerivative x { a, b, c, d } =
3 * a * x ^ 2 + 2 * b * x + c

secondDerivative : Float -> Coefficients -> Float
secondDerivative x { a, b, c, d } =
6 * a * x + 2 * b

evaluate : (Float -> Coefficients -> a) -> Float -> Spline -> a
evaluate f x (Spline spline) =
let
maxIndex =
Array.length spline.coefficients - 1

index =
(x - spline.start)
/ spline.dx
|> floor
|> clamp 0 maxIndex

offset =
x - spline.dx * (toFloat index) - spline.start

coeff =
Array.get index spline.coefficients
|> Maybe.withDefault { a = 0, b = 0, c = 0, d = 0 }
in
f offset coeff

{-| -}
type alias LocalCurve =
{ value : Float
, slope : Float
, concavity : Float
}

{-| Compute a spline, given the minimum and maximum values of `x` and a
list of data for `f(x)`. The data should be evenly spaced and in order of
increasing `x`.

For example, if we have the data

f(2) = 1
f(3) = 5.2
f(4) = 3.2
f(5) = 0.8

Then we would generate a spline by calling

fSpline = withRange 2 6 [ 1, 5.2, 3.2, 0.8 ]

If there is only one data point, then the spline will be a horizontal line
passing through that point. If the data is empty, the spline will be zero
everywhere.
-}
withRange : Float -> Float -> List Float -> Spline
withRange start end heights =
let
n =
List.length heights |> toFloat

dx =
(end - start) / (n - 1)
in
if 1 < n then
withDelta start dx heights
else
{ coefficients = degenerateCoefficients heights
, start = 0
, dx = 0
}
|> Spline

{-| Same as `withRange`, except instead of passing the endpoint as the
second argument, you pass the x-distance between data points.

fSpline = withDelta 2 1 [ 1, 5.2, 3.2, 0.8 ]
-- equivalent to withRange 2 6 [ 1, 5.2, 3.2, 0.8 ]
-}
withDelta : Float -> Float -> List Float -> Spline
withDelta start dx heights =
{ coefficients = findCoefficients dx heights
, start = start
, dx = dx
}
|> Spline

degenerateCoefficients : List Float -> Array Coefficients
degenerateCoefficients heights =
{ a = 0
, b = 0
, c = 0
, d = List.head heights |> Maybe.withDefault 0
}
|> Array.repeat 1

findCoefficients : Float -> List Float -> Array Coefficients
findCoefficients dx heights =
let
concavity y0 y1 y2 =
6 * (y0 + y2 - 2 * y1) / dx ^ 2

sweep y ( uPrev, vPrev ) =
( 1 / (4 - uPrev), (y - vPrev) / (4 - uPrev) )

backSub ( u, v ) mNext =
v - u * mNext

piece ( y0, m0 ) ( y1, m1 ) =
{ a = (m1 - m0) / 6 / dx
, b = m0 / 2
, c = (y1 - y0) / dx - (m1 + 2 * m0) * dx / 6
, d = y0
}
in
mapTriple concavity heights
|> List.scanl sweep ( 0, 0 )
|> scanr backSub 0
|> List.map2 (,) heights
|> mapPair piece
|> Array.fromList

mapTriple : (a -> a -> a -> b) -> List a -> List b
mapTriple f x0 =
let
tail =
List.tail x0

doubleTail =
tail |> Maybe.andThen List.tail
in
case ( tail, doubleTail ) of
( Just x1, Just x2 ) ->
List.map3 f x0 x1 x2

otherwise ->
[]

mapPair : (a -> a -> b) -> List a -> List b
mapPair f x0 =
case (List.tail x0) of
Just x1 ->
List.map2 f x0 x1

otherwise ->
[]

scanr : (a -> b -> b) -> b -> List a -> List b
scanr f init =
List.reverse >> List.scanl f init >> List.reverse

{-| -}
type Spline
= Spline
{ coefficients : Array Coefficients
, start : Float
, dx : Float
}

type alias Coefficients =
{ a : Float
, b : Float
, c : Float
, d : Float
}
```
```