## Modelica_LinearSystems2.Math

Additional functions for Modelica.Math

### Information

Extends from Modelica.Icons.Library (Icon for library).

### Package Content

Name | Description
--- | ---
Complex | Record defining a Complex number
Polynomial | Record defining the data for a polynomial
Matrices | Additional functions for Modelica.Math
Vectors | Functions operating on vectors

## Modelica_LinearSystems2.Math.Complex

Record defining a Complex number

### Information

This record defines a complex number consisting of a real and an imaginary part. Additionally, various utility functions are provided in the record to operate on instances of this record. Example (note: "j" in the comments is defined as j=sqrt(-1)):

       import Modelica.Utilities.Streams;
       import Modelica_LinearSystems2.Math.Complex;
       Complex c1=Complex(re=2, im=3) "= 2 + 3j";
       Complex c2=Complex(3,4) "= 3 + 4j";
     algorithm
       c3 := Complex.'+'(c1, c2) "= c1 + c2";
       Streams.print("c3 = " + Complex.'String'(c3));
       Streams.print("c3 = " + Complex.'String'(c3,"i"));
       // This gives the following print-out:
       c3 = 5 + 7j
       c3 = 5 + 7i

The utility functions are written in such a way that it is convenient to work with them, once operator overloading is provided in Modelica and Modelica tools. Example:

       // Assume operator overloading is available (this is not yet the case):
       Complex j  = Complex.j();
       Complex c4 = -2 + 5*j;
       // A Modelica tool will transform the statement of c4 into
       Complex c4 = Complex.'+'( Complex.(-2,0), Complex.'*'(Complex(5,0),j)));

The utility functions are implemented so that they can be easily inlined by a tool. In such a case, the above statement will not lead to any overhead.

Extends from Modelica.Icons.Record (Icon for a record).

### Parameters

Type | Name | Default | Description
--- | --- | --- | ---
Real | re |  | Real part of complex number
Real | im |  | Imaginary part of complex number

