-- nuderiv.adb  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
--
with real_Arrays;
use  real_Arrays;
with inverse;

procedure nuderiv(order : integer;           -- derivative order
                  npoint : integer;          -- number of points
                  point : integer;           -- where derivative computed
                  X : real_vector;           -- X input values
                  C : in out real_vector) is -- coefficients of Y's
  n : integer := npoint;
  A : real_matrix(1..n,1..n);
  AI : real_matrix(1..n,1..n);
  pwr : real;
  fct : real_vector(1..n);
begin
  for i in 1..n loop -- build matrix that will be inverted
    pwr := 1.0;
    for j in 1..n loop
      A(i,j) := pwr;
      pwr := pwr*X(i-1+X'first);
    end loop;
  end loop;
  AI := inverse(A); -- invert the matrix
  -- form derivative for order and point
  for i in 1..n loop -- compute factors for order
    fct(i) := 1.0;
  end loop;
  for j in 1..order loop
    for i in 1..n loop
      fct(i) := fct(i)*real(i-j);
    end loop;
  end loop;
  for i in 1..n loop
    C(i-1+C'first) := 0.0;
    pwr := 1.0;
    for j in Order+1..n loop
      C(i-1+C'first) := C(i-1+C'first) + fct(j)*AI(j,i)*pwr;
      pwr := pwr * X(point); -- pass actual point
    end loop;
  end loop;
end nuderiv;