package DifferentialOperators1D "Finite Difference differential operators"

  

  function pder "Returns vector of spatial derivative values of a 1D field"

    input Fields.Field1D f;

    input Integer x=1 "Diff order - first or second derivative";

    output Real df[size(f.val;1)];

  algorithm

    df := if x == 1 then SecondOrder.diff1(f.val, f.domain.dx)

    else SecondOrder.diff2(f.val, f.domain.dx);

  end pder;



  package FirstOrder "First order polynomial derivative approximations"

  

    function diff1 "First derivative"

      input Real u[:];

      input Real dx;

      output Real du[size(u,1)];        

    algorithm         // Left, central,  and right differences

      du := cat(1, {(u[2] - u[1])/dx}, (u[3:end] - u[1:end-2])/2/dx,

                   {(u[end] - u[end-1])/dx});

    end diff1;

  end FirstOrder;  

      

  package SecondOrder "Second order polynomial derivative approximations"

  

    function diff1 "First derivative"

      input Real u[:];

      input Real dx;

      output Real du[size(u,1)];

    algorithm

      du := cat(1, {(-3*u[1] + 4*u[2] - u[3])/2/dx },

                   {u[3:end] - u[1:end-2])/2/dx,

                   {(3*u[end] - 4*u[end-1] + u[end-2])/2/dx} );

    end diff1;

    

    function diff2 "Second derivative"

      input Real u[:];

      input Real dx;

      output Real du2[size(u,1)];

    algorithm

      du2 := cat(1, {(2*u[1] - 5*u[2] + 4*u[3] - u[4])/dx/dx },

                    (u[3:end] - 2*u[2:end-1] + u[1:end-2])/dx/dx,

                    {(2*u[end] - 5*u[end-1] + 4*u[end-2] - u[end-3])/dx/dx} );

    end diff2;

  end SecondOrder;

end DifferentialOperators1D;