! nuderivC.f90  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
!

subroutine nuderivC(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 complex, dimension(npoint), intent(in) :: X ! X input values
  double complex, dimension(npoint), intent(inout) :: C ! coefficients of Y's
  integer :: i, j, n 
  double complex, dimension(npoint,npoint) :: A
  double complex, dimension(npoint,npoint) :: AI
  double complex :: pwr
  double complex, dimension(npoint) :: fct

  interface
    subroutine inverseC(n, sz, A, AI)
      integer, intent(in) :: n  ! number of equations
      integer, intent(in) :: sz ! dimension of arrays
      double complex, dimension(sz,sz), intent(in) :: A
      double complex, dimension(sz,sz), intent(inout) :: AI
    end subroutine inverseC 
  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 inverseC(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 nuderivC