Line data    Source code

       1             : """
       2             : Piece of an [`AbstractHermiteSpline`](@ref) that can be evaluated.
       3             : 
       4             : A piece is defined by two knots `t0, t1` and values and derivarives at
       5             : knots.
       6             : 
       7             : See also [`AbstractHermiteSpline`](@ref), [`piece`](@ref),
       8             : [`somepiece`](@ref)
       9             : """
      10             : struct HermitePiece{S, V}
      11             :     x0::V
      12             :     dx0::V
      13             :     x1::V
      14             :     dx1::V
      15             :     t0::S
      16             :     t1::S
      17             : 
      18          15 :     HermitePiece{S, V}(x0, dx0, x1, dx1, t0, t1) where {S, V} =
      19             :         new(x0, dx0, x1, dx1, t0, t1)
      20             : end
      21             : 
      22             : """
      23             :     p = HermitePiece(x0, dx0, x1, dx1, t0, t1)
      24             : 
      25             : Construct [`HermitePiece`](@ref)
      26             : """
      27          15 : HermitePiece(x0::V, dx0::V, x1::V, dx1::V,
      28             :              t0::S, t1::S) where {S<:Real, V} =
      29             :     HermitePiece{S, V}(x0, dx0, x1, dx1, t0, t1)
      30             : 
      31             : """
      32             :     q = similar(p::HermitePiece)
      33             : 
      34             : Get *undefined* [`HermitePiece`](@ref) of same type
      35             : """
      36           1 : function Base.similar(p::HermitePiece{S, V}) where {S, V}
      37           1 :     v = zeros(V)
      38           1 :     HermitePiece(v, v, v, v, S(Inf), S(Inf))
      39             : end
      40             : 
      41             : """
      42             :     cx0, cdx0, cx1, cdx1 = hermitebasis(t)
      43             : 
      44             : Evaluate Hermite basis for computing *values* as weighted sum of
      45             : coefficients `x0, dx0, x1, dx1` for `0 <= t <= 1`.
      46             : """
      47        5970 : function hermitebasis(t)
      48        5970 :     cx0  = @evalpoly(t, 1, 0, -3, 2) # (1 + 2*t) * ((1-t)*(1-t))
      49        5970 :     cdx0 = @evalpoly(t, 0, 1, -2, 1) # t * ((1-t)*(1-t))  ???
      50        5970 :     cx1  = @evalpoly(t, 0, 0, 3, -2) # (t*t) * (3 - 2*t)
      51        5970 :     cdx1 = @evalpoly(t, 0, 0, -1, 1) # (t*t) * (t-1)
      52             : 
      53        5970 :     cx0, cdx0, cx1, cdx1
      54             : end
      55             : 
      56             : """
      57             :     ex0, edx0, ex1, edx1 = hermitederivativebasis(t)
      58             : 
      59             : Evaluate Hermite basis for computing *derivatives* as weighted sum of
      60             : coefficients `x0, dx0, x1, dx1` for `0 <= t <= 1`.
      61             : """
      62        3171 : function hermitederivativebasis(t)
      63        3171 :     ex0  = @evalpoly(t, 0, -6, 6) # t * (t-1) * 6
      64        3171 :     edx0 = @evalpoly(t, 1, -4, 3) # (t-1) * (3*t - 1)
      65        3171 :     ex1  = -ex0
      66        3171 :     edx1 = @evalpoly(t, 0, -2, 3) # (3*t - 2) * t
      67             : 
      68        3171 :     ex0, edx0, ex1, edx1
      69             : end
      70             : 
      71             : """
      72             :     p::HermitePiece
      73             :     x, dx = p(t, ValueDerivative)
      74             : 
      75             : Evaluate [`HermitePiece`](@ref) at `t`.
      76             : """
      77         693 : function (s::HermitePiece{S, V})(t::S, ::TValueDerivative) where {S, V}
      78         693 :     (t == s.t0) && return (s.x0, s.dx0)
      79         693 :     (t == s.t1) && return (s.x1, s.dx1)
      80             : 
      81         693 :     h = s.t1 - s.t0
      82             : 
      83         693 :     t = (t - s.t0) / h
      84         693 :     @assert isfinite(t)
      85             : 
      86         693 :     cx0, cdx0, cx1, cdx1 = hermitebasis(t)
      87         693 :     x = s.x0*cx0 + s.dx0*(cdx0*h) + s.x1*cx1 + s.dx1*(cdx1*h)
      88             : 
      89         693 :     ex0, edx0, ex1, edx1 = hermitederivativebasis(t)
      90         693 :     dx = s.x0*(ex0/h) + s.dx0*edx0 + s.x1*(ex1/h) + s.dx1*edx1
      91             : 
      92         693 :     x, dx
      93             : end
      94             : 
      95             : """
      96             :     p::HermitePiece
      97             :     x = p(t, Value)
      98             : 
      99             : Evaluate [`HermitePiece`](@ref) at `t`.
     100             : """
     101           3 : function (s::HermitePiece{S, V})(t::S, ::TValue) where {S, V}
     102           3 :     (t == s.t0) && return s.x0
     103           3 :     (t == s.t1) && return s.x1
     104             : 
     105           3 :     h = s.t1 - s.t0
     106           3 :     t = (t - s.t0) / h
     107             : 
     108           3 :     @assert isfinite(t)
     109             : 
     110           3 :     cx0, cdx0, cx1, cdx1 = hermitebasis(t)
     111           3 :     x = s.x0*cx0 + s.dx0*(cdx0*h) + s.x1*cx1 + s.dx1*(cdx1*h)
     112             : 
     113           3 :     x
     114             : end
     115             : 
     116             : """
     117             :     p::HermitePiece
     118             :     t in p
     119             : 
     120             : Is `t` in [`knots`](@ref) interval of [`HermitePiece`](@ref) `p`?
     121             : """
     122           8 : function Base.in(t::S, s::HermitePiece{S, V}) where {S, V}
     123           8 :     (s.t0 <= s.t1) && return (s.t0 <= t <= s.t1)
     124           3 :     (s.t1 <= t <=s.t0)
     125             : end