-- laphi.adb  Lagrange Phi functions, orthogonal polynomials
--         function prototype naming convention:
--           phi            basic Lagrange polynomial
--           phip           first derivative "p for prime"
--           phipp          second derivative,  etc
--
-- example  phip(x, i, n, xg()) is
--                     ^__ nth-1 degree
--                  ^__ ith polynomial in set
--             ^___first derivative "prime"
--  xg() is x  grid values, not necessarily uniformly spaced.
--  xg(1) is minimum, left boundary.
--  xg(n) is maximum, right boundary.
--  All grid coordinates must be in increasing order.

with Real_Arrays; -- types real and real_vector and real_matrix
use Real_Arrays;

package body laphi is
  -- phi
  function phi(x:real; j:integer; n:integer; xg:real_vector) return real is
    y : real;
    denom : real := 1.0;
  begin
    for i in 1..n loop
      if i /= j then denom := denom * (xg(j)-xg(i)); end if;
    end loop;
    -- evaluate function, n,j,x,xg,denom
    y := 1.0; -- accumulate product
    for i in 1..n loop
      if i /= j  then y := y * (x-xg(i)); end if;
    end loop;
    y := y/denom;
    return y;
  end phi;

  -- first derivative
  function phip(x:real; j:integer; n:integer; xg:real_vector) return real is
    yp, prod : real;
    denom : real := 1.0;
  begin
    for i in 1..n loop
      if i /= j then denom := denom * (xg(j)-xg(i)); end if;
    end loop;
    -- evaluate first derivative of function, n,j,x,xg,denom
    yp := 0.0; -- accumulate sum of products
    for i in 1..n loop
      if i /= j then
        prod := 1.0;
        for ii in 1..n loop
          if ii/=i and ii/=j then
            prod := prod * (x-xg(ii));
          end if;
        end loop;
        yp := yp + prod;
      end if;
    end loop;
    yp := yp/denom; -- have first derivative
    return yp;
  end phip;


  -- second derivative
  function phipp(x:real; j:integer; n:integer; xg:real_vector) return real is
    ypp, prod : real;
    denom : real := 1.0;
  begin
    for i in 1..n loop
      if i /= j then denom := denom * (xg(j)-xg(i)); end if;
    end loop;
    -- evaluate second derivative of function, n,j,x,xg,denom
    ypp := 0.0; -- accumulate sum of products
    for i in 1..n loop
      if i /= j then
        for ii in 1..n loop
          if ii/=i and ii/=j then
            prod := 1.0;
            for iii in 1..n loop
              if iii/=i and iii/=ii and iii/=j then
                prod := prod * (x-xg(iii));
              end if;
            end loop;
            ypp := ypp + prod;
          end if;
        end loop;
      end if;
    end loop;
    ypp := ypp/denom; -- have second derivative
    return ypp;
  end phipp;


  -- third derivative
  function phippp(x:real; j:integer; n:integer; xg:real_vector) return real is
    yppp, prod : real;
    denom : real := 1.0;
  begin
    for i in 1..n loop
      if i /= j then denom := denom * (xg(j)-xg(i)); end if;
    end loop;
    -- evaluate third derivative of function, n,j,x,xg,denom
    yppp := 0.0; -- accumulate sum of products
    for i in 1..n loop
      if i /= j then
        for ii in 1..n loop
          if ii/=i and ii/=j then
            for iii in 1..n loop
              if iii/=i and iii/=ii and iii/=j then
                prod := 1.0;
                for iiii in 1..n loop
                  if iiii/=i and iiii/=ii and iiii/=iii and iiii/=j then
                    prod := prod * (x-xg(iiii));
                  end if;
                end loop;
                yppp := yppp + prod;
              end if;
            end loop;
          end if;
        end loop;
      end if;
    end loop;
    yppp := yppp/denom;  -- have third derivative
    return yppp;
  end phippp;


  -- fourth derivative
  function phipppp(x:real; j:integer; n:integer; xg:real_vector) return real is
    ypppp, prod : real;
    denom : real := 1.0;
  begin
    for i in 1..n loop
      if i /= j then denom := denom * (xg(j)-xg(i)); end if;
    end loop;
    -- evaluate fourth derivative of function, n,j,x,xg,denom
    ypppp := 0.0; -- accumulate sum of products
    for i in 1..n loop
      if i /= j then
        for ii in 1..n loop
          if ii/=i and ii/=j then
            for iii in 1..n loop
              if iii/=i and iii/=ii and iii/=j then
                for iiii in 1..n loop
                  if iiii/=i and iiii/=ii and iiii/=iii and iiii/=j then
                    prod := 1.0;
                    for iiiii in 1..n loop
                      if iiiii/=i and iiiii/=ii and iiiii/=iii and
                         iiiii/=iiii and iiiii/=j then
                        prod := prod * (x-xg(iiiii));
                      end if;
                    end loop;
                    ypppp := ypppp + prod;
                  end if;
                end loop;
              end if;
            end loop;
          end if;
        end loop;
      end if;
    end loop;
    ypppp := ypppp/denom; -- have fourth derivative
    return ypppp;
  end phipppp;

end Laphi;