// TestDeriv.scala   order 1 through 4, points order+1 through order+3
import scala.math._

object TestDeriv
{
  var np = 7
  var D = Deriv
  def main(args: Array[String])
  {
    println("TestDeriv.scala running")
    println("TestDeriv.scala sine wave, to 4th order to np="+ np)
    println("np is number of points used to compute derivative")
    var b:Int = 0
    var hx:Double = 0.05
    var A = Array.ofDim[Int](np+1) // b is A(npoints)
    var X = Array.ofDim[Double](np)
    var Y = Array.ofDim[Double](np)
    var Yp = Array.ofDim[Double](np)
    var Ypp = Array.ofDim[Double](np)
    var Yppp = Array.ofDim[Double](np)
    var Ypppp = Array.ofDim[Double](np)
    var Cx = Array.ofDim[Double](np)
    var Cxx = Array.ofDim[Double](np)
    var Cxxx = Array.ofDim[Double](np)
    var Cxxxx = Array.ofDim[Double](np)
    var order = 1
    var up = 0.0
    var upp = 0.0
    var uppp = 0.0
    var upppp = 0.0 
    var maxerr = 0.0
    var errd = 0.0
    var errr = 0.0
    var errnu = 0.0

    b = D.gcd(36, 216)
    println("gcd(36,216)="+b)

    // set up data for all cases
    for(i <- 0 until np) // set up data and exact result
    {
      X(i) = hx*i
      Y(i) = f(X(i)) 
      Yp(i) = fp(X(i)) 
      Ypp(i) = fpp(X(i)) 
      Yppp(i) = fppp(X(i)) 
      Ypppp(i) = fpppp(X(i)) 
      println("X("+i+")="+X(i)+",Y="+Y(i)+", Yp="+Yp(i))
      println("Ypp="+Ypp(i)+",Yppp="+Yppp(i)+", Ypppp="+Ypppp(i))
    }

    order = 1
    for(npoint <- order+1 until np)
    {
      maxerr = 0.0
      for(point <- 0 until npoint)
      {
        println("order="+order+", npoint="+npoint+", point="+point)
        A = D.deriv(order, npoint, point)
        b = A(npoint)
        println("D.deriv b="+b)
        for(i <- 0 until npoint)
        {
          println("A("+i+")="+A(i))
        }
        up = 0.0
        for(j <- 0 until npoint)
        {
          up = up + Y(j)*(A(j).toDouble)   // Y(j) = f(X(j)) 
        }
        up = up * 1.0/(hx*(b.toDouble))
        errd = (Yp(point)-up).abs
        println("D.deriv error = "+errd)
        maxerr = max(maxerr, errd)
        println(" ")


        Cx = D.rderiv(order, npoint, point, hx)
        for(i <- 0 until npoint)
        {
          println("rderiv Cx("+i+")="+Cx(i))
        }
        up = 0.0
        for(j <- 0 until npoint)
        {
          up = up + Cx(j)*Y(j) // Y(j) = f(X(j)) 
        }
        errr = (Yp(point)-up).abs
        println("D.rderiv error = "+errr)
        maxerr = max(maxerr, errr)

        Cx = D.nuderiv(order, npoint, point, X)
        up = 0.0
        for(j <- 0 until npoint)
        {
          up = up + Cx(j)*Y(j) // Y(j) = f(X(j)) 
        }
        errnu = (Yp(point)-up).abs
        println("D.nuderiv error = "+errnu)
        maxerr = max(maxerr, errnu)

      }
    }
    println("order="+order+", maxerr="+maxerr)

    println("TestDeriv.scala finished")
  } // end main of TestDeriv.scala

  def f(x:Double):Double= // function
  {
    return sin(x)
  } // end f
  def fp(x:Double):Double= // first derivative 
  {
    return cos(x)
  }  // end fp
  def fpp(x:Double):Double= // second derivative 
  {
    return -sin(x)
  } // end fpp
  def fppp(x:Double):Double= // third derivative 
  {
    return -cos(x) 
  } // end fppp
  def fpppp(x:Double):Double= // fourth derivative 
  {
    return sin(x)
  } // end fpppp

} // end TestDeriv.scala