//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author  Hao Peng
 *  @version 1.6
 *  @date    Sun Nov 12 12:27:00 EST 2017
 *  @see     LICENSE (MIT style license file).
 *
 *  @note    Derivatives for Fourier Basis Functions
 *
 *  @see https://en.wikipedia.org/wiki/Fourier_series
 */

package scalation
package calculus

import scala.Double.NaN
import scala.math.{Pi, cos, sin}

import scalation.mathstat._

//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `DFourier` class provides Fourier basis functions with derivatives.
 *  Such basis functions are useful are useful for fitting periodic data in
 *  Functional Data Analysis.
 *  @see en.wikipedia.org/wiki/Fourier_series
 *  @param w     the fundamental frequency parameter
 *  @param mMax  the number of sin/cos pairs to be used in the basis function
 */
class DFourier (w: Double = 2.0 * Pi, mMax: Int = 4)
      extends Fourier (w, mMax)
         with DBasisFunction:

    //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Obtain the value of nth derivative of the m-th order 'j'-th basis function at time 't'.
     *  Or alternatively, obtain the nth derivative basis function by calling dnbf(n)(m)(j) only.
     *  Ex: val x = dnbf(n)(m)(j)(t) retrieves the nth derivative value of the j-th basis function at 't'.
     *      val f = dnbf(n)(m)(j)    retrieves the nth derivative of the j-th basis function.
     *  @param n  the order of the derivative
     *  @param m  the order of the basis function
     *  @param j  indicates which basis function
     *  @param t  the time parameter
     */
    def dnbf (n: Int)(m: Int)(j: Int)(t: Double): Double =
        val c  = (j+1)/2 * w                               // constant for the inner parts of the sin/cos functions
        val cn = c~^n                                      // c^n
        val cSin = if n/2 % 2 == 0 then     cn else -cn    // multiplicative constant for derivatives of Sines
        val cCos = if (n+1)/2 % 2 == 0 then cn else -cn    // multiplicative constant for derivatives of Cosines

        if j > 2*m then NaN                                // invalid input, j must be <= 2m
        else if j == 0 then if n > 0 then 0.0 else 1.0
        else if j % 2 == 0 then if n % 2 == 0 then cSin * sin (c * t) else cSin * cos (c * t)
        else if n % 2 == 0 then cCos * cos (c * t)
        else cCos * sin (c * t)
    end dnbf

end DFourier


//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `dFourierTest` object is used to test the `DFourier` class.
 *  > runMain scalation.calculus.dFourierTest
 */
@main def dFourierTest (): Unit =

    val m = 1
    val dfour = new DFourier ()
//  val t = VectorD (0.0 to 5 by 0.01)
    val t = VectorD.range (0, 501) * 0.01

    def x (n: Int, tt: Double) = (0 until 2*m+1).map (j => dfour.dnbf(n)(m)(j)(tt)).sum

    val y   = for (tt <- t) yield x(0, tt)
    val dy  = for (tt <- t) yield x(1, tt)
    val d2y = for (tt <- t) yield x(2, tt)

    new Plot (t, dy, y, "dy vs y")
    new Plot (t, d2y, dy, "d2y vs dy")

end dFourierTest