### Modelica definition

    record Complex "Record defining a Complex number"
      extends Modelica.Icons.Record;
      Real re "Real part of complex number";
      Real im "Imaginary part of complex number";

     encapsulated package Examples
        "Library demonstrating the usage of complex numbers"

        import Modelica;
        extends Modelica.Icons.Library;

        function addTwoComplexNumbers "Show how to add 2 complex number"
          import Modelica_LinearSystems2.Math.Complex;
          import Modelica.Utilities.Streams;

          output Boolean ok;
        protected
          Complex j = Complex.j();
          Complex c1=2+3*j;
          Complex c2=3+4*j;
          Complex c3;
        algorithm

          Streams.print("c1 = " + String(c1));
          Streams.print("c2 = " + String(c2));

          c3 := c1 + c2;
          Streams.print("c3 = c1 + c2 = " + String(c3));
          ok := true;
        end addTwoComplexNumbers;

     end Examples;

    encapsulated package Vectors
     function print
          import Modelica.Utilities.Streams.print;
          import Modelica_LinearSystems2.Math.Complex;

      input String name="" "Name of complex vector";
      input Complex c[:] "Complex vector to be printed";

      input String fileName=""
            "Print to fileName; empty fileName prints to the log window";

     algorithm
      print(name + " =", fileName);
      for i in 1:size(c, 1) loop
         print("   " + String(c[i]), fileName);
      end for;
     end print;

    function length "Return length of a complex vector"
          import Modelica_LinearSystems2.Math.Complex;
      input Complex v[:] "Vector";
      output Real result "Length of vector v";

    algorithm
      result := sqrt(sum(v[i].re^2 + v[i].im^2 for i in 1:size(v,1)));
    end length;

    function norm "Returns the norm of a complex vector"
          import Modelica;
          import Modelica_LinearSystems2.Math.Complex;
      input Complex v[:] "Vector";
      input Real p(min=1) = 2
            "Type of p-norm (often used: 1, 2, or Modelica.Constants.inf)";
      output Real result "p-norm of vector v";

    algorithm
      if p == 2 then
        result:= sqrt(sum(v[i].re^2 + v[i].im^2 for i in 1:size(v,1)));
      elseif p == Modelica.Constants.inf then
        result:= Complex.'max'(v);
      elseif p == 1 then
        result:= sum(Complex.'abs'(v[i]) for i in 1:size(v,1));
      else
        result:=(sum(Complex.'abs'(v[i])^p for i in 1:size(v, 1)))^(1/p);
      end if;
    end norm;

    function normalize
          "Return normalized complex vector such that length = 1 and prevent zero-division for zero vector"
          import Modelica;
          import Modelica_LinearSystems2.Math.Complex;

      input Complex v[:] "Vector";
      input Real eps = 100*Modelica.Constants.eps "if |v| < eps then result = v";
      output Complex result[size(v, 1)] "Input vector v normalized to length=1";

        protected
      Real length_v = Complex.Vectors.length(v);
      Complex j = Complex.j();
    algorithm
      if length_v >= eps then
         for i in 1:size(v,1) loop
             result[i] :=v[i].re/length_v + (v[i].im/length_v)*j;
         end for;
      else
         result :=v;
      end if;
    end normalize;

    function sortComplex "Sort elements of complex vector"
          import Modelica_LinearSystems2.Math.Complex;
      input Complex v[:] "Vector to be sorted";
      input Boolean ascending = true
            "= true if ascending order, otherwise descending order";
      input Boolean sortFrequency=true
            "= true, if sorting is first for imaginary then for real value; = false, if sorting is for absolute value";
      output Complex sorted_v[size(v,1)] = v "Sorted vector";
      output Integer indices[size(v,1)] = 1:size(v,1) "sorted_v = v[indices]";

      /* shellsort algorithm; should be improved later */
        protected
      Integer gap;
      Integer i;
      Integer j;
      Complex wv;
      Integer wi;
      Integer nv = size(v,1);
      Boolean swap;
      Integer k1;
      Integer k2;
    algorithm
      gap := div(nv,2);

      while gap > 0 loop
         i := gap;
         while i < nv loop
            j := i-gap;
            if j>=0 then
               k1 := j+1;
               k2 := j + gap + 1;
               if sortFrequency then
                  if ascending then
                     swap := abs(sorted_v[k1].im) >  abs(sorted_v[k2].im) or
                             abs(sorted_v[k1].im) == abs(sorted_v[k2].im) and
                             (sorted_v[k1].re  > sorted_v[k2].re or
                              sorted_v[k1].re  == sorted_v[k2].re and sorted_v[k1].im < sorted_v[k2].im);
                  else
                     swap := abs(sorted_v[k1].im) <  abs(sorted_v[k2].im) or
                             abs(sorted_v[k1].im) == abs(sorted_v[k2].im) and
                             (sorted_v[k1].re  < sorted_v[k2].re or
                              sorted_v[k1].re  == sorted_v[k2].re and sorted_v[k1].im < sorted_v[k2].im);
                  end if;
               else
                  if ascending then
                     swap := Complex.'abs'(sorted_v[k1]) > Complex.'abs'(sorted_v[k2]);
                  else
                     swap := Complex.'abs'(sorted_v[k1]) < Complex.'abs'(sorted_v[k2]);
                  end if;
               end if;
            else
               swap := false;
            end if;

            while swap loop
               wv := sorted_v[j+1];
               wi := indices[j+1];
               sorted_v[j+1] := sorted_v[j+gap+1];
               sorted_v[j+gap+1] := wv;
               indices[j+1] := indices[j+gap+1];
               indices[j+gap+1] := wi;
               j := j - gap;
               if j >= 0 then
                  k1 := j+1;
                  k2 := j + gap + 1;
                  if sortFrequency then
                     if ascending then
                        swap := abs(sorted_v[k1].im) >  abs(sorted_v[k2].im) or
                                abs(sorted_v[k1].im) == abs(sorted_v[k2].im) and
                                (sorted_v[k1].re  > sorted_v[k2].re or
                                 sorted_v[k1].re  == sorted_v[k2].re and sorted_v[k1].im < sorted_v[k2].im);
                     else
                        swap := abs(sorted_v[k1].im) <  abs(sorted_v[k2].im) or
                                abs(sorted_v[k1].im) == abs(sorted_v[k2].im) and
                                (sorted_v[k1].re  < sorted_v[k2].re or
                                 sorted_v[k1].re  == sorted_v[k2].re and sorted_v[k1].im < sorted_v[k2].im);
                     end if;
                  else
                     if ascending then
                        swap := Complex.'abs'(sorted_v[k1]) > Complex.'abs'(sorted_v[k2]);
                     else
                        swap := Complex.'abs'(sorted_v[k1]) < Complex.'abs'(sorted_v[k2]);
                     end if;
                  end if;
               else
                  swap := false;
               end if;
            end while;
            i := i + 1;
         end while;
         gap := div(gap,2);
      end while;
    end sortComplex;

    function multiply "scalar product of two complex vectors"
      extends Modelica.Icons.Function;
          import Modelica_LinearSystems2.Math.Complex;
          import Modelica;

      input Complex v1[:];
      input Complex v2[size(v1,1)];
      output Complex result=Complex(0);

    algorithm
      for i in 1:size(v1,1) loop
        result := result + v1[i]*v2[i];
      end for;

    end multiply;

        function reverse
          "Reverse vector elements (e.g. v[1] becomes last element)"
        extends Modelica.Icons.Function;

          import Modelica;
          import Modelica_LinearSystems2.Math.Complex;

        input Complex v[:] "Vector";
        output Complex result[size(v, 1)] "Elements of vector v in reversed order";

        algorithm
        result := {v[end-i+1] for i in 1:size(v,1)};
        end reverse;

    end Vectors;

    encapsulated package Matrices
        import Modelica;
        import Modelica_LinearSystems2;

    function print "print matrix"
          import Modelica_LinearSystems2.StateSpace;
          import Modelica.Utilities.Strings;
          import Modelica_LinearSystems2.Math.Complex;

      input Complex M[:,:];
      input Integer significantDigits=6
            "Number of significant digits that are shown";
      input String name="M" "Independent variable name used for printing";
      output String s="";
        protected
      String blanks=Strings.repeat(significantDigits);
      String space=Strings.repeat(8);
      String space2=Strings.repeat(3);
      Integer r=size(M, 1);
      Integer c=size(M, 2);

    algorithm
      if r == 0 or c == 0 then
        s := name + " = []";
      else
        s := "\n" + name + " = \n";
        for i in 1:r loop
          s := s + space;
          for j in 1:c loop
            if Complex.'abs'(M[i, j]) >= 0 then
              s := s + " ";
            end if;
            s := s + String(M[i, j], significantDigits=significantDigits) +
              Strings.repeat(significantDigits + 8 - Strings.length(String(Complex.'abs'(M[i,j]))));

          end for;
          s := s + "\n";
        end for;

      end if;
    end print;

    encapsulated function matMatMul "Multiply two complex matrices"
          import Modelica_LinearSystems2.Math.Complex;
          import Re = Modelica_LinearSystems2.Math.Complex.real;
          import Im = Modelica_LinearSystems2.Math.Complex.imag;

      input Complex m1[:,:] "Complex matrix 1";
      input Complex m2[size(m1, 2),:] "Complex matrix 2";
      output Complex m3[size(m1, 1),size(m2, 2)] "= m1*m2";
    //   Complex j=Complex.j();
    //   Real m1Re[size(m1,1),size(m1,2)]=Re(m1);
    //   Real m1Im[size(m1,1),size(m1,2)]=Im(m1);
    //   Real m2Re[size(m2,1),size(m2,2)]=Re(m2);
    //   Real m2Im[size(m2,1),size(m2,2)]=Im(m2);
        protected
      Integer l1;
      Integer l2;
    algorithm
    //  m3 := m1[:, :].re*m2[:, :].re - m1[:, :].im*m2[:, :].im + j*(m1[:, :].re*m2[:,:].im + m1[:, :].im*m2[:, :].re);
    //  m3 :=m1Re*m2Re - m1Im*m2Im + j*(m1Re*m2Im + m1Im*m2Re);

     for l1 in 1:size(m1,1) loop
       for l2 in 1:size(m2,2) loop
         m3[l1,l2] := Complex.Vectors.multiply(m1[l1,:],m2[:,l2]);
       end for;
       end for;

    end matMatMul;

    encapsulated function matVecMul
          "Multiply a complex matrices with a complex vector"
          import Modelica_LinearSystems2.Math.Complex;
    //  import Re = Modelica_LinearSystems2.Math.Complex.real;
    //  import Im = Modelica_LinearSystems2.Math.Complex.imag;

      input Complex m[:,:] "Complex matrix";
      input Complex vi[size(m, 2)] "Complex vector";
      output Complex vo[size(m, 1)] "= m*vi";
    //   Complex j=Complex.j();
    //   Real Rem[size(m, 1),size(m, 2)]=Re(m);
    //   Real Imm[size(m, 1),size(m, 2)]=Im(m);
    //   Real Revi[size(vi, 1)]=Re(vi);
    //   Real Imvi[size(vi, 1)]=Im(vi);
        protected
      Integer l1;

    algorithm
    //  vo := Rem*Revi - Imm*Imvi + j*(Rem*Imvi + Imm*Revi);

    //    for l1 in 1:size(m, 1) loop
    //      vo[l1] := Complex(0);
    //      for l2 in 1:size(m, 2) loop
    //        vo[l1] := vo[l1] + m[l1, l2]*vi[l2];
    //      end for;//l2
    //    end for;  //l1

    for l1 in 1:size(m, 1) loop
      vo[l1] := Complex.Vectors.multiply(m[l1,:],vi);
    end for;

    end matVecMul;

    end Matrices;

      encapsulated operator 'constructor'
        function fromReal
          import Modelica_LinearSystems2.Math.Complex;
          input Real re "Real part of complex number";
          input Real im=0 "Imaginary part of complex number";
          output Complex result(re=re, im=im) "Complex number";
        algorithm

        end fromReal;
      end 'constructor';

      encapsulated operator '-' "Unary and binary minus"
        function negate "Unary minus (multiply complex number by -1)"
          import Modelica_LinearSystems2.Math.Complex; // changed to Modelica_LinearSystems2
          input Complex c1 "Complex number";
          output Complex c2 "= -c1";
        algorithm
          c2 := Complex(-c1.re, -c1.im);
        end negate;

        function subtract "Subtract two complex numbers"
          import Modelica_LinearSystems2.Math.Complex;// changed to Modelica_LinearSystems2
          input Complex c1 "Complex number 1";
          input Complex c2 "Complex number 2";
          output Complex c3 "= c1 - c2";
        algorithm
          c3 := Complex(c1.re - c2.re, c1.im - c2.im);
        end subtract;
      end '-';

      encapsulated operator function '+' "Add two complex numbers"
        import Modelica_LinearSystems2.Math.Complex;

          input Complex c1 "Complex number 1";
          input Complex c2 "Complex number 2";
          output Complex c3 "= c1 + c2";
      algorithm
        c3 := Complex(c1.re + c2.re, c1.im + c2.im);
      end '+';

      encapsulated operator function '*' "Multiply two complex numbers"
        import Modelica_LinearSystems2.Math.Complex;

          input Complex c1 "Complex number 1";
          input Complex c2 "Complex number 2";
          output Complex c3 "= c1*c2";
      algorithm
          c3 := Complex(c1.re*c2.re - c1.im*c2.im, c1.re*c2.im + c1.im*c2.re);
      end '*';

      encapsulated operator function '/' "Divide two complex numbers"
        import Modelica_LinearSystems2.Math.Complex;

          input Complex c1 "Complex number 1";
          input Complex c2 "Complex number 2";
          output Complex c3 "= c1/c2";
      algorithm
          c3 := Complex((c1.re*c2.re + c1.im*c2.im)/(c2.re^2 + c2.im^2), (-c1.re*c2.im + c1.im*c2.re)/(c2.re^2 + c2.im^2));
      end '/';

      encapsulated operator function '=='
        "Test whether two complex numbers are identical"
        import Modelica_LinearSystems2.Math.Complex;
          input Complex c1 "Complex number 1";
          input Complex c2 "Complex number 2";
          output Boolean result "c1 == c2";
      algorithm
          result := c1.re == c2.re and c1.im == c2.im;
      end '==';

      encapsulated operator function 'String'
        "Transform Complex number into a String representation"
        import Modelica_LinearSystems2.Math.Complex;
          input Complex c
          "Complex number to be transformed in a String representation";
          input String name="j"
          "Name of variable representing sqrt(-1) in the string";
          input Integer significantDigits=6
          "Number of significant digits that are shown";
          output String s="";
      algorithm
          s := String(c.re, significantDigits=significantDigits);
          if c.im <> 0 then
            if c.im > 0 then
              s := s + " + ";
            else
              s := s + " - ";
            end if;
            s := s + String(abs(c.im), significantDigits=significantDigits) + "*" + name;
          end if;
        //    end toString;
      end 'String';

      encapsulated function j "Returns sqrt(-1)"
        import Modelica_LinearSystems2.Math.Complex;

        output Complex c "= sqrt(-1)";
      algorithm
        c := Complex(0,1);
      end j;

      encapsulated function 'abs' "Absolute value of complex number"
        import Modelica_LinearSystems2.Math.Complex;

        input Complex c "Complex number";
        output Real result "= abs(c)";
      algorithm
        result := (c.re^2 + c.im^2)^0.5; //changed from sqrt
      end 'abs';

      encapsulated function 'sqrt' "Square root of complex number"
        import Modelica.Math;
        import Modelica_LinearSystems2.Math.Complex;

        input Complex c1 "Complex number";
        output Complex c2 "= sqrt(c1)";
      algorithm
        c2 := Complex(sqrt(Complex.'abs'(c1))*Math.cos(Complex.arg(c1)/2), sqrt(
          Complex.'abs'(c1))*Math.sin(Complex.arg(c1)/2));
      end 'sqrt';

    encapsulated function 'max'
        "Return maximum element of complex vector"
        import Modelica_LinearSystems2.Math.Complex;

      input Complex v[:] "Vector";
      output Real result "Element of v with largest absolute value";
      output Integer index "v[index] has the largest absolute value";

      protected
      Real absv_i;
    algorithm
      if size(v,1) > 0 then
         result := Complex.'abs'(v[1]);
         index  := 1;
         for i in 2:size(v,1) loop
            absv_i :=Complex.'abs'(v[i]);
            if absv_i > result then
               result := absv_i;
               index := i;
            end if;
         end for;
      else
         result := 0;
         index  := 0;
      end if;
    end 'max';

      encapsulated function exp "Exponential of complex number"
        import Modelica_LinearSystems2.Math.Complex;
        import Modelica.Math;

        input Complex c1 "Complex number";
        output Complex c2 "= exp(c1)";
      algorithm
         c2 := Complex(Math.exp(c1.re)*Math.cos(c1.im), Math.exp(c1.re)*Math.sin(c1.im));
      end exp;

      encapsulated function log "Logarithm of complex number"

        import Modelica;
        import Modelica_LinearSystems2.Math.Complex;

        input Complex c1 "Complex number";
        output Complex c2 "= log(c1)";
      algorithm
        c2 := Complex(Modelica.Math.log(Complex.'abs'(c1)), Complex.arg(c1));
      end log;

      encapsulated function sin "Sine of complex number"

        import Modelica_LinearSystems2.Math.Complex;

        input Complex c1 "Complex number";
        output Complex c2 "sin(c1)";

      algorithm
         c2 := (Complex.exp(Complex(-c1.im, +c1.re)) - Complex.exp(Complex(+c1.im, -c1.re)))/
          Complex(0, 2);

      end sin;

      encapsulated function cos "Cosine of complex number"
        import Modelica_LinearSystems2.Math.Complex;
        import Modelica_LinearSystems2.Math.Complex.exp;

        input Complex c1 "Complex number";
        output Complex c2 "= cos(c1)";

      algorithm
       c2 := (exp(Complex(-c1.im, +c1.re)) + exp(Complex(+c1.im, -c1.re)))/2;
      end cos;

      encapsulated function arg "Phase angle of complex number"
        import Modelica;
        import Modelica_LinearSystems2.Math;
        import Modelica_LinearSystems2.Math.Complex;

        input Complex c "Complex number";
        input Modelica.SIunits.Angle phi0=0
          "phase angle phi shall be in the range: -pi < phi-phi0 < pi";
        output Modelica.SIunits.Angle phi "= phase angle of c";

      algorithm
        phi := Modelica.Math.atan3(
            c.im,
            c.re,
            phi0);
      end arg;

      encapsulated function conj "Conjugate of complex number"
        import Modelica_LinearSystems2.Math.Complex;

        input Complex c1 "Complex number";
        output Complex c2 "= c1.re - j*c1.im";
      algorithm
        c2 := Complex(c1.re, -c1.im);
      end conj;

      encapsulated function real "Real part of complex number"
        import Modelica_LinearSystems2.Math.Complex;

        input Complex c "Complex number";
        output Real r "= c.re ";
      algorithm
        r := c.re;
      end real;

      encapsulated function imag "imaginary part of complex number"
        import Modelica_LinearSystems2.Math.Complex;

        input Complex c "Complex number";
        output Real r "= c.im ";
      algorithm
        r := c.im;
      end imag;

      encapsulated function eigenValues
        "compute eingenvalues of a matrix A, using lapack routine dgeevx"

        import Modelica_LinearSystems2.Math.Complex;
        import Modelica_LinearSystems2.Math.Matrices.LAPACK;
        import Modelica_LinearSystems2.Math.Matrices.LAPACK.dgeevx;

        input Real A[:,size(A, 1)] "Real matrix";
        output Complex eigval[size(A, 1)]
          "Finite, invariant zeros of ss; size(Zeros,1) <= size(ss.A,1)";

      protected
        Integer nx=size(A, 1) "Number of states";

        Real alphaReal[nx];
        Real alphaImag[nx];

        Integer info;

      algorithm
       // Compute eigenvalues

        (alphaReal,alphaImag,,,,info) := dgeevx(A);
        assert(info == 0,
          "Failed to compute eigenvalues with function eigenValues_dgeevx(..)");

        for i in 1:nx loop
          eigval[i].re := alphaReal[i];
          eigval[i].im := alphaImag[i];
        end for;

      end eigenValues;

    encapsulated function eigenVectors
        "Calculate the rigth eigenvectors of a linear state space system and write them columnwise in a matrix."
        import Modelica;
        import Modelica.Math.Matrices.LAPACK;
        import Modelica_LinearSystems2.Math.Complex;
        import Re = Modelica_LinearSystems2.Math.Complex.real;
        import Im = Modelica_LinearSystems2.Math.Complex.imag;

      input Real A[:,size(A, 1)] "real square matrix";
      output Complex eigvec[size(A, 1),size(A, 2)] "eigen values of the system";
      output Complex eigval[size(A, 1)]=fill(Complex(0), size(A, 1))
          "eigen values of the system";
      protected
      Integer info;
      Real eigvecRe[size(A, 1),size(A, 2)];
      Real eigvalRe[size(A, 1)]=fill(0, size(A, 1));
      Real eigvalIm[size(A, 1)]=fill(0, size(A, 1));
      Integer n=size(A, 1);
      Integer i;
      Complex j=Complex.j();
    algorithm
      if size(A, 1) > 0 then

        (eigvalRe,eigvalIm,eigvecRe,info) := LAPACK.dgeev(A);
        for i in 1:size(A, 1) loop
          eigval[i].re := eigvalRe[i];
          eigval[i].im := eigvalIm[i];
        end for;

        assert(info == 0, "Calculating the eigen values with function
    \"StateSpace.Analysis.eigenVectors\" is not possible, since the
    numerical algorithm does not converge.");

        i := 1;
        while i < n loop
          if abs(eigvalIm[i]) > 0 then
            for ii in 1:n loop
              eigvec[ii, i] := eigvecRe[ii, i] + j*eigvecRe[ii, i + 1];
              eigvec[ii, i + 1] := eigvecRe[ii, i] - j*eigvecRe[ii, i + 1];
            end for;
            i := i + 2;
          else
            for ii in 1:n loop
              eigvec[ii, i] := Complex(1)*eigvecRe[ii, i];
            end for;
            i := i + 1;
          end if;
        end while;

      end if;

    end eigenVectors;

      encapsulated function frequency
        "Frequency and damping of conjugated complex pole pair"
        import Modelica;
        import Modelica_LinearSystems2.Math;
        import Modelica_LinearSystems2.Math.Complex;
        input Complex c "Complex number";
        output Modelica.SIunits.Frequency f "Frequency of c (= c.im in Hz)";
        output Real damping "Damping of c (= c.re/c.im)";

      protected
        Real abs_ev=(c.re^2 + c.im^2)^0.5;
      algorithm
        f := if abs(c.im) > 10*Modelica.Constants.eps then abs_ev/(2*Modelica.Constants.pi) else 0;
        damping := if abs(c.im) > 10*Modelica.Constants.eps then if abs_ev > Modelica.Constants.eps then -c.re/abs_ev else 0.0 else
          1.0;
      end frequency;

    encapsulated package Internal
        import Modelica;
     function eigenValues_dhseqr
          "compute eingenvalues of a upper Hessenberg matrix using lapack routine DHSEQR"

          import Modelica;
          import Modelica_LinearSystems2.Math.Complex;
          import Modelica_LinearSystems2.Math.Matrices.Internal.eigenvaluesHessenberg;

      input Real H[:,size(H, 1)] "Real upper Hessenberg matrix";
      output Complex zeros[size(H, 1)]
            "Finite, invariant zeros of ss; size(Zeros,1) <= size(ss.A,1)";

        protected
      Integer nx=size(H, 1) "Number of states";
      Real alphaReal[nx];
      Real alphaImag[nx];
      Integer info;

     algorithm
      (alphaReal,alphaImag,info) := eigenvaluesHessenberg(H);
      assert(info == 0,
        "Failed to compute eigenvalues with function Internal.eigenValues_dhseqr(..)");

      for i in 1:nx loop
        zeros[i].re := alphaReal[i];
        zeros[i].im := alphaImag[i];
      end for;

     end eigenValues_dhseqr;

      function C_transpose
          "Computes the transposed matrix of a complex matrix"
        extends Modelica.Icons.Function;

          import Modelica_LinearSystems2.Math.Complex;
          import Re = Modelica_LinearSystems2.Math.Complex.real;
          import Im = Modelica_LinearSystems2.Math.Complex.imag;

        input Complex C[:,:];
        output Complex CT[size(C, 2),size(C, 1)];
        protected
        Complex j=Complex.j();
        Integer l1;
        Integer l2;
      algorithm
      //  CT := Complex(1)*transpose(C[:,:].re)-j*transpose(C[:,:].im);// too slow
        for l1 in 1:size(C, 1) loop
          for l2 in 1:size(C, 2) loop
      //      CT[l2, l1] := Re(C[l1, l2]) - j*Im(C[l1, l2]);
            CT[l2, l1] := Complex.conj(C[l1,l2]);
          end for;
        end for;

      end C_transpose;

    function frobeniusNorm "Return the Frobenius norm of a matrix"
      extends Modelica.Icons.Function;
          import Modelica_LinearSystems2.Math.Complex;
      input Complex A[:,:] "Input matrix";
      output Real result=0.0 "frobenius norm of matrix A";
    algorithm
      for i1 in 1:size(A, 1) loop
        for i2 in 1:size(A, 2) loop
          result := result + Complex.real(A[i1, i2]*Complex.conj(A[i1, i2]));
        end for;
      end for;
      result := sqrt(result);
    end frobeniusNorm;

    end Internal;

    end Complex;

## Modelica_LinearSystems2.Math.Polynomial

Record defining the data for a polynomial

### Information

This record defines a polynomial, e.g., y = 2*x^2 + 3*x + 1. The general form is:

       y = c[1]*x^n + c[2]*x^(n-1) + ... + c[n]*x + c[n+1];

In the record, the coefficients c[i] are stored. Usually, the record is not directly accessed. Instead, a polynomial is generated with the functions provided in the record such as Polynomial.base(..), Polynomial.fitting(..). Several functions are provided that operate on polynomials.

Extends from Modelica.Icons.Record (Icon for a record).

### Parameters

Type | Name | Default | Description
--- | --- | --- | ---
Real | c[:] |  | Polynomial coefficients (c[1]*x^n + ... c[n]*x + c[n+1])

### Modelica definition

    record Polynomial "Record defining the data for a polynomial"
      extends Modelica.Icons.Record;
      Real c[:] "Polynomial coefficients (c[1]*x^n + ... c[n]*x + c[n+1])";

      encapsulated package Examples
        "Examples demonstrating the usage of Polynomials"
        import Modelica;

        function polynomialOperations1
          "Demonstrate basic algebraic operations on polynomials"
          import Modelica.Utilities.Streams.print;
          import Modelica_LinearSystems2.Math.Polynomial;
          import Modelica_LinearSystems2.Math.Polynomial.*;

          output Boolean ok;

        protected
          Polynomial t=Polynomial.x();
          Polynomial p;
        algorithm
          p := t^2 + 8*t + 16;

          print("y1 = " + String(p, name="t"));
          print("y2 = " + String(p*p, name="t"));
          ok := true;
        end polynomialOperations1;

        function polynomialOperations2
          "Demonstrate differentiation and integration of polynomials"
          import Modelica;
          import Modelica_LinearSystems2.Math.Polynomial;
          import Modelica_LinearSystems2.Math.Polynomial.*;
          import Modelica_LinearSystems2.TransferFunction;
          import Modelica.Utilities.Streams.print;

          output Boolean ok;
        protected
          Polynomial p=Polynomial({6,4,-3}) "6*x^2 + 4*x - 3";
          Polynomial int_p;
          Polynomial der_p;
          TransferFunction tf;
          Real x;
          Real int_val1;
          Real int_val2;
          Real der_val1;
          Real der_val2;
          Real der_val3;
          Real der_val4;
        algorithm
          print("Show integration and differentation of y = " + String(p));
          int_p := integral(p);
          der_p := derivative(p);

          print("  Integral of polynomial  : " + String(int_p));
          print("  Derivative of polynomial: " + String(der_p));

          print(
            "Compute derivatives and integral directly and via polynomials above (for p(x=2)):");
          x := 2;
          int_val1 := integralValue(p, x);
          int_val2 := evaluate(int_p, x);
          print("      integralValue(p,x) = " + String(int_val1) +
            ", evaluate(integral(p),x) = " + String(int_val2));

          der_val1 := derivativeValue(p, x);
          der_val2 := evaluate(der_p, x);
          print("    derivativeValue(p,x) = " + String(der_val1) +
            ", evaluate(derivative(p),x) = " + String(der_val2));
          der_val3 := derivativeValue(
                p,
                x,
                2);
          der_val4 := evaluate(derivative(der_p), x);
          print("  derivativeValue(p,x,2) = " + String(der_val3) +
            ", evaluate(derivative(derivative(p)),x) = " + String(der_val4));
          tf := der_p/int_p;
          print("  derivative/integral of polynomial = " + String(tf, name="x"));
          ok := true;
        end polynomialOperations2;

        function PascalTriangle "Generate and print Pascals triangle"
          import Modelica_LinearSystems2.Math.Polynomial;
          import Modelica_LinearSystems2.Math.Polynomial.*;
          import Modelica.Utilities.Streams.print;

          output Boolean ok;
        protected
          Polynomial p=Polynomial.x();
          Polynomial p1;
        algorithm
          p1 := p + 1;
          print("Pascals triangle up to order 20:");
          for i in 0:20 loop
            print("  " + String(p1^i));
          end for;
          ok := true;
        end PascalTriangle;

        function plotPolynomial
          "Demonstrate plotting of polynomial with automatic range and legend determination"
          import Modelica_LinearSystems2;
          import Modelica_LinearSystems2.Math;
          import Modelica_LinearSystems2.Math.Polynomial;
          import Modelica_LinearSystems2.Math.Complex;

          output Boolean ok;

        protected
          Complex j = Modelica_LinearSystems2.Math.Complex.j();
          Complex c[:]={-2+0*j,Complex(0),2+0*j,7+j,7-j};
          Polynomial p=Polynomial(c);
        algorithm
          Polynomial.plot(p);
          ok := true;
        end plotPolynomial;
      end Examples;

      encapsulated operator 'constructor'
        function fromVector
          "Constructor for a polynomial by providing a vector of coefficients"
          import Modelica;
          import Modelica_LinearSystems2.Math.Polynomial;

          input Real c[:] "Coefficients of polynomial in descending order";
          output Polynomial p(redeclare Real c[size(c, 1)]) "Polynomial";
        algorithm
          p.c := c;
        end fromVector;

        function fromReal
          "Generate a Polynomial data record from a Real value"
          import Modelica;
          import Modelica_LinearSystems2.Math.Polynomial;

          input Real r "Value of Real variable";
          output Polynomial p(redeclare Real c[1]);
        algorithm
          p.c := {r};
        end fromReal;

        function fromZeros "Generate a polynomial from given zeros"
          import Modelica;
          import Modelica_LinearSystems2.Math.Polynomial;
          import Modelica.Utilities.Streams;
          import Modelica_LinearSystems2.Math.Complex;

          input Complex roots[:]
            "Zeros of polynomial (must be real or conjugate complex pairs)";
          output Polynomial p(redeclare Real c[size(roots, 1) + 1])
            "Polynomial that corresponds to the zeros";

        protected
          Integer nr=size(roots, 1) "Number of roots";
          Integer nc=nr + 1 "Length of coefficent vector";
          Integer nn;
          Integer i;
        algorithm
          p.c[nc] := 1;
          i := 1;
          nn := 1;
          while i <= nr loop
             // Determine whether zero is real or conjugate complex
            if roots[i].im == 0.0 then
                // real root
              p.c := Polynomial.Internal.mult(p.c, nn, {1,-roots[i].re}, nc);

              i := i + 1;
              nn := nn + 1;
            else
                // complex root; check that it is a conjugate complex pair
              assert(i < nr, "Roots do not define a real valued polynomial\n" +
                "(roots[" + String(nr) + "] is complex without complex conjugate)");
              assert(roots[i].re == roots[i + 1].re,
                "No conjugate complex zero pair\n" + "  roots[" + String(i) +
                "].re = " + String(roots[i].re) + "\n" + "  roots[" + String(i + 1)
                 + "].re = " + String(roots[i].re) + "\n" +
                "and the two values should be identical since conjugate complex pair required.");
              assert(roots[i].im == -roots[i + 1].im,
                "No conjugate complex zero pair\n" + "  roots[" + String(i) +
                "].im = " + String(roots[i].im) + "\n" + "  roots[" + String(i + 1)
                 + "].im = " + String(roots[i + 1].im) + "\n" +
                "and the two values should be identical with opposite sign, since conjugate complex pair required.");
              p.c := Polynomial.Internal.mult(
                    p.c,
                    nn,
                    {1,-2*roots[i].re,roots[i].re^2 + roots[i].im^2},
                    nc);
              i := i + 2;
              nn := nn + 2;
            end if;
          end while;
        end fromZeros;

      end 'constructor';

      encapsulated operator '-'
        encapsulated function negate
          "Unary minus (multiply polynomial by -1)"

          import Modelica_LinearSystems2.Math.Polynomial;
          input Polynomial p;
          output Polynomial result(redeclare Real c[size(p.c, 1)]) "= -p";
        algorithm
          result.c := -p.c;
        end negate;

        encapsulated function subtract
          "Subtract two polynomials (p1 - p2)"
          import Modelica_LinearSystems2.Math.Polynomial;

          input Polynomial p1;
          input Polynomial p2;
          output Polynomial result(redeclare Real c[max(size(p1.c, 1), size(p2.c, 1))])
            "= p1 - p2";
        algorithm
        // Auxiliary variables not used, to enforce function inlining
          result.c := cat(
                1,
                zeros(max(size(p1.c, 1), size(p2.c, 1)) - size(p1.c, 1)),
                p1.c) - cat(
                1,
                zeros(max(size(p1.c, 1), size(p2.c, 1)) - size(p2.c, 1)),
                p2.c);
        end subtract;
      end '-';

      encapsulated operator function '+' "Add two polynomials (p1 + p2)"
        import Modelica_LinearSystems2.Math.Polynomial;

         input Polynomial p1;
         input Polynomial p2;
         output Polynomial result(redeclare Real c[max(size(p1.c, 1), size(p2.c, 1))])
          "= p1 + p2";
      algorithm
        // Auxiliary variables not used, to enforce function inlining
          result.c := cat(
                1,
                zeros(max(size(p1.c, 1), size(p2.c, 1)) - size(p1.c, 1)),
                p1.c) + cat(
                1,
                zeros(max(size(p1.c, 1), size(p2.c, 1)) - size(p2.c, 1)),
                p2.c);
      end '+';

      encapsulated operator function '*'
        "Multiply two polynomials (p1 * p2)"
        import Modelica_LinearSystems2.Math.Polynomial;

          input Polynomial p1;
          input Polynomial p2;
          output Polynomial result(redeclare Real c[size(p1.c, 1) + size(p2.c, 1)
               - 1]) "= p1 * p2";
      protected
          Integer n1=size(p1.c, 1);
          Integer n2=size(p2.c, 1);
          Integer n3=n1 + n2 - 1;
          Real ck;
      algorithm
          for k in 1:n3 loop
            ck := 0.0;
            for j in max(1, k + 1 - n2):min(k, n1) loop
              ck := ck + p1.c[j]*p2.c[k + 1 - j];
            end for;
            result.c[k] := ck;
          end for;
      end '*';

      encapsulated operator function '/'
        import Modelica_LinearSystems2.Math.Polynomial;
        import Modelica_LinearSystems2.TransferFunction;

          input Polynomial p1;
          input Polynomial p2;

          output TransferFunction tf(n=p1.c, d=p2.c);//only for tfpoly
      algorithm

          assert(size(p2.c, 1) > 0,
            "Denominator polynomial p2 must have at least one element, however\n"
             + "denominator is an empty polynomial. This is not allowed for p1/p2.");

      end '/';

      encapsulated operator function '^'
        import Modelica_LinearSystems2.Math.Polynomial;
        import Modelica.Utilities.Streams.print;

          input Polynomial p;
          input Integer n(min=0) = 1 "p^n shall be computed";
          output Polynomial result(redeclare Real c[max((size(p.c, 1) - 1)*n + 1, 1)])
          "= p^n";
      protected
          Integer n_p=size(p.c, 1);
          Integer n_power_p=max((n_p - 1)*n + 1, 1);
      algorithm
          if n == 0 then
            result.c := {1};
          else
            result.c[n_power_p - n_p + 1:n_power_p] := p.c;
            for i in 2:n loop
              result.c := Polynomial.Internal.mult(
                    result.c,
                    (n_p - 1)*(i - 1) + 1,
                    p.c,
                    n_power_p);
            end for;
          end if;
      end '^';

      encapsulated operator function '=='
        "Check whether two polynomials are identical"
        import Modelica_LinearSystems2.Math.Polynomial;

          input Polynomial p1;
          input Polynomial p2;
          input Real eps(min=0) = 0
          "Two coefficients c1 and c2 of the two polynomials are identical if abs(c1-c2) <= eps";
          output Boolean same "=true, if identical";
      protected
          Integer n1=size(p1.c, 1);
          Integer n2=size(p2.c, 1);
      algorithm
          if n1 == n2 then
            same := true;
            for i in 1:n1 loop
              if abs(p1.c[i] - p2.c[i]) > eps then
                same := false;
              end if;
            end for;
          else
            same := false;
          end if;
      end '==';

      encapsulated operator function 'String'
        "Transform Polynomial into a String representation"
        import Modelica_LinearSystems2.Math.Polynomial;

          input Polynomial p
          "Polynomial to be transformed in a String representation";
          input Integer significantDigits=6
          "Number of significant digits that are shown";
          input String name="x" "Independent variable name used for printing";
          output String s="";
      protected
          Boolean outputCoefficient;
          Integer power;
          Real ci;
          String v;
          Integer n=size(p.c, 1);
          Boolean printZero=true;
      algorithm
          if n == 0 then
            s := "0";
          else
            for i in 1:n loop
              if p.c[i] <> 0 and printZero then
                printZero := false;
              end if;

              if p.c[i] <> 0 or printZero then
                power := n - i;
                ci := p.c[i];
                if i > 1 then
                  if ci > 0 then
                    s := s + " + ";
                  elseif ci < 0 then
                    s := s + " - ";
                    ci := abs(ci);
                  end if;
                end if;

                outputCoefficient := power == 0 or ci <> 1;
                if outputCoefficient then
                  s := s + String(ci, significantDigits=significantDigits);
                end if;
                if outputCoefficient and power >= 1 then
                  s := s + "*";
                end if;
                if name == "" then
                  v := "?";
                else
                  v := name;
                end if;
                if power >= 2 then
                  s := s + v + "^" + String(power);
                elseif power == 1 then
                  s := s + v;
                end if;
              end if;
            end for;
          end if;
      end 'String';

      encapsulated function x "Generate a base polynomial y=x"
        import Modelica_LinearSystems2.Math.Polynomial;

        output Polynomial p(redeclare Real c[2]) "Polynomial";
      algorithm
        p.c := {1,0};
      end x;

      encapsulated function fitting
        "Computes a Polynomial that fits a set of data points in a least-squares sense"
        import Modelica;
        import Modelica_LinearSystems2.Math.Polynomial;

        input Real x[:] "Abscissa data values";
        input Real y[size(x, 1)] "Ordinate data values";
        input Integer order(min=1)
          "Order of desired polynomial that fits the data points (x,y)";
        output Polynomial p(redeclare Real c[order + 1])
          "Polynomial that fits the date points in a least squares sense";

      protected
        Real V[size(x, 1),order + 1] "Vandermonde matrix";
      algorithm
          // Construct Vandermonde matrix
        V[:, order + 1] := ones(size(x, 1));
        for j in order:-1:1 loop
          V[:, j] := {x[i]*V[i, j + 1] for i in 1:size(x, 1)};
        end for;

          // Solve least squares problem
        p.c := Modelica.Math.Matrices.leastSquares(V, y);
      end fitting;

      encapsulated function degree "Return degree of polynomial"
        import Modelica_LinearSystems2.Math.Polynomial;

        input Polynomial p;
        output Integer result "Degree of polynomial p";

      algorithm
        result := size(p.c, 1) - 1;
      end degree;

      encapsulated function degree2 "Return degree of polynomial"
        import Modelica_LinearSystems2.Math.Polynomial;

        input Polynomial p;
        output Integer result "Degree of polynomial p";
      protected
        Integer s;
      algorithm
        s := size(p.c, 1);
        for i in 1:s loop
                          // added correct code for degree calculation
          if p.c[i] <> 0 then
            result := s - i;
            break;
          end if;
        end for;
      end degree2;

      encapsulated function plot "Plot polynomial y=p(x)"

        import Modelica;
        import Modelica_LinearSystems2.Math.Polynomial;
        import Modelica.Utilities.Strings;
        import Modelica_LinearSystems2.Math.Complex;

        input Polynomial p "Polynomial to be plotted";
        input Integer nPoints(min=2) = 200 "Number of points";
        input Boolean autoLabel=true "= true, if automatically selected labels";
        input String xLabel="" "Abszissa description, if autoLabel = false";
        input String yLabel="" "Ordinate description, if autoLabel = false";
        input Boolean autoRange=true
          "= true, if abszissa range is automatically determined";
        input Real x_min=-1.0 "Minimum abszissa value, if autoRange = false";
        input Real x_max=1.0 "Maximum abszissa value, if autoRange = false";
      protected
        Real x_min2;
        Real x_max2;
        Real dx;
        Real x[nPoints];
        Real y[nPoints];
        Boolean OK;
        String yLabel2;
        Complex p_zeros[:];
        Complex pd_zeros[:];
        Complex points[:];
        String argument="x";
      algorithm
      /* Determine suitable x_min and x_max:
            Plotted range should contain roots and extrema
         */

        if autoRange then

          x_min2 := -1;
          x_max2 := +1;
          if size(p.c, 1) > 1 then
            p_zeros := Polynomial.roots(p);
            pd_zeros := Polynomial.roots(Polynomial.derivative(p));
            points := cat(
                1,
                p_zeros[:],
                pd_zeros[:]);

            x_min2 := points[1].re;
            x_max2 := points[1].re;
            for i in 2:size(points, 1) loop
              x_min2 := min(x_min2, points[i].re);
              x_max2 := max(x_max2, points[i].re);
            end for;

            if x_max2 > x_min2 then
              x_min2 := x_min2 - 0.1*(x_max2 - x_min2);
              x_max2 := x_max2 + 0.1*(x_max2 - x_min2);
            else
              x_min2 := x_min2 - 1;
              x_max2 := x_max2 + 1;
            end if;
          end if;
        else
          x_min2 := x_min;
          x_max2 := x_max;
        end if;

        // Compute polynomial values
        dx := (x_max2 - x_min2)/(nPoints - 1);
        x := x_min2:dx:x_max2;
        y := Polynomial.evaluate(p, x);

        // Determine labels
        if autoLabel then
            yLabel2 := String(p);
           if Strings.length(yLabel2) >= 100 then
            yLabel2 := Strings.substring(
                yLabel2,
                1,
                96) + " ...";
          end if;
        else
          yLabel2 := yLabel;
        end if;

        // Plot Polynomial
        OK := plotArray(
            x,
            y,
            legend=yLabel2);

      equation

      end plot;

      encapsulated function derivative "Derivative of Polynomial"
        import Modelica_LinearSystems2.Math.Polynomial;

        input Polynomial p;
        output Polynomial der_p(redeclare Real c[size(p.c, 1) - 1])
          "Derivative of Polynomial p";
      protected
        Integer n=size(p.c, 1);
      algorithm
        for j in 1:n - 1 loop
          der_p.c[j] := p.c[j]*(n - j);
        end for;
      end derivative;

      encapsulated function integral "Indefinite integral of Polynomial"
        import Modelica_LinearSystems2.Math.Polynomial;

        input Polynomial p "Polynomial";
        output Polynomial integral_p(redeclare Real c[size(p.c, 1) + 1])
          "Indefinite integral of polynomial p";
      protected
        Integer n=size(p.c, 1) + 1 "Number of coefficients of integral";
      algorithm
        for j in 1:n - 1 loop
          integral_p.c[j] := p.c[j]/(n - j);
        end for;
        integral_p.c[n] := 0;
      end integral;

      encapsulated function evaluate
        "Evaluate a Polynomial at a given (Real) abszissa value"
      annotation(derivative=Polynomial.evaluate_der);
        import Modelica_LinearSystems2.Math.Polynomial;

        input Polynomial p "Polynomial to be evaluated";
        input Real x "Abszissa value";
        output Real y "Value of polynomial at x";

      protected
        Integer n=size(p.c, 1);
      algorithm
        y := p.c[1];
        for j in 2:n loop
          y := p.c[j] + x*y;
        end for;
      end evaluate;

      encapsulated function evaluateMatrix
        "Evaluate a Polynomial with a matrix argument"
        import Modelica_LinearSystems2.Math.Polynomial;

        input Polynomial p "Polynomial to be evaluated";
        input Real X[:,size(X, 1)] "Square matrix argument";
        output Real Y[size(X, 1),size(X, 2)] "Value of polynomial at X";

      protected
        Integer n=size(p.c, 1);
      algorithm
        Y := zeros(size(X, 1), size(X, 2));
        for j in 1:n loop
          Y := X*Y + diagonal(p.c[j]*ones(size(X, 1)));
        end for;
      end evaluateMatrix;

      encapsulated function evaluateComplex
        "Evaluate a Polynomial at a given (Complex) abszissa value"
        import Modelica_LinearSystems2.Math.Complex;
        import Modelica_LinearSystems2.Math.Polynomial;

        input Polynomial p "Polynomial to be evaluated";
        input Complex x "Complex abszissa value";
        output Complex y "Complex value of polynomial at x";
      protected
        Integer n=size(p.c, 1);
      algorithm
        // Horners scheme for complex numbers
        y := Complex(p.c[1], 0);
        for j in 2:n loop
          y := Complex(p.c[j], 0) + x*y;
        end for;
      end evaluateComplex;

      encapsulated function derivativeValue
        "Integral of polynomial p(x) from x_low to x_high"
        import Modelica_LinearSystems2.Math.Polynomial;
        import Modelica.Utilities.Streams.print;

        input Polynomial p "Polynomial";
        input Real x "Abszissa value";
        input Integer i(min=1) = 1
          "i-th derivative to be evaluated; i=1 is first derivative";
        output Real der_y "Value of i-th derivative of polynomial at x";
        // annotation(derivative(zeroDerivative=p, zeroDerivative=i)=Polynomial.Internal.derivativeValue_der);
      protected
        Integer n=size(p.c, 1);
      algorithm
        if i > n - 1 then
          der_y := 0.0;
        else
          der_y := p.c[1]*product(n - k for k in 1:i);
          for j in 2:n - i loop
            der_y := p.c[j]*product(n - j + 1 - k for k in 1:i) + x*der_y;
          end for;
        end if;
      end derivativeValue;

      encapsulated function integralValue
        "Integral of polynomial p(x) from x_low to x_high"
      annotation(derivative=Polynomial.integralValue_der);
        import Modelica_LinearSystems2.Math.Polynomial;

        input Polynomial p "Polynomial to be integrated";
        input Real x_high "High integrand value";
        input Real x_low=0 "Low integrand value, default 0";
        output Real integral=0.0 "Integral of polynomial p from x_low to x_high";

      protected
        Integer n=size(p.c, 1) "Order of integrated polynomial";
        Real y_low=0 "value at lower integrand";
      algorithm
        for j in 1:n loop
          integral := x_high*(p.c[j]/(n - j + 1) + integral);
          y_low := x_low*(p.c[j]/(n - j + 1) + y_low);
        end for;
        integral := integral - y_low;
      end integralValue;

      encapsulated function roots
        "Determine zeros of polynomial, i.e., points x with p(x)=0"
        import Modelica_LinearSystems2.Math.Polynomial;
        import Modelica_LinearSystems2.Math.Complex;

        input Polynomial p "Polynomial";
        input Boolean printRoots=false "= true, if roots shall be pretty printed";
        output Complex result[:]=fill(Complex(0, 0),
            Polynomial.numberOfRoots(p)) "Roots of polynomial";

      algorithm
        result := Polynomial.rootsOfNonZeroHighestCoefficientPolynomial(p,
          Polynomial.numberOfRoots(p));
        if printRoots then
          Complex.Vectors.print("", result);
        end if;
      end roots;

        encapsulated function numberOfRoots
        "Determine number of roots of polynomial"
        import Modelica_LinearSystems2.Math.Polynomial;
           input Polynomial p "Polynomial";
           output Integer result "Number of roots of p";
      protected
           Integer nc=size(p.c, 1);
           Integer i;
        algorithm
        if nc <= 1 then
             // no roots
          result := 0;
        else
             // Remove all leading zero coefficients
          i := 1;
          result := 0;
          while i <= nc - 1 loop
            if p.c[i] <> 0.0 then
              result := nc - i;
              i := nc;
            else
              i := i + 1;
            end if;
          end while;
        end if;
        end numberOfRoots;

        encapsulated function rootsOfNonZeroHighestCoefficientPolynomial
        "Determine zeros of polynomial where highest coefficient of polynomial is not zero"
        import Modelica.Math.Matrices;
        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.Math.Polynomial;
        import Modelica_LinearSystems2.Math.Complex;

          input Polynomial p "Polynomial";
          input Integer numberOfRoots "Number of roots of polynomial";
          output Complex result[:]=fill(Complex(0, 0), numberOfRoots)
          "Roots of polynomial";
      protected
          Integer nc=size(p.c, 1);
          Integer i_start=nc - numberOfRoots;
          Integer n=numberOfRoots;
          Real A[n,n] "Companion matrix";
          Real ev[n,2] "Eigen values";
                                                Complex j = Modelica_LinearSystems2.Math.Complex.j();
        algorithm
          assert(numberOfRoots >= 0 and numberOfRoots < nc,
            "Argument numberOfRoots (= " + String(numberOfRoots) +
            ") is not in the range\n" + "0 <= numberOfRoots <= " + String(nc - 1));
          assert(p.c[i_start] <> 0, "p.c[" + String(i_start) +
            "] = 0. Probably wrong argument numberOfRoots (=" + String(
            numberOfRoots) + ")");

          if numberOfRoots > 0 then
            // companion matrix
            A[1, :] := -p.c[i_start + 1:nc]/p.c[i_start];
            A[2:n, :] := [identity(n - 1),zeros(n - 1)];

            // roots are eigenvalues of companion matrix
            ev := Matrices.eigenValues(A);
            for i in 1:n loop
              result[i] := ev[i, 1]+ j*ev[i, 2];
            end for;
          end if;
        end rootsOfNonZeroHighestCoefficientPolynomial;

        encapsulated function evaluate_der
        "Evaluate derivative of polynomial at a given abszissa value"
        import Modelica_LinearSystems2.Math.Polynomial;

          input Polynomial p "Polynomial";
          input Real x "Abszissa value";
          input Real dx "Derivative of abszissa value, der(x)";
          output Real dy "Derivative value of polynomial at x";
      protected
          Integer n=size(p.c, 1);
        algorithm
          dy := p.c[1]*(n - 1);
          for j in 2:n - 1 loop
            dy := p.c[j]*(n - j) + x*dy;
          end for;
          dy := dy*dx;
        end evaluate_der;

      encapsulated function integralValue_der
        "Evaluate derivative of integral of polynomial p(x) from x_low to x_high, assuming only x_high as time-dependent (Leibnitz rule)"
        import Modelica_LinearSystems2.Math.Polynomial;

        input Polynomial p "Polynomial";
        input Real x_high "High integrand value";
        input Real x_low=0 "Low integrand value, default 0";
        input Real dx_high "High integrand value";
        input Real dx_low=0 "Low integrand value, default 0";
        output Real dintegral=0.0 "Integral of polynomial p from u_low to u_high";
      algorithm
        dintegral := Polynomial.evaluate(p, x_high)*dx_high;
      end integralValue_der;

      encapsulated package Internal
        "Internal library of record Polynomial (should not be directly used by user)"
        import Modelica;
        import Modelica_LinearSystems2.Math.Polynomial;
        extends Modelica.Icons.Library;

      function mult
          "Multiply two polynomials (polynomials are defined by vectors)"
          import Modelica.Utilities.Streams.print;

        input Real p1[:];
        input Integer n1
            "Number of coefficients of p1 to be used, i.e., (end-n1+1:end)";
        input Real p2[:];
        input Integer n3_max "Dimension of output vector";
        output Real p3[n3_max];
        protected
        Integer n1_max = size(p1,1);
        Integer n2=size(p2,1);
        Integer n3=n1+n2-1;
        Real ck;
      algorithm
        for k in 1:n3 loop
           ck := 0.0;
           for j in max(1,k+1-n2):min(k,n1) loop
             ck := ck + p1[n1_max - n1 + j]*p2[k + 1 - j];
           end for;
           p3[n3_max-n3+k] := ck;
        end for;
      end mult;
      end Internal;

    end Polynomial;

## Modelica_LinearSystems2.Math.Complex.real

Real part of complex number

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Complex | c |  | Complex number

### Outputs

Type | Name | Description
--- | --- | ---
Real | r | = c.re

### Modelica definition

    encapsulated function real "Real part of complex number"
      import Modelica_LinearSystems2.Math.Complex;

      input Complex c "Complex number";
      output Real r "= c.re ";
    algorithm
      r := c.re;
    end real;

## Modelica_LinearSystems2.Math.Complex.imag

imaginary part of complex number

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Complex | c |  | Complex number

### Outputs

Type | Name | Description
--- | --- | ---
Real | r | = c.im

### Modelica definition

    encapsulated function imag "imaginary part of complex number"
      import Modelica_LinearSystems2.Math.Complex;

      input Complex c "Complex number";
      output Real r "= c.im ";
    algorithm
      r := c.im;
    end imag;

## Modelica_LinearSystems2.Math.Complex.exp

Exponential of complex number

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Complex | c1 |  | Complex number

### Outputs

Type | Name | Description
--- | --- | ---
Complex | c2 | = exp(c1)

### Modelica definition

    encapsulated function exp "Exponential of complex number"
      import Modelica_LinearSystems2.Math.Complex;
      import Modelica.Math;

      input Complex c1 "Complex number";
      output Complex c2 "= exp(c1)";
    algorithm
       c2 := Complex(Math.exp(c1.re)*Math.cos(c1.im), Math.exp(c1.re)*Math.sin(c1.im));
    end exp;

## Modelica_LinearSystems2.Math.Polynomial.evaluate_der

Evaluate derivative of polynomial at a given abszissa value

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Polynomial | p |  | Polynomial
Real | x |  | Abszissa value
Real | dx |  | Derivative of abszissa value, der(x)

### Outputs

Type | Name | Description
--- | --- | ---
Real | dy | Derivative value of polynomial at x

### Modelica definition

    encapsulated function evaluate_der
      "Evaluate derivative of polynomial at a given abszissa value"
    import Modelica_LinearSystems2.Math.Polynomial;

      input Polynomial p "Polynomial";
      input Real x "Abszissa value";
      input Real dx "Derivative of abszissa value, der(x)";
      output Real dy "Derivative value of polynomial at x";
    protected
      Integer n=size(p.c, 1);
    algorithm
      dy := p.c[1]*(n - 1);
      for j in 2:n - 1 loop
        dy := p.c[j]*(n - j) + x*dy;
      end for;
      dy := dy*dx;
    end evaluate_der;

## Modelica_LinearSystems2.Math.Polynomial.integralValue_der

Evaluate derivative of integral of polynomial p(x) from x_low to x_high, assuming only x_high as time-dependent (Leibnitz rule)

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Polynomial | p |  | Polynomial
Real | x_high |  | High integrand value
Real | x_low | 0 | Low integrand value, default 0
Real | dx_high |  | High integrand value
Real | dx_low | 0 | Low integrand value, default 0

### Outputs

Type | Name | Description
--- | --- | ---
Real | dintegral | Integral of polynomial p from u_low to u_high

### Modelica definition

    encapsulated function integralValue_der
      "Evaluate derivative of integral of polynomial p(x) from x_low to x_high, assuming only x_high as time-dependent (Leibnitz rule)"
      import Modelica_LinearSystems2.Math.Polynomial;

      input Polynomial p "Polynomial";
      input Real x_high "High integrand value";
      input Real x_low=0 "Low integrand value, default 0";
      input Real dx_high "High integrand value";
      input Real dx_low=0 "Low integrand value, default 0";
      output Real dintegral=0.0 "Integral of polynomial p from u_low to u_high";
    algorithm
      dintegral := Polynomial.evaluate(p, x_high)*dx_high;
    end integralValue_der;