# nuderiv.jl  non uniformly spaced derivative coefficients
#    given  x1, x2, x3, x4  non uniformly spaced
#    find the coefficients of y1, y2, y3, y4
#    in order to compute the derivatives:
#      y'(x1) := c11*y1 + c12*y2 + c13*y3 + c14*y4
#      y'(x2) := c21*y1 + c22*y2 + c23*y3 + c24*y4
#      y'(x3) := c31*y1 + c32*y2 + c33*y3 + c34*y4
#      y'(x4) := c41*y1 + c42*y2 + c43*y3 + c44*y4
#
#  Method:  fit data to  y(x)    := a + b*x + c*x^2 + d*x^3
#                        y'(x)   :=     b   + 2*c*x + 3*d*x^2
#                        y''(x)  :=           2*c   + 6*d*x
#                        y'''(x) :=                   6*d
#
#           with y1, y2, y3 and y4 variables:
#
#                              y1  y2  y3  y4
#  form:  | 1  x1  x1^2  x1^3   1   0   0   0 | := a
#         | 1  x2  x2^2  x2^3   0   1   0   0 | := b
#         | 1  x3  x3^2  x3^3   0   0   1   0 | := c
#         | 1  x4  x4^2  x4^3   0   0   0   1 | := d
#
# reduce  | 1   0     0     0  a1  a2  a3  a4 | := a
#         | 0   1     0     0  b1  b2  b3  b4 | := b
#         | 0   0     1     0  c1  c2  c3  c4 | := c
#         | 0   0     0     1  d1  d2  d3  d4 | := d
#
#                              this is just the inverse
#
#   y'(x1) :=  b1*y1        + b2*y2        + b3*y3        + b4*y4        +
#              2*c1*y1*x1   + 2*c2*y2*x1   + 2*c3*y3*x1   + 2*c4*y4*x1   +
#              3*d1*y1*x1^2 + 3*d2*y2*x1^2 + 3*d2*y2*x1^2 + 3*d4*y4*x1^2
#
#   or  c11 := b1 + 2*c1*x1 + 3*d1*x1^2
#       c12 := b2 + 2*c2*x1 + 3*d2*x1^2
#       c13 := b3 + 2*c3*x1 + 3*d3*x1^2
#       c14 := b4 + 2*c4*x1 + 3*d4*x1^2
#
#   y'(x1) := c11*y1 + c12*y2 + c13*y3 +c14*y4
#
#   y'(x2) :=  b1*y1        + b2*y2        + b3*y3        + b4*y4        +
#              2*c1*y1*x2   + 2*c2*y2*x2   + 2*c3*y3*x2   + 2*c4*y4*x2   +
#              3*d1*y1*x2^2 + 3*d2*y2*x2^2 + 3*d3*y3*x2^2 + 3*d4*y4*x2^2
#
#   or  c21 := b1 + 2*c1*x2 + 3*d1*x2^2
#       c22 := b2 + 2*c2*x2 + 3*d2*x2^2
#       c23 := b3 + 2*c3*x2 + 3*d3*x2^2
#       c24 := b4 + 2*c4*x2 + 3*d4*x2^2
#
#       cij := bj + 2*cj*xi + 3*dj*xi^2
#
#
#   y'(x2) := c21*y1 + c22*y2 + c23*y3 +c24*y4
#
#   y''(x1) := 2*c1*y1    + 2*c2*y2    + 2*c3*y3    + 2*c4*y4    +
#              6*d1*y1*x1 + 6*d2*y2*x1 + 6*d3*y3*x1 + 6*d4*y4*x1
#
#   or  c11 := 2*c1 + 6*d1*x1
#       c12 := 2*c1 + 6*d2*x1
#       c13 := 2*c2 + 6*d3*x1
#       c14 := 2*c3 + 6*d4*x1
#
#   y'''(x1) := 6*d1*y1 + 6*d2*y2 + 6*d3*y3 + 6*d4*y4
#
#   or  c11 := 6*d1    independent of x
#       c12 := 6*d2
#       c13 := 6*d3
#       c14 := 6*d4
#
using LinearAlgebra   # needs  inv
using Printf

function nuderiv(order, npoint, point, X, C)
  A = Array{Float64,2}(undef, npoint, npoint)
  fct = Array{Float64,1}(undef, npoint)
  n = npoint
  
  for i in 1:n    # build matrix that will be inverted
    pwr = 1.0
    for j in 1:n
      A[i,j] = pwr
      pwr = pwr*X[i]
    end  # j
  end  # i
  AI = inv(A)
  # form derivative for order and point
  for i in 1:n         # compute factors for order
    fct[i] = 1.0
  end  # i
  for j in 1:order
    for i in 1:n
      fct[i] = fct[i]*(i-j)
    end  # i
  end  # j
  for i in 1:n
    C[i] = 0.0
    pwr = 1.0
    for j in order+1:n
      C[i] = C[i] + fct[j]*AI[j,i]*pwr
      pwr = pwr * X[point] # pass actual point
    end  # j
  end  # i
end # function test_nuderiv

function vecprt(name,V) # works for row and column vector
  n = size(V,1)
  if n>1
    for i in 1:n
      @printf "%s[%d]=%f \n" name i  V[i]
    end
    println(" ")
  else
    n = size(V,2)
    for i in 1:n
      @printf "%s[1,%d]=%f \n" name i  V[1,i]
    end
    println(" ")
  end # end if else
  return nothing
end # end vecprt

function matprt(name,A) # works for 2D matrix
  nrow = size(A,1)
  ncol = size(A,2)
  for i in 1:nrow    # default subscripts start at 1 like Fortran and Matlab
    for j in 1:ncol
      @printf "%s[%1d][%1d]=%8.5f \n" name i  j A[i,j]
    end
  end
  println(" ")
  return nothing
end # end matprt

function f(x)
  return x*x*x + 2.0*x*x + 3.0*x + 4.0
end # end f(x)

function fx(x)
  return 3.0*x*x + 4.0*x + 3.0
end # end f(x)

function fxx(x)
  return 6.0*x + 4.0
end # end f(x)

println("test_nuderiv.jl running")
println("f(x) = x*x*x + 2.0*x*x + 3.0*x + 4.0 ")
println("f'(x) = fx(x) = 3.0*x*x + 4.0*x + 3.0 ")
println("f''(x) = fxx(x) = 6.0*x + 4.0 ")
println(" ")


order = 1
npoint = 9
point = 5
@printf "order=%d  npoint=%d  point=%d \n" order npoint point
X = [1.0, 1.1, 1.4, 1.6, 1.8, 1.9, 2.1, 2.2, 2.3] 
C = zeros(Float64, npoint)
vecprt("X",X)
nuderiv(order, npoint, point, X, C)

vecprt("C",C)
deriv = 0.0
for i in 1:npoint
  global deriv
  deriv += C[i]*f(X[i])
end
@printf "deriv=%f at x=%f \n" deriv X[point]
@printf "true =%f \n" fx(X[point])
@printf " \n"

order = 2
npoint = 9
point = 4
@printf "order=%d  npoint=%d  point=%d \n" order npoint point
X = [1.0, 1.1, 1.4, 1.6, 1.8, 1.9, 2.1, 2.2, 2.3] 
C = zeros(Float64, npoint)
nuderiv(order, npoint, point, X, C)

vecprt("C",C)
deriv = 0.0
for i in 1:npoint
  global deriv
  deriv += C[i]*f(X[i])
end
@printf "deriv=%f at x=%f \n" deriv X[point]
@printf "true =%f \n" fxx(X[point])
@printf " \n"

println("test_nuderiv.jl finished")
# end test_nuderiv.jl