! nuderiv.f90  non uniformly spaced derivative coefficients
!    given  x0, x1, x2, x3  non uniformly spaced
!    find the coefficients of y0, y1, y2, y3
!    in order to compute the derivatives:
!      y'(x0) = c00*y0 + c01*y1 + c02*y2 + c03*y3
!      y'(x1) = c10*y0 + c11*y1 + c12*y2 + c13*y3
!      y'(x2) = c20*y0 + c21*y1 + c22*y2 + c23*y3
!      y'(x3) = c30*y0 + c31*y1 + c32*y2 + c33*y3
!
!  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 y0, y1, y2 and y3 variables:
!
!                              y0  y1  y2  y3
!  form:  | 1  x0  x0^2  x0^3   1   0   0   0 | = a
!         | 1  x1  x1^2  x1^3   0   1   0   0 | = b
!         | 1  x2  x2^2  x2^3   0   0   1   0 | = c
!         | 1  x3  x3^2  x3^3   0   0   0   1 | = d
!
! reduce  | 1   0     0     0  a0  a1  a2  a3 | = a
!         | 0   1     0     0  b0  b1  b2  b3 | = b
!         | 0   0     1     0  c0  c1  c2  c3 | = c
!         | 0   0     0     1  d0  d1  d2  d3 | = d
!
!                              this is just the inverse
!
!   y'(x0) =  b0*y0        + b1*y1        + b2*y2        + b3*y3        +
!              2*c0*y0*x0   + 2*c1*y1*x0   + 2*c2*y2*x0   + 2*c3*y3*x0   +
!              3*d0*y0*x0^2 + 3*d1*y1*x0^2 + 3*d2*y2*x0^2 + 3*d3*y3*x0^2
!
!   or  c00 = b0 + 2*c0*x0 + 3*d0*x0^2
!       c01 = b1 + 2*c1*x0 + 3*d1*x0^2
!       c02 = b2 + 2*c2*x0 + 3*d2*x0^2
!       c03 = b3 + 2*c3*x0 + 3*d3*x0^2
!
!   y'(x0) = c00*y0 + c01*y1 + c02*y2 +c03*y3
!
!   y'(x1) =  b0*y0        + b1*y1        + b2*y2        + b3*y3        +
!              2*c0*y0*x1   + 2*c1*y1*x1   + 2*c2*y2*x1   + 2*c3*y3*x1   +
!              3*d0*y0*x1^2 + 3*d1*y1*x1^2 + 3*d2*y2*x1^2 + 3*d3*y3*x1^2
!
!   or  c10 = b0 + 2*c0*x1 + 3*d0*x1^2
!       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
!
!       cij = bj + 2*cj*xi + 3*dj*xi^2
!
!
!   y'(x1) = c10*y0 + c11*y1 + c12*y2 +c13*y3
!
!   y''(x0) = 2*c0*y0    + 2*c1*y1    + 2*c2*y2    + 2*c3*y3    +
!              6*d0*y0*x0 + 6*d1*y1*x0 + 6*d2*y2*x0 + 6*d3*y3*x0
!
!   or  c00 = 2*c0 + 6*d0*x0
!       c01 = 2*c1 + 6*d1*x0
!       c02 = 2*c2 + 6*d2*x0
!       c03 = 2*c3 + 6*d3*x0
!
!   y'''(x0) = 6*d0*y0 + 6*d1*y1 + 6*d2*y2 + 6*d3*y3
!
!   or  c00 = 6*d0    independent of x
!       c01 = 6*d1
!       c02 = 6*d2
!       c03 = 6*d3
!

subroutine nuderiv(order, npoint, point, X, C)
  implicit none
  integer, intent(in) :: order           ! derivative order
  integer, intent(in) :: npoint          ! number of points
  integer, intent(in) :: point           ! where derivative computed
  double precision, dimension(npoint), intent(in) :: X ! X input values
  double precision, dimension(npoint), intent(inout) :: C ! coefficients of Y's
  integer :: i, j, n 
  double precision, dimension(npoint,npoint) :: A
  double precision, dimension(npoint,npoint) :: AI
  double precision :: pwr
  double precision, dimension(npoint) :: fct

  interface
    subroutine inverse(n, sz, A, AI)
      integer, intent(in) :: n  ! number of equations
      integer, intent(in) :: sz ! dimension of arrays
      double precision, dimension(sz,sz), intent(in) :: A
      double precision, dimension(sz,sz), intent(inout) :: AI
    end subroutine inverse 
  end interface

  n = npoint
  do i=1,n    ! build matrix that will be inverted
    pwr = 1.0D0
    do j=1,n
      A(i,j) = pwr
      pwr = pwr*X(i)
    end do ! j
  end do ! i
  call inverse(n, npoint, A, AI) ! invert the matrix
  ! form derivative for order and point
  do i=1,n         ! compute factors for order
    fct(i) = 1.0D0
  end do ! i
  do j=1,order
    do i=1,n
      fct(i) = fct(i)*dble(i-j)
    end do ! i
  end do ! j
  do i=1,n
    C(i) = 0.0D0
    pwr = 1.0D0
    do j=order+1,n
      C(i) = C(i) + fct(j)*AI(j,i)*pwr
      pwr = pwr * X(point) ! pass actual point
    end do ! j
  end do ! i
end subroutine nuderiv