## Modelica_LinearSystems2.Internal

Internal library of library LinearSystems (should not be directly used by user)

### Information

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

### Package Content

Name | Description
--- | ---
Streams |
Eigenvalue | Record containing a eigen value or a pair of conjugated complex pair, repectively and additionally characteristics of the eigenvalue(s)
PolesAndZeros | Record containing poles and zeros of a system in two real vectors containing the real parts and the imaginary parts respctively
StateSpaceR | Continuous state space description of a linear, time invariant differential equation system (data + operations)
AnalyseOptions | Defines the characteristics of the eigenvalues to be print or to be plot
AnalyseOptions2 | Defines the characteristics of the eigenvalues to be print or to be plot
DefaultDiagramPolesAndZeros | Default diagram for polesAndZeros plot
DefaultDiagramBodePlot | Default diagram for Bode plot
DefaultDiagramTimeResponse | Default diagram for a time response plot
PartialPlotFunction | Interface of a plot function
PartialPlotFunctionMIMO | Interface of a plot function for MIMO systems
defaultColumnLabels | Compute default column labels for plot
defaultImpulseResponseHeadings | Compute default headings for impulse response plot
defaultInitialResponseHeadings | Compute default headings for initial response plot
defaultRampResponseHeadings | Compute default headings for ramp response plot
defaultStepResponseHeadings | Compute default headings for step response plot
extractRealZeros | Extract real and conjugate complex elements from a complex vector
frequencyRangeBode | Determine min. and max. frequencies for Bode plot
frequencyRangeZeros | Determine min. and max. frequencies for a vector of zeros (numerator or denominator zeros)
frequencyVector | Determine frequency vector for Bode plot
numberOfRealZeros | Determine the number of elements of a Complex vector where the imaginary part is zero
printComplexVector |
reorderEigenvalue | Extract real and conjugate complex elements from a complex vector as part of a record Eigenvalue
reorderZeros | Extract real and conjugate complex elements from a complex vector
sortEigenvalue | Sort elements of Eigenvalue-record depending on the eigenvalues
timeResponseMask | Declares the common structure for the set of response functions
timeResponseMask2 | Declares the common structure for the set of response functions
timeResponseMask2_tf | Declares the common structure for the set of response functions
timeResponseMask2_zp | Declares the common structure for the set of response functions
timeResponseMask_tf | Declares the common structure for the set of response functions for a transfer function
timeResponseMask_zp | Declares the common structure for the set of response functions for a zeros and poles transfer function
timeResponseSamples | Estimate reasonable discretisation sample time and simulation time span for time response plot
from_Hz | Convert from Hz to rad/s
to_Hz | Convert from rad/s to Hz

## Modelica_LinearSystems2.Internal.Eigenvalue

Record containing a eigen value or a pair of conjugated complex pair, repectively and additionally characteristics of the eigenvalue(s)

### Information

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

### Modelica definition

    record Eigenvalue
      "Record containing a eigen value or a pair of conjugated complex pair, repectively and additionally characteristics of the eigenvalue(s)"
      import Modelica_LinearSystems2.Math.Complex;
      extends Modelica.Icons.Record;

      Complex ev;
      Boolean imag = false;
      Boolean isStable = false;
      Boolean isControllable = false;
      Boolean isStabilizable = false;
      Boolean isObservable = false;
      Boolean isDetectable = false;
      Real frequency;
      Real damping;
      Real timeConstant;

    encapsulated function constructor
        "Default constructor for eigenvalue"
        import Modelica;
        import Modelica_LinearSystems2.Math.Complex;
        import Modelica_LinearSystems2.Internal;
        import Modelica_LinearSystems2.Internal.Eigenvalue;

      input Complex ev_in=Complex(0);
      input Boolean isControllable=false;
      input Boolean isObservable = false;
      input Integer maxIndex1=0;
      input Integer maxIndex2=0;
      input Real Teps = 1e6
          "maximum time constant before regarded as infinity, i.e. a real eigenvalue is zero";

      output Eigenvalue ev;

      protected
      Boolean isStable = ev_in.re<0;
      Boolean isImag = abs(ev_in.im)>Modelica.Constants.eps;
      Real abs_ev = (ev_in.re^2+ev_in.im^2)^0.5;

    algorithm
      ev.ev := ev_in;
      ev.imag := isImag;
      ev.isStable := isStable;
      ev.isControllable := isControllable;
      ev.isStabilizable := isStable or isControllable;
      ev.isObservable := isObservable;
      ev.isDetectable := isStable or isObservable;
      ev.frequency := if ev.imag then abs_ev/(2*Modelica.Constants.pi) else 0.0;
      ev.damping := if ev.imag then if abs_ev>Modelica.Constants.eps then -ev.ev.re/abs_ev else 0.0 else 1.0;
      ev.timeConstant := if ev.imag then 0.0 else if abs(ev.ev.re) > 1/Teps then 1/abs(ev.ev.re) else Teps;
    end constructor;

    end Eigenvalue;

## Modelica_LinearSystems2.Internal.PolesAndZeros

Record containing poles and zeros of a system in two real vectors containing the real parts and the imaginary parts respctively

### Information

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

### Modelica definition

    record PolesAndZeros
      "Record containing poles and zeros of a system in two real vectors containing the real parts and the imaginary parts respctively"
      import Modelica_LinearSystems2.Math.Complex;
      extends Modelica.Icons.Record;

      Real p_real[:];
      Real p_im[:];
      Real z_real[:];
      Real z_im[:];

      Integer norz_p "number of real zeros of p";
      Integer norz_z "number of real zeros of z";
      encapsulated function constructor
        "Default constructor for poles and zeros"
        import Modelica;
        import Modelica_LinearSystems2.Math.Complex;
        import PolesAndZeros2 = Modelica_LinearSystems2.Internal.PolesAndZeros;
        import Modelica_LinearSystems2.Internal;

        input Complex p[:]={Complex(0, 0)} "Complex zeros";
        input Complex z[:]={Complex(0, 0)} "Complex zeros";
        output PolesAndZeros2 pz(
          redeclare Real p_real[size(p, 1)],
          redeclare Real p_im[size(p, 1)],
          redeclare Real z_real[size(z, 1)],
          redeclare Real z_im[size(z, 1)]);
      algorithm
        pz.p_real := p[:].re;
        pz.p_im := p[:].im;
        pz.z_real := z[:].re;
        pz.z_im := z[:].im;
        pz.norz_p := Internal.numberOfRealZeros(p);
        pz.norz_z := Internal.numberOfRealZeros(z);

      end constructor;

    end PolesAndZeros;

## Modelica_LinearSystems2.Internal.StateSpaceR

Continuous state space description of a linear, time invariant differential equation system (data + operations)

### Information

This record defines a linear time invariant differential equation system in state space form:

        der(x) = A * x + B * u
            y  = C * x + D * u

with

- u - the input vector
- y - the output vector
- x - the state vector
- A,B,C,D - matrices of appropriate dimensions

Compared to record StateSpace StateSpaceR cantains an additional Integer value r

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

### Parameters

Type | Name | Default | Description
--- | --- | --- | ---
der(x) = A*x + B*u; y = C*x + D*u
Real | A[:, :] |  |
Real | B[size(A, 1), :] |  |
Real | C[:, size(A, 1)] |  |
Real | D[size(C, 1), size(B, 2)] |  |

### Modelica definition

    record StateSpaceR
      "Continuous state space description of a linear, time invariant differential equation system (data + operations)"

      extends Modelica.Icons.Record;

      Real A[:,:];
      Real B[size(A, 1),:];
      Real C[:,size(A, 1)];
      Real D[size(C, 1),size(B, 2)];
      Integer r=0;

      encapsulated function constructor
        "Default constructor for a StateSpace record"
        import Modelica;
        import Modelica_LinearSystems2.Internal.StateSpaceR;

        input Real A[:,size(A, 1)] "der(x) = A*x + B*u";
        input Real B[size(A, 1),:] "der(x) = A*x + B*u";
        input Real C[:,size(A, 1)] "y = C*x + D*u";
        input Real D[size(C, 1),size(B, 2)] "y = C*x + D*u";
        input Integer r;
        output StateSpaceR result(
          redeclare Real A[size(A, 1),size(A, 2)],
          redeclare Real B[size(B, 1),size(B, 2)],
          redeclare Real C[size(C, 1),size(C, 2)],
          redeclare Real D[size(D, 1),size(D, 2)]);
      algorithm
        result.A := A;
        result.B := B;
        result.C := C;
        result.D := D;
        result.r := r;
      end constructor;
    end StateSpaceR;

## Modelica_LinearSystems2.Internal.AnalyseOptions

Defines the characteristics of the eigenvalues to be print or to be plot

### Information

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

### Parameters

Type | Name | Default | Description
--- | --- | --- | ---
analyse options
Boolean | plotEigenValues | true | plot eigenvalues
Boolean | plotInvariantZeros | true | plot invariant zeros
Boolean | plotStepResponse | true | plot step respones. Only for SISO system
Boolean | plotFrequencyResponse | true | plot bode diagram. Only for SISO system
Boolean | printEigenValues | true | write eigenvalues into the report
Boolean | printEigenValueProperties | true | write eigenvalues with properties
Boolean | printInvariantZeros | true | write invariant zreos into the report
Boolean | printControllability | true | indicates controllability of every single pole
Boolean | printObservability | true | indicates observability of every single pole

### Modelica definition

    record AnalyseOptions
      "Defines the characteristics of the eigenvalues to be print or to be plot"
      extends Modelica.Icons.Record;
      Boolean plotEigenValues=true "plot eigenvalues";
      Boolean plotInvariantZeros=true "plot invariant zeros";
      Boolean plotStepResponse=true "plot step respones. Only for SISO system";
      Boolean plotFrequencyResponse=true "plot bode diagram. Only for SISO system";
      Boolean printEigenValues=true "write eigenvalues into the report";
      Boolean printEigenValueProperties=true "write eigenvalues with properties";

      Boolean printInvariantZeros=true "write invariant zreos into the report";
      Boolean printControllability=true
        "indicates controllability of every single pole";
      Boolean printObservability=true
        "indicates observability of every single pole";
      String headingEigenValues="Eigenvalues";
      String headingInvariantzeros="Invariant zeros";
      String headingStepResponse="Step response";
      String headingFrequencyResponse="Frequency response";

    end AnalyseOptions;

## Modelica_LinearSystems2.Internal.AnalyseOptions2

Defines the characteristics of the eigenvalues to be print or to be plot

### Information

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

### Parameters

Type | Name | Default | Description
--- | --- | --- | ---
analyse options
Boolean | plotEigenValues | true | plot eigenvalues
Boolean | plotInvariantZeros | true | plot invariant zeros
Boolean | plotStepResponse | true | plot step respones. Only for SISO system
Boolean | plotFrequencyResponse | true | plot bode diagram. Only for SISO system
Boolean | printEigenValues | true | write eigenvalues into the report
Boolean | printInvariantZeros | true | write invariant zreos into the report
Boolean | printControllability | false | indicates controllability of every single pole
Boolean | printObservability | false | indicates observability of every single pole

### Modelica definition

    record AnalyseOptions2
      "Defines the characteristics of the eigenvalues to be print or to be plot"
      extends Modelica.Icons.Record;
      Boolean plotEigenValues=true "plot eigenvalues";
      Boolean plotInvariantZeros=true "plot invariant zeros";
      Boolean plotStepResponse=true "plot step respones. Only for SISO system";
      Boolean plotFrequencyResponse=true "plot bode diagram. Only for SISO system";
      Boolean printEigenValues=true "write eigenvalues into the report";

      Boolean printInvariantZeros=true "write invariant zreos into the report";
      Boolean printControllability=false
        "indicates controllability of every single pole";
      Boolean printObservability=false
        "indicates observability of every single pole";
      String headingEigenValues="Eigenvalues";
      String headingInvariantzeros="Invariant zeros";
      String headingStepResponse="Step response";
      String headingFrequencyResponse="Frequency response";

    final Boolean printEigenValueProperties = false
        "write eigenvalues with properties";

    end AnalyseOptions2;

## Modelica_LinearSystems2.Internal.DefaultDiagramPolesAndZeros

Default diagram for polesAndZeros plot

### Information

This record contains the default diagram options for pole/zero plots.

Extends from Modelica.Icons.Record (Icon for a record), Modelica_LinearSystems2.Utilities.Plot.Records.DefaultDiagram (Diagram with no curves (might be used as default diagram)).

### Parameters

Type | Name | Default | Description
--- | --- | --- | ---
Curve | curve[:] | fill(Curve(x=fill(0.0, 0), y... | Properties of the curves in one diagram of the figure
String | heading | "Poles (x) and invariant zer... | Heading displayed above diagram
Real | heightRatio | 0.6 | Height of diagram = heightRatio*diagramWidth
Boolean | grid | true | = true, if grid is shown
Axes
String | xLabel | "Real part" | String displayed at horizontal axis
String | yLabel | "Imaginary part" | String displayed at vertical axis
Boolean | logX | false | = true, if logarithmic scale of x-axis
Boolean | logY | false | = true, if logarithmic scale of y-axis
Boolean | uniformScaling | false | = true, if same vertical and horizontal axis increment
Legend
Boolean | legend | false | = true, if legend is shown
Boolean | legendFrame | false | = true, if frame around legend
Boolean | legendHorizontal | true | = true, if horizontal legend (provided it is meaningful)
LegendLocation | legendLocation | Modelica_LinearSystems2.Util... | Legend placement

### Modelica definition

    record DefaultDiagramPolesAndZeros
      "Default diagram for polesAndZeros plot"
      import Modelica;

       extends Modelica.Icons.Record;
       extends Modelica_LinearSystems2.Utilities.Plot.Records.DefaultDiagram(
                          heading="Poles (x) and invariant zeros (o)",
                          heightRatio=0.6,
                          legend=false,
                          xLabel="Real part",
                          yLabel="Imaginary part");
    end DefaultDiagramPolesAndZeros;

## Modelica_LinearSystems2.Internal.DefaultDiagramBodePlot

Default diagram for Bode plot

### Information

This record contains the default diagram options for magnitude in bode plots.

Extends from Modelica.Icons.Record (Icon for a record), Modelica_LinearSystems2.Utilities.Plot.Records.DefaultDiagram (Diagram with no curves (might be used as default diagram)).

### Parameters

Type | Name | Default | Description
--- | --- | --- | ---
Curve | curve[:] | fill(Curve(x=fill(0.0, 0), y... | Properties of the curves in one diagram of the figure
String | heading | "Bode-Diagram" | Heading displayed above diagram
Real | heightRatio | 0.4 | Height of diagram = heightRatio*diagramWidth
Boolean | grid | true | = true, if grid is shown
Axes
String | xLabel | "Frequency [Hz]" | String displayed at horizontal axis
String | yLabel | "" | String displayed at vertical axis
Boolean | logX | true | = true, if logarithmic scale of x-axis
Boolean | logY | true | = true, if logarithmic scale of y-axis
Boolean | uniformScaling | false | = true, if same vertical and horizontal axis increment
Legend
Boolean | legend | false | = true, if legend is shown
Boolean | legendFrame | false | = true, if frame around legend
Boolean | legendHorizontal | true | = true, if horizontal legend (provided it is meaningful)
LegendLocation | legendLocation | Modelica_LinearSystems2.Util... | Legend placement

### Modelica definition

    record DefaultDiagramBodePlot "Default diagram for Bode plot"
      import Modelica;

       extends Modelica.Icons.Record;
       extends Modelica_LinearSystems2.Utilities.Plot.Records.DefaultDiagram(
                          heading="Bode-Diagram",
                          heightRatio=0.4,
                          legend=false,
                          xLabel="Frequency [Hz]",
                          yLabel="",
                          logX=true,
                          logY=true);
    end DefaultDiagramBodePlot;

## Modelica_LinearSystems2.Internal.DefaultDiagramTimeResponse

Default diagram for a time response plot

### Information

This record contains the default diagram options for pole/zero plots.

Extends from Modelica.Icons.Record (Icon for a record), Modelica_LinearSystems2.Utilities.Plot.Records.DefaultDiagram (Diagram with no curves (might be used as default diagram)).

### Parameters

Type | Name | Default | Description
--- | --- | --- | ---
Curve | curve[:] | fill(Curve(x=fill(0.0, 0), y... | Properties of the curves in one diagram of the figure
String | heading | "time response" | Heading displayed above diagram
Real | heightRatio | 0.4 | Height of diagram = heightRatio*diagramWidth
Boolean | grid | true | = true, if grid is shown
Axes
String | xLabel | "time [s]" | String displayed at horizontal axis
String | yLabel | "system output" | String displayed at vertical axis
Boolean | logX | false | = true, if logarithmic scale of x-axis
Boolean | logY | false | = true, if logarithmic scale of y-axis
Boolean | uniformScaling | false | = true, if same vertical and horizontal axis increment
Legend
Boolean | legend | false | = true, if legend is shown
Boolean | legendFrame | false | = true, if frame around legend
Boolean | legendHorizontal | true | = true, if horizontal legend (provided it is meaningful)
LegendLocation | legendLocation | Modelica_LinearSystems2.Util... | Legend placement

### Modelica definition

    record DefaultDiagramTimeResponse
      "Default diagram for a time response plot"
      import Modelica;

       extends Modelica.Icons.Record;
       extends Modelica_LinearSystems2.Utilities.Plot.Records.DefaultDiagram(
                          heading="time response",
                          heightRatio=0.4,
                          legend=false,
                          xLabel="time [s]",
                          yLabel="system output");
    end DefaultDiagramTimeResponse;

## Modelica_LinearSystems2.Internal.PartialPlotFunction

Interface of a plot function

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Diagram | defaultDiagram |  | Default diagram layout
Device | device | Modelica_LinearSystems2.Util... | Properties of device where figure is shown

### Modelica definition

    partial function PartialPlotFunction "Interface of a plot function"

      import Modelica_LinearSystems2.Utilities.Plot;

      input Plot.Records.Diagram defaultDiagram "Default diagram layout";
      input Plot.Records.Device device=Modelica_LinearSystems2.Utilities.Plot.Records.Device()
        "Properties of device where figure is shown";
    end PartialPlotFunction;

## Modelica_LinearSystems2.Internal.PartialPlotFunctionMIMO

Interface of a plot function for MIMO systems

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Diagram | defaultDiagram |  | Default diagram layout
Device | device | Modelica_LinearSystems2.Util... | Properties of device where figure is shown

### Modelica definition

    partial function PartialPlotFunctionMIMO
      "Interface of a plot function for MIMO systems"

      import Modelica_LinearSystems2.Utilities.Plot;

      input Plot.Records.Diagram defaultDiagram "Default diagram layout";
      input Plot.Records.Device device=Modelica_LinearSystems2.Utilities.Plot.Records.Device()
        "Properties of device where figure is shown";
    end PartialPlotFunctionMIMO;

## Modelica_LinearSystems2.Internal.defaultColumnLabels

Compute default column labels for plot

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Integer | sizeC1 |  |

### Outputs

Type | Name | Description
--- | --- | ---
String | columnLabels[sizeC1 + 1] |

### Modelica definition

    function defaultColumnLabels "Compute default column labels for plot"
      import Modelica_LinearSystems2.StateSpace;

      input Integer sizeC1;
      output String columnLabels[sizeC1 + 1];

    algorithm
      columnLabels := cat(
        1,
        {"time [t]"},
        {"Out[" + String(i) + "]" for i in 1:sizeC1}) "Column labels";
    end defaultColumnLabels;

## Modelica_LinearSystems2.Internal.defaultImpulseResponseHeadings

Compute default headings for impulse response plot

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Integer | sizeB2 |  |
Integer | sizeC1 |  |

### Outputs

Type | Name | Description
--- | --- | ---
String | heading[sizeC1, sizeB2] |

### Modelica definition

    function defaultImpulseResponseHeadings
      "Compute default headings for impulse response plot"

      input Integer sizeB2;
      input Integer sizeC1;
      output String heading[sizeC1,sizeB2];
    algorithm
      heading[:, :] := {{"Impulse response[" + String(i) + "," + String(j) + "]"
        for j in 1:sizeB2} for i in 1:sizeC1};
    end defaultImpulseResponseHeadings;

## Modelica_LinearSystems2.Internal.defaultInitialResponseHeadings

Compute default headings for initial response plot

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Integer | sizeC1 |  |

### Outputs

Type | Name | Description
--- | --- | ---
String | heading[sizeC1, 1] |

### Modelica definition

    function defaultInitialResponseHeadings
      "Compute default headings for initial response plot"

      input Integer sizeC1;
      output String heading[sizeC1,1];
    algorithm
      heading[:, :] := matrix({"Initial response (Out " + String(i) + ")" for i in
        1:sizeC1});
    end defaultInitialResponseHeadings;

## Modelica_LinearSystems2.Internal.defaultRampResponseHeadings

Compute default headings for ramp response plot

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Integer | sizeB2 |  |
Integer | sizeC1 |  |

### Outputs

Type | Name | Description
--- | --- | ---
String | heading[sizeC1, sizeB2] |

### Modelica definition

    function defaultRampResponseHeadings
      "Compute default headings for ramp response plot"

      input Integer sizeB2;
      input Integer sizeC1;
      output String heading[sizeC1,sizeB2];
    algorithm
      heading[:, :] := {{"Ramp response[" + String(i) + "," + String(j) + "]" for j in
            1:sizeB2} for i in 1:sizeC1};
    end defaultRampResponseHeadings;

## Modelica_LinearSystems2.Internal.defaultStepResponseHeadings

Compute default headings for step response plot

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Integer | sizeB2 |  |
Integer | sizeC1 |  |

### Outputs

Type | Name | Description
--- | --- | ---
String | heading[sizeC1, sizeB2] |

### Modelica definition

    function defaultStepResponseHeadings
      "Compute default headings for step response plot"

      input Integer sizeB2;
      input Integer sizeC1;
      output String heading[sizeC1,sizeB2];
    algorithm
      heading[:, :] := {{"Step response[" + String(i) + "," + String(j) + "]" for j in
            1:sizeB2} for i in 1:sizeC1};
    end defaultStepResponseHeadings;

## Modelica_LinearSystems2.Internal.extractRealZeros

Extract real and conjugate complex elements from a complex vector

### Information

#### Syntax

>                     realZeros = Matrices.extractRealZeros(complexVector, numberOfRealRoots);
>     (realZeros, complexZeros) = Matrices.extractRealZeros(complexVector, numberOfRealRoots,
>                                                           name="complexVector");

#### Description

Function extractRealZeros(..) extracts the real zeros from the Complex vector "complexVector". It is required that all elements of complexVector define either a real zero (complexVector[i].im=0) or a conjugate complex zero pair (complexVector[i].re == complexVector[i+1].re and complexVector[i].im == -complexVector[i+1].im). The second argument "numberOfRealZeros" is determined by a function call of Internal.numberOfRealZeros(). The optional input argument "name" is used as name of "complexVector" in error messages.

The function returns the real elements of complexVector in vector "realZeros" and the real and imaginary part of a conjugate complex zero pair in matrix complexZeros[:]".

#### Example

>       // c = {0; 1+2j; 1-2j; 2; -3; -1-j; -1+j};
>         Integer n = numberOfRealZeros(c);
>         Real realZeros[n];
>         Real complexZeros[:] = fill(Complex(0), integer((size(c,1)-n)/2));
>       algorithm
>       (realZeros, complexZeros) := extractRealZeros(c, n);
>                -> realZeros    = {0, 2, (-3)};
>                   complexZeros = { 1+2j,
>                                   -1+j}

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Complex | complexVector[:] |  | Zeros of a polynomial with real coefficients, e.g., eigen values of a real matrix
Integer | numberOfRealZeros |  | Number of real zeros of ComplexVector determined with function numberOfRealZeros
String | name | "complexVector" | Name of complexVector to be used in error message

### Outputs

Type | Name | Description
--- | --- | ---
Real | realZeros[numberOfRealZeros] | Real zeros of complexVector
Complex | complexZeros[:] | Complex zeros without the corresponding conjugate complex pair element

### Modelica definition

    function extractRealZeros
      "Extract real and conjugate complex elements from a complex vector"
      import Modelica_LinearSystems2.Math.Complex;

      input Complex complexVector[:]
        "Zeros of a polynomial with real coefficients, e.g., eigen values of a real matrix";
      input Integer numberOfRealZeros
        "Number of real zeros of ComplexVector determined with function numberOfRealZeros";
      input String name="complexVector"
        "Name of complexVector to be used in error message";
      output Real realZeros[numberOfRealZeros] "Real zeros of complexVector";
      output Complex complexZeros[:]=fill(Complex(0), integer((size(complexVector, 1) -
          numberOfRealZeros)/2))
        "Complex zeros without the corresponding conjugate complex pair element";

    protected
      Integer n=size(complexVector, 1);
      Integer i;
      Integer i1;
      Integer i2;
    algorithm
      i1 := 0;
      i2 := 0;
      i := 1;
      while i <= n loop
        if complexVector[i].im == 0.0 then
          i1 := i1 + 1;
          realZeros[i1] := complexVector[i].re;
          i := i + 1;
        else
            // check that the next two zeros are a conjugate complex pair
          assert(i < n, "Argument " + name +
            " does not define a real valued polynomial\n" + name + "[" + String(n)
             + "] is complex without complex conjugate.");

          assert(complexVector[i].re == complexVector[i + 1].re,
            "No conjugate complex pair (checked real parts)\n" + "  " + name + "["
             + String(i) + "] = " + String(complexVector[i]) + "\n" + "  " + name
             + "[" + String(i + 1) + "] = " + String(complexVector[i]) + "\n" +
            "and the real parts of these two complex numbers should be identical\n"
             + "since conjugate complex pairs required.");
          assert(complexVector[i].im == -complexVector[i + 1].im,
            "No conjugate complex pair (checked the imaginary parts)\n" + "  " +
            name + "[" + String(i) + "] = " + String(complexVector[i]) + "\n" +
            "  " + name + "[" + String(i + 1) + "] = " + String(complexVector[i + 1])
             + "\n" +
            "and the imaginary parts of these two complex numbers should be identical\n"
             + "with opposite sign, since a conjugate complex pair is required.");
          i2 := i2 + 1;
            // Store the zero with the positive imaginary part
          if complexVector[i].im >= 0 then
            complexZeros[i2] := complexVector[i];
          else
            complexZeros[i2] := complexVector[i + 1];
          end if;
          i := i + 2;
        end if;
      end while;
    end extractRealZeros;

## Modelica_LinearSystems2.Internal.frequencyRangeBode

Determine min. and max. frequencies for Bode plot

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Complex | numZeros[:] |  | Zeros of numerator
Complex | denZeros[:] |  | Zeros of denominator

### Outputs

Type | Name | Description
--- | --- | ---
AngularVelocity | w_min | [rad/s]
AngularVelocity | w_max | [rad/s]

### Modelica definition

    function frequencyRangeBode
      "Determine min. and max. frequencies for Bode plot"
      import Modelica;
      import Modelica_LinearSystems2.Math.Complex;
      import Modelica_LinearSystems2.Internal;
      import Modelica_LinearSystems2.TransferFunction;
      import Modelica.Utilities.Streams;

      input Complex numZeros[:] "Zeros of numerator";
      input Complex denZeros[:] "Zeros of denominator";
      output Modelica.SIunits.AngularVelocity w_min;
      output Modelica.SIunits.AngularVelocity w_max;
    protected
      Real phi_min=Modelica.SIunits.Conversions.from_deg(3);
      Real real_min=1.0e-4;
      Real pi=Modelica.Constants.pi;
      Integer n_num;
      Integer n_den;
      Real w_min1;
      Real w_min2;
      Real w_max1;
      Real w_max2;
      Real f_min;
      Real f_max;
    algorithm
      // Compute frequencies for numerator
      n_num := size(numZeros, 1);
      if n_num > 0 then
        (w_min1,w_max1) := Internal.frequencyRangeZeros(
          numZeros,
          phi_min,
          real_min);
      end if;

      // Compute frequencies for denominator
      n_den := size(denZeros, 1);
      if n_den > 0 then
        (w_min2,w_max2) := Internal.frequencyRangeZeros(
          denZeros,
          phi_min,
          real_min);
      end if;

      // Use largest range
      if n_num == 0 and n_den == 0 then
        w_min := 0.1;
        w_max := 10;
      elseif n_num == 0 then
        w_min := w_min2;
        w_max := w_max2;
      elseif n_den == 0 then
        w_min := w_min1;
        w_max := w_max1;
      else
        w_min := min(w_min1, w_min2);
        w_max := max(w_max1, w_max2);
      end if;
      f_min := Internal.to_Hz(w_min);
      f_max := Internal.to_Hz(w_max);
    end frequencyRangeBode;

## Modelica_LinearSystems2.Internal.frequencyRangeZeros

Determine min. and max. frequencies for a vector of zeros (numerator or denominator zeros)

### Information

This function estimates a useful frequency range for the Bode plot of a vector of zeros (numerator or denominator zeros). This frequency range is estimated such that the phase angle of one zero is in the range:

       phi_min/n_zeros ≤ |phase angle| ≤ pi/2 - phi_min/n_zeros

where n_zeros is the number of zeros. Note, the phase angle of one zero for a frequency of 0 up to infinity is in the range:

       0 ≤ |phase angle| ≤ pi/2

Therefore, the frequency range is estimated such that the essential part of the phase angle (defined by phi_min) is present.

If the real part of a complex zero vanishes (i.e., the zero is located on the imaginary axis), the maximum value of the bode plot magnitude of the zero is infinity. In order to avoid difficulties, zeros close to the imaginary axis are shifted by the input argument real_min along the real axis.

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Complex | z[:] |  | Vector of zeros
Angle | phi_min | Modelica.SIunits.Conversions... | Minimum phase angle [rad]
Real | real_min | 1.e-4 | [r| < real_min are treated as |real_min|

### Outputs

Type | Name | Description
--- | --- | ---
AngularVelocity | w_min | Minimum frequency [rad/s]
AngularVelocity | w_max | Maximum frequency [rad/s]

### Modelica definition

    function frequencyRangeZeros
      "Determine min. and max. frequencies for a vector of zeros (numerator or denominator zeros)"
      import Modelica;
      import Modelica.Math;
      import SI = Modelica.SIunits;
      import Modelica_LinearSystems2.Math.Complex;

      input Complex z[:] "Vector of zeros";
      input SI.Angle phi_min(min=10*Modelica.Constants.eps)=
        Modelica.SIunits.Conversions.from_deg(5) "Minimum phase angle";
      input Real real_min(min=0) = 1.e-4 "[r| < real_min are treated as |real_min|";
      output Modelica.SIunits.AngularVelocity w_min "Minimum frequency";
      output Modelica.SIunits.AngularVelocity w_max "Maximum frequency";

    protected
      Integer nz=size(z, 1);
      Real tan_min;
      Real tan_max1;
      Real tan_max2;
      Real z_re;
      Real w_min1;
      Real w_max1;
      Real z_abs2;
      Real k1;
      Real k2;
    algorithm
      /* - Real zero:
           tan(phi_desired) = w/re -> w = re*tan(phi_desired)

         - Conjugate complex zero:
           tan(phi1) = (w + im)/re
           tan(phi2) = (w - im)/re
           phi1 + phi2 = phi_desired -> tan(phi1 + phi2) = (tan(phi1) + tan(phi2))/(1-tan(phi1)*tan(phi2))
           -> after a longer derivation it follows:
              w = sqrt( re^2/tan(phi_desired)^2 + (re^2 + im^2) ) - re/tan(phi_desired)
      */
      assert(nz > 0, "Vector z of zeros has dimension 0, This is not allowed");
      tan_min := Math.tan(phi_min);
      tan_max1 := Math.tan(Modelica.Constants.pi/2 - phi_min);
      tan_max2 := Math.tan(Modelica.Constants.pi - phi_min);
      for i in 1:size(z, 1) loop
        if z[i].im >= 0.0 then
          z_re := max(abs(z[i].re), real_min);
          if z[i].im > 0.0 then
            z_abs2 := z_re^2 + z[i].im^2;
            k1 := z_re/tan_min;
            w_min1 := sqrt(k1^2 + z_abs2) - k1;
            k2 := z_re/tan_max2;
            w_max1 := sqrt(k2^2 + z_abs2) - k2;
          else
            w_min1 := z_re*tan_min;
            w_max1 := z_re*tan_max1;
          end if;

          if i == 1 then
            w_min := w_min1;
            w_max := w_max1;
          else
            w_min := min(w_min, w_min1);
            w_max := max(w_max, w_max1);
          end if;
        end if;
      end for;
    end frequencyRangeZeros;

## Modelica_LinearSystems2.Internal.frequencyVector

Determine frequency vector for Bode plot

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Integer | nPoints | 200 | Number of points
Boolean | autoRange | true | = true, if abszissa range is automatically determined
Frequency | f_min | 0.1 | Minimum frequency value, if autoRange = false [Hz]
Frequency | f_max | 10 | Maximum frequency value, if autoRange = false [Hz]
Complex | numZeros[:] | fill(Complex(0), 0) | Zeros of numerator
Complex | denZeros[:] | fill(Complex(0), 0) | Zeros of denominator

### Outputs

Type | Name | Description
--- | --- | ---
Frequency | f[nPoints] | Frequency vector (automatic or manual) [Hz]

### Modelica definition

    function frequencyVector "Determine frequency vector for Bode plot"
      import Modelica;
      import Modelica_LinearSystems2.Math.Complex;
      import LinearSystems = Modelica_LinearSystems2;
      import Modelica_LinearSystems2.Internal;
      import SI = Modelica.SIunits;

      input Integer nPoints(min=2) = 200 "Number of points";
      input Boolean autoRange=true
        "= true, if abszissa range is automatically determined";
      input Modelica.SIunits.Frequency f_min(min=0) = 0.1
        "Minimum frequency value, if autoRange = false";
      input Modelica.SIunits.Frequency f_max(min=0) = 10
        "Maximum frequency value, if autoRange = false";
      input Complex numZeros[:]=fill(Complex(0), 0) "Zeros of numerator";
      input Complex denZeros[:]=fill(Complex(0), 0) "Zeros of denominator";
      output SI.Frequency f[nPoints] "Frequency vector (automatic or manual)";
    protected
      Real w_min;
      Real w_max;
      Real f_min2;
      Real f_max2;
      Real f_log[nPoints];
    algorithm
      // Determine f_min2, f_max2 (auto or manual)
      if autoRange then
        (w_min,w_max) := Internal.frequencyRangeBode(numZeros, denZeros);
        f_min2 := Internal.to_Hz(w_min);
        f_max2 := Internal.to_Hz(w_max);

      else
        f_min2 := f_min;
        f_max2 := f_max;
      end if;
      w_min := Internal.from_Hz(f_min2);
      w_max := Internal.from_Hz(f_max2);

      // Compute logarithmic vector of frequency points
      f_log := linspace(
        Modelica.Math.log10(f_min2),
        Modelica.Math.log10(f_max2),
        nPoints);
      for i in 1:nPoints loop
        f[i] := 10^f_log[i];
      end for;
    end frequencyVector;

## Modelica_LinearSystems2.Internal.numberOfRealZeros

Determine the number of elements of a Complex vector where the imaginary part is zero

### Information

#### Syntax

>          Matrices.numberOfRealZeros(complexVector);

#### Description

Function numberOfRealZeros(..) determines the number of elements of vector "complexVector" with vanishing imaginary part, i.e., complexVector[i].im = 0.

#### Example

>       // c = {0; 1+2*j; 1-2*j; 2; -3; -1-j; -1+j};
>       result = Matrices.numberOfRealZeros(c);
>                -> result = 3;

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Complex | complexVector[:] |  | Complex vector

### Outputs

Type | Name | Description
--- | --- | ---
Integer | result | Number of elements of v with v.im = 0

### Modelica definition

    function numberOfRealZeros
      "Determine the number of elements of a Complex vector where the imaginary part is zero"
      import Modelica_LinearSystems2.Math.Complex;

      input Complex complexVector[:] "Complex vector";
      output Integer result "Number of elements of v with v.im = 0";

    algorithm
      result := 0;
      for i in 1:size(complexVector, 1) loop
        if complexVector[i].im == 0 then
          result := result + 1;
        end if;
      end for;
    end numberOfRealZeros;

## Modelica_LinearSystems2.Internal.printComplexVector

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
String | name |  |
Complex | c[:] |  |

### Modelica definition

    function printComplexVector
      import Modelica.Utilities.Streams.print;
      import Modelica_LinearSystems2.Math.Complex;

      input String name;
      input Complex c[:];

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

## Modelica_LinearSystems2.Internal.reorderEigenvalue

Extract real and conjugate complex elements from a complex vector as part of a record Eigenvalue

### Information

#### Syntax

>                    reorderedZero = Matrices.reorderZeros(complexVector);
>     (reorderedZeros, nRealZeros) = Matrices.reorderZeros(complexVector,
>                                                          name="complexVector");

#### Description

Function reorderZeros(..) reorders the zeros from the Complex vector "complexVector" such that the returned Complex vector reorderedZeros contains first all real Zeros and afterwards the conjugate complex zero pairs. It is required that all elements of complexVector define either a real zero (complexVector[i].im=0) or a conjugate complex zero pair (complexVector[i].re == complexVector[i+1].re and complexVector[i].im == -complexVector[i+1].im). The optional input argument "name" is used as name of "complexVector" in error messages.

The function returns the vector element reordered, as well as the number of real zeros (nRealZeros).

#### Example

>       // c = {0; 1+2j; 1-2j; 2; -3; -1-j; -1+j};
>         Real complexZeros[:] = fill(Complex(0), integer((size(c,1)-n)/2));
>       algorithm
>       (reorderedZeros, nRealZeros) := reorderZeros(c);
>           -> reorderedZeros = {0, 2, (-3), 1+2j, 1-2j, -1+j, -1-j}
>              nRealZeros     = 3

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Eigenvalue | EigenvalueVector[:] |  | Zeros of a polynomial with real coefficients, e.g., eigen values of a real matrix
String | name | "EigenvalueVector" | Name of complexVector to be used in error message

### Outputs

Type | Name | Description
--- | --- | ---
Eigenvalue | reorderedEigenvalues[size(EigenvalueVector, 1)] | Reordered zeros
Integer | nRealEigenvalues | Number of real zeros (EigenvalueVector[1:nRealEigenvalues] are the real zeros)

### Modelica definition

    function reorderEigenvalue
      "Extract real and conjugate complex elements from a complex vector as part of a record Eigenvalue"
      import Modelica_LinearSystems2.Math.Complex;
      import Modelica_LinearSystems2.Internal.Eigenvalue;
      input Eigenvalue EigenvalueVector[:]
        "Zeros of a polynomial with real coefficients, e.g., eigen values of a real matrix";
      input String name="EigenvalueVector"
        "Name of complexVector to be used in error message";
      output Eigenvalue reorderedEigenvalues[size(EigenvalueVector, 1)]
        "Reordered zeros";
      output Integer nRealEigenvalues
        "Number of real zeros (EigenvalueVector[1:nRealEigenvalues] are the real zeros)";

    protected
      Integer n=size(EigenvalueVector, 1);
      Integer i;
      Integer jr;
      Integer jc;
      Complex complexVector[size(EigenvalueVector, 1)];
    algorithm
    nRealEigenvalues:=numberOfRealZeros(complexVector);

      i := 1;
      jr := 1;
      jc := nRealEigenvalues + 1;
      while i <= n loop
        if EigenvalueVector[i].ev.im == 0.0 then

          reorderedEigenvalues[jr] := EigenvalueVector[i];
          reorderedEigenvalues[jr].ev.im := 0.0;
          i := i + 1;
          jr := jr + 1;
        else
    // check that the next two zeros are a conjugate complex pair
          assert(i < n, "Argument " + name +
            " does not define a real valued polynomial\n" + name + "[" + String(n)
             + "] is complex without complex conjugate.");
          assert(abs(EigenvalueVector[i].ev.re - EigenvalueVector[i+1].ev.re) < max(Modelica.Constants.eps,
            abs(EigenvalueVector[i+1].ev.re)*100*Modelica.Constants.eps),
            "No conjugate complex pair (checked the real parts)\n" + "  " + name +
            "[" + String(i) + "] = " + String(EigenvalueVector[i].ev) + "\n" + "  " +
            name + "[" + String(i + 1) + "] = " + String(EigenvalueVector[i].ev) + "\n"
             +
            "and the real parts of these two complex numbers should be identical\n"
             + "since conjugate complex pairs required.");
           assert(abs(EigenvalueVector[i].ev.im + EigenvalueVector[i+1].ev.im) < max(Modelica.Constants.eps,
            abs(EigenvalueVector[i+1].ev.im)*100*Modelica.Constants.eps),
            "No conjugate complex pair (checked the imaginary parts)\n" + "  " +
            name + "[" + String(i) + "] = " + String(EigenvalueVector[i].ev) + "\n" +
            "  " + name + "[" + String(i + 1) + "] = " + String(EigenvalueVector[i + 1].ev)
             + "\n" +
            "and the imaginary parts of these two complex numbers should be identical\n"
             + "with opposite sign, since a conjugate complex pair is required.");

            // Store the zero with the positive imaginary part

          if EigenvalueVector[i].ev.im >= 0 then
            reorderedEigenvalues[jc] := EigenvalueVector[i];
            reorderedEigenvalues[jc + 1] := EigenvalueVector[i+1];
          else
            reorderedEigenvalues[jc] := EigenvalueVector[i+1];
            reorderedEigenvalues[jc + 1] := EigenvalueVector[i];
          end if;
          i := i + 2;
          jc := jc + 2;
        end if;
      end while;
    end reorderEigenvalue;

## Modelica_LinearSystems2.Internal.reorderZeros

Extract real and conjugate complex elements from a complex vector

### Information

#### Syntax

>                    reorderedZero = Internal.reorderZeros(complexVector);
>     (reorderedZeros, nRealZeros) = Internal.reorderZeros(complexVector,
>                                                          name="complexVector");

#### Description

Function reorderZeros(..) reorders the zeros from the Complex vector "complexVector" such that the returned Complex vector reorderedZeros contains first all real Zeros and afterwards the conjugate complex zero pairs. It is required that all elements of complexVector define either a real zero (complexVector[i].im=0) or a conjugate complex zero pair (complexVector[i].re == complexVector[i+1].re and complexVector[i].im == -complexVector[i+1].im). The optional input argument "name" is used as name of "complexVector" in error messages.

The function returns the vector element reordered, as well as the number of real zeros (nRealZeros).

#### Example

>       // c = {0; 1+2j; 1-2j; 2; -3; -1-j; -1+j};
>         Real complexZeros[:] = fill(Complex(0), integer((size(c,1)-n)/2));
>       algorithm
>       (reorderedZeros, nRealZeros) := reorderZeros(c);
>           -> reorderedZeros = {0, 2, (-3), 1+2j, 1-2j, -1+j, -1-j}
>              nRealZeros     = 3

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Complex | complexVector[:] |  | Zeros of a polynomial with real coefficients, e.g., eigen values of a real matrix
String | name | "complexVector" | Name of complexVector to be used in error message

### Outputs

Type | Name | Description
--- | --- | ---
Complex | reorderedZeros[size(complexVector, 1)] | Reordered zeros
Integer | nRealZeros | Number of real zeros (reorderedZeros[1:nRealZeros] are the real zeros)

### Modelica definition

    function reorderZeros
      "Extract real and conjugate complex elements from a complex vector"
      import Modelica_LinearSystems2.Math.Complex;

      input Complex complexVector[:]
        "Zeros of a polynomial with real coefficients, e.g., eigen values of a real matrix";
      input String name="complexVector"
        "Name of complexVector to be used in error message";
      output Complex reorderedZeros[size(complexVector, 1)] "Reordered zeros";
      output Integer nRealZeros=numberOfRealZeros(complexVector)
        "Number of real zeros (reorderedZeros[1:nRealZeros] are the real zeros)";

    protected
      Integer n=size(complexVector, 1);
      Integer i;
      Integer jr;
      Integer jc;
    algorithm
      i := 1;
      jr := 1;
      jc := nRealZeros + 1;
      while i <= n loop
        if abs(complexVector[i].im) < Modelica.Constants.eps then
          reorderedZeros[jr].re := complexVector[i].re;
          reorderedZeros[jr].im := 0.0;
          i := i + 1;
          jr := jr + 1;
        else
    // check that the next two zeros are a conjugate complex pair
          assert(i < n, "Argument " + name +
            " does not define a real valued polynomial\n" + name + "[" + String(n)
             + "] is complex without complex conjugate.");
          assert(abs(complexVector[i].re - complexVector[i + 1].re) < max(Modelica.Constants.eps,
            abs(complexVector[i + 1].re)*100*Modelica.Constants.eps),
            "No conjugate complex pair (checked the real parts)\n" + "  " + name +
            "[" + String(i) + "] = " + String(complexVector[i]) + "\n" + "  " +
            name + "[" + String(i + 1) + "] = " + String(complexVector[i]) + "\n"
             +
            "and the real parts of these two complex numbers should be identical\n"
             + "since conjugate complex pairs required.");
          assert(abs(complexVector[i].im + complexVector[i + 1].im) < max(Modelica.Constants.eps,
            abs(complexVector[i + 1].im)*100*Modelica.Constants.eps),
            "No conjugate complex pair (checked the imaginary parts)\n" + "  " +
            name + "[" + String(i) + "] = " + String(complexVector[i]) + "\n" +
            "  " + name + "[" + String(i + 1) + "] = " + String(complexVector[i + 1])
             + "\n" +
            "and the imaginary parts of these two complex numbers should be identical\n"
             + "with opposite sign, since a conjugate complex pair is required.");

            // Store the zero with the positive imaginary part
          if complexVector[i].im >= 0 then
            reorderedZeros[jc] := complexVector[i];
            reorderedZeros[jc + 1] := complexVector[i + 1];
          else
            reorderedZeros[jc] := complexVector[i + 1];
            reorderedZeros[jc + 1] := complexVector[i];
          end if;
          i := i + 2;
          jc := jc + 2;
        end if;
      end while;
    end reorderZeros;

## Modelica_LinearSystems2.Internal.sortEigenvalue

Sort elements of Eigenvalue-record depending on the eigenvalues

### Information

#### Syntax

>                sorted_v = Vectors.sort(v);
>     (sorted_v, indices) = Vectors.sort(v, ascending=true);

#### Description

Function sort(..) sorts a Real vector v in ascending order and returns the result in sorted_v. If the optional argument "ascending" is false, the vector is sorted in descending order. In the optional second output argument the indices of the sorted vector with respect to the original vector are given, such that sorted_v = v[indices].

#### Example

>       (v2, i2) := Vectors.sort({-1, 8, 3, 6, 2});
>            -> v2 = {-1, 2, 3, 6, 8}
>               i2 = {1, 5, 3, 4, 2}

Extends from Modelica.Icons.Function (Icon for a function).

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Eigenvalue | ev[:] |  | Vector to be sorted
Boolean | ascending | true | = true if ascending order, otherwise descending order
Boolean | sortFrequency | true | = true, if sorting is first for imaginary then for real value; = false, if sorting is for absolute value

### Outputs

Type | Name | Description
--- | --- | ---
Eigenvalue | sorted_ev[size(ev, 1)] | Sorted vector
Integer | indices[size(ev, 1)] | sorted_ev = ev[indices]

### Modelica definition

    function sortEigenvalue
      "Sort elements of Eigenvalue-record depending on the eigenvalues"
      extends Modelica.Icons.Function;

      import Modelica_LinearSystems2.Internal.Eigenvalue;
      import Modelica_LinearSystems2.Math.Complex;
      input Eigenvalue ev[:] "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 Eigenvalue sorted_ev[size(ev,1)] = ev "Sorted vector";
      output Integer indices[size(ev,1)] = 1:size(ev,1) "sorted_ev = ev[indices]";

      /* shellsort algorithm; should be improved later */
    protected
      Integer gap;
      Integer i;
      Integer j;
      Eigenvalue wev;
      Integer wi;
      Integer nev = size(ev,1);
      Boolean swap;
      Integer k1;
      Integer k2;
    algorithm
      gap := div(nev,2);

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

            while swap loop
               wev := sorted_ev[j+1];
               wi := indices[j+1];
               sorted_ev[j+1] := sorted_ev[j+gap+1];
               sorted_ev[j+gap+1] := wev;
               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_ev[k1].ev.im) >  abs(sorted_ev[k2].ev.im) or
                                abs(sorted_ev[k1].ev.im) == abs(sorted_ev[k2].ev.im) and
                                (sorted_ev[k1].ev.re  > sorted_ev[k2].ev.re or
                                 sorted_ev[k1].ev.re  == sorted_ev[k2].ev.re and sorted_ev[k1].ev.im < sorted_ev[k2].ev.im);
                     else
                        swap := abs(sorted_ev[k1].ev.im) <  abs(sorted_ev[k2].ev.im) or
                                abs(sorted_ev[k1].ev.im) == abs(sorted_ev[k2].ev.im) and
                                (sorted_ev[k1].ev.re  < sorted_ev[k2].ev.re or
                                 sorted_ev[k1].ev.re  == sorted_ev[k2].ev.re and sorted_ev[k1].ev.im < sorted_ev[k2].ev.im);
                     end if;
                  else
                     if ascending then
                        swap := Complex.'abs'(sorted_ev[k1].ev) > Complex.'abs'(sorted_ev[k2].ev);
                     else
                        swap := Complex.'abs'(sorted_ev[k1].ev) < Complex.'abs'(sorted_ev[k2].ev);
                     end if;
                  end if;
               else
                  swap := false;
               end if;
            end while;
            i := i + 1;
         end while;
         gap := div(gap,2);
      end while;
    end sortEigenvalue;

## Modelica_LinearSystems2.Internal.timeResponseMask

Declares the common structure for the set of response functions

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | sc |  |
Real | dt | 0 | Sample time [s]
Real | tSpan | 0 | Simulation time span [s]
Integer | toWindow[:, :] | fill(-1, size(sc.C, 1), size... | -1/0/>0 plot in new window/last window/window ID
Boolean | clearWindow[:, :] | fill(false, size(sc.C, 1), s... | = true, if previous window content is removed
String | heading[:, :] | fill("Time response", size(s... | Heading of the response diagram
String | columnLabels[size(sc.C, 1) + 1] | Modelica_LinearSystems2.Inte... |

### Outputs

Type | Name | Description
--- | --- | ---
Real | y[:, size(sc.C, 1), size(sc.B, 2)] | Output response: (number of samples) x (number of outputs) x (number of inputs)
Real | t[:] | Time vector: (number of samples)
Real | x_continuous[:, size(sc.A, 1), size(sc.B, 2)] | State trajectories: (number of samples) x (number of states) x (number of inputs)

### Modelica definition

    partial function timeResponseMask
      "Declares the common structure for the set of response functions"
      input StateSpace sc;
      input Real dt=0 "Sample time [s]";
      input Real tSpan=0 "Simulation time span [s]";
      input Integer toWindow[:,:]=fill(
          -1,
          size(sc.C, 1),
          size(sc.B, 2)) "-1/0/>0 plot in new window/last window/window ID";
      /* the Dialog annotation seems not to work properly. why? */
      input Boolean clearWindow[:,:]=fill(
          false,
          size(sc.C, 1),
          size(sc.B, 2)) "= true, if previous window content is removed";
      input String heading[:,:]=fill(
          "Time response",
          size(sc.C, 1),
          size(sc.B, 2)) "Heading of the response diagram";

      input String columnLabels[size(sc.C, 1) + 1]=
          Modelica_LinearSystems2.Internal.defaultColumnLabels(size(sc.C, 1));

      replaceable output Real y[:,size(sc.C, 1),size(sc.B, 2)]
        "Output response: (number of samples) x (number of outputs) x (number of inputs)";
      output Real t[:] "Time vector: (number of samples)";
      replaceable output Real x_continuous[:,size(sc.A, 1),size(sc.B, 2)]
        "State trajectories: (number of samples) x (number of states) x (number of inputs)";

    end timeResponseMask;

## Modelica_LinearSystems2.Internal.timeResponseMask2

Declares the common structure for the set of response functions

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | sc |  |
Real | dt | 0 | Sample time [s]
Real | tSpan | 0 | Simulation time span [s]

### Outputs

Type | Name | Description
--- | --- | ---
Real | y[:, size(sc.C, 1), size(sc.B, 2)] | Output response: (number of samples) x (number of outputs) x (number of inputs)
Real | t[:] | Time vector: (number of samples)
Real | x_continuous[:, size(sc.A, 1), size(sc.B, 2)] | State trajectories: (number of samples) x (number of states) x (number of inputs)

### Modelica definition

    partial function timeResponseMask2
      "Declares the common structure for the set of response functions"
      input StateSpace sc;
      input Real dt=0 "Sample time [s]";
      input Real tSpan=0 "Simulation time span [s]";

      replaceable output Real y[:,size(sc.C, 1),size(sc.B, 2)]
        "Output response: (number of samples) x (number of outputs) x (number of inputs)";
      output Real t[:] "Time vector: (number of samples)";
      replaceable output Real x_continuous[:,size(sc.A, 1),size(sc.B, 2)]
        "State trajectories: (number of samples) x (number of states) x (number of inputs)";

    end timeResponseMask2;

## Modelica_LinearSystems2.Internal.timeResponseMask2_tf

Declares the common structure for the set of response functions

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
TransferFunction | tf |  |
Real | dt | 0 | Sample time [s]
Real | tSpan | 0 | Simulation time span [s]

### Outputs

Type | Name | Description
--- | --- | ---
Real | y[:, 1, 1] | Output response: (number of samples) x (number of outputs) x (number of inputs)
Real | t[:] | Time vector: (number of samples)
Real | x_continuous[:, Modelica_LinearSystems2.TransferFunction.Analysis.denominatorDegree(tf), 1] | State trajectories: (number of samples) x (number of states) x (number of inputs)

### Modelica definition

    partial function timeResponseMask2_tf
      "Declares the common structure for the set of response functions"
      input Modelica_LinearSystems2.TransferFunction tf;
      input Real dt=0 "Sample time [s]";
      input Real tSpan=0 "Simulation time span [s]";

      replaceable output Real y[:,1,1]
        "Output response: (number of samples) x (number of outputs) x (number of inputs)";
      output Real t[:] "Time vector: (number of samples)";
      replaceable output Real x_continuous[:,Modelica_LinearSystems2.TransferFunction.Analysis.denominatorDegree(tf), 1]
        "State trajectories: (number of samples) x (number of states) x (number of inputs)";

    end timeResponseMask2_tf;

## Modelica_LinearSystems2.Internal.timeResponseMask2_zp

Declares the common structure for the set of response functions

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
ZerosAndPoles | zp |  |
Real | dt | 0 | Sample time [s]
Real | tSpan | 0 | Simulation time span [s]

### Outputs

Type | Name | Description
--- | --- | ---
Real | y[:, 1, 1] | Output response: (number of samples) x (number of outputs) x (number of inputs)
Real | t[:] | Time vector: (number of samples)
Real | x_continuous[:, Modelica_LinearSystems2.ZerosAndPoles.Analysis.denominatorDegree(zp), 1] | State trajectories: (number of samples) x (number of states) x (number of inputs)

### Modelica definition

    partial function timeResponseMask2_zp
      "Declares the common structure for the set of response functions"
      import Modelica_LinearSystems2;

      input Modelica_LinearSystems2.ZerosAndPoles zp;
      input Real dt=0 "Sample time [s]";
      input Real tSpan=0 "Simulation time span [s]";

      replaceable output Real y[:,1,1]
        "Output response: (number of samples) x (number of outputs) x (number of inputs)";
      output Real t[:] "Time vector: (number of samples)";
      replaceable output Real x_continuous[:,Modelica_LinearSystems2.ZerosAndPoles.Analysis.denominatorDegree(zp), 1]
        "State trajectories: (number of samples) x (number of states) x (number of inputs)";

    end timeResponseMask2_zp;

## Modelica_LinearSystems2.Internal.timeResponseMask_tf

Declares the common structure for the set of response functions for a transfer function

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
TransferFunction | tf |  |
Real | dt | 0 | Sample time [s]
Real | tSpan | 0 | Simulation time span [s]
Integer | toWindow[:, :] | fill(-1, 1, 1) | -1/0/>0 plot in new window/last window/window ID
Boolean | clearWindow[:, :] | fill(false, 1, 1) | = true, if previous window content is removed
String | heading[:, :] | fill("Time response", 1, 1) | Heading of the response diagram
String | columnLabels[2] | Modelica_LinearSystems2.Inte... |

### Outputs

Type | Name | Description
--- | --- | ---
Real | y[:, 1, 1] | Output response: (number of samples) x (number of outputs) x (number of inputs)
Real | t[:] | Time vector: (number of samples)
Real | x_continuous[:, Modelica_LinearSystems2.TransferFunction.Analysis.denominatorDegree(tf), 1] | State trajectories: (number of samples) x (number of states) x (number of inputs)

### Modelica definition

    partial function timeResponseMask_tf
      "Declares the common structure for the set of response functions for a transfer function"
      input TransferFunction tf;
      input Real dt=0 "Sample time [s]";
      input Real tSpan=0 "Simulation time span [s]";
      input Integer toWindow[:,:]=fill(
          -1,
          1,
          1) "-1/0/>0 plot in new window/last window/window ID";
      /* the Dialog annotation seems not to work properly. why? */
      input Boolean clearWindow[:,:]=fill(
          false,
          1,
          1) "= true, if previous window content is removed";
      input String heading[:,:]=fill(
          "Time response",
          1,
          1) "Heading of the response diagram";

      input String columnLabels[2]=
          Modelica_LinearSystems2.Internal.defaultColumnLabels(1);

      replaceable output Real y[:,1,1]
        "Output response: (number of samples) x (number of outputs) x (number of inputs)";
      output Real t[:] "Time vector: (number of samples)";
      replaceable output Real x_continuous[:, Modelica_LinearSystems2.TransferFunction.Analysis.denominatorDegree(tf), 1]
        "State trajectories: (number of samples) x (number of states) x (number of inputs)";

    end timeResponseMask_tf;

## Modelica_LinearSystems2.Internal.timeResponseMask_zp

Declares the common structure for the set of response functions for a zeros and poles transfer function

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
ZerosAndPoles | zp |  |
Real | dt | 0 | Sample time [s]
Real | tSpan | 0 | Simulation time span [s]
Integer | toWindow[:, :] | fill(-1, 1, 1) | -1/0/>0 plot in new window/last window/window ID
Boolean | clearWindow[:, :] | fill(false, 1, 1) | = true, if previous window content is removed
String | heading[:, :] | fill("Time response", 1, 1) | Heading of the response diagram
String | columnLabels[2] | Modelica_LinearSystems2.Inte... |

### Outputs

Type | Name | Description
--- | --- | ---
Real | y[:, 1, 1] | Output response: (number of samples) x (number of outputs) x (number of inputs)
Real | t[:] | Time vector: (number of samples)
Real | x_continuous[:, Modelica_LinearSystems2.ZerosAndPoles.Analysis.denominatorDegree(zp), 1] | State trajectories: (number of samples) x (number of states) x (number of inputs)

### Modelica definition

    partial function timeResponseMask_zp
      "Declares the common structure for the set of response functions for a zeros and poles transfer function"
      input ZerosAndPoles zp;
      input Real dt=0 "Sample time [s]";
      input Real tSpan=0 "Simulation time span [s]";
      input Integer toWindow[:,:]=fill(
          -1,
          1,
          1) "-1/0/>0 plot in new window/last window/window ID";
      /* the Dialog annotation seems not to work properly. why? */
      input Boolean clearWindow[:,:]=fill(
          false,
          1,
          1) "= true, if previous window content is removed";
      input String heading[:,:]=fill(
          "Time response",
          1,
          1) "Heading of the response diagram";

      input String columnLabels[2]=
          Modelica_LinearSystems2.Internal.defaultColumnLabels(1);

      replaceable output Real y[:,1,1]
        "Output response: (number of samples) x (number of outputs) x (number of inputs)";
      output Real t[:] "Time vector: (number of samples)";
      replaceable output Real x_continuous[:, Modelica_LinearSystems2.ZerosAndPoles.Analysis.denominatorDegree(zp), 1]
        "State trajectories: (number of samples) x (number of states) x (number of inputs)";

    end timeResponseMask_zp;

## Modelica_LinearSystems2.Internal.timeResponseSamples

Estimate reasonable discretisation sample time and simulation time span for time response plot

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | sc |  |

### Outputs

Type | Name | Description
--- | --- | ---
Real | dt | Sample time
Real | tSpan | Time span

### Modelica definition

    function timeResponseSamples
      "Estimate reasonable discretisation sample time and simulation time span for time response plot"
      import Modelica_LinearSystems2;

      input Modelica_LinearSystems2.StateSpace sc;
      output Real dt "Sample time";
      output Real tSpan "Time span";
    protected
      Modelica_LinearSystems2.Math.Complex eig[size(sc.A, 1)];
      Real realp[size(sc.A, 1)];
      Real sorted[size(sc.A, 1)];
      Real indices[size(sc.A, 1)];
      Integer i;
    algorithm
      eig := Modelica_LinearSystems2.StateSpace.Analysis.eigenValues(sc);
      for i in 1:size(eig, 1) loop
        realp[i] := eig[i].re;
      end for;
      (sorted,indices) := Modelica.Math.Vectors.sort(realp);

      // Estimate simulation time span
      if sorted[end] < 0 then
        tSpan := -5*sorted[end];
      elseif sorted[end] > 0 then
        tSpan := 15/sorted[end];
      elseif sorted[end] == 0 then
        tSpan := 15000;
      end if;
      // Curb simulation time span to maximal 15000s
      if tSpan > 15000 then
        tSpan := 15000;
      end if;

      // set sample time:
      dt := tSpan/1000;

    end timeResponseSamples;

## Modelica_LinearSystems2.Internal.from_Hz

Convert from Hz to rad/s

### Information

Extends from Modelica.SIunits.Conversions.ConversionIcon (Base icon for conversion functions).

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Frequency | f |  | frequency [Hz]

### Outputs

Type | Name | Description
--- | --- | ---
AngularVelocity | w | angular velocity [rad/s]

### Modelica definition

    function from_Hz "Convert from Hz to rad/s"
      extends Modelica.SIunits.Conversions.ConversionIcon;
      input Modelica.SIunits.Frequency f "frequency";
      output Modelica.SIunits.AngularVelocity w "angular velocity";

    algorithm
      w := 2*Modelica.Constants.pi*f;
    end from_Hz;

## Modelica_LinearSystems2.Internal.to_Hz

Convert from rad/s to Hz

### Information

Extends from Modelica.SIunits.Conversions.ConversionIcon (Base icon for conversion functions).

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
AngularVelocity | w |  | angular velocity [rad/s]

### Outputs

Type | Name | Description
--- | --- | ---
Frequency | f | frequency [Hz]

### Modelica definition

    function to_Hz "Convert from rad/s to Hz"
      extends Modelica.SIunits.Conversions.ConversionIcon;
      input Modelica.SIunits.AngularVelocity w "angular velocity";
      output Modelica.SIunits.Frequency f "frequency";
    algorithm
      f := w/(2*Modelica.Constants.pi);
    end to_Hz;

## Modelica_LinearSystems2.StateSpace.Internal

Internal library of record StateSpace (should not be directly used by user)

### Information

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

### Package Content

Name | Description
--- | ---
isSISO | To check a state space system to be SISO (or not)
invariantZeros2 | Compute invariant zeros of linear SISO state space system with a generalized system matrix [A, B, C, D] which is of upper Hessenberg form
characterizeEigenvalue | Check stability, stabilizability, controllability, observability nad detectability of the single poles
isStabilizableSISO | To check wether a SISO system is stabliziable
isStabilizableMIMO | To check wether a MIMO system is stabliziable
isDetectableSISO | To check wether a SISO system is detectable
isDetectableMIMO | To check wether a MIMO system is detectable
isObservableSISO | To check wether a SISO system is observable
isControllableSISO | To check a SISO system wether it is controllable
isControllableMIMO | To check a MIMO system wether it is controllable
isObservableMIMO | To check a MIMO system wether it is observable
isControllableAndObservableSISO | To check whether a SISO system is controllable and observable
readLength_nx | Read the order nx of a state space system from a file
staircaseQR | Staircase algorithm to put a state space system to controller Hessenberg form
staircaseSVD | Staircase algorithm based on singular value decomposition
partialGain |
assignOneOrTwoPoles_alpha | Algorithm to assign p (p = 1 or 2) eigenvalues
assignOneOrTwoPoles | Algorithm to assign p (p = 1 or 2) eigenvalues
readSystemDimension | Read the order nx of state matrix and the numbers nu and ny of inputs and outputs
readLength_nu | Read the number of inputs nu of a state space system from a file
readLength_ny | Read the number of outputs ny of a state space system from a file
dgreeOfRedSys | Calculate the controllable and observable part of a state space system
householder |
numberOfZeros | Calculate the number of zeros of the related transfer function
complexZeros | Calculate the zeros of the related transfer function
numberOfPoles | Calculate the number of poles of the related transfer function
numberOfPolesAndZeros | Calculate the number poles and of zeros of the related transfer function
complexPoles | Generate a zeros-and-poles representation from state space representation
trianUpperHess | Triangulize an upper Hessenberg matrix by repeatedly applicated householder reflexion
reducedCtrSystem2 | calculate the controllable part of a SISO system
controllablePoles | Compute the controllable and uncontrollable poles of a state space system
polesAndZeros | Generate poles and zeros from state space representation
scaleFactor1 | Return scale factor for first order block
scaleFactor2 | Return scale factor for second order block
invariantZerosHessenberg | Fast version to calculate the system zeros of a SISO system with D=[0] and A has upper Hessenberg form, delivered by StateSpace.reduceSystem
cntrHessenberg | calculate the controllable part of a SISO system
transposeStateSpace | Return the transposed state space system
reduceRosenbrock | Algorithm to compress the generalized system matrix [A, B; C, D] to calculate the invariant zeros of a system
reducedCtrSystem | calculate the controllable part of a SISO system
reducedCtrSystemX | calculate the controllable part of a SISO system
assignPolesMI2 | Pole assigment design algorithm for multi input systems

## Modelica_LinearSystems2.StateSpace.Internal.isSISO

To check a state space system to be SISO (or not)

### Information

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  |

### Outputs

Type | Name | Description
--- | --- | ---
Boolean | isSISO |

### Modelica definition

    encapsulated function isSISO
      "To check a state space system to be SISO (or not)"

        import Modelica_LinearSystems2.StateSpace;

      input StateSpace ss;

      output Boolean isSISO;
    algorithm
      isSISO := size(ss.B, 2) == 1 and size(ss.C, 1) == 1;
    end isSISO;

## Modelica_LinearSystems2.StateSpace.Internal.invariantZeros2

Compute invariant zeros of linear SISO state space system with a generalized system matrix [A, B, C, D] which is of upper Hessenberg form

### Information

Computes the invariant zeros of a system in state space form:

       der(x) = A*x + B*u
            y = C*x + D*u

The invariant zeros of this system are defined as the variables z that make the following matrix singular:

        | A B |     | I 0 |
        |     | - z*|     |
        | C D |     | 0 0 |

where I is the identity matrix of the same size as A and 0 are zero matrices of appropriate dimensions.

Unlike to function StateSpace.Analysis.invariantZeros for general systems, it is assumned in StateSpace.Analysis.invariantZeros that the generalized system matrix [A, B; C, D] hast upper Hessenberg form. Especially for SISO system this is achieved when A is of upper Hessenberg form and [1, n] matrix C is of form C = k*[0, 0, ..., 0, 1].

The function uses the LAPACK routine DHGEQZ. Look at Modelica_LinearSystems2.Math.Matrices.LAPACK.dhgeqz for details.

The advantage of this function in comparision to the general invariantZeros function is the lower computatioal effort bacause systems with arbitrary system functions are first transformed into an upper Hessenberg form system.

This function is used in fromStateSpace transformation functions which use Hessenberg form systems anyway.

Currently, there is the restriction that the number of inputs and the number of outputs must be identical. Other systems have to be treated like p*q SISO systems where p is the number of putputs and q the number of inputs of the MIMO system.

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  | Linear system in state space form

### Outputs

Type | Name | Description
--- | --- | ---
Complex | Zeros[:] | Finite, invariant zeros of ss; size(Zeros,1) <= size(ss.A,1)

### Modelica definition

    encapsulated function invariantZeros2
      "Compute invariant zeros of linear SISO state space system with a generalized system matrix [A, B, C, D] which is of upper Hessenberg form"
          import Modelica;
          import Modelica_LinearSystems2.StateSpace;
          import Modelica_LinearSystems2.Math.Complex;
          import Modelica_LinearSystems2.Math.Matrices;
          import Modelica_LinearSystems2.Math.Matrices.LAPACK;

      input StateSpace ss "Linear system in state space form";
      output Complex Zeros[:]
        "Finite, invariant zeros of ss; size(Zeros,1) <= size(ss.A,1)";

    protected
      Integer nx=size(ss.A, 1) "Number of states";
      Integer nu=size(ss.B, 2) "Number of inputs";
      Integer ny=size(ss.C, 1) "Number of outputs";
      Integer na=nx + nu;
      Real A[nx + ny,nx + nu]=[ss.A,ss.B; ss.C,ss.D];
      Real B[nx + ny,nx + nu]=[identity(nx),zeros(nx, nu); zeros(ny, nx + nu)];
      Real alphaReal[na];
      Real alphaImag[na];
      Real beta[na];
      Integer info;
      Real beta_small=100*Modelica.Constants.eps;
      Integer nZeros;
      Complex z[size(ss.A, 1)];
      Integer j;
      Real normB=max(beta_small,Modelica.Math.Matrices.norm(ss.B, p=1));
    algorithm
      assert(nu == ny, "Function invariantZeros requires currently that the number of
    inputs (= "
              + String(nu) + ") = number of outputs (= " + String(ny) + ")
    This condition is however not fulfilled");

    // Compute zeros

      (alphaReal,alphaImag,beta,info) :=
        Matrices.generalizedEigenvaluesTriangular(A, B);

      assert(info == 0,
        "Failed to compute invariant zeros with function invariantZeros(..)");

    // If beta[i] is zero, then zero i is infinite.
      j := 1;
      for i in 1:na loop

       if beta[i] >= normB*1e-6 then
         // finite eigenvalue
          z[j].re := if abs(alphaReal[i]) >= normB*1e-12 then alphaReal[i]/beta[i] else
                  0;
          z[j].im := if abs(alphaImag[i]) >= normB*1e-12 then alphaImag[i]/beta[i] else
                  0;
          j := j + 1;
        end if;
      end for;
      nZeros := j - 1;
      Zeros := z[1:nZeros];
    end invariantZeros2;

## Modelica_LinearSystems2.StateSpace.Internal.characterizeEigenvalue

Check stability, stabilizability, controllability, observability nad detectability of the single poles

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss | Modelica_LinearSystems2.Stat... |
Eigenvalue | evin[:] |  |

### Outputs

Type | Name | Description
--- | --- | ---
Eigenvalue | ev[size(ss.A, 1)] |

### Modelica definition

    function characterizeEigenvalue
      "Check stability, stabilizability, controllability, observability nad detectability of the single poles"

          import Modelica_LinearSystems2;
          import Modelica_LinearSystems2.StateSpace;
          import Modelica_LinearSystems2.Internal.Eigenvalue;
          import Modelica_LinearSystems2.Math.Complex;
          import Modelica_LinearSystems2.Internal;

      input StateSpace ss=Modelica_LinearSystems2.StateSpace(
          A=fill(
            0,
            0,
            0),
          B=fill(
            0,
            0,
            0),
          C=fill(
            0,
            0,
            0),
          D=fill(
            0,
            0,
            0));
      input Eigenvalue evin[:];
      output Eigenvalue ev[size(ss.A, 1)];

    protected
      Real cPoles[:,2] "controllable poles";
      Real ncPoles[:,2] "uncontrollable poles";
      Real poles[size(ss.A, 1),2] "controllable and uncontrollable poles";
      Integer n_c;
      Integer nx=size(ss.A, 1);

      Eigenvalue cdummy[size(ss.A, 1)];
      Eigenvalue odummy[size(ss.A, 1)];
      StateSpace sst=StateSpace(
          A=transpose(ss.A),
          B=transpose(ss.C),
          C=transpose(ss.B),
          D=transpose(ss.D));
      Boolean equal;
      Integer ii;

      Real eps=1e-16*Modelica.Math.Matrices.norm(A=ss.A, p=1);
      Real factor_eps=1;
      Real absVector[nx];
      Integer indices[:];
      Integer indexMin;
      Real indexVector[:];
      Integer vv[:];
      Complex j = Modelica_LinearSystems2.Math.Complex.j();

    algorithm
      for i in 1:nx loop
        ev[i].ev := evin[i].ev;
      end for;

      (cPoles,ncPoles,poles) := StateSpace.Internal.controllablePoles(ss);
      n_c := size(cPoles, 1);
      for i in 1:n_c loop
        cdummy[i].ev := cPoles[i, 1]+cPoles[i, 2]*j;
        cdummy[i].isControllable := true;
      end for;
      for i in 1:size(ncPoles, 1) loop
        cdummy[n_c + i].ev := ncPoles[i, 1]+ncPoles[i, 2]*j;
        cdummy[n_c + i].isControllable := false;
      end for;

    // controllable poles of the trnasposed system are the observable poles
      (cPoles,ncPoles,poles) := StateSpace.Internal.controllablePoles(sst);
      n_c := size(cPoles, 1);
      for i in 1:n_c loop
        odummy[i].ev := cPoles[i, 1]+cPoles[i, 2]*j;
        odummy[i].isObservable := true;
      end for;
      for i in 1:size(ncPoles, 1) loop
        odummy[n_c + i].ev := ncPoles[i, 1]+ncPoles[i, 2]*j;
        odummy[n_c + i].isObservable := false;
      end for;

    // using ev.imag as an flag to mark the eigenvalues
      for i in 1:nx loop
        absVector := fill(1e50, nx);
        for ii in 1:nx loop
          if not odummy[ii].imag then
            absVector[ii] := Complex.'abs'(ev[i].ev - odummy[ii].ev);
          end if;
        end for;
        (absVector,indices) := Modelica.Math.Vectors.sort(absVector);
        indexMin := indices[1];
        ev[i].isObservable := odummy[indexMin].isObservable;
        odummy[indexMin].imag := true;
      end for;

      for i in 1:nx loop
        absVector := fill(1e50, nx);
        for ii in 1:nx loop
          if not cdummy[ii].imag then
            absVector[ii] := Complex.'abs'(ev[i].ev - cdummy[ii].ev);
          end if;
        end for;
        (absVector,indices) := Modelica.Math.Vectors.sort(absVector);
        indexMin := indices[1];
        ev[i].isControllable := cdummy[indexMin].isControllable;
        cdummy[indexMin].imag := true;
      end for;

      for i in 1:nx loop
        ev[i] := Eigenvalue(
          ev[i].ev,
          ev[i].isControllable,
          ev[i].isObservable);
      end for;

    end characterizeEigenvalue;

## Modelica_LinearSystems2.StateSpace.Internal.isStabilizableSISO

To check wether a SISO system is stabliziable

### Information

This function checks whether a SISO state space system is stabilizable or not.

A system is stabilizable for the continuous-time case if all of the uncontrollable eigenvalues have neagtive real part or for the discrete-time case if all of the uncontrollable eigenvalues are in the complex unit circle respectively. Hence, a controllable system is always stabilizable of course.

To check stabilizability, ths system is transformed to to upper controller Hessenberg form

>                    | *   *   ...   ...    * |               | * |
>                    | *   *   ...   ...    * |               | 0 |
>      Q*A*Q ' = H = | 0   *   ...   ...    * |,    Q*b = q = | . |,   c*Q = ( *, ..., * )
>                    | .   .    .     .     . |               | . |
>                    | 0  ...   0     *     * |               | 0 |

The system can be partitioned to

>     H=[H11,H12; H21, H22], q=[q1;0],

H11, q1) contains the controllable part of the system, that is, rank(H) = rank(H11). For stabilizability the H22 has to be stable.

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  |

### Outputs

Type | Name | Description
--- | --- | ---
Boolean | stabilizable |

### Modelica definition

    encapsulated function isStabilizableSISO
      "To check wether a SISO system is stabliziable"

        import Modelica;
        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.StateSpace;
        import Modelica_LinearSystems2.Math.Complex;

      input StateSpace ss;

    protected
      input Modelica_LinearSystems2.Internal.StateSpaceR ssm=
          StateSpace.Internal.cntrHessenberg(ss);
    public
      output Boolean stabilizable;

    protected
      Complex evd[:]=fill(Complex(0), size(ss.A, 1) - ssm.r);

    algorithm
      if size(ss.C, 1) <> 1 or size(ss.B, 2) <> 1 then
        assert(size(ss.B, 2) == 1,
          "A SISO-system is expected as input\n but the number of inputs is " +
          String(size(ss.B, 2)) + " instead of 1");
        assert(size(ss.C, 1) == 1,
          " A SISO-system is expected as input\n but the number of outputs is " +
          String(size(ss.C, 1)) + " instead of 1");
      end if;
      evd := Complex.eigenValues(ssm.A[ssm.r + 1:size(ss.A, 1), ssm.r + 1:size(ss.A,
        1)]);
      stabilizable := true;

      if size(ss.A, 1) > ssm.r then
         for i1 in 1:size(evd, 1) loop
          stabilizable := stabilizable and evd[i1].re < 0;
        end for;
      end if;

    end isStabilizableSISO;

## Modelica_LinearSystems2.StateSpace.Internal.isStabilizableMIMO

To check wether a MIMO system is stabliziable

### Information

This function checks whether a MIMO state space system is stabilizable or not.

A system is stabilizable for the continuous-time case if all of the uncontrollable eigenvalues have neagtive real part or for the discrete-time case if all of the uncontrollable eigenvalues are in the complex unit circle respectively. Hence, a controllable system is always stabilizable of course.

To check stabilizability, staircase algorithm is used to separate the controllable subspace from the uncontrollable subspace. The uncontrollable poles are checked to to stable.

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  |

### Outputs

Type | Name | Description
--- | --- | ---
Boolean | stabilizable |

### Modelica definition

    encapsulated function isStabilizableMIMO
      "To check wether a MIMO system is stabliziable"

        import Modelica;
        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.StateSpace;
        import Modelica_LinearSystems2.Math.Complex;

        input StateSpace ss;

        output Boolean stabilizable;

    protected
        Complex evnc[:] "complex vector of uncontrollable poles";
        Real cPoles[:,2] "controllable poles";
        Real ncPoles[:,2] "uncontrollable poles";
        Real poles[size(ss.A, 1),2] "controllable and uncontrollable poles";
        Complex j=Modelica_LinearSystems2.Math.Complex.j();

    algorithm
        (cPoles,ncPoles,poles) := StateSpace.Internal.controllablePoles(ss);
        evnc := fill(Complex(0), size(ncPoles, 1));
        for i1 in 1:size(ncPoles, 1) loop
          evnc[i1] := ncPoles[i1, 1] + j*ncPoles[i1, 2];
        end for;

        stabilizable := true;

        if size(ss.A, 1) > size(cPoles, 1) then
          for i1 in 1:size(ncPoles, 1) loop
            stabilizable := stabilizable and ncPoles[i1, 1] < 0;
          end for;
        end if;

    end isStabilizableMIMO;

## Modelica_LinearSystems2.StateSpace.Internal.isDetectableSISO

To check wether a SISO system is detectable

### Information

This function checks whether a SISO state space system is detectable or not.

A system is detectable for the continuous-time case if all of the unobservable eigenvalues have neagtive real part or for the discrete-time case if all of the unobservable eigenvalues are in the complex unit circle respectively. Hence, a oberservable system is always detectable of course.

As observability is a dual concept of controllability, the concept of detectability is dual to stabilizability, that is, a system is detectable if the pair (A', C') is stabilizable. Therefore, the same algorithm to check stabilizability are applied to the dual pair (A', C') of the system:

To check stabilizability (see Modelica_LinearSystems2.StateSpace.Analysis.isStabilizable) , ths system is transformed to to upper controller Hessenberg form

>                    | *   *   ...   ...    * |               | * |
>                    | *   *   ...   ...    * |               | 0 |
>      Q*A*Q ' = H = | 0   *   ...   ...    * |,    Q*b = q = | . |,   c*Q = ( *, ..., * )
>                    | .   .    .     .     . |               | . |
>                    | 0  ...   0     *     * |               | 0 |

The system can be partitioned to

>     H=[H11,H12; H21, H22], q=[q1;0],

H11, q1) contains the controllable part of the system, that is, rank(H) = rank(H11). For stabilizability the H22 has to be stable.

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  |

### Outputs

Type | Name | Description
--- | --- | ---
Boolean | detectable |

### Modelica definition

    encapsulated function isDetectableSISO
      "To check wether a SISO system is detectable"

        import Modelica;
        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.StateSpace;
        import Modelica_LinearSystems2.Math.Complex;

      input StateSpace ss;

    protected
      input Modelica_LinearSystems2.Internal.StateSpaceR ssm=
          StateSpace.Internal.cntrHessenberg(
          StateSpace.Internal.transposeStateSpace(ss));
    public
      output Boolean detectable;
    protected
      Complex evd[:]=fill(Complex(0), size(ss.A, 1) - ssm.r);

    algorithm
      if size(ss.C, 1) <> 1 or size(ss.B, 2) <> 1 then
        assert(size(ss.B, 2) == 1,
          "A SISO-system is expected as input\n but the number of inputs is " +
          String(size(ss.B, 2)) + " instead of 1");
        assert(size(ss.C, 1) == 1,
          " A SISO-system is expected as input\n but the number of outputs is " +
          String(size(ss.C, 1)) + " instead of 1");
      end if;
      evd := Complex.eigenValues(ssm.A[ssm.r + 1:size(ss.A, 1), ssm.r + 1:size(ss.A,1)]);

      if size(ss.A, 1) == ssm.r then
        detectable := true;

      else
        detectable := true;
        for i1 in 1:size(evd, 1) loop
          detectable := detectable and evd[i1].re < 0;
        end for;
      end if;

    end isDetectableSISO;

## Modelica_LinearSystems2.StateSpace.Internal.isDetectableMIMO

To check wether a MIMO system is detectable

### Information

This function checks whether a MIMO state space system is detectable or not.

A system is detectable for the continuous-time case if all of the unobservable eigenvalues have negative real part or for the discrete-time case if all of the unobservable eigenvalues are in the complex unit circle respectively. Hence, a observable system is always detectable of course.

To check detectability, staircase algorithm is used to separate the observable subspace from the unobservable subspace. The unobservable poles are checked to be stable.

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  |

### Outputs

Type | Name | Description
--- | --- | ---
Boolean | detectable |

### Modelica definition

    encapsulated function isDetectableMIMO
      "To check wether a MIMO system is detectable"

        import Modelica;
        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.StateSpace;
        import Modelica_LinearSystems2.Math.Complex;

        input StateSpace ss;

        output Boolean detectable;

    protected
        StateSpace sst=StateSpace.Internal.transposeStateSpace(ss);

        Complex evnd[:] "complex vector of uncontrollable poles";
        Real dPoles[:,2] "controllable poles";
        Real ndPoles[:,2] "uncontrollable poles";
        Real poles[size(ss.A, 1),2] "controllable and uncontrollable poles";
        Complex j=Modelica_LinearSystems2.Math.Complex.j();

    algorithm
        (dPoles,ndPoles,poles) := StateSpace.Internal.controllablePoles(sst);
        evnd := fill(Complex(0), size(ndPoles, 1));
        for i1 in 1:size(ndPoles, 1) loop
          evnd[i1] := ndPoles[i1, 1] + j*ndPoles[i1, 2];
        end for;

        detectable := true;

        if size(sst.A, 1) == size(dPoles, 1) then
          detectable := true;
        else
          for i1 in 1:size(ndPoles, 1) loop
            detectable := detectable and ndPoles[i1, 1] < 0;
          end for;
        end if;

    end isDetectableMIMO;

## Modelica_LinearSystems2.StateSpace.Internal.isObservableSISO

To check wether a SISO system is observable

### Information

This function is to calculate whether a SISO state space system is observable or not. Therefore, the dual System (A', c', b', d') it is transformed to upper observer Hessenberg form

>                    | *   *   ...   ...    * |             | * |
>                    | *   *   ...   ...    * |             | . |
>      Q*A'*Q' = H = | 0   *   ...   ...    * |,    Q*c' =  | . |,   b'*Q = q = ( 0, ..., 0, * )
>                    | .   .    .     .     . |             | * |
>                    | 0  ...   0     *     * |             | * |

Note, that

>                              n-1                          n-1
>     rank(c'; c'*A'; ...; c'*A'   ) = rank(q; q*H; ...; q*H   )

and that

>                       n-1
>      (q; H*q; ...; q*H  )

is a lower triangular matrix and has full rank if and only if none of the elements in the diagonal is zero. That is, that neither qn or hi,i-1, i = 2,..., n may be zero.

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  |

### Outputs

Type | Name | Description
--- | --- | ---
Boolean | observable |

### Modelica definition

    encapsulated function isObservableSISO
      "To check wether a SISO system is observable"

        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.StateSpace;

       input StateSpace ss;

    protected
       StateSpace ss2=StateSpace.Internal.transposeStateSpace(ss);

       Modelica_LinearSystems2.Internal.StateSpaceR ssm2=
           StateSpace.Internal.reducedCtrSystem(ss2);
    public
       output Boolean observable;
    algorithm
     if size(ss.C, 1) <> 1 or size(ss.B, 2) <> 1 then
           assert(size(ss.B, 2) == 1,
             "A SISO-system is expected as input\n but the number of inputs is "
              + String(size(ss.B, 2)) + " instead of 1");
           assert(size(ss.C, 1) == 1,
             " A SISO-system is expected as input\n but the number of outputs is "
              + String(size(ss.C, 1)) + " instead of 1");
         end if;
       observable := size(ss.A, 1) == ssm2.r;

    end isObservableSISO;

## Modelica_LinearSystems2.StateSpace.Internal.isControllableSISO

To check a SISO system wether it is controllable

### Information

This function is to calculate whether a SISO state space system is controllable or not. Therefore, it is transformed to lower controller Hessenberg form

>                    | *    *     0   ...  0 |               | 0 |
>                    | .    .     .    .   . |               | . |
>      Q*A*Q ' = H = | *   ...   ...   *   0 |,    Q*b = q = | . |,   c*Q = ( *, ..., * )
>                    | *   ...   ...   *   * |               | 0 |
>                    | *   ...   ...   *   * |               | * |

Note, that

>                        n-1                        n-1
>     rank(b, A*b, ..., A  *b) = rank(q, H*q, ..., H  *q)

and that

>                      n-1
>      (q, H*q, ..., H  *q)

is a lower triangular matrix and has full rank if and only if none of the elements in the diagonal is zero. That is, that neither qn or hi,i+1, i = 1,..., n-1 may be zero.

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  |

### Outputs

Type | Name | Description
--- | --- | ---
Boolean | controllable |

### Modelica definition

    encapsulated function isControllableSISO
      "To check a SISO system wether it is controllable"

        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.StateSpace;

      input StateSpace ss;

    protected
      Modelica_LinearSystems2.Internal.StateSpaceR ssm=
          StateSpace.Internal.reducedCtrSystem(ss);

    public
      output Boolean controllable;
    algorithm
      if size(ss.C, 1) <> 1 or size(ss.B, 2) <> 1 then
        assert(size(ss.B, 2) == 1,
          "A SISO-system is expected as input\n but the number of inputs is "
           + String(size(ss.B, 2)) + " instead of 1");
        assert(size(ss.C, 1) == 1,
          " A SISO-system is expected as input\n but the number of outputs is "
           + String(size(ss.C, 1)) + " instead of 1");
      end if;
      controllable := size(ss.A, 1) == ssm.r;

    end isControllableSISO;

## Modelica_LinearSystems2.StateSpace.Internal.isControllableMIMO

To check a MIMO system wether it is controllable

### Information

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  |
StaircaseMethod | method | Modelica_LinearSystems2.Type... |

### Outputs

Type | Name | Description
--- | --- | ---
Boolean | controllable |

### Modelica definition

    encapsulated function isControllableMIMO
      "To check a MIMO system wether it is controllable"

        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.StateSpace;

        input StateSpace ss;
        input Modelica_LinearSystems2.Types.StaircaseMethod method=
            Modelica_LinearSystems2.Types.StaircaseMethod.SVD;

        output Boolean controllable;
    algorithm
        assert(method == Modelica_LinearSystems2.Types.StaircaseMethod.SVD or
          method == Modelica_LinearSystems2.Types.StaircaseMethod.QR, "\nMethods for staircase algorithm are QR factorization or singular value decomposition. Therefore,
    the variable \"method\" in \"Modelica_LinearSystems2.StateSpace.Internal.isControllableMIMO\" has to be qr or svd but is method = "
           + String(method));
        if min(size(ss.B)) == 0 then
          controllable := false;
        else
          if method == Modelica_LinearSystems2.Types.StaircaseMethod.QR then
            controllable := StateSpace.Internal.staircaseQR(ss);
          else
            controllable := StateSpace.Internal.staircaseSVD(ss);
          end if;
        end if;
    end isControllableMIMO;

## Modelica_LinearSystems2.StateSpace.Internal.isObservableMIMO

To check a MIMO system wether it is observable

### Information

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  |
StaircaseMethod | method | Modelica_LinearSystems2.Type... |

### Outputs

Type | Name | Description
--- | --- | ---
Boolean | observable |

### Modelica definition

    encapsulated function isObservableMIMO
      "To check a MIMO system wether it is observable"

        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.StateSpace;

        input StateSpace ss;
        input Modelica_LinearSystems2.Types.StaircaseMethod method=
            Modelica_LinearSystems2.Types.StaircaseMethod.SVD;

    protected
        StateSpace ss2=StateSpace.Internal.transposeStateSpace(ss);

    public
        output Boolean observable;
    algorithm
        assert(method == Modelica_LinearSystems2.Types.StaircaseMethod.SVD or
          method == Modelica_LinearSystems2.Types.StaircaseMethod.QR, "\nMethods for staircase algorithm are QR factorization or singular value decomposition. Therefore,
    the variable \"method\" in \"Modelica_LinearSystems2.StateSpace.Internal.isControllableMIMO\" has to be qr or svd but is method = "
           + String(method));
        if min(size(ss.C)) == 0 then
          observable := false;
        else
          if method == Modelica_LinearSystems2.Types.StaircaseMethod.QR then
            observable := StateSpace.Internal.staircaseQR(ss2);
          else
            observable := StateSpace.Internal.staircaseSVD(ss2);
          end if;

        end if;

    end isObservableMIMO;

## Modelica_LinearSystems2.StateSpace.Internal.isControllableAndObservableSISO

To check whether a SISO system is controllable and observable

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  |

### Outputs

Type | Name | Description
--- | --- | ---
Boolean | controllableAndObservable |

### Modelica definition

    encapsulated function isControllableAndObservableSISO
      "To check whether a SISO system is controllable and observable"

         import Modelica_LinearSystems2;
         import Modelica_LinearSystems2.StateSpace;

       input StateSpace ss;

    protected
       Modelica_LinearSystems2.Internal.StateSpaceR ssm1=
           StateSpace.Internal.reducedCtrSystem(ss);
       Integer nx=size(ss.A, 1);
       Integer rankQ=ssm1.r;
       StateSpace ss2=StateSpace(
             A=transpose(ssm1.A[nx - rankQ + 1:nx, nx - rankQ + 1:nx]),
             B=transpose(ssm1.C[:, nx - rankQ + 1:nx]),
             C=transpose(ssm1.B[nx - rankQ + 1:nx, :]),
             D=ssm1.D);
       Integer nx2=ssm1.r;
       Modelica_LinearSystems2.Internal.StateSpaceR ssm2=
           StateSpace.Internal.reducedCtrSystem(ss2);
    public
       output Boolean controllableAndObservable;
    algorithm
       if size(ss.C, 1) <> 1 or size(ss.B, 2) <> 1 then
         assert(size(ss.B, 2) == 1,
           "A SISO-system is expected as input\n but the number of inputs is " +
           String(size(ss.B, 2)) + " instead of 1");
         assert(size(ss.C, 1) == 1,
           " A SISO-system is expected as input\n but the number of outputs is "
            + String(size(ss.C, 1)) + " instead of 1");
       end if;
       controllableAndObservable := size(ss.A, 1) == ssm2.r;
    equation

    end isControllableAndObservableSISO;

## Modelica_LinearSystems2.StateSpace.Internal.readLength_nx

Read the order nx of a state space system from a file

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
String | fileName | "ss_siso.mat" | Name of the state space system data file

### Outputs

Type | Name | Description
--- | --- | ---
Integer | nx |

### Modelica definition

    encapsulated function readLength_nx
      "Read the order nx of a state space system from a file"

      input String fileName="ss_siso.mat"
        "Name of the state space system data file";
      output Integer nx;
    protected
      Real nxMat[1,1]=readMatrix(
              fileName,
              "nx",
              1,
              1);
    algorithm
      nx := integer(nxMat[1, 1]);
    end readLength_nx;

## Modelica_LinearSystems2.StateSpace.Internal.staircaseQR

Staircase algorithm to put a state space system to controller Hessenberg form

### Information

This algorithm usues QR factorization to generate staircase form i.e. block upper Hessenberg form of the pair (A,B). Due to the well known problem to determine numerically reliable the rank of a matrix, this algorithm should only be used to well conditioned systems. The best way for rank decision would be singular value decomposition, that is used in staicasSVD.

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  |

### Outputs

Type | Name | Description
--- | --- | ---
Boolean | isControllable |
StateSpaceR | ssm1 | controllable state space system
Real | PP[:, :] |

### Modelica definition

    encapsulated function staircaseQR
      "Staircase algorithm to put a state space system to controller Hessenberg form"

        import Modelica;
        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.StateSpace;
        import Modelica_LinearSystems2.Math.Matrices;
        import Modelica_LinearSystems2.Math.Vectors;

      input StateSpace ss;

      output Boolean isControllable;
      output Modelica_LinearSystems2.Internal.StateSpaceR ssm1(
        redeclare Real A[size(ss.A, 1),size(ss.A, 2)],
        redeclare Real B[size(ss.B, 1),size(ss.B, 2)],
        redeclare Real C[size(ss.C, 1),size(ss.C, 2)],
        redeclare Real D[size(ss.D, 1),size(ss.D, 2)])
        "controllable state space system";
      output Real PP[:,:];
    protected
      Real A[:,:];
      Real B[size(ss.B, 1),size(ss.B, 2)];
      Real C[:,:];
      Real Q[:,:];
      Real Q2[:,:];
      Real R[:,:];

      Real P[:,:];
      Real tau[:];

      Integer nx=size(ss.A, 1);
      Integer nu=size(ss.B, 2);
      Integer n1;
      Boolean stop;
      Real normA=Modelica.Math.Matrices.norm(A=ss.A, p=1);
      Real eps=normA*1e-10;
      Integer stairStepinSys;
      Integer info;
      Integer nn;
      Integer stairStep;
      Integer rankR;

    algorithm
      if nu > 0 then
      if nx > 1 then

         //#####  first step of staircase
           // transform b->Q'b = {*,0,...,0} and c->cQ, A->Q'AQ

        (Q,R,tau,Q2) := Matrices.QR( ss.B);
        B := [R; zeros(nx - nu, nu)];  // should be the same as transopse(Q2)*ss.B

        A := transpose(Q2)*ss.A;
        A := A*Q2;
        C := ss.C*Q2;
        PP := transpose(Q2);

        stairStep := 0;
        rankR := 0;
         // for i in 1:size(R, 1) loop
         //   if abs(R[i, i]) > eps then
         //     rankR := rankR + 1;
         //   end if;
         // end for;

    //  !!!! rank has to be determined. In the case of ill conditioned systems svd should be used
        for i in 1:min(size(R, 1), size(R, 2)) loop
          if abs(R[i, size(R, 2) - min(size(R, 1), size(R, 2)) + i]) > eps then
            rankR := rankR + 1;
          end if;
        end for;
        stairStep := stairStep + rankR;
        n1 := nx - stairStep;
        stop := false;

         // #######  buildig rest of staircase
        while not stop loop

          (Q,R,tau,Q2) := Matrices.QR( A[stairStep + 1:nx, stairStep - rankR +
            1:stairStep]);
          P := [identity(nx - n1),zeros(nx - n1, n1); zeros(n1, nx - n1),Q2];
          PP := transpose(P)*PP;
          A := [A[1:stairStep, 1:stairStep],A[1:stairStep, stairStep + 1:nx]*Q2;
            transpose(Q2)*A[stairStep + 1:nx, 1:stairStep],transpose(Q2)*A[
            stairStep + 1:nx, stairStep + 1:nx]*Q2];
           //=transpose(P)*A*P = [A11, A12*Q2; transpose(Q2)*A21, transpose(Q2)*A22*Q2]
          C[:, nx - n1 + 1:nx] := C[:, nx - n1 + 1:nx]*Q2;

          rankR := 0;
           // for i in 1:size(R, 1) loop
           //   if abs(R[i, i]) > eps then
           //     rankR := rankR + 1;
           //   end if;
           // end for;

    //  !!!! rank has to be determined. In the case of ill conditioned systems svd should be used
          for i in 1:min(size(R, 1), size(R, 2)) loop
            if abs(R[i, size(R, 2) - min(size(R, 1), size(R, 2)) + i]) > eps then
              rankR := rankR + 1;
            end if;
          end for;
          stairStep := stairStep + rankR;
          n1 := if rankR < 1 then -1 else n1 - rankR;
          stop := n1 <= 0;

        end while;
      else
        stairStep := if Modelica.Math.Matrices.isEqual(ss.B, zeros(size(ss.B, 1), size(ss.B, 2))) then 0 else
                1;
        A := ss.A;
        B := ss.B;
        C := ss.C;
      end if;

      ssm1 := Modelica_LinearSystems2.Internal.StateSpaceR(
            A=A,
            B=B,
            C=C,
            D=ss.D,
            r=stairStep);

      isControllable := stairStep == nx;
      else // no inputs, nu==0
      isControllable := false;
      ssm1 :=  Modelica_LinearSystems2.Internal.StateSpaceR(
            A=ss.A,
            B=ss.B,
            C=ss.C,
            D=ss.D,
            r=0);
      P := identity(nu);
      end if;

    end staircaseQR;

## Modelica_LinearSystems2.StateSpace.Internal.staircaseSVD

Staircase algorithm based on singular value decomposition

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  |

### Outputs

Type | Name | Description
--- | --- | ---
Boolean | isControllable |
StateSpaceR | ssm1 | upper block Hessenberg form state space system
Real | P[:, :] |

### Modelica definition

    encapsulated function staircaseSVD
      "Staircase algorithm based on singular value decomposition"

        import Modelica;
        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.StateSpace;
        import Modelica_LinearSystems2.Math.Vectors;
        import Modelica_LinearSystems2.Math.Complex;

      input StateSpace ss;

      output Boolean isControllable;
      output Modelica_LinearSystems2.Internal.StateSpaceR ssm1(
        redeclare Real A[size(ss.A, 1),size(ss.A, 2)],
        redeclare Real B[size(ss.B, 1),size(ss.B, 2)],
        redeclare Real C[size(ss.C, 1),size(ss.C, 2)],
        redeclare Real D[size(ss.D, 1),size(ss.D, 2)])
        "upper block Hessenberg form state space system";
      output Real P[:,:];

    protected
      Real A[:,:];
      Real B[size(ss.B, 1),size(ss.B, 2)];
      Real C[:,:];
      Real U[:,:];
      Real VT[:,:];

      Real Q[:,:];
      Real Pi[:,:];
      Real tau[:];
      Real sigma[:];

      Integer nx=size(ss.A, 1);
      Integer nu=size(ss.B, 2);
      Integer ni;
      Boolean stop;
      Real normA=Modelica.Math.Matrices.norm(A=ss.A, p=1);
      Real eps=normA*1e-10;
      Integer stairStepinSys;
      Integer info;
      Integer nn;
      Integer stairStep;
      Integer rankS;
      Integer maxA=max(size(ss.A));
    //   Real evc[:,2];
    //   Real evnc[:,2];

    algorithm
      if nu > 0 then
      if nx > 1 then
    //#####  first step of staircase
            // transform b->Q'b = {*,0,...,0} and c->cQ, A->Q'AQ
        (sigma,U,VT) := Modelica.Math.Matrices.singularValues(ss.B);

        rankS := 0;
        for i in 1:size(sigma, 1) loop
          if sigma[i] > maxA*sigma[1]*Modelica.Constants.eps then
            rankS := rankS + 1;
          end if;
        end for;

        B := [diagonal(sigma[1:rankS]),zeros(rankS, nu - rankS); zeros(nx - rankS, nu)];
        P := transpose(U);
        Q := transpose(VT);
        A := transpose(U)*ss.A*U;
        C := ss.C*U;
        (sigma,U,VT) := Modelica.Math.Matrices.singularValues(A[rankS + 1:nx, 1:rankS]);

        stairStep := rankS;
        rankS := 0;
        if size(sigma, 1) > 1 then
          for i in 1:size(sigma, 1) loop
            if sigma[i] > maxA*sigma[1]*Modelica.Constants.eps then
              rankS := rankS + 1;
            end if;
          end for;
        else
          rankS := if size(sigma, 1) > 0 then if sigma[1] > maxA*Modelica.Constants.eps then 1 else 0 else 0;
        end if;

        stairStep := stairStep + rankS;
        Pi := [VT,zeros(size(VT, 1), size(U, 1)); zeros(size(U, 2), size(VT, 2)),transpose(U)];
        B := Pi*B;
        P := Pi*P;

     // P could be ambigious according to the sign
        if transpose(P)*B[:,1]*ss.B[:,1]<0 then
          Pi:=-Pi;
          P:=-P;
        end if;

        A := Pi*A*transpose(Pi);

       // should be made better because of many zeros in B

        while stairStep < nx and rankS > 0 and not
            Modelica.Math.Matrices.isEqual(
                A[stairStep + 1:nx, stairStep - rankS + 1:stairStep],
                zeros(stairStep - rankS, rankS),
                eps) loop

          (sigma,U,VT) := Modelica.Math.Matrices.singularValues(A[stairStep + 1:nx, stairStep
             - rankS + 1:stairStep]);

          Pi := [identity(stairStep - rankS),zeros(stairStep - rankS, nx -
            stairStep + rankS); zeros(nx - stairStep + rankS, stairStep - rankS),
            [VT,zeros(rankS, nx - stairStep); zeros(nx - stairStep, rankS),
            transpose(U)]];
          P := Pi*P;
          A := Pi*A*transpose(Pi);

    //new implenmentation necessary because of many zeros in Pi
          C := C*transpose(Pi);
          rankS := 0;
          if size(sigma, 1) > 1 then
            for i in 1:size(sigma, 1) loop
              if sigma[i] > maxA*sigma[1]*Modelica.Constants.eps then
                rankS := rankS + 1;
              end if;
            end for;
          else
            rankS := if size(sigma, 1) > 0 then if sigma[1] > maxA*Modelica.Constants.eps then
                    1 else 0 else 0;
          end if;
          stairStep := stairStep + rankS;
        end while;

      else
        stairStep := if Modelica.Math.Matrices.isEqual(ss.B, zeros(size(ss.B, 1),
          size(ss.B, 2))) then 0 else 1;
        A := ss.A;
        B := ss.B;
        C := ss.C;
      end if;

      ssm1 := Modelica_LinearSystems2.Internal.StateSpaceR(
            A=A,
            B=B,
            C=C,
            D=ss.D,
            r=stairStep);

    //   evc := Modelica.Math.Matrices.eigenValues(A[1:stairStep, 1:stairStep]);
    //   evnc := Modelica.Math.Matrices.eigenValues(A[stairStep + 1:nx, stairStep + 1:nx]);

      isControllable := stairStep == nx;

      else // no inputs, nu==0
      isControllable := false;
      ssm1 :=  Modelica_LinearSystems2.Internal.StateSpaceR(
            A=ss.A,
            B=ss.B,
            C=ss.C,
            D=ss.D,
            r=0);
      P := identity(nu);
      end if;

    end staircaseSVD;

## Modelica_LinearSystems2.StateSpace.Internal.partialGain

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | H[:, size(H, 1)] |  |
Real | b[size(H, 1)] |  |

### Outputs

Type | Name | Description
--- | --- | ---
Real | result |

### Modelica definition

    encapsulated function partialGain
        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.StateSpace;
        import Modelica_LinearSystems2.StateSpace.Internal;

      input Real H[:,size(H, 1)];
      input Real b[size(H, 1)];
      output Real result;
    protected
      Real Hh[:,:]=H;
      Real bh[:]=b;
      Integer q=size(H, 1);
    algorithm

      (Hh,bh) := Internal.trianUpperHess(Hh, bh);
      result := bh[q]/Hh[q, q];

    end partialGain;

## Modelica_LinearSystems2.StateSpace.Internal.assignOneOrTwoPoles_alpha

Algorithm to assign p (p = 1 or 2) eigenvalues

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | F[:, size(F, 1)] |  | system matrix of order p=1 or p=2
Real | G[size(F, 1), :] |  | control input matrix p rows
Complex | gamma[size(F, 1)] |  |
Real | tolerance | Modelica.Constants.eps |

### Outputs

Type | Name | Description
--- | --- | ---
Real | K[:, size(F, 1)] | feedback matrix p columns

### Modelica definition

    encapsulated function assignOneOrTwoPoles_alpha
      "Algorithm to assign p (p = 1 or 2) eigenvalues"

        import Modelica;
        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.Math.Complex;
        import Modelica_LinearSystems2.Math.Vectors;

      input Real F[:,size(F, 1)] "system matrix of order p=1 or p=2";
      input Real G[size(F, 1),:] "control input matrix p rows";
      input Complex gamma[size(F, 1)];
      input Real tolerance=Modelica.Constants.eps;
      output Real K[:,size(F, 1)] "feedback matrix p columns";

    protected
      Real Gamma[:,:];
      Integer rankGs;
      Real Fs[size(F, 1),size(F, 2)];
      Real Gs[size(G, 1),size(G, 2)];
      Real Gst[:,:]=transpose(G);
      Real Ks[:,size(F, 1)];
      Real c;
      Real s;
      Real r;

      Real V1[size(G, 2),size(G, 2)];
      Real V2[size(G, 2),size(G, 2)];
      Real V[size(G, 2),size(G, 2)];
      Real U[size(F, 1),size(F, 2)];

      Real u1[:];
      Real u2[:];

    algorithm
      assert(size(F, 1) >= size(gamma, 1),
        "\n In function StateSpace.Internal.assignOneOrTwoPoles() matrix F is of size ["
         + String(size(F, 1)) + "," + String(size(F, 1)) + "] and " + String(
        size(F, 1)) + " demanded assigned poles are expected. However, " +
        String(size(gamma, 1)) + " poles are given");
    //assert(not Modelica.Math.Matrices.isEqual(G,zeros(size(G,1),size(G,2)),tolerance),"A subsystem (F, G) in StateSpace.Internal.assignOneOrTwoPoles() is not controllable, since G is equal to zero matrix ");
      if size(gamma, 1) == 1 then
        assert(gamma[1].im == 0, "\n In function StateSpace.Internal.assignOneOrTwoPoles() matrix F has size [" + String(size(F, 1)) + "," + String(size(F, 1)) +
          "], therefore, the demanded assigned pole must be real. However, the imaginary part is "
           + String(gamma[1].im));
      elseif abs(gamma[1].im) > 0 or abs(gamma[2].im) > 0 then
        assert(gamma[1].re == gamma[2].re and gamma[1].im == -gamma[2].im,
          "\nThe assigned pole pair given in function StateSpace.Internal.assignOneOrTwoPoles() must be conjungated complex. However, the poles are\npole1 = "
           + String(gamma[1]) + "\npole2 = " + String(gamma[2]) +
          ". \nTry\npole1 = " + String(gamma[1]) + "\npole2 = " + String(
          Complex.conj(gamma[1])) + "\ninstead");
      end if;

      if not Modelica.Math.Matrices.isEqual(
              G,
              zeros(size(G, 1), size(G, 2)),
              tolerance) then
        if size(G, 2) == 1 then
          V := [1];
          if size(G, 1) == 1 then
            U := [1];
          else
             // Givens
            r := sqrt(G[1, 1]^2 + G[2, 1]^2);
            c := G[1, 1]/r;
            s := G[2, 1]/r;
            U := [c,s; -s,c];
          end if;
          Gs := U*G;

          rankGs := if abs(Gs[1, 1]) > tolerance then 1 else 0;
        else
         // size(G, 2)>1

          if size(G, 1) == 1 then // U=I, compute V by just one Householder transformation
            U := [1];
            u1 := cat(1, Vectors.householderVector(Gst[:, 1], cat(1, {1}, zeros(size(G, 2) - 1))));
                                     // Householder vector
            Gst := Modelica_LinearSystems2.Math.Matrices.householderReflexion(
              Gst, u1);

            V := identity(size(G, 2)) - 2*matrix(u1)*transpose(matrix(u1))/(u1*u1);
            Gs := transpose(Gst);
            rankGs := if abs(Gs[1, 1]) > tolerance then 1 else 0;

          else
            //2xHH + Givens
            u1 := cat(1, Vectors.householderVector(Gst[:, 1], cat(
                  1,
                  {1},
                  zeros(size(G, 2) - 1))));
                                     // Householder vector1
            Gst := Modelica_LinearSystems2.Math.Matrices.householderReflexion(
              Gst, u1);
            V1 := identity(size(G, 2)) - 2*matrix(u1)*transpose(matrix(u1))/(u1
              *u1);

    // if rank of G of a multi input system is equal to 1
            if Modelica.Math.Vectors.isEqual(Gst[:, 2], zeros(size(G, 2)), tolerance) or
              Modelica.Math.Matrices.isEqual(Gst[2:size(Gst, 1), :], zeros(size(Gst, 1) - 1, size(Gst, 2)), tolerance) then
              V := V1;
              rankGs := if abs(Gs[1, 1]) > tolerance then 1 else 0;
            else

              u2 := cat(
                    1,
                    zeros(1),
                    Vectors.householderVector(Gst[2:size(G, 2), 2], cat(
                      1,
                      {1},
                      zeros(size(G, 2) - 2))));
                                       // Householder vector2
              Gst := Modelica_LinearSystems2.Math.Matrices.householderReflexion(
                Gst, u2);

              V2 := identity(size(G, 2)) - 2*matrix(u2)*transpose(matrix(u2))/(
                u1*u1);
              V := V2*V1;

            end if;

            Gs := transpose(Gst);

            rankGs := 0;
            for i in 1:2 loop
              if abs(Gs[i, i]) > tolerance then
                rankGs := rankGs + 1;
              end if;
            end for;
         // Givens rotation to transfotm Gs[1:2,1:2] to right upper triangle
            r := sqrt(Gs[1, 1]^2 + Gs[2, 1]^2);
            c := Gs[1, 1]/r;
            s := Gs[2, 1]/r;
            U := [c,s; -s,c];
            Gs := U*Gs;

          end if;
              end if;

    // check controllability
        assert(not Modelica.Math.Matrices.isEqual(
              Gs,
              zeros(size(Gs, 1), size(Gs, 2)),
              tolerance),
          "A subsystem in StateSpace.Internal.assignOneOrTwoPoles() is not controllable");

        Ks := fill(
              0,
              rankGs,
              size(F, 1));
        Fs := U*F*transpose(U);

        if size(F, 1) == 1 then
          Ks := matrix((Fs[1, 1] - gamma[1].re)/Gs[1, 1]);
        else
          if rankGs == size(F, 1) then

            // Gamma:= if size(F,1)==1 then [gamma[1].re] else [gamma[1].re, -(gamma[1].im)^2;1, gamma[2].re];
            //  Ks :=  Modelica_LinearSystems2.Math.Matrices.solve2(Gs, Fs - Gamma);
            Ks := [(Fs[1, 1] - gamma[1].re)/Gs[1, 1] - Gs[1, 2]*(Fs[2, 1] - 1)/Gs[1, 1]/Gs[2,
              2],(Fs[1, 2] + (gamma[1].im)^2)/Gs[1, 1] - Gs[1, 2]*(Fs[2, 2] -
              gamma[2].re)/Gs[1, 1]/Gs[2, 2]; (Fs[2, 1] - 1)/Gs[2, 2],(Fs[2, 2]
               - gamma[2].re)/Gs[2, 2]];
          else

            Ks[1, 1] := (gamma[1].re + gamma[2].re - Fs[1, 1] - Fs[2, 2])/Gs[1,
              1];
            Ks[1, 2] := Ks[1, 1]*Fs[2, 2]/Fs[2, 1] + (Fs[1, 1]*Fs[2, 2] - Fs[1,
              2]*Fs[2, 1] - (gamma[1].re*gamma[2].re - gamma[1].im*gamma[2].im))
              /Fs[2, 1]/Gs[1, 1];
            Ks := -Ks;
          end if;
        end if;

    //    K := transpose(V)*[Ks; zeros(size(G, 2) - rankGs, size(Ks, 2))]*(U);
        K := transpose(V[1:size(Ks, 1), :])*Ks*U;

      else
        Modelica.Utilities.Streams.print(
          "\n A subsystem (F, G) in StateSpace.Internal.assignOneOrTwoPoles() is not controllable, since G is equal to zero matrix. Therefore, K is set to zero matrix and the eigenvalues are retained");
        K := zeros(size(G, 2), size(F, 1));
      end if;

    end assignOneOrTwoPoles_alpha;

## Modelica_LinearSystems2.StateSpace.Internal.assignOneOrTwoPoles

Algorithm to assign p (p = 1 or 2) eigenvalues

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | F[:, size(F, 1)] |  | system matrix of order p=1 or p=2
Real | G[size(F, 1), :] |  | control input matrix p rows
Complex | gamma[size(F, 1)] |  |
Real | tolerance | Modelica.Constants.eps |

### Outputs

Type | Name | Description
--- | --- | ---
Real | K[:, size(F, 1)] | feedback matrix p columns

### Modelica definition

    encapsulated function assignOneOrTwoPoles
      "Algorithm to assign p (p = 1 or 2) eigenvalues"

        import Modelica;
        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.Math.Complex;
        import Modelica_LinearSystems2.Math.Vectors;

      input Real F[:,size(F, 1)] "system matrix of order p=1 or p=2";
      input Real G[size(F, 1),:] "control input matrix p rows";
      input Complex gamma[size(F, 1)];
      input Real tolerance=Modelica.Constants.eps;
      output Real K[:,size(F, 1)] "feedback matrix p columns";

    protected
      Real Gamma[:,:];
      Integer rankGs;
      Real Fs[size(F, 1),size(F, 2)];
      Real Gs[size(G, 1),size(G, 2)];
      Real Gst[:,:]=transpose(G);
      Real Ks[:,size(F, 1)];
      Real c;
      Real s;
      Real r;
      Integer p=size(F,1);
      Real sigmaG[:];

      Real V[size(G, 2),size(G, 2)];
      Real U[size(F, 1),size(F, 2)];

      Real u1[:];
      Real u2[:];
      Integer i;
      Complex system_ev[:];

    algorithm
      assert(size(F, 1) >= size(gamma, 1),
        "\n In function StateSpace.Internal.assignOneOrTwoPoles() matrix F is of size ["
         + String(size(F, 1)) + "," + String(size(F, 1)) + "] and " + String(
        size(F, 1)) + " demanded assigned poles are expected. However, " +
        String(size(gamma, 1)) + " poles are given");
    //assert(not Modelica.Math.Matrices.isEqual(G,zeros(size(G,1),size(G,2)),tolerance),"A subsystem (F, G) in StateSpace.Internal.assignOneOrTwoPoles() is not controllable, since G is equal to zero matrix ");
      if size(gamma, 1) == 1 then
        assert(gamma[1].im == 0, "\n In function StateSpace.Internal.assignOneOrTwoPoles() matrix F has size [" + String(size(F, 1)) + "," + String(size(F, 1)) +
          "], therefore, the demanded assigned pole must be real. However, the imaginary part is "
           + String(gamma[1].im));
      elseif abs(gamma[1].im) > 0 or abs(gamma[2].im) > 0 then
        assert(gamma[1].re == gamma[2].re and gamma[1].im == -gamma[2].im,
          "\nThe assigned pole pair given in function StateSpace.Internal.assignOneOrTwoPoles() must be conjungated complex. However, the poles are\npole1 = "
           + String(gamma[1]) + "\npole2 = " + String(gamma[2]) +
          ". \nTry\npole1 = " + String(gamma[1]) + "\npole2 = " + String(
          Complex.conj(gamma[1])) + "\ninstead");
      end if;

      if not Modelica.Math.Matrices.isEqual(
          G,
          zeros(size(G, 1), size(G, 2)),
          tolerance) then
        if size(G, 2) == 1 then
          V := [1];
          if size(G, 1) == 1 then
            U := [1];
          else
             // Givens
            r := sqrt(G[1, 1]^2 + G[2, 1]^2);
            c := G[1, 1]/r;
            s := G[2, 1]/r;
            U := [c,s; -s,c];
          end if;
          Gs := U*G;

          rankGs := if abs(Gs[1, 1]) > tolerance then 1 else 0;
        else
         // size(G, 2)>1

          if size(G, 1) == 1 then // U=I, compute V by just one Householder transformation
            U := [1];
            u1 := cat(1, Vectors.householderVector(Gst[:, 1],
                         cat(1, {1}, zeros(size(G, 2) - 1))));// Householder vector
            Gst := Modelica_LinearSystems2.Math.Matrices.householderReflexion(Gst, u1);

            V := identity(size(G, 2)) - 2*matrix(u1)*transpose(matrix(u1))/(u1*u1);
            Gs := transpose(Gst);
            rankGs := if abs(Gs[1, 1]) > tolerance then 1 else 0;

          else
    // systems with p==2 and m>1 are transformed by svd
            (sigmaG,U,V) := Modelica.Math.Matrices.singularValues(G);
            rankGs := 0;
            i := size(sigmaG, 1);
            while i > 0 loop
              if sigmaG[i] > 1e-10 then
                rankGs := i;
                i := 0;
              end if;
              i := i - 1;
            end while;
            Gs := zeros(p, size(G, 2));
            for i in 1:rankGs loop
              Gs[i, i] := sigmaG[i];
            end for;

          end if;
          V := transpose(V);
        end if;

    // check controllability
        assert(not Modelica.Math.Matrices.isEqual(
          Gs,
          zeros(size(Gs, 1), size(Gs, 2)),
          tolerance), "A subsystem in StateSpace.Internal.assignOneOrTwoPoles() is not controllable");

        Ks := fill(
          0,
          rankGs,
          size(F, 1));
        Fs := U*F*transpose(U);

        if size(F, 1) == 1 then
          Ks := matrix((Fs[1, 1] - gamma[1].re)/Gs[1, 1]);
        else
          if rankGs == size(F, 1) then

            // Gamma:= if size(F,1)==1 then [gamma[1].re] else [gamma[1].re, -(gamma[1].im)^2;1, gamma[2].re];
            //  Ks :=  Modelica_LinearSystems2.Math.Matrices.solve2(Gs, Fs - Gamma);
            Ks := [(Fs[1, 1] - gamma[1].re)/Gs[1, 1] - Gs[1, 2]*(Fs[2, 1] - 1)/Gs[1, 1]/Gs[2,2],
            (Fs[1, 2] + (gamma[1].im)^2)/Gs[1, 1] - Gs[1, 2]*(Fs[2, 2] - gamma[2].re)/Gs[1, 1]/Gs[2, 2];
            (Fs[2, 1] - 1)/Gs[2, 2],(Fs[2, 2] - gamma[2].re)/Gs[2, 2]];
          else

            Ks[1, 1] := (gamma[1].re + gamma[2].re - Fs[1, 1] - Fs[2, 2])/Gs[1, 1];
            Ks[1, 2] := Ks[1, 1]*Fs[2, 2]/Fs[2, 1] + (Fs[1, 1]*Fs[2, 2] - Fs[1, 2]*
              Fs[2, 1] - (gamma[1].re*gamma[2].re - gamma[1].im*gamma[2].im))/Fs[2,1]/Gs[1, 1];
            Ks := -Ks;
          end if;
        end if;

        K := V[:, 1:size(Ks, 1)]*Ks*U;

      else
        if p == 1 then
          Modelica.Utilities.Streams.print("\n A subsystem (F, G) in StateSpace.Internal.assignOneOrTwoPoles() is not controllable, since G is equal to zero matrix. Therefore, K is set to zero matrix and the eigenvalues are retained.\n
          That is, " + String(F[1, 1]) + " remains and " + String(gamma[1].re) + " cannot be realized");
        else
          system_ev := Complex.eigenValues(F);
          Modelica.Utilities.Streams.print("\n A subsystem (F, G) in StateSpace.Internal.assignOneOrTwoPoles() is not controllable, since G is equal to zero matrix. Therefore, K is set to zero matrix and the eigenvalues are retained.\n
          That is, " + String(system_ev[1].re) + (if abs(system_ev[1].im) > 0 then " + " else
                  " - ") + String(system_ev[1].im) + "j and " + String(system_ev[2].re)
             + (if abs(system_ev[2].im) > 0 then " + " else " - ") + String(
            system_ev[2].im) + "j remain and " + String(gamma[1].re) + (if abs(
            gamma[1].im) > 0 then (if gamma[1].im > 0 then " + " else " - " +
            String(gamma[1].im) + "j") else "" + " and ") + String(gamma[2].re) + (
            if abs(gamma[2].im) > 0 then (if gamma[2].im > 0 then " + " else " - " +
            String(gamma[2].im) + "j") else "") + " cannot be realized");
        end if;
        K := zeros(size(G, 2), size(F, 1));
      end if;

    end assignOneOrTwoPoles;

## Modelica_LinearSystems2.StateSpace.Internal.readSystemDimension

Read the order nx of state matrix and the numbers nu and ny of inputs and outputs

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
String | fileName | "stateSpace.mat" |
String | matrixName | "ABCD" | Name of the generalized state space system matrix

### Outputs

Type | Name | Description
--- | --- | ---
Integer | xuy[3] |

### Modelica definition

    encapsulated function readSystemDimension
      "Read the order nx of state matrix and the numbers nu and ny of inputs and outputs"
        import Modelica_LinearSystems2;
      input String fileName="stateSpace.mat";
      input String matrixName="ABCD"
        "Name of the generalized state space system matrix";
      output Integer xuy[3];

    protected
      Real sizeA[1,1]=readMatrix(
              fileName,
              "nx",
              1,
              1);

      Integer ABCDsizes[2]=readMatrixSize(fileName, matrixName);

    algorithm
      xuy[1] := integer(sizeA[1, 1]);
      xuy[2] := ABCDsizes[2] - xuy[1];
      xuy[3] := ABCDsizes[1] - xuy[1];

    end readSystemDimension;

## Modelica_LinearSystems2.StateSpace.Internal.readLength_nu

Read the number of inputs nu of a state space system from a file

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
String | fileName | "ss_siso.mat" | Name of the state space system data file
String | matrixName | "ABCD" | Name of the state space system matrix

### Outputs

Type | Name | Description
--- | --- | ---
Integer | nu |

### Modelica definition

    encapsulated function readLength_nu
      "Read the number of inputs nu of a state space system from a file"

      input String fileName="ss_siso.mat"
        "Name of the state space system data file";
      input String matrixName="ABCD" "Name of the state space system matrix";

      output Integer nu;
    protected
      Real nxMat[1,1]=readMatrix(
              fileName,
              "nx",
              1,
              1);
      Integer ABCDsizes[2]=readMatrixSize(fileName, matrixName);
      Integer nx=integer(nxMat[1, 1]);

    algorithm
      nu := ABCDsizes[2] - nx;
    end readLength_nu;

## Modelica_LinearSystems2.StateSpace.Internal.readLength_ny

Read the number of outputs ny of a state space system from a file

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
String | fileName | "ss_siso.mat" | Name of the state space system data file
String | matrixName | "ABCD" | Name of the state space system matrix

### Outputs

Type | Name | Description
--- | --- | ---
Integer | ny |

### Modelica definition

    encapsulated function readLength_ny
      "Read the number of outputs ny of a state space system from a file"

      input String fileName="ss_siso.mat"
        "Name of the state space system data file";
      input String matrixName="ABCD" "Name of the state space system matrix";

      output Integer ny;
    protected
      Real nxMat[1,1]=readMatrix(
              fileName,
              "nx",
              1,
              1);
      Integer ABCDsizes[2]=readMatrixSize(fileName, matrixName);
      Integer nx=integer(nxMat[1, 1]);

    algorithm
      ny := ABCDsizes[1] - nx;
    end readLength_ny;

## Modelica_LinearSystems2.StateSpace.Internal.dgreeOfRedSys

Calculate the controllable and observable part of a state space system

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  |

### Outputs

Type | Name | Description
--- | --- | ---
Integer | degree_rs |

### Modelica definition

    encapsulated function dgreeOfRedSys
      "Calculate the controllable and observable part of a state space system"

        import Modelica_LinearSystems2.StateSpace;
        import Modelica_LinearSystems2.Math.Complex;
        import Modelica_LinearSystems2;

      input StateSpace ss;
      output Integer degree_rs;

    protected
      StateSpace ssm1=
          Modelica_LinearSystems2.StateSpace.Internal.reducedCtrSystem(ss);

      StateSpace ss2=StateSpace(
              A=transpose(ssm1.A),
              B=transpose(ssm1.C),
              C=transpose(ssm1.B),
              D=ssm1.D);

      StateSpace ssm2=
          Modelica_LinearSystems2.StateSpace.Internal.reducedCtrSystem(ss2);

    algorithm
      degree_rs := size(ssm2.A, 1);

    end dgreeOfRedSys;

## Modelica_LinearSystems2.StateSpace.Internal.householder

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  |
Real | v[size(ss.A, 1)] |  |

### Outputs

Type | Name | Description
--- | --- | ---
StateSpace | ssh |

### Modelica definition

    encapsulated function householder
        import Modelica_LinearSystems2.StateSpace;
        import Modelica_LinearSystems2.Math.Matrices;
        import Modelica_LinearSystems2.Math.Vectors;

      input StateSpace ss;
      input Real v[size(ss.A, 1)];

      output StateSpace ssh;

    protected
      Real Ah[size(ss.A, 1),size(ss.A, 2)];
      Real Bh[size(ss.B, 1),size(ss.B, 2)];
      Real Ch[size(ss.C, 1),size(ss.C, 2)];
    algorithm
      Ah := Matrices.householderSimilarityTransformation(ss.A, v);
      Bh := Matrices.householderReflexion(ss.B, v);
      Ch := transpose(Matrices.householderReflexion(transpose(ss.C), v));

      ssh := StateSpace(
            A=Ah,
            B=Bh,
            C=Ch,
            D=ss.D);

    end householder;

## Modelica_LinearSystems2.StateSpace.Internal.numberOfZeros

Calculate the number of zeros of the related transfer function

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  | StateSpace object

### Outputs

Type | Name | Description
--- | --- | ---
Integer | numberOfZeros |

### Modelica definition

    encapsulated function numberOfZeros
      "Calculate the number of zeros of the related transfer function"

        import Modelica_LinearSystems2.StateSpace;

      input StateSpace ss "StateSpace object";

      output Integer numberOfZeros;

    protected
      Integer nop;
      Integer noz;

    algorithm
      (nop,noz) := StateSpace.Internal.numberOfPolesAndZeros(ss);
      numberOfZeros := noz;
    end numberOfZeros;

## Modelica_LinearSystems2.StateSpace.Internal.complexZeros

Calculate the zeros of the related transfer function

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  | StateSpace object

### Outputs

Type | Name | Description
--- | --- | ---
Complex | zeros2[:] |

### Modelica definition

    encapsulated function complexZeros
      "Calculate the zeros of the related transfer function"

        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.Math.Complex;
        import Modelica_LinearSystems2.StateSpace;

      input StateSpace ss "StateSpace object";
      output Complex zeros2[:]=fill(Complex(0, 0),
          StateSpace.Internal.numberOfZeros(ss));

    protected
      Integer nx=size(ss.A, 2);

      Complex zeros[:];
      Complex poles[:];
      Real eval[nx,2];
      Real evec[nx,nx];

      Integer index[:]=fill(0, nx) "indices of zeros which are equal to poles";
      Integer i;
      Integer j;
      Integer k;
      Boolean h;
      Integer nzero;

    algorithm
      if size(ss.C, 1) <> 1 or size(ss.B, 2) <> 1 then
        assert(size(ss.B, 2) == 1,
          " function fromStateSpaceSISO expects a SISO-system as input\n but the number of inputs is "
           + String(size(ss.B, 2)) + " instead of 1");
        assert(size(ss.C, 1) == 1,
          " function fromStateSpaceSISO expects a SISO-system as input\n but the number of outputs is "
           + String(size(ss.C, 1)) + " instead of 1");
      end if;

      zeros := StateSpace.Analysis.invariantZeros( ss);

      poles := Modelica_LinearSystems2.Math.Complex.eigenValues(ss.A);

      for i in 1:size(zeros, 1) loop
        for j in 1:size(poles, 1) loop
           if zeros[i] == poles[j] then
            h := false;
            k := 1;
            while ((k < i) and (not h)) loop
              h := if (j == index[k]) then true else false;
              k := k + 1;
            end while;
            index[i] := if h then 0 else j;
          end if;
        end for;
      end for;

      j := 0;
      for i in 1:size(zeros, 1) loop
        if index[i] == 0 then
          j := j + 1;
          zeros2[j] := zeros[i];
        end if;
      end for;

    end complexZeros;

## Modelica_LinearSystems2.StateSpace.Internal.numberOfPoles

Calculate the number of poles of the related transfer function

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  | StateSpace object

### Outputs

Type | Name | Description
--- | --- | ---
Integer | numberOfPoles |

### Modelica definition

    encapsulated function numberOfPoles
      "Calculate the number of poles of the related transfer function"

        import Modelica_LinearSystems2.StateSpace;

      input StateSpace ss "StateSpace object";

      output Integer numberOfPoles=StateSpace.Internal.numberOfPolesAndZeros(ss);

    algorithm
    end numberOfPoles;

## Modelica_LinearSystems2.StateSpace.Internal.numberOfPolesAndZeros

Calculate the number poles and of zeros of the related transfer function

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  | StateSpace object

### Outputs

Type | Name | Description
--- | --- | ---
Integer | numberOfPoles |
Integer | numberOfZeros |

### Modelica definition

    encapsulated function numberOfPolesAndZeros
      "Calculate the number poles and of zeros of the related transfer function"

        import Modelica;
        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.Math.Complex;
        import Modelica_LinearSystems2.ZerosAndPoles;
        import Modelica_LinearSystems2.StateSpace;
        import Modelica_LinearSystems2.Math.Matrices.LAPACK;

      input StateSpace ss "StateSpace object";
      output Integer numberOfPoles;
      output Integer numberOfZeros;

    protected
      Integer nx=size(ss.A, 2);
      Complex zeros[:];
      Complex zeros2[:]
        "Finite, invariant zeros of ss; size(Zeros,1) <= size(ss.A,1)";
      Complex poles[:];
      Complex poles2[:] "eigenvalues of ss";
      Real eval[nx,2];
      Real evec[nx,nx];

      Integer index[:]=fill(0, nx) "indices of zeros which are equal to poles";
      Integer i;
      Integer j;
      Integer k;
      Boolean h;
      Integer nzero;

    algorithm
      if size(ss.C, 1) <> 1 or size(ss.B, 2) <> 1 then
        assert(size(ss.B, 2) == 1,
          " function fromStateSpaceSISO expects a SISO-system as input\n but the number of inputs is "
           + String(size(ss.B, 2)) + " instead of 1");
        assert(size(ss.C, 1) == 1,
          " function fromStateSpaceSISO expects a SISO-system as input\n but the number of outputs is "
           + String(size(ss.C, 1)) + " instead of 1");
      end if;

      zeros := StateSpace.Analysis.invariantZeros(ss);
      zeros2 := zeros;

      poles := Modelica_LinearSystems2.Math.Complex.eigenValues(ss.A);
      poles2 := poles;

    //Reduce terms which are in nominator as well as in denominator
      for i in 1:size(zeros, 1) loop
        for j in 1:size(poles, 1) loop
          if zeros[i] == poles[j] then
            h := false;
            k := 1;
            while ((k < i) and (not h)) loop
              h := if (j == index[k]) then true else false;
              k := k + 1;
            end while;
            index[i] := if h then 0 else j;
          end if;
        end for;
      end for;

      j := 0;
      for i in 1:size(zeros, 1) loop
        if index[i] == 0 then
          j := j + 1;
          zeros2[j] := zeros[i];
        end if;
      end for;
      nzero := j;
      j := 0;
      for i in 1:size(poles, 1) loop
        h := false;
        k := 1;
        while (k <= size(zeros, 1) and (not h)) loop
          h := if i == index[k] then true else false;
          k := k + 1;
        end while;
        if not h then
          j := j + 1;
          poles2[j] := poles[i];

        end if;
      end for;

      numberOfPoles := nx - size(zeros, 1) + nzero;
      numberOfZeros := nzero;

    end numberOfPolesAndZeros;

## Modelica_LinearSystems2.StateSpace.Internal.complexPoles

Generate a zeros-and-poles representation from state space representation

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  | StateSpace object

### Outputs

Type | Name | Description
--- | --- | ---
Complex | poles2[:] |

### Modelica definition

    encapsulated function complexPoles
      "Generate a zeros-and-poles representation from state space representation"

        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.Math.Complex;
        import Modelica_LinearSystems2.ZerosAndPoles;
        import Modelica_LinearSystems2.StateSpace;

      input StateSpace ss "StateSpace object";

      output Complex poles2[:]=fill(Complex(0, 0),
          StateSpace.Internal.numberOfPoles(ss));

    protected
      Integer nx=size(ss.A, 2);

      Complex zeros[:]
        "Finite, invariant zeros of ss; size(Zeros,1) <= size(ss.A,1)";
      Complex zeros2[:];
      Complex poles[:] "eigenvalues of ss";
      Real eval[nx,2];
      Real evec[nx,nx];

      Integer index[:]=fill(0, nx) "indices of zeros which are equal to poles";
      Integer i;
      Integer j;
      Integer k;
      Boolean h;
      Integer nzero;

    algorithm
      if size(ss.C, 1) <> 1 or size(ss.B, 2) <> 1 then
        assert(size(ss.B, 2) == 1,
          " function fromStateSpaceSISO expects a SISO-system as input\n but the number of inputs is "
           + String(size(ss.B, 2)) + " instead of 1");
        assert(size(ss.C, 1) == 1,
          " function fromStateSpaceSISO expects a SISO-system as input\n but the number of outputs is "
           + String(size(ss.C, 1)) + " instead of 1");
      end if;

      h := false;
      zeros := StateSpace.Analysis.invariantZeros( ss);
      zeros2 := zeros;

      poles := Modelica_LinearSystems2.Math.Complex.eigenValues(ss.A);

      for i in 1:size(zeros, 1) loop
        for j in 1:size(poles, 1) loop
          if zeros[i] == poles[j] then
            h := false;
            k := 1;
            while ((k < i) and (not h)) loop
              h := if (j == index[k]) then true else false;
              k := k + 1;
            end while;
            index[i] := if h then 0 else j;
          end if;
        end for;
      end for;

      j := 0;
      for i in 1:size(zeros, 1) loop
        if index[i] == 0 then
          j := j + 1;
          zeros2[j] := zeros[i];
        end if;
      end for;
      nzero := j;
      j := 0;
      for i in 1:size(poles, 1) loop
        h := false;
        k := 1;
        while (k <= size(zeros, 1) and (not h)) loop
          h := if i == index[k] then true else false;
          k := k + 1;
        end while;
        if not h then
          j := j + 1;
          poles2[j] := poles[i];

        end if;
      end for;

    end complexPoles;

## Modelica_LinearSystems2.StateSpace.Internal.trianUpperHess

Triangulize an upper Hessenberg matrix by repeatedly applicated householder reflexion

### Information

This function computes an triangular matrix from an upper Hessenberg matrix by stepwise annihilation of the subdiagonal elements.

>     A -> QA = T

It is assumend that the original matrix has upper hessenberg form. Additionally the vector b is transformed in the same way

>     b -> Qb = q

The function is primarily used to calculate the transfer function gain from a SISO state space system in observer Hessenberg form

>          ( *   *   ...   ...    * )          ( * )
>          ( *   *   ...   ...    * )          ( . )
>      A = ( 0   *   ...   ...    * ),    b =  ( . ),   c = ( 0, ..., 0, * )
>          ( .   .    .     .     . )          ( * )
>          ( 0  ...   0     *     * )          ( * )

If A is upper Hessenberg and T = Q*A is triangular then obviously H(s) = Q*(s*I -A) = s*I - T.

Further on, if T is triangular then also H = sI - T is and the element l_nn of L = inv(H) is given by 1/h_nn. The frequency response G(s0)for a given s0 that is neither zero nor pole of the system can be calculated by

>                         -1               -1     -1         -1           -1  -1        -1
>     G(s0)  = c*(s0*I -A)  *b = c*(s0*I -A)  *Q*Q  *b = c*(Q  *(s0*I -A))  *Q  *b = c*H(s0)*q

and because only the n'th element of c is different to zero the gain k is given by

>         q_nn*c_nn     product(s0 - poles_i)
>     k = ---------- * ----------------------
>            h_nn       product(s0 - zeros_i)

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | H[:, :] |  |
Real | b[size(H, 1)] |  |

### Outputs

Type | Name | Description
--- | --- | ---
Real | Ht[size(H, 1), size(H, 2)] |
Real | bt[size(b, 1)] |

### Modelica definition

    encapsulated function trianUpperHess
      "Triangulize an upper Hessenberg matrix by repeatedly applicated householder reflexion"
        import Modelica;
        import Modelica_LinearSystems2.Math.Matrices;
        import Modelica_LinearSystems2.Math.Vectors;

      input Real H[:,:];
      input Real b[size(H, 1)];

      output Real Ht[size(H, 1),size(H, 2)];
      output Real bt[size(b, 1)];

    protected
      Integer q=size(H, 1);
      Real u[:] "householder vector";
      Integer ll;

    algorithm
      Ht := H;
      bt := b;

      for ll in 1:q - 1 loop
        u := cat(
              1,
              zeros(ll - 1),
              cat(
                1,
                Vectors.householderVector(vector(Ht[ll:ll + 1, ll]), {1,0}),
                zeros(q - ll - 1)));
        Ht := Matrices.householderReflexion(Ht, u);
        bt := Vectors.householderReflexion(bt, u);
      end for;

    end trianUpperHess;

## Modelica_LinearSystems2.StateSpace.Internal.reducedCtrSystem2

calculate the controllable part of a SISO system

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  |
Real | eps | 0 |

### Outputs

Type | Name | Description
--- | --- | ---
StateSpace | ssm1 | controllable state space system

### Modelica definition

    encapsulated function reducedCtrSystem2
      "calculate the controllable part of a SISO system"

        import Modelica;
        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.StateSpace;
        import Modelica_LinearSystems2.Math.Matrices;
        import Modelica_LinearSystems2.Math.Vectors;

      input StateSpace ss;
      input Real eps=0;

    protected
      StateSpace sst(
        redeclare Real A[size(ss.A, 1),size(ss.A, 2)],
        redeclare Real B[size(ss.B, 1),size(ss.B, 2)],
        redeclare Real C[size(ss.C, 1),size(ss.C, 2)],
        redeclare Real D[size(ss.D, 1),size(ss.D, 2)])
        "tranformed state space system";

      Integer nx=size(ss.A, 1);
      Real Ah1[size(ss.A, 1),size(ss.A, 2)];
      Real bh1[size(ss.A, 1)];
      Real ch1[size(ss.A, 1)];

      Real u[:] "householder vector";
      Real cpoles[:,2]=Modelica_LinearSystems2.StateSpace.Internal.controllablePoles(ss);
      Integer rankQc=size(cpoles,1);

      Integer rankQc2;
      Real Qc2[nx,nx];
      Real sigma[:];
      Real eps2;

      Real Ah2[rankQc,rankQc];
      Real bh2[rankQc];
      Real ch2[rankQc];

      Integer ll;
      Integer r;

      Boolean h;

    public
      output StateSpace ssm1(
        redeclare Real A[rankQc,rankQc],
        redeclare Real B[rankQc,1],
        redeclare Real C[1,rankQc],
        redeclare Real D[size(ss.D, 1),size(ss.D, 2)])
        "controllable state space system";

    algorithm
      if size(ss.C, 1) <> 1 or size(ss.B, 2) <> 1 then
        assert(size(ss.B, 2) == 1,
          "A SISO-system is expected as input\n but the number of inputs is "
           + String(size(ss.B, 2)) + " instead of 1");
        assert(size(ss.C, 1) == 1,
          " A SISO-system is expected as input\n but the number of outputs is "
           + String(size(ss.C, 1)) + " instead of 1");
      end if;

      Ah1 := ss.A;
      bh1 := ss.B[:, 1];
      ch1 := ss.C[1, :];

      if nx > 1 then

        u := Vectors.householderVector(bh1, cat(
              1,
              fill(0, nx - 1),
              {1}));

        Ah1 := Matrices.householderSimilarityTransformation(Ah1, u);

        bh1 := Vectors.householderReflexion(bh1, u);
        ch1 := Vectors.householderReflexion(ch1, u);

        ll := nx;

        h := true;
        r := 1;

        while r <= nx - 2 and h loop
          if max(Ah1[1:ll - 1, ll]) <= 1e-8 then

            u := cat(
                  1,
                  Vectors.householderVector(Ah1[1:ll - 1, ll], cat(
                    1,
                    fill(0, ll - 2),
                    {1})),
                  fill(0, nx - ll + 1));
            Ah1 := Matrices.householderSimilarityTransformation(Ah1, u);
            ch1 := Vectors.householderReflexion(ch1, u);

            Ah1[1:ll - 2, ll] := fill(0, ll - 2);

            ll := ll - 1;

          else
            h := false;

          end if;
          r := r + 1;
        end while;
      end if;

      Qc2 := cat(
            2,
            Ah1[:, 2:nx],
            matrix(bh1));
      sigma := Modelica.Math.Matrices.singularValues(Qc2);
      eps2 := if eps > 0 then eps else 1000*sigma[1]*Modelica.Constants.eps;
      rankQc2 := 0;

      for i in 1:nx loop
        Modelica.Utilities.Streams.print(" s[" + String(i) + "] = " + String(
          sigma[i]));
        if sigma[i] > eps2 then
          rankQc2 := rankQc2 + 1;
        end if;
      end for;

      Modelica.Utilities.Streams.print("rankQc = " + String(rankQc) +
        "     rankQc2 = " + String(rankQc2) + "     eps2 = " + String(eps2));

      sst := StateSpace(
            A=Ah1,
            B=matrix(bh1),
            C=transpose(matrix(ch1)),
            D=ss.D);

      Ah2 := Ah1[nx - rankQc2 + 1:nx, nx - rankQc2 + 1:nx];
      bh2 := bh1[nx - rankQc2 + 1:nx];
      ch2 := ch1[nx - rankQc2 + 1:nx];
      ssm1 := StateSpace(
            A=Ah2,
            B=matrix(bh2),
            C=transpose(matrix(ch2)),
            D=ss.D);

    end reducedCtrSystem2;

## Modelica_LinearSystems2.StateSpace.Internal.controllablePoles

Compute the controllable and uncontrollable poles of a state space system

### Information

The function uses the SVD based staircase algorithm to transform the state space representation into a similar state space to separate the uncontrollable poles from the controllable poles.

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss | StateSpace(A=[-1], B=[1], C=... |

### Outputs

Type | Name | Description
--- | --- | ---
Real | cPoles[:, 2] | controllable poles
Real | ncPoles[:, 2] | uncontrollable poles
Real | poles[size(ss.A, 1), 2] | controllable and uncontrollable poles

### Modelica definition

    encapsulated function controllablePoles
      "Compute the controllable and uncontrollable poles of a state space system"

        import Modelica;
        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.StateSpace;
        import Modelica_LinearSystems2.Internal;

      input StateSpace ss=StateSpace(
          A=[-1],
          B=[1],
          C=[0],
          D=[0]);

      output Real cPoles[:,2] "controllable poles";
      output Real ncPoles[:,2] "uncontrollable poles";
      output Real poles[size(ss.A, 1),2] "controllable and uncontrollable poles";
    protected
      Modelica_LinearSystems2.Internal.StateSpaceR ssch(
        redeclare Real A[size(ss.A, 1),size(ss.A, 2)],
        redeclare Real B[size(ss.B, 1),size(ss.B, 2)],
        redeclare Real C[size(ss.C, 1),size(ss.C, 2)],
        redeclare Real D[size(ss.D, 1),size(ss.D, 2)])
        "upper block controller Hessenberg form state space system";
      Boolean isControllable;

    algorithm
      if size(ss.B, 2) == 0 then
        poles := Modelica.Math.Matrices.eigenValues(ss.A);
        ncPoles := poles;
        cPoles := fill(0, 0, 2);
      else
    // build upper Hessenberg staircase to decomposite controllable/uncontrollable subspaces
    // The controllable part of A is in A[1:ssch.r, 1:ssch.r]
        (isControllable,ssch) := StateSpace.Internal.staircaseSVD(ss);
        if isControllable then
          poles := Modelica.Math.Matrices.eigenValues(ss.A);
          cPoles := poles;
          ncPoles := fill(0, 0, 2);
        else
          cPoles := Modelica.Math.Matrices.eigenValues(ssch.A[1:ssch.r, 1:ssch.r]);
          ncPoles := Modelica.Math.Matrices.eigenValues(ssch.A[ssch.r + 1:size(ss.A,
            1), ssch.r + 1:size(ss.A, 1)]);
          poles := [cPoles; ncPoles];
        end if;
      end if;

    end controllablePoles;

## Modelica_LinearSystems2.StateSpace.Internal.polesAndZeros

Generate poles and zeros from state space representation

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  | StateSpace object
StateSpace | ssm | Modelica_LinearSystems2.Stat... |

### Outputs

Type | Name | Description
--- | --- | ---
PolesAndZeros | pz |

### Modelica definition

    encapsulated function polesAndZeros
      "Generate poles and zeros from state space representation"

        import Modelica;
        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.StateSpace;
        import Modelica_LinearSystems2.Internal.PolesAndZeros;
        import Modelica_LinearSystems2.Math.Complex;
        import Modelica_LinearSystems2.Internal;

      input StateSpace ss "StateSpace object";
      input StateSpace ssm=Modelica_LinearSystems2.StateSpace.Transformation.toIrreducibleForm(
                                                                           ss);
      output Internal.PolesAndZeros pz(
        redeclare Real p_real[size(ssm.A, 1)],
        redeclare Real p_im[size(ssm.A, 1)],
        redeclare Real z_real[size(StateSpace.Analysis.invariantZeros(ssm), 1)],
        redeclare Real z_im[size(StateSpace.Analysis.invariantZeros(ssm), 1)]);
    protected
      Complex poles[:]=Complex.eigenValues(ssm.A);
      Complex zeros[:]=StateSpace.Analysis.invariantZeros( ssm);

    algorithm
      pz.p_real := poles[:].re;
      pz.p_im := poles[:].im;
      pz.z_real := zeros[:].re;
      pz.z_im := zeros[:].im;
      pz.norz_p := Internal.numberOfRealZeros(poles);
      pz.norz_z := Internal.numberOfRealZeros(zeros);

    end polesAndZeros;

## Modelica_LinearSystems2.StateSpace.Internal.scaleFactor1

Return scale factor for first order block

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | n |  | (s+n)/(s+d)
Real | d |  | (s+n)/(s+d)
Real | small | 100*Modelica.Constants.eps |

### Outputs

Type | Name | Description
--- | --- | ---
Real | k | = d/n, if d,n are not zero, otherwise special cases

### Modelica definition

    encapsulated function scaleFactor1
      "Return scale factor for first order block"
        import Modelica;
      input Real n "(s+n)/(s+d)";
      input Real d "(s+n)/(s+d)";
      input Real small=100*Modelica.Constants.eps;
      output Real k "= d/n, if d,n are not zero, otherwise special cases";
    algorithm
      k := (if abs(d) > small then abs(d) else 1)/(if abs(n) > small then abs(n) else
              1);
    end scaleFactor1;

## Modelica_LinearSystems2.StateSpace.Internal.scaleFactor2

Return scale factor for second order block

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | n1 |  | (s^2 + n1*s + n2)/(s^2 + d1*s + d2)
Real | n2 |  | (s^2 + n1*s + n2)/(s^2 + d1*s + d2)
Real | d1 |  | (s^2 + n1*s + n2)/(s^2 + d1*s + d2)
Real | d2 |  | (s^2 + n1*s + n2)/(s^2 + d1*s + d2)
Real | small | 100*Modelica.Constants.eps |

### Outputs

Type | Name | Description
--- | --- | ---
Real | k | = d2/n2, if d2,n2 are not zero, otherwise special cases

### Modelica definition

    function scaleFactor2 "Return scale factor for second order block"
        import Modelica;
      input Real n1 "(s^2 + n1*s + n2)/(s^2 + d1*s + d2)";
      input Real n2 "(s^2 + n1*s + n2)/(s^2 + d1*s + d2)";
      input Real d1 "(s^2 + n1*s + n2)/(s^2 + d1*s + d2)";
      input Real d2 "(s^2 + n1*s + n2)/(s^2 + d1*s + d2)";
      input Real small=100*Modelica.Constants.eps;
      output Real k "= d2/n2, if d2,n2 are not zero, otherwise special cases";
    algorithm
      k := (if abs(d2) > small then abs(d2) else (if abs(d1) > small then abs(
        d1) else 1))/(if abs(n2) > small then abs(n2) else (if abs(n1) > small then
              abs(n1) else 1));
    end scaleFactor2;

## Modelica_LinearSystems2.StateSpace.Internal.invariantZerosHessenberg

Fast version to calculate the system zeros of a SISO system with D=[0] and A has upper Hessenberg form, delivered by StateSpace.reduceSystem

### Information

Computes the invariant zeros of a system in state space form:

       der(x) = A*x + B*u
            y = C*x + D*u

The invariant zeros of this system are defined as the variables z that make the following matrix singular:

        | A B |     | I 0 |
        |     | - z*|     |
        | C D |     | 0 0 |

where I is the identity matrix of the same size as A and 0 are zero matrices of appropriate dimensions.

Currently, there is the restriction that the number of inputs and the number of outputs must be identical.

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  | Linear system in state space form

### Outputs

Type | Name | Description
--- | --- | ---
Complex | Zeros[:] | Finite, invariant zeros of ss; size(Zeros,1) <= size(ss.A,1)

### Modelica definition

    encapsulated function invariantZerosHessenberg
      "Fast version to calculate the system zeros of a SISO system with D=[0] and A has upper Hessenberg form, delivered by StateSpace.reduceSystem"
        import Modelica;
        import Modelica_LinearSystems2.StateSpace;
        import Modelica_LinearSystems2.Math.Complex;
        import Modelica_LinearSystems2.Math.Matrices.LAPACK;

      input StateSpace ss "Linear system in state space form";
      output Complex Zeros[:]
        "Finite, invariant zeros of ss; size(Zeros,1) <= size(ss.A,1)";

    protected
      Integer nx=size(ss.A, 1) "Number of states";
      Integer nu=size(ss.B, 2) "Number of inputs";
      Integer ny=size(ss.C, 1) "Number of outputs";
      Real Ah[nx,nx]=ss.A;
      Real eps=100*Modelica.Constants.eps;
      Integer k;
      Boolean h;

    algorithm
      assert(nu == 1, "Function invariantZeros2 requires currently a SISO-input system.\n
    This condition is however not fulfilled because the number of inputs is nu = "
         + String(nu));
      assert(ny == 1, "Function invariantZeros2 requires currently a SISO-input system.\n
    This condition is however not fulfilled because the number of outputs is ny = "
         + String(ny));

      h := true;
      k := nx + 1;

      if size(ss.B, 2) > 0 then
        while k >= 1 and h loop
          k := k - 1;
          if abs(ss.B[k, 1]) >= eps then

            h := false;
          end if;

        end while;

        Zeros := fill(Complex(0), k - 1);

        if k > 1 then
          Ah[:, k - 1] := ss.A[:, k - 1] - ss.A[k, k - 1]/ss.B[k, 1]*ss.B[:, 1];

    //    Zeros := Complex.eigenValues(Ah[1:k - 1, 1:k - 1]);
          Zeros := Complex.Internal.eigenValues_dhseqr(Ah[1:k - 1, 1:k - 1]);

          for i in 1:k - 1 loop
            if Complex.'abs'(Zeros[i]) < Modelica.Math.Matrices.norm(Ah[1:k - 1, 1:
                k - 1], p=1)*1e-12 then
              Zeros[i] := Complex(0);
            end if;

          end for;

        end if;
      else
        Zeros := fill(Complex(0, 0), 0);
      end if;

    end invariantZerosHessenberg;

## Modelica_LinearSystems2.StateSpace.Internal.cntrHessenberg

calculate the controllable part of a SISO system

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  |

### Outputs

Type | Name | Description
--- | --- | ---
StateSpaceR | ssm1 | controllable state space system

### Modelica definition

    encapsulated function cntrHessenberg
      "calculate the controllable part of a SISO system"

        import Modelica;
        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.StateSpace;
        import Modelica_LinearSystems2.Math.Matrices;
        import Modelica_LinearSystems2.Math.Vectors;
        import Modelica_LinearSystems2.Math.Complex;

      input StateSpace ss;

      output Modelica_LinearSystems2.Internal.StateSpaceR ssm1(
        redeclare Real A[size(ss.A, 1),size(ss.A, 2)],
        redeclare Real B[size(ss.B, 1),1],
        redeclare Real C[1,size(ss.C, 2)],
        redeclare Real D[size(ss.D, 1),size(ss.D, 2)])
        "controllable state space system";

    protected
      Integer nx=size(ss.A, 1);
      Real Ah1[nx,nx];
      Real bh1[nx];
      Real ch1[nx];
      Real u[:] "householder vector";
      Real Q[nx,nx];
      Real V[size(ss.A, 1),size(ss.A, 2)];
      Real tau[nx - 1];
      Real Qc[:,:];
      Real svd[:];
      Real normA=Modelica.Math.Matrices.norm(A=ss.A, p=1);
      Integer rankMinSys;
      Boolean isZero=false;

    algorithm
      if size(ss.C, 1) <> 1 or size(ss.B, 2) <> 1 then
        assert(size(ss.B, 2) == 1,
          "A SISO-system is expected as input\n but the number of inputs is "
           + String(size(ss.B, 2)) + " instead of 1");
        assert(size(ss.C, 1) == 1,
          " A SISO-system is expected as input\n but the number of outputs is "
           + String(size(ss.C, 1)) + " instead of 1");
      end if;

      Ah1 := ss.A;
      bh1 := ss.B[:, 1];
      ch1 := ss.C[1, :];

      if Modelica.Math.Vectors.length(bh1) > 0 then

        if nx > 1 then

            // transform b->Qb = {*,0,...,0} and c->cQ, A->QAQ
          u := Vectors.householderVector(bh1, cat(
                1,
                {1},
                fill(0, nx - 1)));  //householder vector to compute a housholder reflector S = I - 2*u*u'/u'*u
          Ah1 := Matrices.householderSimilarityTransformation(Ah1, u);
          bh1 := Vectors.householderReflexion_e1(bh1, u);
          ch1 := Vectors.householderReflexion(ch1, u);

          (Ah1,V,tau) := Matrices.toUpperHessenberg(
                Ah1,
                1,
                nx);
          Q := Matrices.orthogonalQ(
                V,
                tau,
                1,
                nx);
          ch1 := ch1*Q;

        end if;

        rankMinSys := 1;
        while rankMinSys < nx and not isZero loop
          isZero := abs(Ah1[rankMinSys + 1, rankMinSys]) < normA*1e-10;
          rankMinSys := rankMinSys + 1;
        end while;

        ssm1 := Modelica_LinearSystems2.Internal.StateSpaceR(
              A=Ah1,
              B=matrix(bh1),
              C=transpose(matrix(ch1)),
              D=ss.D,
              r=if isZero then rankMinSys - 1 else rankMinSys);

      end if;

    //equation

    //algorithm
    end cntrHessenberg;

## Modelica_LinearSystems2.StateSpace.Internal.transposeStateSpace

Return the transposed state space system

### Information

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  |

### Outputs

Type | Name | Description
--- | --- | ---
StateSpace | sst |

### Modelica definition

    encapsulated function transposeStateSpace
      "Return the transposed state space system"

        import Modelica_LinearSystems2.StateSpace;

      input StateSpace ss;

      output StateSpace sst=StateSpace(
              A=transpose(ss.A),
              B=transpose(ss.C),
              C=transpose(ss.B),
              D=transpose(ss.D),
              uNames=ss.yNames,
              yNames=ss.uNames);
    algorithm

    end transposeStateSpace;

## Modelica_LinearSystems2.StateSpace.Internal.reduceRosenbrock

Algorithm to compress the generalized system matrix [A, B; C, D] to calculate the invariant zeros of a system

### Information

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | A[:, :] |  |
Real | B[:, :] |  |
Real | C[:, :] |  |
Real | D[:, :] |  |

### Outputs

Type | Name | Description
--- | --- | ---
Real | Ar[:, :] |
Real | Br[:, :] |
Real | Cr[:, :] |
Real | Dr[:, :] |
Integer | n |
Integer | m |
Integer | p |

### Modelica definition

    encapsulated function reduceRosenbrock
      "Algorithm to compress the generalized system matrix [A, B; C, D] to calculate the invariant zeros of a system"

        import Modelica;
        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.StateSpace;
        import Modelica_LinearSystems2.Math.Matrices;
        import Modelica_LinearSystems2.Math.Vectors;

        import Modelica.Utilities.Streams.print;

      input Real A[:,:];
      input Real B[:,:];
      input Real C[:,:];
      input Real D[:,:];

      output Real Ar[:,:];
      output Real Br[:,:];
      output Real Cr[:,:];
      output Real Dr[:,:];
      output Integer n;
      output Integer m;
      output Integer p;

    protected
      Real A2[:,:];
      Real B2[:,:];
      Real C2[:,:];
      Real CC[:,:];
      Real Co[:,:];
      Real Cu[:,:];
      Real D2[:,:];
      Real DD[:,:];
      Real Mr[:,:];
      Real Vf[:,:];
      Real V[:,:];
      Real V2[:,:];
      Real R[:,:];
      Real tau[:];

      Integer nx=size(A, 1);
      Integer nu=size(B, 2);
      Integer ny=size(C, 1);

      Integer nue;
      Integer delta;
      Integer rho;
      Integer mue;
      Integer sigma;
      Integer j;
      Boolean stop;
      Boolean stop1 "reduction finished";
      Boolean stop2 "system has no zeros";
      Integer rankR;
      Real normA=Modelica.Math.Matrices.norm(A=A, p=1);
      Real eps=normA*1e-10;

    algorithm
      if nx > 0 then

        A2 := A;
        B2 := B;
        C2 := C;
        D2 := D;
        stop := false;
        stop1 := false;
        stop2 := false;
        nue := nx;
        delta := 0;
        mue := ny;
        sigma := ny;
        j := 1;

        while not stop loop
          (V,R,tau,V2) := Matrices.QR( D2);

          rankR := 0;
    //  !!!! rank has to be determined. In the case of ill conditioned systems svd should be used
          for i in 1:min(size(R, 1), size(R, 2)) loop
            if abs(R[i, size(R, 2) - min(size(R, 1), size(R, 2)) + i]) > eps then
              rankR := rankR + 1;
            end if;
          end for;

    //rankR:=Modelica.Math.Matrices.rank(R);

          DD := R[1:rankR, :];

          CC := transpose(V2)*C2;

          sigma := rankR;
          stop1 := size(CC,1) == rankR;

          if not stop1 then
            Cu := CC[sigma + 1:end, :];
            Co := CC[1:sigma, :];

            (V,R,tau,V2) := Matrices.QR( Matrices.fliplr(transpose(Cu)));
             Vf:=Matrices.fliplr(V2);

             rankR := 0;
    //  !!!! rank determination
            for i in 1:min(size(R, 1), size(R, 2)) loop
              if abs(R[i, size(R, 2) - min(size(R, 1), size(R, 2)) + i]) > eps then
                rankR := rankR + 1;
              end if;
            end for;
    //rankR:=Modelica.Math.Matrices.rank(R);

            rho := rankR;
            stop1 := rho == 0;
            stop2 := size(Cu,2) == rankR;

            if not stop1 and not stop2 then
              nue := size(Cu, 2) - rankR;
              mue := rho + sigma;
              delta := delta + rho;

              if sigma == 0 then
                Mr := [transpose(Vf)*A2*Vf,transpose(Vf)*B2];
              else
                Mr := [transpose(Vf)*A2*Vf,transpose(Vf)*B2; Co*Vf,DD];
              end if;

              A2 := Mr[1:nue, 1:nue];
              B2 := Mr[1:nue, nue + rho + 1:nue + rho + nu];
              C2 := Mr[nue + 1:nue + mue,1:nue];
              D2 := Mr[nue + 1:nue + mue,nue + rho + 1:nue + rho + nu];

              j := j + 1;
            end if;
           //not stop1 or not stop2

          end if;
           //if not stop1

          stop := stop1 or stop2 or j>3*nx;

        end while;

        if stop1 then
          Ar := A2;
          Br := B2;
          Cr := C2;
          Dr := D2;
          n := nue;
          p := sigma;
          m := nu;
        else
          n := 0;
          p := 0;
          m := 0;
          Ar := fill(0,0,0);
          Br := fill(0,0,0);
          Cr := fill(0,0,0);
          Dr := fill(0,0,0);
        end if;

      else
        n := 0;
        p := 0;
        m := 0;
        A2 := A;
        B2 := B;
        C2 := C;
        D2 := D;
      end if;

    end reduceRosenbrock;

## Modelica_LinearSystems2.StateSpace.Internal.reducedCtrSystem

calculate the controllable part of a SISO system

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  |

### Outputs

Type | Name | Description
--- | --- | ---
StateSpaceR | ssm1 | controllable state space system

### Modelica definition

    encapsulated function reducedCtrSystem
      "calculate the controllable part of a SISO system"

        import Modelica;
        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.StateSpace;
        import Modelica_LinearSystems2.Math.Matrices;
        import Modelica_LinearSystems2.Math.Vectors;
        import Modelica_LinearSystems2.Math.Complex;

      input StateSpace ss;

      output Modelica_LinearSystems2.Internal.StateSpaceR ssm1(
        redeclare Real A[size(ss.A, 1),size(ss.A, 2)],
        redeclare Real B[size(ss.B, 1),1],
        redeclare Real C[1,size(ss.C, 2)],
        redeclare Real D[size(ss.D, 1),size(ss.D, 2)])
        "controllable state space system";

    protected
      Integer nx=size(ss.A, 1);
      Real Ah1[nx,size(ss.A, 2)];
      Real bh1[nx];
      Real ch1[nx];
      Real u[:] "householder vector";
      Integer ll;
      Integer r=1;
      Real maxa;
      Real normA=Modelica.Math.Matrices.norm(A=ss.A, p=1);

    algorithm
      if size(ss.C, 1) <> 1 or size(ss.B, 2) <> 1 then
        assert(size(ss.B, 2) == 1,
          "A SISO-system is expected as input\n but the number of inputs is "
           + String(size(ss.B, 2)) + " instead of 1");
        assert(size(ss.C, 1) == 1,
          " A SISO-system is expected as input\n but the number of outputs is "
           + String(size(ss.C, 1)) + " instead of 1");
      end if;

      Ah1 := ss.A;
      bh1 := ss.B[:, 1];
      ch1 := ss.C[1, :];

      if Modelica.Math.Vectors.length(bh1) > 0 then
        r := 1;
        if nx > 1 then

          u := Vectors.householderVector(bh1, cat(
                1,
                fill(0, nx - 1),
                {1}));

          Ah1 := Matrices.householderSimilarityTransformation(Ah1, u);
          bh1 := Vectors.householderReflexion_en(bh1, u);
          ch1 := Vectors.householderReflexion(ch1, u);
          bh1[1:nx - 1] := fill(0, nx - 1);

          ll := nx;
          maxa := max(abs(Ah1[1:ll - 1, ll]));

          while r <= nx - 1 and maxa > nx*normA*1e-10 loop
            u := cat(
                  1,
                  Vectors.householderVector(Ah1[1:ll - 1, ll], cat(
                    1,
                    fill(0, ll - 2),
                    {1})),
                  fill(0, nx - ll + 1));

            Ah1 := Matrices.Internal.hohoTrafoLowerHess(
                  Ah1,
                  u,
                  r);

            ch1 := Vectors.householderReflexion(ch1, u);
            ll := ll - 1;
            maxa := max(abs(Ah1[1:ll - 1, ll]));

            r := r + 1;
          end while;

        end if;

        ssm1 := Modelica_LinearSystems2.Internal.StateSpaceR(
              A=Ah1,
              B=matrix(bh1),
              C=transpose(matrix(ch1)),
              D=ss.D,
              r=r);
        ssm1.r := r;
      end if;
    end reducedCtrSystem;

## Modelica_LinearSystems2.StateSpace.Internal.reducedCtrSystemX

calculate the controllable part of a SISO system

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  |

### Outputs

Type | Name | Description
--- | --- | ---
StateSpaceR | ssm1 | controllable state space system

### Modelica definition

    encapsulated function reducedCtrSystemX
      "calculate the controllable part of a SISO system"

        import Modelica;
        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.StateSpace;
        import Modelica_LinearSystems2.Math.Matrices;
        import Modelica_LinearSystems2.Math.Vectors;
        import Modelica_LinearSystems2.Math.Complex;

      input StateSpace ss;

      output Modelica_LinearSystems2.Internal.StateSpaceR ssm1(
        redeclare Real A[size(ss.A, 1),size(ss.A, 2)],
        redeclare Real B[size(ss.B, 1),1],
        redeclare Real C[1,size(ss.C, 2)],
        redeclare Real D[size(ss.D, 1),size(ss.D, 2)])
        "controllable state space system";

    protected
      Integer nx=size(ss.A, 1);
      Real Ah1[nx,size(ss.A, 2)];
      Real bh1[nx];
      Real ch1[nx];
      Real u[:] "householder vector";
      Integer ll;
      Integer r=1;
      Real maxa;
      Real normA=Modelica.Math.Matrices.norm(A=ss.A, p=1);

    algorithm
      if size(ss.C, 1) <> 1 or size(ss.B, 2) <> 1 then
        assert(size(ss.B, 2) == 1,
          "A SISO-system is expected as input\n but the number of inputs is "
           + String(size(ss.B, 2)) + " instead of 1");
        assert(size(ss.C, 1) == 1,
          " A SISO-system is expected as input\n but the number of outputs is "
           + String(size(ss.C, 1)) + " instead of 1");
      end if;

      Ah1 := ss.A;
      bh1 := ss.B[:, 1];
      ch1 := ss.C[1, :];

      if Modelica.Math.Vectors.length(bh1) > 0 then
        r := 1;
        if nx > 1 then

          u := Vectors.householderVector(bh1, cat(
                1,
                fill(0, nx - 1),
                {1}));

          Ah1 := Matrices.householderSimilarityTransformation(Ah1, u);
          bh1 := Vectors.householderReflexion_en(bh1, u);
          ch1 := Vectors.householderReflexion(ch1, u);
          bh1[1:nx - 1] := fill(0, nx - 1);

          ll := nx;
          maxa := max(abs(Ah1[1:ll - 1, ll]));

          while r <= nx - 1 and maxa > nx*normA*1e-5 loop
            u := cat(
                  1,
                  Vectors.householderVector(Ah1[1:ll - 1, ll], cat(
                    1,
                    fill(0, ll - 2),
                    {1})),
                  fill(0, nx - ll + 1));

            Ah1 := Matrices.Internal.hohoTrafoLowerHess(
                  Ah1,
                  u,
                  r);

            ch1 := Vectors.householderReflexion(ch1, u);
            ll := ll - 1;
            maxa := max(abs(Ah1[1:ll - 1, ll]));

            r := r + 1;
          end while;

        end if;

        ssm1 := Modelica_LinearSystems2.Internal.StateSpaceR(
              A=Ah1,
              B=matrix(bh1),
              C=transpose(matrix(ch1)),
              D=ss.D,
              r=r);
        ssm1.r := r;
      end if;
    end reducedCtrSystemX;

## Modelica_LinearSystems2.StateSpace.Internal.assignPolesMI2

Pole assigment design algorithm for multi input systems

### Information

#### Syntax

(K, S, po, nfp, nap, nup) | = | StateSpace.Design.assignPolesMI(ss, gamma, np, alpha, tol)

#### Description

The purpose of this function is to determine the state feedback matrix K for a given time invariant multi input state system (A,B) such that the closed-loop state matrix A-B*K has specified eigenvalues. The feedback matrix K is calculated by factorization following [1]. The algorithm modifies the eigenvalues sequentially and also allows partial eigenvalue assignment.
At the beginning of the algorithm, the feedback matrix K is set to zero (K = 0) and the matrix A is reduced to an ordered real Schur form by separating its spectrum in two parts

>                   | F1  F3|
>      F = Q*A*Q' = |       |
>                   | 0   F2|

in such a way, that F1 contains the eigenvalues that will be retained and F3 contains the eigenvalues going to be modified. On the suggestion of [1] the eigenvalues evr to be retained are chosen as

>       evr = {s in C: Re(s) < -alpha, alpha > =0}

but other specification are conceivable of course.
Let

>       G = [G1;G2] = Q*B

with an appropriate partition according to F2. (F2, G2) has to be controllable.
If the feedback matrix K is taken in a form

>      K = [0, K2]

the special structure of F and K results in a closed loop state matrix

>               |F1 F3 - G1*K2|
>     F - G*K = |             |
>               |0  F2 - G2*K2|

with only the eigenvalues of F2 are modified. This approach to modify separated eigenvalues is used to sequentially shift one real eigenvalue ore two complex conjugated eigenvalues stepwise until all assigned eigenvalues are placed. Therefore, at each step i always the (two) lower right eingenvalue(s) are modified by an appropriate feedback matrix Ki. The matrix F - G*Ki remains in real Schur form. The assigned eigenvalue(s) is (are) then moved to another diagonal position of the real Schur form using reordering techniques F < -- Qi*F*Qi' and a new block is transferred to the lower right diagonal position. The transformations are accumulated in Qi and are also applicated to the matrices

>      G < - Qi*G Q < - Qi*Q

The eigenvalue(s) to be assigned at each step is (are) chosen such that the norm of each Ki is minimized [1].

#### Example

>        Modelica_LinearSystems2.StateSpace ss=Modelica_LinearSystems2.StateSpace(
>           A=[-1, 1, 1;0, 1, 1;0, 0, 1],
>           B=[0; 0; 1],
>           C=[0, 1, 0],
>           D=[0]);

>        Real Q[3,3];

>     algorithm
>       Q := Modelica_LinearSystems2.StateSpace.Analysis.observabilityMatrix(ss);
>     // Q = [0, 1, 0; 0, 1, 1; 1, 1, 2]

#### References

[1] | Varga A. | "A Schur method for pole assignment" | IEEE Trans. Autom. Control, Vol. AC-26, pp. 517-519, 1981

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
StateSpace | ss |  | state space system
Complex | gamma[:] | fill(Complex(0), size(ss.A, ... | Designed Poles
Integer | np | size(gamma, 1) | number of given eigenvalues to assign
Boolean | giveEigenvalues | false | Boolean parameter to display the order of the eigenvalues
Real | alpha | -1e10 | maximum admissible value for real parts(continuous) or for moduli (discrete) of the eigenvalues of A which will not be modified by the eigenvalue assignment algorithm
Real | tolerance | Modelica.Math.Matrices.norm(... | The tolerance to be used in determining the controllability of (A,B)

### Outputs

Type | Name | Description
--- | --- | ---
Real | K[size(ss.B, 2), size(ss.A, 1)] | State feedback matrix assigning the desired poles
Real | S[:, :] | Closed loop System matrix
Complex | po[size(ss.A, 1)] | poles of the closed loop system
Integer | nfp | number of eigenvalues that are not modified with respect to alpha
Integer | nap | number of assigned eigenvalues
Integer | nup | number of uncontrollable eigenvalues
Complex | X[size(ss.A, 1), size(ss.A, 1)] | eigenvectors of the closed loop system

### Modelica definition

    encapsulated function assignPolesMI2
      "Pole assigment design algorithm for multi input systems"

          import Modelica_LinearSystems2;
          import Modelica_LinearSystems2.Math.Complex;
          import Modelica_LinearSystems2.StateSpace;
          import Modelica;
          import Modelica.Utilities.Streams.print;
          import Modelica_LinearSystems2.TransferFunction;
          import Modelica_LinearSystems2.Math.Matrices;

      input StateSpace ss "state space system";

      input Complex gamma[:]=fill(Complex(0), size(ss.A,1)) "Designed Poles";
      input Integer np=size(gamma, 1) "number of given eigenvalues to assign";
      input Boolean giveEigenvalues=false
        "Boolean parameter to display the order of the eigenvalues";
      input Real alpha=-1e10
        "maximum admissible value for real parts(continuous) or for moduli (discrete) of the eigenvalues of A which will not be modified by the eigenvalue assignment algorithm";
      input Real tolerance=Modelica.Math.Matrices.norm(ss.A, 1)*1e-12
        "The tolerance to be used in determining the controllability of (A,B)";
      output Real K[size(ss.B, 2),size(ss.A, 1)]
        "State feedback matrix assigning the desired poles";
      output Real S[:,:] "Closed loop System matrix";
      output Complex po[size(ss.A, 1)] "poles of the closed loop system";
      output Integer nfp
        "number of eigenvalues that are not modified with respect to alpha";
      output Integer nap "number of assigned eigenvalues";
      output Integer nup "number of uncontrollable eigenvalues";
      output Complex X[size(ss.A, 1),size(ss.A, 1)]
        "eigenvectors of the closed loop system";

    protected
      Real A_rsf[size(ss.A, 1),size(ss.A, 2)];
      Real B_rsf[size(ss.B, 1),size(ss.B, 2)];
      Real Q[size(ss.A, 1),size(ss.A, 1)];
      Real Ks1[:,:];
      Real Ks2[:,:];
      Real Q2[:,:];
      Real A_rsf_1[:,:];
      Real Q1[:,:];
      Boolean select[:];
      Boolean rselectA[:];
      Real Z[:,:] "orthogonal transformation matrix";
      Real ZT[:,:] "orthogonal transformation matrix";
      Complex pf[:];
      Complex gammaReordered[:]=gamma;
      Integer info;
      Real wr[size(gamma, 1)];
      Real wi[size(gamma, 1)];
      Boolean imag=false;
      Integer i;
      Integer ii;
      Integer counter;
      Integer n=size(ss.A, 1);
      Integer ipf;
      Integer nccA "number of conjugated complex pole pairs of unmodified system";
      Integer nccg "number of conjugated complex pole pairs of gamma";
      Integer ncc;

      Complex h;
      Integer h2[2];
      Real alphaReal[size(ss.A, 1)] "Real part of eigenvalue=alphaReal+i*alphaImag";
      Real alphaImag[size(ss.A, 1)]
        "Imaginary part of eigenvalue=(alphaReal+i*alphaImag";

      Complex pi[:]=Complex.eigenValues(ss.A);

      Boolean complex_assignedPoles=false
        "true if there is at least one conjugated comples pole pair in the set of the assigned poles";
      Boolean complex_originalPoles=false
        "true if there is at least one conjugated comples pole pair in the set of unmodified system poles";
      Boolean consistency;

      Complex SS[:,:];
      Complex Xj[:,:];

      Integer markA[n]=fill(1,n);
      Integer markg[n]=fill(1,n);

      Complex ev[:];

    algorithm
      if giveEigenvalues then
        (,,alphaReal,alphaImag,info) :=
          Modelica_LinearSystems2.Math.Matrices.LAPACK.dgees(ss.A);
        assert(info == 0, "The output info of LAPACK.dgees should be zero, else if\n
         info < 0:  if info = -i, the i-th argument of dgees had an illegal value\n
         info > 0:  if INFO = i, and i is
                   <= N: the QR algorithm failed to compute all the
                         eigenvalues; elements 1:ILO-1 and i+1:N of WR and WI
                         contain those eigenvalues which have converged; if
                         JOBVS = 'V', VS contains the matrix which reduces A
                         to its partially converged Schur form.\n");
        for i in 1:n loop
          po[i].re := alphaReal[i];
          po[i].im := alphaImag[i];
        end for;
        Complex.Vectors.print("The eigenvalues of the open loop system are sorted to\n eigenvalues", pi);
      else
        assert(size(gamma, 1) <= size(ss.A, 1), "At most n (order of ss) eigenvalues can be assigned");

     /* Extraction of Poles (Variable conversation) and pole sequence check */
      for i in 1:size(gamma, 1) loop
        wr[i] := gamma[i].re;
        wi[i] := gamma[i].im;
        if imag then
          assert(wi[i - 1] == -wi[i] and wr[i - 1] == wr[i],
            "Poles are in wrong sequence");
          imag := false;
        elseif wi[i] <> 0 then
          imag := true;
        end if;
      end for;

      // put matrix ss.A to real Schur form A <- QAQ' and compute B <- QB
      (A_rsf,Z,alphaReal,alphaImag) := Matrices.rsf2(ss.A);
      ZT := transpose(Z);

      // determine number of poles not to be assigned according to alpha
         nfp := 0;
         nccA := 0;
         nccg := 0;
         for i in 1:n loop
           if alphaReal[i] < alpha then
             nfp := nfp + 1;
             markA[i] := 0;
             markg[i] := 0;
           end if;
           if abs(alphaImag[i]) > 0 then
             markA[i] := 2;
             nccA := nccA + 1;
           end if;
           if abs(gammaReordered[i].im) > 0 then
             markg[i] := 2;
             nccg := nccg + 1;
           end if;
         end for;
         nap := n - nfp;
         nccA := div(nccA, 2);
         nccg := div(nccg, 2);
         ncc := max(nccA, nccg);

       // reorder real Schur form according to alpha
       (A_rsf,Z,alphaReal,alphaImag) := Matrices.Internal.reorderRSF2(
           A_rsf,
           identity(size(A_rsf, 1)),
           alphaReal,
           alphaImag,
           alpha);
       ZT := transpose(Z)*ZT;
       B_rsf := ZT*ss.B;

        Modelica_LinearSystems2.Math.Vectors.printVector(alphaReal,6,"alphaReal");
        Modelica_LinearSystems2.Math.Vectors.printVector(alphaImag,6,"alphaImag");
       Modelica_LinearSystems2.Math.Complex.Vectors.print("gammaReordered1",gammaReordered);

      // Reorder gammaReordered according to alpha
        ii := 1;
        for i in 1:n loop
          if markg[i]==0 then
            h := gammaReordered[ii];
            gammaReordered[ii] := gammaReordered[i];
            gammaReordered[i] := h;
            ii := ii+1;
          end if;
        end for;

       // check consistency of poles assignment, i.e. complex pole pairs in gammaReordered and alphaReal and alphaImag may not be separated
        consistency:=true;
        i:=1;
        while i<n loop
          if markA[i]==2 then
            consistency :=markA[i + 1] == 2 and ((markg[i] == 2 and markg[i + 1] == 2)
              or (markg[i] == 1 and markg[i + 1] == 1));
            i := i+2;
          elseif markg[i]==2 then
            consistency :=markg[i + 1] == 2 and (markA[i] == 1 and markA[i + 1] == 1);
            i := i+2;
          else
            i := i+1;
          end if;
        end while;
       assert(consistency,"System poles and assigned poles have to be assigned consistently, i.e. complex pole pairs may not be separated");

      // main algorithm
      K := zeros(size(ss.B, 2), size(ss.A, 1));
      counter := nfp + 1;
       while counter<=n loop
        if markA[n + nfp + 1 - counter] == 2 or markg[n + nfp + 1 - counter] == 2 then

    //#############################
    Modelica_LinearSystems2.Math.Complex.Vectors.print("g2",gammaReordered[n + nfp - counter:n + nfp + 1 - counter]);
    Modelica_LinearSystems2.Math.Vectors.printVector(alphaReal[n + nfp - counter:n + nfp + 1 - counter],6,"ar2");
    ev:=Modelica_LinearSystems2.Math.Complex.eigenValues( A_rsf[n - 1:n, n - 1:n]);
    Modelica_LinearSystems2.Math.Complex.Vectors.print("ev2",ev);
    //#############################

          Ks2 := StateSpace.Internal.assignOneOrTwoPoles(
            A_rsf[n - 1:n, n - 1:n],
            matrix(B_rsf[n - 1:n, :]),
            gammaReordered[n + nfp - counter:n + nfp +
            1 - counter],
            tolerance);
          K := K + [zeros(size(Ks2, 1), size(K, 2) - 2),Ks2]*ZT;
          A_rsf := A_rsf - B_rsf*[zeros(size(Ks2, 1), size(K, 2) - 2),Ks2];
          select := fill(false, n - counter + 1);
          select[n - counter:n - counter + 1] := {true,true};

          (A_rsf[counter:n, counter:n],Q2) := Matrices.LAPACK.dtrsen(
            "E",
            "V",
            select,
            A_rsf[counter:n, counter:n],
            identity(n - counter + 1));  //The Schur vector matrix is identity, since A_rsf already has Schur form

          A_rsf[1:counter - 1, counter:n] := A_rsf[1:counter - 1, counter:n]*Q2;
          B_rsf[counter:n, :] := transpose(Q2)*B_rsf[counter:n, :];
          ZT[counter:n, :] := transpose(Q2)*ZT[counter:n, :];
    //       h2 := markA[counter:counter + 1];
    //       markA[counter:counter + 1] := markA[n - 1:n];
    //       markA[n - 1:n] := h2;
    //       h2 := markg[counter:counter + 1];
    //       markg[counter:counter + 1] := markg[n - 1:n];
    //       markg[n - 1:n] := h2;
    //       h := gammaReordered[n - 1];
    //       gammaReordered[n - 1] := gammaReordered[counter];
    //       gammaReordered[counter] := gammaReordered[n - 1];
    //       h := gammaReordered[n];
    //       gammaReordered[n] := gammaReordered[counter + 1];
    //       gammaReordered[counter + 1] := gammaReordered[n];
          counter := counter + 2;
        else
    //#############################
    print("g1 = "+String(gammaReordered[n + nfp + 1 - counter]));
    print("ar1 = "+String(alphaReal[n + nfp + 1 - counter]));
    print("ev1 = "+String(A_rsf[n,n])+"\n");
    //#############################

          Ks1 := StateSpace.Internal.assignOneOrTwoPoles(
            matrix(A_rsf[n, n]),
            transpose(matrix(B_rsf[n, :])),
            {gammaReordered[n + nfp + 1 - counter]},
            tolerance);

          K := K + [zeros(size(Ks1, 1), size(K, 2) - 1),Ks1]*ZT;
          A_rsf := A_rsf - B_rsf*[zeros(size(Ks1, 1), size(K, 2) - 1),Ks1];
          select := fill(false, n - counter + 1);
          select[n - counter + 1] := true;

          (A_rsf[counter:n, counter:n],Q1) := Matrices.LAPACK.dtrsen(
            "E",
            "V",
            select,
            A_rsf[counter:n, counter:n],
            identity(n - counter + 1)); //The Schur vector matrix is identity, since A_rsf already has Schur form

          A_rsf[1:counter - 1, counter:n] := A_rsf[1:counter - 1, counter:n]*Q1;
          B_rsf[counter:n, :] := transpose(Q1)*B_rsf[counter:n, :];
          ZT[counter:n, :] := transpose(Q1)*ZT[counter:n, :];
    //       h2[1] := markA[counter];
    //       markA[counter] := markA[n];
    //       markA[n] := h2[1];
    //       h2[1] := markg[counter];
    //       markg[counter] := markg[n];
    //       markg[n] := h2[1];
    //       h := gammaReordered[n];
    //       gammaReordered[n] := gammaReordered[counter];
    //       gammaReordered[counter] := gammaReordered[n];
          counter := counter + 1;
        end if;

        end while;

      S := ss.A - ss.B*K;
      po := Complex.eigenValues(S);
       X := Complex.eigenVectors(S);

    //    X := fill(Complex(0),n,n);
    //    for i in 1:n loop
    //      SS:=Complex(1)*S;
    //      for ii in 1:n loop
    //        SS[ii,ii] := SS[ii,ii]-po[i];
    //      end for;
    //      Xj := Modelica_LinearSystems2.WorkInProgress.Math.Matrices.C_nullspace(
    //                                 SS);
    //      // Matrices.printMatrix(Complex.real(Xj),6,"ReXj");
    //      // Matrices.printMatrix(Complex.imag(Xj),6,"ImXj");
    //      for ii in 1:n loop
    //        X[ii,i] := Xj[ii,1];
    //      end for;
    //    end for;

    end if;

    end assignPolesMI2;

## Modelica_LinearSystems2.ZerosAndPoles.Internal

Internal library of record Filter (should not be directly used by user)

### Information

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

### Package Content

Name | Description
--- | ---
ZerosAndPoles | Continuous zeros and poles description of a single input, single output system (data + operations)
filter | Generate the data record of a ZerosAndPoles transfer function from a filter description
findMatrixName |
fromFile_pc | Generate a zeros and poles data record by reading the polynomial coefficients from a file (default file name is pc.mat)
fromFile_zp | Generate a zeros and poles data record by reading poles and zeros from a file (default file name is zp.mat)
checkRepresentation | checks whether the system on file is represented by zeros and poles (z, p) or first and second order polynomials (n1, n2, d1, d2)
sameVectorElements | Determine identical elements of a vector
fromRealAndImag | Generate a complex vector from a real part vector and imaginary part vector
sameMatrixRows | Determine identical rows of a [:,2] matrix
firstOrderToString | Transform vector of coefficients of normalized first order polynomials to a string representation
secondOrderToString | Transform vector of coefficients of normalized second order polynomials to a string representation
normalizationFactor | Compute correction factor of low pass filter such that amplitude at cut-off frequency is 1/sqrt(2) (= 3db)
BesselCoefficients | Return coefficients of normalized low pass Bessel filter (= gain at cut-off frequency 1 rad/s is decreased 3dB
'p+a' |
'p^2+k[1]*p+k[2]' |
roots | Determine zeros of factorized polynomial
frequencyRange | Determine min. and max. resonance frequencies
frequencyRangeBode | Determine min. and max. frequencies for Bode plot
filterToNormalized | Given [p^2,p] and [p] coefficients, transform to normalized form with highest power of p equal 1
numberOfRealZeros | Calculate number of real zeros
numberOfRealZeros2 | Calculate number of real zeros
numberOfRealPoles | Calculate number of real poles
isRoot | check if frequency is an element of the complex vector zeros
numberOfRealZerosAndPoles_zp | Get the number of first oder polynomials (n1, d1) and second order polynomials (n2, d2) of zeros and poles from zeros and poles written in a MAT-file
numberOfRealZerosAndPoles_pc | Generate a zeros and poles data record by reading the polynomial coefficients from a file (default file name is zp.mat)
isControllableAndObservableSISO | To check whether a SISO system is controllable and observable

## Modelica_LinearSystems2.ZerosAndPoles.Internal.ZerosAndPoles

Continuous zeros and poles description of a single input, single output system (data + operations)

### Information

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

### Parameters

Type | Name | Default | Description
--- | --- | --- | ---
y = k*(product(p+n1[i]) * product(p^2+n2[i,1]*p+n2[i,2])) / (product(p+d1[i])*product(p^2+d2[i,1]*p+d2[i,2])) *u
Real | k | 1.0 | Multiplicative factor of transfer function
Real | n1[:] |  | [p^0] coefficients of 1st order numerator polynomials
Real | n2[:, 2] |  | [p,p^0] coefficients of 2nd order numerator polynomials
Real | d1[:] |  | [p^0] coefficients of 1st order denominator polynomials
Real | d2[:, 2] |  | [p,p^0] coefficients of 2nd order denominator polynomials

### Modelica definition

    record ZerosAndPoles
      "Continuous zeros and poles description of a single input, single output system (data + operations)"
      extends Modelica.Icons.Record;

      Real k=1.0 "Multiplicative factor of transfer function";
      Real n1[:] "[p^0] coefficients of 1st order numerator polynomials";
      Real n2[:,2] "[p,p^0] coefficients of 2nd order numerator polynomials";
      Real d1[:] "[p^0] coefficients of 1st order denominator polynomials";
      Real d2[:,2] "[p,p^0] coefficients of 2nd order denominator polynomials";

    end ZerosAndPoles;

## Modelica_LinearSystems2.ZerosAndPoles.Internal.filter

Generate the data record of a ZerosAndPoles transfer function from a filter description

### Information

#### Syntax

filterFunction | = | Modelica_LinearSystems2.ZerosAndPoles.Design.filter(analogFilter, filterType, order, f_cut, gain, A_ripple, normalized)

#### Description

This function constructs a ZerosAndPoles transfer function description of low and high pass filters. Typical frequency responses for the 4 supported low pass filter types are shown in the next figure (this figure was generated with function ZerosAndPoles.Examples.bodePlot_filter2):

The step responses of the same low pass filters are shown in the next figure, starting from a steady state initial filter with initial input = 0.2:

Obviously, the frequency responses give a somewhat wrong impression of the filter characteristics: Although Butterworth and Chebyshev filters have a significantly steeper magnitude as the CriticalDamping and Bessel filters, the step responses of the latter ones are much better. This means for example, that a CriticalDamping or a Bessel filter should be selected, if a filter is mainly used to make a non-linear inverse model realizable.

Typical frequency responses for the 4 supported high pass filter types are shown in the next figure:

The corresponding step responses of these high pass filters are shown in the next figure:

All filters are available in normalized (default) and non-normalized form. In the normalized form, the amplitude of the filter transfer function at the cutoff frequency is 1/sqrt(2) (= 3 dB). Note, when comparing the filters of this function with other software systems, the setting of "normalized" has to be selected appropriately. For example, the signal processing toolbox of Matlab provides the filters in non-normalized form and therefore a comparision makes only sense, if normalized = false is set.

#### Example

>        Types.AnalogFilter analogFilter=Types.AnalogFilter.CriticalDamping;
>        Integer order=2;
>        Modelica.SIunits.Frequency f_cut=10;

>        ZerosAndPoles zp_filter;

>     algorithm
>         zp_filter=Modelica_LinearSystems2.ZerosAndPoles.Design.filter(
>           order=order,
>           f_cut=f_cut,
>           analogFilter=analogFilter);

>     // zp_filter = 9530.93/( (p + 97.6265)^2 )

#### References

[1] | Tietze U., and Schenk Ch. | "Halbleiter-Schaltungstechnik" | Springer Verlag, 12. Auflage, pp. 815-852, 2002.

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
AnalogFilter | analogFilter | Types.AnalogFilter.CriticalD... | Analog filter characteristics (CriticalDamping/Bessel/Butterworth/Chebyshev)
FilterType | filterType | Types.FilterType.LowPass | Type of filter (LowPass/HighPass)
Integer | order | 2 | Order of filter
Frequency | f_cut | 1/(2*Modelica.Constants.pi) | Cut-off frequency (default is w_cut = 1 rad/s) [Hz]
Real | gain | 1.0 | Gain (= amplitude of frequency response at zero frequency)
Real | A_ripple | 0.5 | Pass band ripple for Chebyshev filter (otherwise not used) [dB]
Boolean | normalized | true | = true, if amplitude of low pass filter at f_cut where the magnitude of the filter sagged for 3dB, otherwise unmodified filter

### Outputs

Type | Name | Description
--- | --- | ---
ZerosAndPoles | filter | Filter transfer function

### Modelica definition

    encapsulated function filter
      "Generate the data record of a ZerosAndPoles transfer function from a filter description"

        import Modelica;
        import Modelica.Utilities.Streams;
        import Modelica.Utilities.Streams.print;
        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.Types;
        import Modelica_LinearSystems2.ZerosAndPoles;

        import Modelica_LinearSystems2.ZerosAndPoles.Internal;
        import Modelica.Math.*;
        import Modelica_LinearSystems2.Math.*;

      input Modelica_LinearSystems2.Types.AnalogFilter analogFilter=Types.AnalogFilter.CriticalDamping
        "Analog filter characteristics (CriticalDamping/Bessel/Butterworth/Chebyshev)";
      input Modelica_LinearSystems2.Types.FilterType filterType=Types.FilterType.LowPass
        "Type of filter (LowPass/HighPass)";
      input Integer order(min=1) = 2 "Order of filter";
      input Modelica.SIunits.Frequency f_cut=1/(2*Modelica.Constants.pi)
        "Cut-off frequency (default is w_cut = 1 rad/s)";
      input Real gain=1.0
        "Gain (= amplitude of frequency response at zero frequency)";
      input Real A_ripple(unit="dB") = 0.5
        "Pass band ripple for Chebyshev filter (otherwise not used)";
      input Boolean normalized=true
        "= true, if amplitude of low pass filter at f_cut where the magnitude of the filter sagged for 3dB, otherwise unmodified filter";
      output Internal.ZerosAndPoles filter(
        redeclare Real n1[if filterType == Types.FilterType.LowPass then 0 else (
          if analogFilter == Types.AnalogFilter.CriticalDamping then order else
          mod(order, 2))],
        redeclare Real n2[if filterType == Types.FilterType.LowPass then 0 else (
          if analogFilter == Types.AnalogFilter.CriticalDamping then 0 else
          integer((order + 1)/2)),2],
        redeclare Real d1[if analogFilter == Types.AnalogFilter.CriticalDamping then
                order else mod(order, 2)],
        redeclare Real d2[if analogFilter == Types.AnalogFilter.CriticalDamping then
                0 else integer(order/2),2]) "Filter transfer function";

    protected
      Integer n_num1=size(filter.n1, 1);
      Integer n_num2=size(filter.n2, 1);
      Integer n_den1=size(filter.d1, 1);
      Integer n_den2=size(filter.d2, 1);
      Integer n_num=n_num1 + 2*n_num2;
      Integer n_den=n_den1 + 2*n_den2;
      Real pi=Modelica.Constants.pi;
      Boolean evenOrder=mod(order, 2) == 0
        "= true, if even filter order (otherwise uneven)";
      Modelica.SIunits.AngularVelocity w_cut=2*pi*f_cut "Cut-off angular frequency";
      Real w_cut2 "= w_cut*w_cut";
      Real alpha=1.0 "Frequency correction factor";
      Real alpha2 "= alpha*alpha";

      Real alphax;

      Real poleReal "Real part of pole";
      Real poleImag "Imaginary part of pole";
      Real epsilon "Ripple size";
      Real fac "arsinh(epsilon)";
      Real A2 "poleReal^2 + poleImag^2";
      Real A "Amplitude at w_cut";
      Real num1[n_num1]
        "[p] coefficients of numerator first order polynomials (a*p + 1)";
      Real num2[n_num2,2]
        "[p^2, p] coefficients of numerator second order polynomials (b*p^2 + a*p + 1I)";
      Real den1[n_den1]
        "[p] coefficients of denominator first order polynomials (a*p + 1)";
      Real den2[n_den2,2]
        "[p^2, p] coefficients of denominator second order polynomials (b*p^2 + a*p + 1I)";
      Real aux;
      Real k;
    algorithm
      // Set properties that are common for all filters
      filter.k := gain;

      /* Compute filter coefficients of prototype low pass filter. If another filter
         characteristics is desired (e.g. high pass filter), it is derived
         from the low pass filter coefficients below
      */
      if analogFilter == Types.AnalogFilter.CriticalDamping then
        if normalized then
          alpha := sqrt(2^(1/order) - 1);
    //alpha := sqrt(10^(3/10/order)-1)
        else
          alpha := 1;
        end if;
        for i in 1:n_den1 loop
          den1[i] := alpha;
        end for;

      elseif analogFilter == Types.AnalogFilter.Bessel then
        (den1,den2,alpha) := Internal.BesselCoefficients(order);
        if not normalized then
          alpha2 := alpha*alpha;
          for i in 1:n_den2 loop
            den2[i, 1] := den2[i, 1]*alpha2;
            den2[i, 2] := den2[i, 2]*alpha;
          end for;
          if not evenOrder then
            den1[1] := den1[1]*alpha;
          end if;
        end if;

      elseif analogFilter == Types.AnalogFilter.Butterworth then
         // Original filter is already normalized
        for i in 1:n_den2 loop
          den2[i, 1] := 1.0;
          den2[i, 2] := -2*cos(pi*(0.5 + (i - 0.5)/order));
        end for;
        if not evenOrder then
          den1[1] := 1.0;
        end if;

      elseif analogFilter == Types.AnalogFilter.Chebyshev then
        epsilon := sqrt(10^(A_ripple/10) - 1);
        fac := asinh(1/epsilon)/order;

        for i in 1:n_den2 loop
           // poles of prototype lowpass
          poleReal := cos(pi/2 + pi*(2*i - 1)/(2*order))*sinh(fac);
          poleImag := -sin(pi/2 + pi*(2*i - 1)/(2*order))*cosh(fac);

           // second order filter coefficients
          A2 := poleReal^2 + poleImag^2;
          den2[i, 1] := 1/A2;
          den2[i, 2] := -2*poleReal/A2;
        end for;
        if not evenOrder then
          den1[1] := -1.0/(cos(pi/2 + pi*(2*n_den - 1)/(2*order))*sinh(fac));
        end if;

         /* Transformation of filter transfer function with "new(p) = alpha*p"
            in order that the filter transfer function has an amplitude of
            1/sqrt(2) at the cutoff frequency
         */
        if normalized then
          alpha := Internal.normalizationFactor(den1, den2);
          alpha2 := alpha*alpha;
          for i in 1:n_den2 loop
            den2[i, 1] := den2[i, 1]*alpha2;
            den2[i, 2] := den2[i, 2]*alpha;
          end for;
          if not evenOrder then
            den1[1] := den1[1]*alpha;
          end if;
        end if;

      else
        Streams.error("analogFilter (= " + String(analogFilter) +
          ") is not supported");
      end if;

      // Compute amplitude at w=1
    /*
      A := 1.0;
      for i in 1:n_den2 loop
         A := A*(1 + den2[i,2]^2 - 2*den2[i,1]
                   + den2[i,1]^2);
      end for;
      for i in 1:n_den1 loop
         A := A*(1 + den1[i]^2);
      end for;
      A := 1/sqrt(A);
      print("A = " + String(A));
    */

      // Determine normalized denominator polynomials with highest power of p equal to one
      filter.n1 := zeros(n_num1);
      filter.n2 := zeros(n_num2, 2);
      (filter.d1,filter.d2,k) := Internal.filterToNormalized(den1, den2);
      filter.k := filter.k/k;

      // Compute desired filter characteristics from low pass filter coefficients
      if filterType == Types.FilterType.HighPass then
         /* The high pass filter is derived from the low pass filter by
            the transformation new(p) = 1/p
            1/(p^2 + a*p + b) -> 1/((1/p)^2 + a*(1/p) + b) = (1/b)*p^2 / (p^2 + (a/b)*p + 1/b)
            1/(p + a)         -> 1/((1/p) + a) = (1/a)*p / (p + (1/a))
         */
        assert(n_num1 == n_den1 and n_num2 == n_den2,
          "Internal error 1, should not occur");
        filter.n1 := zeros(n_num1);
        filter.n2 := zeros(n_num2, 2);
        for i in 1:n_num1 loop
          filter.k := filter.k/filter.d1[i];
          filter.d1[i] := 1/filter.d1[i];
        end for;
        for i in 1:n_num2 loop
          filter.k := filter.k/filter.d2[i, 2];
          filter.d2[i, 1] := filter.d2[i, 1]/filter.d2[i, 2];
          filter.d2[i, 2] := 1/filter.d2[i, 2];
        end for;
      end if;

      /* Change filter coefficients according to transformation new(p) = p/w_cut
         Numerator  :     (p/w)^2 + a*(p/w) + b = (1/w^2)*(p^2 + (a*w)*p + b*w^2)
                                      (p/w) + a = (1/w)*(p + w*a)
         Denominator: 1/((p/w)^2 + a*(p/w) + b) = w^2/(p^2 + (a*w)*p + w^2/b)
                                  1/((p/w) + a) = w/(p + w*a)
      */
      w_cut2 := w_cut*w_cut;
      filter.k := filter.k*w_cut^(n_den1 + 2*n_den2 - n_num1 - 2*n_num2);
      filter.n1 := w_cut*filter.n1;
      filter.d1 := w_cut*filter.d1;
      filter.n2 := [w_cut*filter.n2[:, 1],w_cut2*filter.n2[:, 2]];
      filter.d2 := [w_cut*filter.d2[:, 1],w_cut2*filter.d2[:, 2]];

    end filter;

## Modelica_LinearSystems2.ZerosAndPoles.Internal.findMatrixName

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
String | filename |  |
String | matName | "z" |

### Outputs

Type | Name | Description
--- | --- | ---
Integer | m |

### Modelica definition

    function findMatrixName
      input String filename;
      input String matName="z";

      output Integer m;

    external "C" m = findMatrixName(
            filename,
            matName,
            "NoClass");

    end findMatrixName;

## Modelica_LinearSystems2.ZerosAndPoles.Internal.fromFile_pc

Generate a zeros and poles data record by reading the polynomial coefficients from a file (default file name is pc.mat)

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
String | fileName | "pc.mat" | Name of the zeros and poles data file

### Outputs

Type | Name | Description
--- | --- | ---
ZerosAndPoles | zp |

### Modelica definition

    encapsulated function fromFile_pc
      "Generate a zeros and poles data record by reading the polynomial coefficients from a file (default file name is pc.mat)"
        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.ZerosAndPoles;
        import Modelica_LinearSystems2.Math.Complex;
        import Modelica;

      input String fileName="pc.mat" "Name of the zeros and poles data file";

    protected
      input Integer n1n2d1d2[4]=
          ZerosAndPoles.Internal.numberOfRealZerosAndPoles_pc(fileName);
      input Integer n1=n1n2d1d2[1];
      input Integer n2=n1n2d1d2[2];
      input Integer d1=n1n2d1d2[3];
      input Integer d2=n1n2d1d2[4];
      input Integer zSize=n1n2d1d2[1] + 2*n1n2d1d2[2];
      input Integer pSize=n1n2d1d2[3] + 2*n1n2d1d2[4];
    public
      output ZerosAndPoles zp(
        n1=fill(0, n1),
        n2=fill(
              0,
              n2,
              2),
        d1=fill(0, d1),
        d2=fill(
              0,
              d2,
              2));

    protected
      Integer n1_2=if n1 > 0 then 1 else 0 "second dimension of n1-matrix";
      Integer n2_2=if n2 > 0 then 2 else 0 "second dimension of n2-matrix";
      Integer d1_2=if d1 > 0 then 1 else 0 "second dimension of d1-matrix";
      Integer d2_2=if d2 > 0 then 2 else 0 "second dimension of d2-matrix";

      Real k=scalar(readMatrix(
            fileName,
            "k",
            1,
            1));
      Real n1Vector[n1]=vector(readMatrix(
            fileName,
            "n1",
            n1,
            n1_2)) "coefficients of first order numenator polynomials";
      Real n2Matrix[n2,n2_2]=readMatrix(
            fileName,
            "n2",
            n2,
            n2_2) "coefficients of second order denominator polynomials";
      Real d1Vector[d1]=vector(readMatrix(
            fileName,
            "d1",
            d2,
            d1_2)) "coefficients of first order denominator polynomials";
      Real d2Matrix[d2,d2_2]=readMatrix(
            fileName,
            "d2",
            d2,
            d2_2) "coefficients of second order numenator polynomials";

    algorithm
      zp.k := k;
      zp.n1 := if n1 > 0 then n1Vector else fill(0, 0);
      zp.n2 := if n2 > 0 then n2Matrix else fill(
          0,
          0,
          2);
      zp.d1 := if d1 > 0 then d1Vector else fill(0, 0);
      zp.d2 := if d2 > 0 then d2Matrix else fill(
          0,
          0,
          2);

    end fromFile_pc;

## Modelica_LinearSystems2.ZerosAndPoles.Internal.fromFile_zp

Generate a zeros and poles data record by reading poles and zeros from a file (default file name is zp.mat)

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
String | fileName | "zp.mat" | Name of the zeros and poles data file

### Outputs

Type | Name | Description
--- | --- | ---
ZerosAndPoles | zp |

### Modelica definition

    encapsulated function fromFile_zp
      "Generate a zeros and poles data record by reading poles and zeros from a file (default file name is zp.mat)"

        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.ZerosAndPoles;
        import Modelica_LinearSystems2.Math.Complex;

      input String fileName="zp.mat" "Name of the zeros and poles data file";
    protected
      input Integer n1n2d1d2[4]=
          ZerosAndPoles.Internal.numberOfRealZerosAndPoles_zp(fileName);
      input Integer n1=n1n2d1d2[1];
      input Integer n2=n1n2d1d2[2];
      input Integer d1=n1n2d1d2[3];
      input Integer d2=n1n2d1d2[4];
      input Integer zSize=n1n2d1d2[1] + 2*n1n2d1d2[2];
      input Integer pSize=n1n2d1d2[3] + 2*n1n2d1d2[4];
    public
      output ZerosAndPoles zp(
        n1=fill(0, n1),
        n2=fill(
              0,
              n2,
              2),
        d1=fill(0, d1),
        d2=fill(
              0,
              d2,
              2));

    protected
      Integer z_2=if zSize > 0 then 2 else 0 "second dimension of zeros-matrix";
      Integer p_2=if pSize > 0 then 2 else 0 "second dimension of poles-matrix";

      Real k=scalar(readMatrix(
            fileName,
            "k",
            1,
            1));
      Real zerosMatrix[zSize,z_2]=readMatrix(
            fileName,
            "z",
            zSize,
            z_2) "zeros in rows of real parts and imaginary parts";
      Real polesMatrix[pSize,p_2]=readMatrix(
            fileName,
            "p",
            pSize,
            p_2) "poles in rows of real parts and imaginary parts";
      Complex zeros[:]=if zSize > 0 then ZerosAndPoles.Internal.fromRealAndImag(
          zerosMatrix[:, 1], zerosMatrix[:, z_2]) else fill(Complex(0), 0);
      Complex poles[:]=if pSize > 0 then ZerosAndPoles.Internal.fromRealAndImag(
          polesMatrix[:, 1], polesMatrix[:, p_2]) else fill(Complex(0), 0);

    algorithm
      zp := ZerosAndPoles(
          k=k,
          z=zeros,
          p=poles);
    end fromFile_zp;

## Modelica_LinearSystems2.ZerosAndPoles.Internal.checkRepresentation

checks whether the system on file is represented by zeros and poles (z, p) or first and second order polynomials (n1, n2, d1, d2)

### Information

The function output is true if the system is given in zeros and poles representation. Therefore, it is assumend that the used array names are "z" and "p" or "n1, n2, d1" and "d2" respectively.

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
String | fileName | "zp.mat" | Name of the zeros and poles data file

### Outputs

Type | Name | Description
--- | --- | ---
Boolean | iszp |

### Modelica definition

    function checkRepresentation
      "checks whether the system on file is represented by zeros and poles (z, p) or first and second order polynomials (n1, n2, d1, d2)"
      import Modelica_LinearSystems2.ZerosAndPoles.Internal;
      import Modelica_LinearSystems2;
      input String fileName="zp.mat" "Name of the zeros and poles data file";
      output Boolean iszp=true;
    protected
      Integer m=0;

    algorithm
      m := Internal.findMatrixName(fileName, "z");
      m := m + Internal.findMatrixName(fileName, "p");
      iszp := m == 2;

    end checkRepresentation;

## Modelica_LinearSystems2.ZerosAndPoles.Internal.sameVectorElements

Determine identical elements of a vector

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | v[:] |  | Vector
Integer | startIndex | 1 | Start index

### Outputs

Type | Name | Description
--- | --- | ---
Integer | endIndex | startIndex:endIndex are identical elements in v

### Modelica definition

    function sameVectorElements
      "Determine identical elements of a vector"

      input Real v[:] "Vector";
      input Integer startIndex=1 "Start index";
      output Integer endIndex=startIndex
        "startIndex:endIndex are identical elements in v";
    protected
      Integer nv=size(v, 1);
      Integer i=startIndex + 1;
      Real v0=v[startIndex];
    algorithm
      while i <= nv loop
        if v[i] == v0 then
          endIndex := i;
          i := i + 1;
        else
          i := nv + 1;
        end if;
      end while;
    end sameVectorElements;

## Modelica_LinearSystems2.ZerosAndPoles.Internal.fromRealAndImag

Generate a complex vector from a real part vector and imaginary part vector

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | real[:] |  |
Real | imag[size(real, 1)] |  |

### Outputs

Type | Name | Description
--- | --- | ---
Complex | result[size(real, 1)] | Number of real zeros

### Modelica definition

    function fromRealAndImag
      "Generate a complex vector from a real part vector and imaginary part vector "

      import Modelica_LinearSystems2.Math.Complex;

      input Real real[:];
      input Real imag[size(real, 1)];
      output Complex result[size(real, 1)] "Number of real zeros";
    algorithm
      for i in 1:size(real, 1) loop
        result[i].re := real[i];
        result[i].im := imag[i];
      end for;
    end fromRealAndImag;

## Modelica_LinearSystems2.ZerosAndPoles.Internal.sameMatrixRows

Determine identical rows of a [:,2] matrix

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | M[:, 2] |  | Matrix
Integer | startIndex | 1 | Start index

### Outputs

Type | Name | Description
--- | --- | ---
Integer | endIndex | startIndex:endIndex are identical rows in M

### Modelica definition

    function sameMatrixRows "Determine identical rows of a [:,2] matrix"

      input Real M[:,2] "Matrix";
      input Integer startIndex=1 "Start index";
      output Integer endIndex=startIndex
        "startIndex:endIndex are identical rows in M";
    protected
      Integer n=size(M, 1);
      Integer i=startIndex + 1;
      Real v0[size(M, 2)]=M[startIndex, :];
    algorithm
      while i <= n loop
        if M[i, 1] == v0[1] and M[i, 2] == v0[2] then
          endIndex := i;
          i := i + 1;
        else
          i := n + 1;
        end if;
      end while;
    end sameMatrixRows;

## Modelica_LinearSystems2.ZerosAndPoles.Internal.firstOrderToString

Transform vector of coefficients of normalized first order polynomials to a string representation

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | c[:] |  | p^0 coefficients of normalized first order polynomials
Integer | significantDigits | 6 | Number of significant digits that are shown
String | name | "p" | Independent variable name used for printing

### Outputs

Type | Name | Description
--- | --- | ---
String | s |

### Modelica definition

    function firstOrderToString
      "Transform vector of coefficients of normalized first order polynomials to a string representation"
      import Modelica_LinearSystems2.Math.Vectors;

      input Real c[:] "p^0 coefficients of normalized first order polynomials";
      input Integer significantDigits=6
        "Number of significant digits that are shown";
      input String name="p" "Independent variable name used for printing";
      output String s="";
    protected
      Integer nc=size(c, 1);
      Real cs[nc]=Modelica.Math.Vectors.sort(c);
      Real cc[nc];
      Integer i;
      Integer j;
      Integer nj;
      Integer k;
      Boolean found=false;
    algorithm
      // Move coefficients with zero to the beginning
      i := 1;
      while i <= nc loop
        if cs[i] == 0 then
          found := true;
          k := i;
          i := nc + 1;
        else
          i := i + 1;
        end if;
      end while;

      if found then
        j := sameVectorElements(cs, k);
        cc := cat(
              1,
              cs[k:j],
              cs[1:k - 1],
              cs[j + 1:end]);
      else
        cc := cs;
      end if;

      // Transform coefficients to string
      i := 1;
      while i <= nc loop
        if i <> 1 then
          s := s + "*";
        end if;
        if cc[i] == 0 then
          s := s + name;
        else
          if cc[i] > 0 then
            s := s + "(" + name + " + ";
          else
            s := s + "(" + name;
          end if;
          s := s + String(cc[i], significantDigits=significantDigits) + ")";
        end if;
        j := sameVectorElements(cc, i);
        nj := j - i + 1;
        if nj > 1 then
          s := s + "^" + String(nj);
        end if;
        i := j + 1;
      end while;
    end firstOrderToString;

## Modelica_LinearSystems2.ZerosAndPoles.Internal.secondOrderToString

Transform vector of coefficients of normalized second order polynomials to a string representation

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | c[:, 2] |  | [p,p^0] coefficients of normalized second order polynomials
Integer | significantDigits | 6 | Number of significant digits that are shown
String | name | "p" | Independent variable name used for printing

### Outputs

Type | Name | Description
--- | --- | ---
String | s |

### Modelica definition

    function secondOrderToString
      "Transform vector of coefficients of normalized second order polynomials to a string representation"
      import Modelica_LinearSystems2.Math.Matrices;

      input Real c[:,2]
        "[p,p^0] coefficients of normalized second order polynomials";
      input Integer significantDigits=6
        "Number of significant digits that are shown";
      input String name="p" "Independent variable name used for printing";
      output String s="";
    protected
      Integer nc=size(c, 1);
      Real cc[nc,2]=Modelica.Math.Matrices.sort(c);
      Integer i=1;
      Integer j;
      Integer nj;
    algorithm
      while i <= nc loop
        if i <> 1 then
          s := s + "*";
        end if;
        j := sameMatrixRows(cc, i);
        nj := j - i + 1;
        if cc[i, 1] == 0 and cc[i, 2] == 0 then
            // case p^2
          s := s + name + "^" + String(2*nj);
        else
          s := s + "(" + name + "^2";

            // b*p term
          if cc[i, 1] == 1 then
            s := s + "+" + name;
          elseif cc[i, 1] == -1 then
            s := s + " - " + name;
          elseif cc[i, 1] <> 0 then
            if cc[i, 1] > 0 then
              s := s + " + ";
            end if;
            s := s + String(cc[i, 1], significantDigits=significantDigits) +
              "*" + name;
          end if;

            // a*p^0 term
          if cc[i, 2] > 0 then
            s := s + " + ";
          end if;
          s := s + String(cc[i, 2], significantDigits=significantDigits) + ")";
          if nj > 1 then
            s := s + "^" + String(nj);
          end if;
        end if;
        i := j + 1;
      end while;
    end secondOrderToString;

## Modelica_LinearSystems2.ZerosAndPoles.Internal.normalizationFactor

Compute correction factor of low pass filter such that amplitude at cut-off frequency is 1/sqrt(2) (= 3db)

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | c1[:] |  | [p] coefficients of denominator polynomials (a*p + 1)
Real | c2[:, 2] |  | [p^2, p] coefficients of denominator polynomials (b*p^2 + a*p + 1)

### Outputs

Type | Name | Description
--- | --- | ---
Real | alpha | Correction factor (replace p by alpha*p)

### Modelica definition

    function normalizationFactor
      "Compute correction factor of low pass filter such that amplitude at cut-off frequency is 1/sqrt(2) (= 3db)"
      import Modelica.Utilities.Streams;

      input Real c1[:] "[p] coefficients of denominator polynomials (a*p + 1)";
      input Real c2[:,2]
        "[p^2, p] coefficients of denominator polynomials (b*p^2 + a*p + 1)";
      output Real alpha "Correction factor (replace p by alpha*p)";
    protected
      Real oldAlpha=1.0;
      Integer counter=0;
      Integer counterMax=100;
      Boolean iterating=true;
      Real residue;
      Real der_residue;

      encapsulated function normalizationFactorResidues
        "Residue and residue derivative for Newton iteration of correction factor computation"
        input Real c1[:] "[p] coefficients of denominator polynomials (a*p + 1)";
        input Real c2[:,2]
          "[p^2, p] coefficients of denominator polynomials (b*p^2 + a*p + 1)";
        input Real alpha;
        output Real residue;
        output Real der_residue;
      protected
        Real beta=2;  //beta = 10^(3/10);
        Real cc1;
        Real cc2;
        Real p;
        Real pd;
        Real alpha2=alpha*alpha;
        Real alpha4=alpha2*alpha2;
        Real A2=1.0;
        Real fac=0.0;
      algorithm
        if size(c1, 1) == 1 then
          cc1 := c1[1]*c1[1];
          p := 1 + cc1*alpha2;
          pd := cc1;
          fac := fac + pd/p;
          A2 := A2*p;
        end if;
        for i in 1:size(c2, 1) loop
          cc1 := c2[i, 2]*c2[i, 2] - 2*c2[i, 1];
          cc2 := c2[i, 1]*c2[i, 1];
          p := 1 + cc1*alpha2 + cc2*alpha4;
          pd := cc1 + 2*cc2*alpha2;
          fac := fac + pd/p;
          A2 := A2*p;
        end for;
        residue := A2 - beta;
        der_residue := 2*alpha*A2*fac;
      end normalizationFactorResidues;
    algorithm
       // Simple Newton iteration
      while iterating loop
        counter := counter + 1;
        (residue,der_residue) := normalizationFactorResidues(
              c1,
              c2,
              oldAlpha);
        alpha := oldAlpha - residue/der_residue;
        if abs(alpha - oldAlpha) <= 100*Modelica.Constants.eps then
          iterating := false;
        elseif counter > counterMax then
          Streams.error("Numerical problems when computing cutoff frequency\n"
             + "correction factor of filter (use option normalized=false)");
        else
          oldAlpha := alpha;
        end if;
      end while;
      assert(alpha > 0.0,
        "Numerical problems when computing cutoff frequency\n" +
        "correction factor of filter (alpha = " + String(alpha) + ")\n" +
        "(use option normalized=false)");
    end normalizationFactor;

## Modelica_LinearSystems2.ZerosAndPoles.Internal.BesselCoefficients

Return coefficients of normalized low pass Bessel filter (= gain at cut-off frequency 1 rad/s is decreased 3dB

### Information

The transfer function H(p) of a n 'th order Bessel filter is given by

>              Bn(0)
>      H(p) = -------
>              Bn(p)

with the denominator polynomial

>               n             n  (2n - k)!       p^k
>      Bn(p) = sum c_k*p^k = sum ----------- * -------   (1)
>              k=0           k=0 (n - k)!k!    2^(n-k)

and the numerator

>                    (2n)!     1
>     Bn(0) = c_0 = ------- * ---- .                     (2)
>                     n!      2^n

Although the coefficients c_k are integer numbers, it is not advisable to use the polynomials in an unfactorized form because the coefficients are fast growing with order n (c_0 is approximately 0.3e24 and 0.8e59 for order n=20 and order n=40 respectively).
Therefore, the polynomial Bn(p) is factorized to first and second order polynomials with real coefficients corresponding to zeros and poles representation that is used in this library.

The function returns the coefficients which resulted from factorization of the normalized transfer function

>     H'(p') = H(p),  p' = p/w0

as well as

>     alpha = 1/w0

the reciprocal of the cut of frequency w0 where the gain of the transfer function is decreased 3dB.

Both, coefficients and cut off frequency were calculated symbolically and were eventually evaluated with high precision calculation. The results were stored in this function as real numbers.

Calculation of normalized Bessel filter coefficients
Equation

>        abs(H(j*w0)) = abs(Bn(0)/Bn(j*w0)) = 10^(-3/20)

which must be fulfilled for cut off frequency w = w0 leads to

>        [Re(Bn(j*w0))]^2 + [Im(Bn(j*w0))]^2 - (Bn(0)^2)*10^(3/10) = 0

which has exactly one real solution w0 for each order n. This solutions of w0 are calculated symbolically first and evaluated by using high precise values of the coefficients c_k calculated by following (1) and (2).
With w0, the coefficients of the factorized polynomial can be computed by calculating the zeros of the denominator polynomial

>              n
>      Bn(p) = sum w0^k*c_k*(p/w0)^k
>              k=0

of the normalized transfer function H'(p'). There exist n/2 of conjugate complex pairs of zeros (beta +-j*gamma) if n is even and one additional real zero (alpha) if n is odd. Finally, the coefficients a, b1_k, b2_k of the polynomials

>      a*p + 1,  n is odd

and

>      b2_k*p^2 + b1_k*p + 1,   k = 1,... div(n,2)

results from

>      a = -1/alpha

and

>      b2_k = 1/(beta_k^2 + gamma_k^2) b1_k = -2*beta_k/(beta_k^2 + gamma_k^2)

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Integer | order |  | Order of filter in the range 1..41

### Outputs

Type | Name | Description
--- | --- | ---
Real | c1[mod(order, 2)] | [p] coefficients of Bessel denominator polynomials (a*p + 1)
Real | c2[integer(order/2), 2] | [p^2, p] coefficients of Bessel denominator polynomials (b2*p^2 + b1*p + 1)
Real | alpha | Normalization factor

### Modelica definition

    function BesselCoefficients
      "Return coefficients of normalized low pass Bessel filter (= gain at cut-off frequency 1 rad/s is decreased 3dB"

      import Modelica.Utilities.Streams;
      import Modelica_LinearSystems2;
      input Integer order "Order of filter in the range 1..41";
      output Real c1[mod(order, 2)]
        "[p] coefficients of Bessel denominator polynomials (a*p + 1)";
      output Real c2[integer(order/2),2]
        "[p^2, p] coefficients of Bessel denominator polynomials (b2*p^2 + b1*p + 1)";
      output Real alpha "Normalization factor";
    algorithm
      if order == 1 then
        alpha := 1.002377293007601;
        c1[1] := 0.9976283451109835;
      elseif order == 2 then
        alpha := 0.7356641785819585;
        c2[1, 1] := 0.6159132201783791;
        c2[1, 2] := 1.359315879600889;
      elseif order == 3 then
        alpha := 0.5704770156982642;
        c1[1] := 0.7548574865985343;
        c2[1, 1] := 0.4756958028827457;
        c2[1, 2] := 0.9980615136104388;
      elseif order == 4 then
        alpha := 0.4737978580281427;
        c2[1, 1] := 0.4873729247240677;
        c2[1, 2] := 1.337564170455762;
        c2[2, 1] := 0.3877724315741958;
        c2[2, 2] := 0.7730405590839861;
      elseif order == 5 then
        alpha := 0.4126226974763408;
        c1[1] := 0.6645723262620757;
        c2[1, 1] := 0.4115231900614016;
        c2[1, 2] := 1.138349926728708;
        c2[2, 1] := 0.3234938702877912;
        c2[2, 2] := 0.6205992985771313;
      elseif order == 6 then
        alpha := 0.3705098000736233;
        c2[1, 1] := 0.3874508649098960;
        c2[1, 2] := 1.219740879520741;
        c2[2, 1] := 0.3493298843155746;
        c2[2, 2] := 0.9670265529381365;
        c2[3, 1] := 0.2747419229514599;
        c2[3, 2] := 0.5122165075105700;
      elseif order == 7 then
        alpha := 0.3393452623586350;
        c1[1] := 0.5927147125821412;
        c2[1, 1] := 0.3383379423919174;
        c2[1, 2] := 1.092630816438030;
        c2[2, 1] := 0.3001025788696046;
        c2[2, 2] := 0.8289928256598656;
        c2[3, 1] := 0.2372867471539579;
        c2[3, 2] := 0.4325128641920154;
      elseif order == 8 then
        alpha := 0.3150267393795002;
        c2[1, 1] := 0.3151115975207653;
        c2[1, 2] := 1.109403015460190;
        c2[2, 1] := 0.2969344839572762;
        c2[2, 2] := 0.9737455812222699;
        c2[3, 1] := 0.2612545921889538;
        c2[3, 2] := 0.7190394712068573;
        c2[4, 1] := 0.2080523342974281;
        c2[4, 2] := 0.3721456473047434;
      elseif order == 9 then
        alpha := 0.2953310177184124;
        c1[1] := 0.5377196679501422;
        c2[1, 1] := 0.2824689124281034;
        c2[1, 2] := 1.022646191567475;
        c2[2, 1] := 0.2626824161383468;
        c2[2, 2] := 0.8695626454762596;
        c2[3, 1] := 0.2302781917677917;
        c2[3, 2] := 0.6309047553448520;
        c2[4, 1] := 0.1847991729757028;
        c2[4, 2] := 0.3251978031287202;
      elseif order == 10 then
        alpha := 0.2789426890619463;
        c2[1, 1] := 0.2640769908255582;
        c2[1, 2] := 1.019788132875305;
        c2[2, 1] := 0.2540802639216947;
        c2[2, 2] := 0.9377020417760623;
        c2[3, 1] := 0.2343577229427963;
        c2[3, 2] := 0.7802229808216112;
        c2[4, 1] := 0.2052193139338624;
        c2[4, 2] := 0.5594176813008133;
        c2[5, 1] := 0.1659546953748916;
        c2[5, 2] := 0.2878349616233292;
      elseif order == 11 then
        alpha := 0.2650227766037203;
        c1[1] := 0.4950265498954191;
        c2[1, 1] := 0.2411858478546218;
        c2[1, 2] := 0.9567800996387417;
        c2[2, 1] := 0.2296849355380925;
        c2[2, 2] := 0.8592523717113126;
        c2[3, 1] := 0.2107851705677406;
        c2[3, 2] := 0.7040216048898129;
        c2[4, 1] := 0.1846461385164021;
        c2[4, 2] := 0.5006729207276717;
        c2[5, 1] := 0.1504217970817433;
        c2[5, 2] := 0.2575070491320295;
      elseif order == 12 then
        alpha := 0.2530051198547209;
        c2[1, 1] := 0.2268294941204543;
        c2[1, 2] := 0.9473116570034053;
        c2[2, 1] := 0.2207657387793729;
        c2[2, 2] := 0.8933728946287606;
        c2[3, 1] := 0.2087600700376653;
        c2[3, 2] := 0.7886236252756229;
        c2[4, 1] := 0.1909959101492760;
        c2[4, 2] := 0.6389263649257017;
        c2[5, 1] := 0.1675208146048472;
        c2[5, 2] := 0.4517847275162215;
        c2[6, 1] := 0.1374257286372761;
        c2[6, 2] := 0.2324699157474680;
      elseif order == 13 then
        alpha := 0.2424910397561007;
        c1[1] := 0.4608848369928040;
        c2[1, 1] := 0.2099813050274780;
        c2[1, 2] := 0.8992478823790660;
        c2[2, 1] := 0.2027250423101359;
        c2[2, 2] := 0.8328117484224146;
        c2[3, 1] := 0.1907635894058731;
        c2[3, 2] := 0.7257379204691213;
        c2[4, 1] := 0.1742280397887686;
        c2[4, 2] := 0.5830640944868014;
        c2[5, 1] := 0.1530858190490478;
        c2[5, 2] := 0.4106192089751885;
        c2[6, 1] := 0.1264090712880446;
        c2[6, 2] := 0.2114980230156001;
      elseif order == 14 then
        alpha := 0.2331902368695848;
        c2[1, 1] := 0.1986162311411235;
        c2[1, 2] := 0.8876961808055535;
        c2[2, 1] := 0.1946683341271615;
        c2[2, 2] := 0.8500754229171967;
        c2[3, 1] := 0.1868331332895056;
        c2[3, 2] := 0.7764629313723603;
        c2[4, 1] := 0.1752118757862992;
        c2[4, 2] := 0.6699720402924552;
        c2[5, 1] := 0.1598906457908402;
        c2[5, 2] := 0.5348446712848934;
        c2[6, 1] := 0.1407810153019944;
        c2[6, 2] := 0.3755841316563539;
        c2[7, 1] := 0.1169627966707339;
        c2[7, 2] := 0.1937088226304455;
      elseif order == 15 then
        alpha := 0.2248854870552422;
        c1[1] := 0.4328492272335646;
        c2[1, 1] := 0.1857292591004588;
        c2[1, 2] := 0.8496337061962563;
        c2[2, 1] := 0.1808644178280136;
        c2[2, 2] := 0.8020517898136011;
        c2[3, 1] := 0.1728264404199081;
        c2[3, 2] := 0.7247449729331105;
        c2[4, 1] := 0.1616970125901954;
        c2[4, 2] := 0.6205369315943097;
        c2[5, 1] := 0.1475257264578426;
        c2[5, 2] := 0.4929612162355906;
        c2[6, 1] := 0.1301861023357119;
        c2[6, 2] := 0.3454770708040735;
        c2[7, 1] := 0.1087810777120188;
        c2[7, 2] := 0.1784526655428406;
      elseif order == 16 then
        alpha := 0.2174105053474761;
        c2[1, 1] := 0.1765637967473151;
        c2[1, 2] := 0.8377453068635511;
        c2[2, 1] := 0.1738525357503125;
        c2[2, 2] := 0.8102988957433199;
        c2[3, 1] := 0.1684627004613343;
        c2[3, 2] := 0.7563265923413258;
        c2[4, 1] := 0.1604519074815815;
        c2[4, 2] := 0.6776082294687619;
        c2[5, 1] := 0.1498828607802206;
        c2[5, 2] := 0.5766417034027680;
        c2[6, 1] := 0.1367764717792823;
        c2[6, 2] := 0.4563528264410489;
        c2[7, 1] := 0.1209810465419295;
        c2[7, 2] := 0.3193782657322374;
        c2[8, 1] := 0.1016312648007554;
        c2[8, 2] := 0.1652419227369036;
      elseif order == 17 then
        alpha := 0.2106355148193306;
        c1[1] := 0.4093223608497299;
        c2[1, 1] := 0.1664014345826274;
        c2[1, 2] := 0.8067173752345952;
        c2[2, 1] := 0.1629839591538256;
        c2[2, 2] := 0.7712924931447541;
        c2[3, 1] := 0.1573277802512491;
        c2[3, 2] := 0.7134213666303411;
        c2[4, 1] := 0.1494828185148637;
        c2[4, 2] := 0.6347841731714884;
        c2[5, 1] := 0.1394948812681826;
        c2[5, 2] := 0.5375594414619047;
        c2[6, 1] := 0.1273627583380806;
        c2[6, 2] := 0.4241608926375478;
        c2[7, 1] := 0.1129187258461290;
        c2[7, 2] := 0.2965752009703245;
        c2[8, 1] := 0.9533357359908857e-1;
        c2[8, 2] := 0.1537041700889585;
      elseif order == 18 then
        alpha := 0.2044575288651841;
        c2[1, 1] := 0.1588768571976356;
        c2[1, 2] := 0.7951914263212913;
        c2[2, 1] := 0.1569357024981854;
        c2[2, 2] := 0.7744529690772538;
        c2[3, 1] := 0.1530722206358810;
        c2[3, 2] := 0.7335304425992080;
        c2[4, 1] := 0.1473206710524167;
        c2[4, 2] := 0.6735038935387268;
        c2[5, 1] := 0.1397225420331520;
        c2[5, 2] := 0.5959151542621590;
        c2[6, 1] := 0.1303092459809849;
        c2[6, 2] := 0.5026483447894845;
        c2[7, 1] := 0.1190627367060072;
        c2[7, 2] := 0.3956893824587150;
        c2[8, 1] := 0.1058058030798994;
        c2[8, 2] := 0.2765091830730650;
        c2[9, 1] := 0.8974708108800873e-1;
        c2[9, 2] := 0.1435505288284833;
      elseif order == 19 then
        alpha := 0.1987936248083529;
        c1[1] := 0.3892259966869526;
        c2[1, 1] := 0.1506640012172225;
        c2[1, 2] := 0.7693121733774260;
        c2[2, 1] := 0.1481728062796673;
        c2[2, 2] := 0.7421133586741549;
        c2[3, 1] := 0.1440444668388838;
        c2[3, 2] := 0.6975075386214800;
        c2[4, 1] := 0.1383101628540374;
        c2[4, 2] := 0.6365464378910025;
        c2[5, 1] := 0.1310032283190998;
        c2[5, 2] := 0.5606211948462122;
        c2[6, 1] := 0.1221431166405330;
        c2[6, 2] := 0.4713530424221445;
        c2[7, 1] := 0.1116991161103884;
        c2[7, 2] := 0.3703717538617073;
        c2[8, 1] := 0.9948917351196349e-1;
        c2[8, 2] := 0.2587371155559744;
        c2[9, 1] := 0.8475989238107367e-1;
        c2[9, 2] := 0.1345537894555993;
      elseif order == 20 then
        alpha := 0.1935761760416219;
        c2[1, 1] := 0.1443871348337404;
        c2[1, 2] := 0.7584165598446141;
        c2[2, 1] := 0.1429501891353184;
        c2[2, 2] := 0.7423000962318863;
        c2[3, 1] := 0.1400877384920004;
        c2[3, 2] := 0.7104185332215555;
        c2[4, 1] := 0.1358210369491446;
        c2[4, 2] := 0.6634599783272630;
        c2[5, 1] := 0.1301773703034290;
        c2[5, 2] := 0.6024175491895959;
        c2[6, 1] := 0.1231826501439148;
        c2[6, 2] := 0.5285332736326852;
        c2[7, 1] := 0.1148465498575254;
        c2[7, 2] := 0.4431977385498628;
        c2[8, 1] := 0.1051289462376788;
        c2[8, 2] := 0.3477444062821162;
        c2[9, 1] := 0.9384622797485121e-1;
        c2[9, 2] := 0.2429038300327729;
        c2[10, 1] := 0.8028211612831444e-1;
        c2[10, 2] := 0.1265329974009533;
      elseif order == 21 then
        alpha := 0.1887494014766075;
        c1[1] := 0.3718070668941645;
        c2[1, 1] := 0.1376151928386445;
        c2[1, 2] := 0.7364290859445481;
        c2[2, 1] := 0.1357438914390695;
        c2[2, 2] := 0.7150167318935022;
        c2[3, 1] := 0.1326398453462415;
        c2[3, 2] := 0.6798001808470175;
        c2[4, 1] := 0.1283231214897678;
        c2[4, 2] := 0.6314663440439816;
        c2[5, 1] := 0.1228169159777534;
        c2[5, 2] := 0.5709353626166905;
        c2[6, 1] := 0.1161406100773184;
        c2[6, 2] := 0.4993087153571335;
        c2[7, 1] := 0.1082959649233524;
        c2[7, 2] := 0.4177766148584385;
        c2[8, 1] := 0.9923596957485723e-1;
        c2[8, 2] := 0.3274257287232124;
        c2[9, 1] := 0.8877776108724853e-1;
        c2[9, 2] := 0.2287218166767916;
        c2[10, 1] := 0.7624076527736326e-1;
        c2[10, 2] := 0.1193423971506988;
      elseif order == 22 then
        alpha := 0.1842668221199706;
        c2[1, 1] := 0.1323053462701543;
        c2[1, 2] := 0.7262446126765204;
        c2[2, 1] := 0.1312121721769772;
        c2[2, 2] := 0.7134286088450949;
        c2[3, 1] := 0.1290330911166814;
        c2[3, 2] := 0.6880287870435514;
        c2[4, 1] := 0.1257817990372067;
        c2[4, 2] := 0.6505015800059301;
        c2[5, 1] := 0.1214765261983008;
        c2[5, 2] := 0.6015107185211451;
        c2[6, 1] := 0.1161365140967959;
        c2[6, 2] := 0.5418983553698413;
        c2[7, 1] := 0.1097755171533100;
        c2[7, 2] := 0.4726370779831614;
        c2[8, 1] := 0.1023889478519956;
        c2[8, 2] := 0.3947439506537486;
        c2[9, 1] := 0.9392485861253800e-1;
        c2[9, 2] := 0.3090996703083202;
        c2[10, 1] := 0.8420273775456455e-1;
        c2[10, 2] := 0.2159561978556017;
        c2[11, 1] := 0.7257600023938262e-1;
        c2[11, 2] := 0.1128633732721116;
      elseif order == 23 then
        alpha := 0.1800893554453722;
        c1[1] := 0.3565232673929280;
        c2[1, 1] := 0.1266275171652706;
        c2[1, 2] := 0.7072778066734162;
        c2[2, 1] := 0.1251865227648538;
        c2[2, 2] := 0.6900676345785905;
        c2[3, 1] := 0.1227944815236645;
        c2[3, 2] := 0.6617011100576023;
        c2[4, 1] := 0.1194647013077667;
        c2[4, 2] := 0.6226432315773119;
        c2[5, 1] := 0.1152132989252356;
        c2[5, 2] := 0.5735222810625359;
        c2[6, 1] := 0.1100558598478487;
        c2[6, 2] := 0.5151027978024605;
        c2[7, 1] := 0.1040013558214886;
        c2[7, 2] := 0.4482410942032739;
        c2[8, 1] := 0.9704014176512626e-1;
        c2[8, 2] := 0.3738049984631116;
        c2[9, 1] := 0.8911683905758054e-1;
        c2[9, 2] := 0.2925028692588410;
        c2[10, 1] := 0.8005438265072295e-1;
        c2[10, 2] := 0.2044134600278901;
        c2[11, 1] := 0.6923832296800832e-1;
        c2[11, 2] := 0.1069984887283394;
      elseif order == 24 then
        alpha := 0.1761838665838427;
        c2[1, 1] := 0.1220804912720132;
        c2[1, 2] := 0.6978026874156063;
        c2[2, 1] := 0.1212296762358897;
        c2[2, 2] := 0.6874139794926736;
        c2[3, 1] := 0.1195328372961027;
        c2[3, 2] := 0.6667954259551859;
        c2[4, 1] := 0.1169990987333593;
        c2[4, 2] := 0.6362602049901176;
        c2[5, 1] := 0.1136409040480130;
        c2[5, 2] := 0.5962662188435553;
        c2[6, 1] := 0.1094722001757955;
        c2[6, 2] := 0.5474001634109253;
        c2[7, 1] := 0.1045052832229087;
        c2[7, 2] := 0.4903523180249535;
        c2[8, 1] := 0.9874509806025907e-1;
        c2[8, 2] := 0.4258751523524645;
        c2[9, 1] := 0.9217799943472177e-1;
        c2[9, 2] := 0.3547079765396403;
        c2[10, 1] := 0.8474633796250476e-1;
        c2[10, 2] := 0.2774145482392767;
        c2[11, 1] := 0.7627722381240495e-1;
        c2[11, 2] := 0.1939329108084139;
        c2[12, 1] := 0.6618645465422745e-1;
        c2[12, 2] := 0.1016670147947242;
      elseif order == 25 then
        alpha := 0.1725220521949266;
        c1[1] := 0.3429735385896000;
        c2[1, 1] := 0.1172525033170618;
        c2[1, 2] := 0.6812327932576614;
        c2[2, 1] := 0.1161194585333535;
        c2[2, 2] := 0.6671566071153211;
        c2[3, 1] := 0.1142375145794466;
        c2[3, 2] := 0.6439167855053158;
        c2[4, 1] := 0.1116157454252308;
        c2[4, 2] := 0.6118378416180135;
        c2[5, 1] := 0.1082654809459177;
        c2[5, 2] := 0.5713609763370088;
        c2[6, 1] := 0.1041985674230918;
        c2[6, 2] := 0.5230289949762722;
        c2[7, 1] := 0.9942439308123559e-1;
        c2[7, 2] := 0.4674627926041906;
        c2[8, 1] := 0.9394453593830893e-1;
        c2[8, 2] := 0.4053226688298811;
        c2[9, 1] := 0.8774221237222533e-1;
        c2[9, 2] := 0.3372372276379071;
        c2[10, 1] := 0.8075839512216483e-1;
        c2[10, 2] := 0.2636485508005428;
        c2[11, 1] := 0.7282483286646764e-1;
        c2[11, 2] := 0.1843801345273085;
        c2[12, 1] := 0.6338571166846652e-1;
        c2[12, 2] := 0.9680153764737715e-1;
      elseif order == 26 then
        alpha := 0.1690795702796737;
        c2[1, 1] := 0.1133168695796030;
        c2[1, 2] := 0.6724297955493932;
        c2[2, 1] := 0.1126417845769961;
        c2[2, 2] := 0.6638709519790540;
        c2[3, 1] := 0.1112948749545606;
        c2[3, 2] := 0.6468652038763624;
        c2[4, 1] := 0.1092823986944244;
        c2[4, 2] := 0.6216337070799265;
        c2[5, 1] := 0.1066130386697976;
        c2[5, 2] := 0.5885011413992190;
        c2[6, 1] := 0.1032969057045413;
        c2[6, 2] := 0.5478864278297548;
        c2[7, 1] := 0.9934388184210715e-1;
        c2[7, 2] := 0.5002885306054287;
        c2[8, 1] := 0.9476081523436283e-1;
        c2[8, 2] := 0.4462644847551711;
        c2[9, 1] := 0.8954648464575577e-1;
        c2[9, 2] := 0.3863930785049522;
        c2[10, 1] := 0.8368166847159917e-1;
        c2[10, 2] := 0.3212074592527143;
        c2[11, 1] := 0.7710664731701103e-1;
        c2[11, 2] := 0.2510470347119383;
        c2[12, 1] := 0.6965807988411425e-1;
        c2[12, 2] := 0.1756419294111342;
        c2[13, 1] := 0.6080674930548766e-1;
        c2[13, 2] := 0.9234535279274277e-1;
      elseif order == 27 then
        alpha := 0.1658353543067995;
        c1[1] := 0.3308543720638957;
        c2[1, 1] := 0.1091618578712746;
        c2[1, 2] := 0.6577977071169651;
        c2[2, 1] := 0.1082549561495043;
        c2[2, 2] := 0.6461121666520275;
        c2[3, 1] := 0.1067479247890451;
        c2[3, 2] := 0.6267937760991321;
        c2[4, 1] := 0.1046471079537577;
        c2[4, 2] := 0.6000750116745808;
        c2[5, 1] := 0.1019605976654259;
        c2[5, 2] := 0.5662734183049320;
        c2[6, 1] := 0.9869726954433709e-1;
        c2[6, 2] := 0.5257827234948534;
        c2[7, 1] := 0.9486520934132483e-1;
        c2[7, 2] := 0.4790595019077763;
        c2[8, 1] := 0.9046906518775348e-1;
        c2[8, 2] := 0.4266025862147336;
        c2[9, 1] := 0.8550529998276152e-1;
        c2[9, 2] := 0.3689188223512328;
        c2[10, 1] := 0.7995282239306020e-1;
        c2[10, 2] := 0.3064589322702932;
        c2[11, 1] := 0.7375174596252882e-1;
        c2[11, 2] := 0.2394754504667310;
        c2[12, 1] := 0.6674377263329041e-1;
        c2[12, 2] := 0.1676223546666024;
        c2[13, 1] := 0.5842458027529246e-1;
        c2[13, 2] := 0.8825044329219431e-1;
      elseif order == 28 then
        alpha := 0.1627710671942929;
        c2[1, 1] := 0.1057232656113488;
        c2[1, 2] := 0.6496161226860832;
        c2[2, 1] := 0.1051786825724864;
        c2[2, 2] := 0.6424661279909941;
        c2[3, 1] := 0.1040917964935006;
        c2[3, 2] := 0.6282470268918791;
        c2[4, 1] := 0.1024670101953951;
        c2[4, 2] := 0.6071189030701136;
        c2[5, 1] := 0.1003105109519892;
        c2[5, 2] := 0.5793175191747016;
        c2[6, 1] := 0.9762969425430802e-1;
        c2[6, 2] := 0.5451486608855443;
        c2[7, 1] := 0.9443223803058400e-1;
        c2[7, 2] := 0.5049796971628137;
        c2[8, 1] := 0.9072460982036488e-1;
        c2[8, 2] := 0.4592270546572523;
        c2[9, 1] := 0.8650956423253280e-1;
        c2[9, 2] := 0.4083368605952977;
        c2[10, 1] := 0.8178165740374893e-1;
        c2[10, 2] := 0.3527525188880655;
        c2[11, 1] := 0.7651838885868020e-1;
        c2[11, 2] := 0.2928534570013572;
        c2[12, 1] := 0.7066010532447490e-1;
        c2[12, 2] := 0.2288185204390681;
        c2[13, 1] := 0.6405358596145789e-1;
        c2[13, 2] := 0.1602396172588190;
        c2[14, 1] := 0.5621780070227172e-1;
        c2[14, 2] := 0.8447589564915071e-1;
      elseif order == 29 then
        alpha := 0.1598706626277596;
        c1[1] := 0.3199314513011623;
        c2[1, 1] := 0.1021101032532951;
        c2[1, 2] := 0.6365758882240111;
        c2[2, 1] := 0.1013729819392774;
        c2[2, 2] := 0.6267495975736321;
        c2[3, 1] := 0.1001476175660628;
        c2[3, 2] := 0.6104876178266819;
        c2[4, 1] := 0.9843854640428316e-1;
        c2[4, 2] := 0.5879603139195113;
        c2[5, 1] := 0.9625164534591696e-1;
        c2[5, 2] := 0.5594012291050210;
        c2[6, 1] := 0.9359356960417668e-1;
        c2[6, 2] := 0.5251016150410664;
        c2[7, 1] := 0.9047086748649986e-1;
        c2[7, 2] := 0.4854024475590397;
        c2[8, 1] := 0.8688856407189167e-1;
        c2[8, 2] := 0.4406826457109709;
        c2[9, 1] := 0.8284779224069856e-1;
        c2[9, 2] := 0.3913408089298914;
        c2[10, 1] := 0.7834154620997181e-1;
        c2[10, 2] := 0.3377643999400627;
        c2[11, 1] := 0.7334628941928766e-1;
        c2[11, 2] := 0.2802710651919946;
        c2[12, 1] := 0.6780290487362146e-1;
        c2[12, 2] := 0.2189770008083379;
        c2[13, 1] := 0.6156321231528423e-1;
        c2[13, 2] := 0.1534235999306070;
        c2[14, 1] := 0.5416797446761512e-1;
        c2[14, 2] := 0.8098664736760292e-1;
      elseif order == 30 then
        alpha := 0.1571200296252450;
        c2[1, 1] := 0.9908074847842124e-1;
        c2[1, 2] := 0.6289618807831557;
        c2[2, 1] := 0.9863509708328196e-1;
        c2[2, 2] := 0.6229164525571278;
        c2[3, 1] := 0.9774542692037148e-1;
        c2[3, 2] := 0.6108853364240036;
        c2[4, 1] := 0.9641490581986484e-1;
        c2[4, 2] := 0.5929869253412513;
        c2[5, 1] := 0.9464802912225441e-1;
        c2[5, 2] := 0.5693960175547550;
        c2[6, 1] := 0.9245027206218041e-1;
        c2[6, 2] := 0.5403402396359503;
        c2[7, 1] := 0.8982754584112941e-1;
        c2[7, 2] := 0.5060948065875106;
        c2[8, 1] := 0.8678535291732599e-1;
        c2[8, 2] := 0.4669749797983789;
        c2[9, 1] := 0.8332744242052199e-1;
        c2[9, 2] := 0.4233249626334694;
        c2[10, 1] := 0.7945356393775309e-1;
        c2[10, 2] := 0.3755006094498054;
        c2[11, 1] := 0.7515543969833788e-1;
        c2[11, 2] := 0.3238400339292700;
        c2[12, 1] := 0.7040879901685638e-1;
        c2[12, 2] := 0.2686072427439079;
        c2[13, 1] := 0.6515528854010540e-1;
        c2[13, 2] := 0.2098650589782619;
        c2[14, 1] := 0.5925168237177876e-1;
        c2[14, 2] := 0.1471138832654873;
        c2[15, 1] := 0.5225913954211672e-1;
        c2[15, 2] := 0.7775248839507864e-1;
      elseif order == 31 then
        alpha := 0.1545067022920929;
        c1[1] := 0.3100206996451866;
        c2[1, 1] := 0.9591020358831668e-1;
        c2[1, 2] := 0.6172474793293396;
        c2[2, 1] := 0.9530301275601203e-1;
        c2[2, 2] := 0.6088916323460413;
        c2[3, 1] := 0.9429332655402368e-1;
        c2[3, 2] := 0.5950511595503025;
        c2[4, 1] := 0.9288445429894548e-1;
        c2[4, 2] := 0.5758534119053522;
        c2[5, 1] := 0.9108073420087422e-1;
        c2[5, 2] := 0.5514734636081183;
        c2[6, 1] := 0.8888719137536870e-1;
        c2[6, 2] := 0.5221306199481831;
        c2[7, 1] := 0.8630901440239650e-1;
        c2[7, 2] := 0.4880834248148061;
        c2[8, 1] := 0.8335074993373294e-1;
        c2[8, 2] := 0.4496225358496770;
        c2[9, 1] := 0.8001502494376102e-1;
        c2[9, 2] := 0.4070602306679052;
        c2[10, 1] := 0.7630041338037624e-1;
        c2[10, 2] := 0.3607139804818122;
        c2[11, 1] := 0.7219760885744920e-1;
        c2[11, 2] := 0.3108783301229550;
        c2[12, 1] := 0.6768185077153345e-1;
        c2[12, 2] := 0.2577706252514497;
        c2[13, 1] := 0.6269571766328638e-1;
        c2[13, 2] := 0.2014081375889921;
        c2[14, 1] := 0.5710081766945065e-1;
        c2[14, 2] := 0.1412581515841926;
        c2[15, 1] := 0.5047740914807019e-1;
        c2[15, 2] := 0.7474725873250158e-1;
      elseif order == 32 then
        alpha := 0.1520196210848210;
        c2[1, 1] := 0.9322163554339406e-1;
        c2[1, 2] := 0.6101488690506050;
        c2[2, 1] := 0.9285233997694042e-1;
        c2[2, 2] := 0.6049832320721264;
        c2[3, 1] := 0.9211494244473163e-1;
        c2[3, 2] := 0.5946969295569034;
        c2[4, 1] := 0.9101176786042449e-1;
        c2[4, 2] := 0.5793791854364477;
        c2[5, 1] := 0.8954614071360517e-1;
        c2[5, 2] := 0.5591619969234026;
        c2[6, 1] := 0.8772216763680164e-1;
        c2[6, 2] := 0.5342177994699602;
        c2[7, 1] := 0.8554440426912734e-1;
        c2[7, 2] := 0.5047560942986598;
        c2[8, 1] := 0.8301735302045588e-1;
        c2[8, 2] := 0.4710187048140929;
        c2[9, 1] := 0.8014469519188161e-1;
        c2[9, 2] := 0.4332730387207936;
        c2[10, 1] := 0.7692807528893225e-1;
        c2[10, 2] := 0.3918021436411035;
        c2[11, 1] := 0.7336507157284898e-1;
        c2[11, 2] := 0.3468890521471250;
        c2[12, 1] := 0.6944555312763458e-1;
        c2[12, 2] := 0.2987898029050460;
        c2[13, 1] := 0.6514446669420571e-1;
        c2[13, 2] := 0.2476810747407199;
        c2[14, 1] := 0.6040544477732702e-1;
        c2[14, 2] := 0.1935412053397663;
        c2[15, 1] := 0.5509478650672775e-1;
        c2[15, 2] := 0.1358108994174911;
        c2[16, 1] := 0.4881064725720192e-1;
        c2[16, 2] := 0.7194819894416505e-1;
      elseif order == 33 then
        alpha := 0.1496489351138032;
        c1[1] := 0.3009752799176432;
        c2[1, 1] := 0.9041725460994505e-1;
        c2[1, 2] := 0.5995521047364046;
        c2[2, 1] := 0.8991117804113002e-1;
        c2[2, 2] := 0.5923764112099496;
        c2[3, 1] := 0.8906941547422532e-1;
        c2[3, 2] := 0.5804822013853129;
        c2[4, 1] := 0.8789442491445575e-1;
        c2[4, 2] := 0.5639663528946501;
        c2[5, 1] := 0.8638945831033775e-1;
        c2[5, 2] := 0.5429623519607796;
        c2[6, 1] := 0.8455834602616358e-1;
        c2[6, 2] := 0.5176379938389326;
        c2[7, 1] := 0.8240517431382334e-1;
        c2[7, 2] := 0.4881921474066189;
        c2[8, 1] := 0.7993380417355076e-1;
        c2[8, 2] := 0.4548502528082586;
        c2[9, 1] := 0.7714713890732801e-1;
        c2[9, 2] := 0.4178579388038483;
        c2[10, 1] := 0.7404596598181127e-1;
        c2[10, 2] := 0.3774715722484659;
        c2[11, 1] := 0.7062702339160462e-1;
        c2[11, 2] := 0.3339432938810453;
        c2[12, 1] := 0.6687952672391507e-1;
        c2[12, 2] := 0.2874950693388235;
        c2[13, 1] := 0.6277828912909767e-1;
        c2[13, 2] := 0.2382680702894708;
        c2[14, 1] := 0.5826808305383988e-1;
        c2[14, 2] := 0.1862073169968455;
        c2[15, 1] := 0.5321974125363517e-1;
        c2[15, 2] := 0.1307323751236313;
        c2[16, 1] := 0.4724820282032780e-1;
        c2[16, 2] := 0.6933542082177094e-1;
      elseif order == 34 then
        alpha := 0.1473858373968463;
        c2[1, 1] := 0.8801537152275983e-1;
        c2[1, 2] := 0.5929204288972172;
        c2[2, 1] := 0.8770594341007476e-1;
        c2[2, 2] := 0.5884653382247518;
        c2[3, 1] := 0.8708797598072095e-1;
        c2[3, 2] := 0.5795895850253119;
        c2[4, 1] := 0.8616320590689187e-1;
        c2[4, 2] := 0.5663615383647170;
        c2[5, 1] := 0.8493413175570858e-1;
        c2[5, 2] := 0.5488825092350877;
        c2[6, 1] := 0.8340387368687513e-1;
        c2[6, 2] := 0.5272851839324592;
        c2[7, 1] := 0.8157596213131521e-1;
        c2[7, 2] := 0.5017313864372913;
        c2[8, 1] := 0.7945402670834270e-1;
        c2[8, 2] := 0.4724089864574216;
        c2[9, 1] := 0.7704133559556429e-1;
        c2[9, 2] := 0.4395276256463053;
        c2[10, 1] := 0.7434009635219704e-1;
        c2[10, 2] := 0.4033126590648964;
        c2[11, 1] := 0.7135035113853376e-1;
        c2[11, 2] := 0.3639961488919042;
        c2[12, 1] := 0.6806813160738834e-1;
        c2[12, 2] := 0.3218025212900124;
        c2[13, 1] := 0.6448214312000864e-1;
        c2[13, 2] := 0.2769235521088158;
        c2[14, 1] := 0.6056719318430530e-1;
        c2[14, 2] := 0.2294693573271038;
        c2[15, 1] := 0.5626925196925040e-1;
        c2[15, 2] := 0.1793564218840015;
        c2[16, 1] := 0.5146352031547277e-1;
        c2[16, 2] := 0.1259877129326412;
        c2[17, 1] := 0.4578069074410591e-1;
        c2[17, 2] := 0.6689147319568768e-1;
      elseif order == 35 then
        alpha := 0.1452224267615486;
        c1[1] := 0.2926764667564367;
        c2[1, 1] := 0.8551731299267280e-1;
        c2[1, 2] := 0.5832758214629523;
        c2[2, 1] := 0.8509109732853060e-1;
        c2[2, 2] := 0.5770596582643844;
        c2[3, 1] := 0.8438201446671953e-1;
        c2[3, 2] := 0.5667497616665494;
        c2[4, 1] := 0.8339191981579831e-1;
        c2[4, 2] := 0.5524209816238369;
        c2[5, 1] := 0.8212328610083385e-1;
        c2[5, 2] := 0.5341766459916322;
        c2[6, 1] := 0.8057906332198853e-1;
        c2[6, 2] := 0.5121470053512750;
        c2[7, 1] := 0.7876247299954955e-1;
        c2[7, 2] := 0.4864870722254752;
        c2[8, 1] := 0.7667670879950268e-1;
        c2[8, 2] := 0.4573736721705665;
        c2[9, 1] := 0.7432449556218945e-1;
        c2[9, 2] := 0.4250013835198991;
        c2[10, 1] := 0.7170742126011575e-1;
        c2[10, 2] := 0.3895767735915445;
        c2[11, 1] := 0.6882488171701314e-1;
        c2[11, 2] := 0.3513097926737368;
        c2[12, 1] := 0.6567231746957568e-1;
        c2[12, 2] := 0.3103999917596611;
        c2[13, 1] := 0.6223804362223595e-1;
        c2[13, 2] := 0.2670123611280899;
        c2[14, 1] := 0.5849696460782910e-1;
        c2[14, 2] := 0.2212298104867592;
        c2[15, 1] := 0.5439628409499822e-1;
        c2[15, 2] := 0.1729443731341637;
        c2[16, 1] := 0.4981540179136920e-1;
        c2[16, 2] := 0.1215462157134930;
        c2[17, 1] := 0.4439981033536435e-1;
        c2[17, 2] := 0.6460098363520967e-1;
      elseif order == 36 then
        alpha := 0.1431515914458580;
        c2[1, 1] := 0.8335881847130301e-1;
        c2[1, 2] := 0.5770670512160201;
        c2[2, 1] := 0.8309698922852212e-1;
        c2[2, 2] := 0.5731929100172432;
        c2[3, 1] := 0.8257400347039723e-1;
        c2[3, 2] := 0.5654713811993058;
        c2[4, 1] := 0.8179117911600136e-1;
        c2[4, 2] := 0.5539556343603020;
        c2[5, 1] := 0.8075042173126963e-1;
        c2[5, 2] := 0.5387245649546684;
        c2[6, 1] := 0.7945413151258206e-1;
        c2[6, 2] := 0.5198817177723069;
        c2[7, 1] := 0.7790506514288866e-1;
        c2[7, 2] := 0.4975537629595409;
        c2[8, 1] := 0.7610613635339480e-1;
        c2[8, 2] := 0.4718884193866789;
        c2[9, 1] := 0.7406012816626425e-1;
        c2[9, 2] := 0.4430516443136726;
        c2[10, 1] := 0.7176927060205631e-1;
        c2[10, 2] := 0.4112237708115829;
        c2[11, 1] := 0.6923460172504251e-1;
        c2[11, 2] := 0.3765940116389730;
        c2[12, 1] := 0.6645495833489556e-1;
        c2[12, 2] := 0.3393522147815403;
        c2[13, 1] := 0.6342528888937094e-1;
        c2[13, 2] := 0.2996755899575573;
        c2[14, 1] := 0.6013361864949449e-1;
        c2[14, 2] := 0.2577053294053830;
        c2[15, 1] := 0.5655503081322404e-1;
        c2[15, 2] := 0.2135004731531631;
        c2[16, 1] := 0.5263798119559069e-1;
        c2[16, 2] := 0.1669320999865636;
        c2[17, 1] := 0.4826589873626196e-1;
        c2[17, 2] := 0.1173807590715484;
        c2[18, 1] := 0.4309819397289806e-1;
        c2[18, 2] := 0.6245036108880222e-1;
      elseif order == 37 then
        alpha := 0.1411669104782917;
        c1[1] := 0.2850271036215707;
        c2[1, 1] := 0.8111958235023328e-1;
        c2[1, 2] := 0.5682412610563970;
        c2[2, 1] := 0.8075727567979578e-1;
        c2[2, 2] := 0.5628142923227016;
        c2[3, 1] := 0.8015440554413301e-1;
        c2[3, 2] := 0.5538087696879930;
        c2[4, 1] := 0.7931239302677386e-1;
        c2[4, 2] := 0.5412833323304460;
        c2[5, 1] := 0.7823314328639347e-1;
        c2[5, 2] := 0.5253190555393968;
        c2[6, 1] := 0.7691895211595101e-1;
        c2[6, 2] := 0.5060183741977191;
        c2[7, 1] := 0.7537237072011853e-1;
        c2[7, 2] := 0.4835036020049034;
        c2[8, 1] := 0.7359601294804538e-1;
        c2[8, 2] := 0.4579149413954837;
        c2[9, 1] := 0.7159227884849299e-1;
        c2[9, 2] := 0.4294078049978829;
        c2[10, 1] := 0.6936295002846032e-1;
        c2[10, 2] := 0.3981491350382047;
        c2[11, 1] := 0.6690857785828917e-1;
        c2[11, 2] := 0.3643121502867948;
        c2[12, 1] := 0.6422751692085542e-1;
        c2[12, 2] := 0.3280684291406284;
        c2[13, 1] := 0.6131430866206096e-1;
        c2[13, 2] := 0.2895750997170303;
        c2[14, 1] := 0.5815677249570920e-1;
        c2[14, 2] := 0.2489521814805720;
        c2[15, 1] := 0.5473023527947980e-1;
        c2[15, 2] := 0.2062377435955363;
        c2[16, 1] := 0.5098441033167034e-1;
        c2[16, 2] := 0.1612849131645336;
        c2[17, 1] := 0.4680658811093562e-1;
        c2[17, 2] := 0.1134672937045305;
        c2[18, 1] := 0.4186928031694695e-1;
        c2[18, 2] := 0.6042754777339966e-1;
      elseif order == 38 then
        alpha := 0.1392625697140030;
        c2[1, 1] := 0.7916943373658329e-1;
        c2[1, 2] := 0.5624158631591745;
        c2[2, 1] := 0.7894592250257840e-1;
        c2[2, 2] := 0.5590219398777304;
        c2[3, 1] := 0.7849941672384930e-1;
        c2[3, 2] := 0.5522551628416841;
        c2[4, 1] := 0.7783093084875645e-1;
        c2[4, 2] := 0.5421574325808380;
        c2[5, 1] := 0.7694193770482690e-1;
        c2[5, 2] := 0.5287909941093643;
        c2[6, 1] := 0.7583430534712885e-1;
        c2[6, 2] := 0.5122376814029880;
        c2[7, 1] := 0.7451020436122948e-1;
        c2[7, 2] := 0.4925978555548549;
        c2[8, 1] := 0.7297197617673508e-1;
        c2[8, 2] := 0.4699889739625235;
        c2[9, 1] := 0.7122194706992953e-1;
        c2[9, 2] := 0.4445436860615774;
        c2[10, 1] := 0.6926216260386816e-1;
        c2[10, 2] := 0.4164072786327193;
        c2[11, 1] := 0.6709399961255503e-1;
        c2[11, 2] := 0.3857341621868851;
        c2[12, 1] := 0.6471757977022456e-1;
        c2[12, 2] := 0.3526828388476838;
        c2[13, 1] := 0.6213084287116965e-1;
        c2[13, 2] := 0.3174082831364342;
        c2[14, 1] := 0.5932799638550641e-1;
        c2[14, 2] := 0.2800495563550299;
        c2[15, 1] := 0.5629672408524944e-1;
        c2[15, 2] := 0.2407078154782509;
        c2[16, 1] := 0.5301264751544952e-1;
        c2[16, 2] := 0.1994026830553859;
        c2[17, 1] := 0.4942673259817896e-1;
        c2[17, 2] := 0.1559719194038917;
        c2[18, 1] := 0.4542996716979947e-1;
        c2[18, 2] := 0.1097844277878470;
        c2[19, 1] := 0.4070720755433961e-1;
        c2[19, 2] := 0.5852181110523043e-1;
      elseif order == 39 then
        alpha := 0.1374332900196804;
        c1[1] := 0.2779468246419593;
        c2[1, 1] := 0.7715084161825772e-1;
        c2[1, 2] := 0.5543001331300056;
        c2[2, 1] := 0.7684028301163326e-1;
        c2[2, 2] := 0.5495289890712267;
        c2[3, 1] := 0.7632343924866024e-1;
        c2[3, 2] := 0.5416083298429741;
        c2[4, 1] := 0.7560141319808483e-1;
        c2[4, 2] := 0.5305846713929198;
        c2[5, 1] := 0.7467569064745969e-1;
        c2[5, 2] := 0.5165224112570647;
        c2[6, 1] := 0.7354807648551346e-1;
        c2[6, 2] := 0.4995030679271456;
        c2[7, 1] := 0.7222060351121389e-1;
        c2[7, 2] := 0.4796242430956156;
        c2[8, 1] := 0.7069540462458585e-1;
        c2[8, 2] := 0.4569982440368368;
        c2[9, 1] := 0.6897453353492381e-1;
        c2[9, 2] := 0.4317502624832354;
        c2[10, 1] := 0.6705970959388781e-1;
        c2[10, 2] := 0.4040159353969854;
        c2[11, 1] := 0.6495194541066725e-1;
        c2[11, 2] := 0.3739379843169939;
        c2[12, 1] := 0.6265098412417610e-1;
        c2[12, 2] := 0.3416613843816217;
        c2[13, 1] := 0.6015440984955930e-1;
        c2[13, 2] := 0.3073260166338746;
        c2[14, 1] := 0.5745615876877304e-1;
        c2[14, 2] := 0.2710546723961181;
        c2[15, 1] := 0.5454383762391338e-1;
        c2[15, 2] := 0.2329316824061170;
        c2[16, 1] := 0.5139340231935751e-1;
        c2[16, 2] := 0.1929604256043231;
        c2[17, 1] := 0.4795705862458131e-1;
        c2[17, 2] := 0.1509655259246037;
        c2[18, 1] := 0.4412933231935506e-1;
        c2[18, 2] := 0.1063130748962878;
        c2[19, 1] := 0.3960672309405603e-1;
        c2[19, 2] := 0.5672356837211527e-1;
      elseif order == 40 then
        alpha := 0.1356742655825434;
        c2[1, 1] := 0.7538038374294594e-1;
        c2[1, 2] := 0.5488228264329617;
        c2[2, 1] := 0.7518806529402738e-1;
        c2[2, 2] := 0.5458297722483311;
        c2[3, 1] := 0.7480383050347119e-1;
        c2[3, 2] := 0.5398604576730540;
        c2[4, 1] := 0.7422847031965465e-1;
        c2[4, 2] := 0.5309482987446206;
        c2[5, 1] := 0.7346313704205006e-1;
        c2[5, 2] := 0.5191429845322307;
        c2[6, 1] := 0.7250930053201402e-1;
        c2[6, 2] := 0.5045099368431007;
        c2[7, 1] := 0.7136868456879621e-1;
        c2[7, 2] := 0.4871295553902607;
        c2[8, 1] := 0.7004317764946634e-1;
        c2[8, 2] := 0.4670962098860498;
        c2[9, 1] := 0.6853470921527828e-1;
        c2[9, 2] := 0.4445169164956202;
        c2[10, 1] := 0.6684507689945471e-1;
        c2[10, 2] := 0.4195095960479698;
        c2[11, 1] := 0.6497570123412630e-1;
        c2[11, 2] := 0.3922007419030645;
        c2[12, 1] := 0.6292726794917847e-1;
        c2[12, 2] := 0.3627221993494397;
        c2[13, 1] := 0.6069918741663154e-1;
        c2[13, 2] := 0.3312065181294388;
        c2[14, 1] := 0.5828873983769410e-1;
        c2[14, 2] := 0.2977798532686911;
        c2[15, 1] := 0.5568964389813015e-1;
        c2[15, 2] := 0.2625503293999835;
        c2[16, 1] := 0.5288947816690705e-1;
        c2[16, 2] := 0.2255872486520188;
        c2[17, 1] := 0.4986456327645859e-1;
        c2[17, 2] := 0.1868796731919594;
        c2[18, 1] := 0.4656832613054458e-1;
        c2[18, 2] := 0.1462410193532463;
        c2[19, 1] := 0.4289867647614935e-1;
        c2[19, 2] := 0.1030361558710747;
        c2[20, 1] := 0.3856310684054106e-1;
        c2[20, 2] := 0.5502423832293889e-1;
      elseif order == 41 then
        alpha := 0.1339811106984253;
        c1[1] := 0.2713685065531391;
        c2[1, 1] := 0.7355140275160984e-1;
        c2[1, 2] := 0.5413274778282860;
        c2[2, 1] := 0.7328319082267173e-1;
        c2[2, 2] := 0.5371064088294270;
        c2[3, 1] := 0.7283676160772547e-1;
        c2[3, 2] := 0.5300963437270770;
        c2[4, 1] := 0.7221298133014343e-1;
        c2[4, 2] := 0.5203345998371490;
        c2[5, 1] := 0.7141302173623395e-1;
        c2[5, 2] := 0.5078728971879841;
        c2[6, 1] := 0.7043831559982149e-1;
        c2[6, 2] := 0.4927768111819803;
        c2[7, 1] := 0.6929049381827268e-1;
        c2[7, 2] := 0.4751250308594139;
        c2[8, 1] := 0.6797129849758392e-1;
        c2[8, 2] := 0.4550083840638406;
        c2[9, 1] := 0.6648246325101609e-1;
        c2[9, 2] := 0.4325285673076087;
        c2[10, 1] := 0.6482554675958526e-1;
        c2[10, 2] := 0.4077964789091151;
        c2[11, 1] := 0.6300169683004558e-1;
        c2[11, 2] := 0.3809299858742483;
        c2[12, 1] := 0.6101130648543355e-1;
        c2[12, 2] := 0.3520508315700898;
        c2[13, 1] := 0.5885349417435808e-1;
        c2[13, 2] := 0.3212801560701271;
        c2[14, 1] := 0.5652528148656809e-1;
        c2[14, 2] := 0.2887316252774887;
        c2[15, 1] := 0.5402021575818373e-1;
        c2[15, 2] := 0.2545001287790888;
        c2[16, 1] := 0.5132588802608274e-1;
        c2[16, 2] := 0.2186415296842951;
        c2[17, 1] := 0.4841900639702602e-1;
        c2[17, 2] := 0.1811322622296060;
        c2[18, 1] := 0.4525419574485134e-1;
        c2[18, 2] := 0.1417762065404688;
        c2[19, 1] := 0.4173260173087802e-1;
        c2[19, 2] := 0.9993834530966510e-1;
        c2[20, 1] := 0.3757210572966463e-1;
        c2[20, 2] := 0.5341611499960143e-1;
      else
        Streams.error("Input argument order (= " + String(order) +
          ") of Bessel filter is not in the range 1..41");
      end if;

    end BesselCoefficients;

## Modelica_LinearSystems2.ZerosAndPoles.Internal.'p+a'

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Complex | p |  |
Real | a |  |

### Outputs

Type | Name | Description
--- | --- | ---
Complex | c |

### Modelica definition

    function 'p+a'
      import Modelica;
      import Modelica_LinearSystems2.Math.Complex;

      input Complex p;
      input Real a;
      output Complex c;
    protected
      Complex j = Modelica_LinearSystems2.Math.Complex.j();
    algorithm
      c := p.re + p.im*j +a;
    end 'p+a';

## Modelica_LinearSystems2.ZerosAndPoles.Internal.'p^2+k[1]*p+k[2]'

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Complex | p |  |
Real | k[2] |  |

### Outputs

Type | Name | Description
--- | --- | ---
Complex | c |

### Modelica definition

    function 'p^2+k[1]*p+k[2]'
      import Modelica;
      import Modelica.Utilities.Streams.print;
      import Modelica_LinearSystems2.Math.Complex;

      input Complex p;
      input Real k[2];
      output Complex c;
    protected
      Complex c1;
      Complex c2;
      Complex c3;
      Complex c4;
      Complex j = Modelica_LinearSystems2.Math.Complex.j();
    algorithm
      c := p.re^2 - p.im^2 + k[1]*p.re + k[2]+p.im*(2*p.re + k[1])*j;

    end 'p^2+k[1]*p+k[2]';

## Modelica_LinearSystems2.ZerosAndPoles.Internal.roots

Determine zeros of factorized polynomial

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | poly1[:] |  | [p^0] coefficients of first order polynomials
Real | poly2[:, 2] |  | [p, p^0] coefficients of second order polynomials
Integer | n_real |  | Number of real zeros computed with Internal.numberOfRealZeros

### Outputs

Type | Name | Description
--- | --- | ---
Real | realZeros[n_real] | All real zeros of poly1 and poly2
Complex | complexZeros[:] | All complex zeros of poly1 and poly2; for a complex conjugate pair, only one zero is stored

### Modelica definition

    function roots "Determine zeros of factorized polynomial"
      import Modelica;
      import Modelica_LinearSystems2.Math.Complex;

      input Real poly1[:] "[p^0] coefficients of first order polynomials";
      input Real poly2[:,2] "[p, p^0] coefficients of second order polynomials";
      input Integer n_real
        "Number of real zeros computed with Internal.numberOfRealZeros";
      output Real realZeros[n_real] "All real zeros of poly1 and poly2";
      output Complex complexZeros[:]=fill(Complex(0, 0), integer((size(poly1, 1)
           + 2*size(poly2, 1) - n_real)/2))
        "All complex zeros of poly1 and poly2; for a complex conjugate pair, only one zero is stored";
    protected
      Integer np1=size(poly1, 1);
      Integer np2=size(poly2, 1);
      Real D;
      Real D2;
      Real b;
      Integer j1;
      Integer j2;
      Complex j = Modelica_LinearSystems2.Math.Complex.j();

    algorithm
      assert(np1 <= n_real, "Size of poly1 = " + String(np1) + " > n_real " +
        " (= " + String(n_real) + ").");
      for i in 1:np1 loop
        realZeros[i] := -poly1[i];
      end for;

      j1 := np1 + 1;
      j2 := 1;
      for i in 1:np2 loop
        b := poly2[i, 1]/2;
        D := b*b - poly2[i, 2];
        D2 := sqrt(abs(D));
        if D >= 0 then
          realZeros[j1] := -b + D2;
          realZeros[j1 + 1] := -b - D2;
          j1 := j1 + 2;
        else
          complexZeros[j2] := -b+D2*j;
          j2 := j2 + 1;
        end if;
      end for;
    end roots;

## Modelica_LinearSystems2.ZerosAndPoles.Internal.frequencyRange

Determine min. and max. resonance frequencies

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | poly1[:] |  |
Real | poly2[:, 2] |  |

### Outputs

Type | Name | Description
--- | --- | ---
Boolean | w_found |
AngularVelocity | w_min | [rad/s]
AngularVelocity | w_max | [rad/s]

### Modelica definition

    function frequencyRange
      "Determine min. and max. resonance frequencies"
      import Modelica;
      import Modelica_LinearSystems2.Math.Complex;

      input Real poly1[:];
      input Real poly2[:,2];
      output Boolean w_found=false;
      output Modelica.SIunits.AngularVelocity w_min;
      output Modelica.SIunits.AngularVelocity w_max;
    protected
      Integer order=size(poly1, 1) + 2*size(poly2, 1);
      Integer n_real=numberOfRealZeros(poly1, poly2);
      Real zeros1[n_real];
      Complex zeros2[:]=fill(Complex(0, 0), integer((order - n_real)/2));
      Real w;
    algorithm
      // Compute zeros
      (zeros1,zeros2) := roots(
            poly1,
            poly2,
            n_real);

      // Compute resonance frequencies
      w_min := Modelica.Constants.inf;
      w_max := -Modelica.Constants.inf;
      for i in 1:size(zeros1, 1) loop
        if zeros1[i] <> 0 then
          w := abs(zeros1[i]);
          w_min := min(w_min, w);
          w_max := max(w_max, w);
          w_found := true;
        end if;
      end for;

      for i in 1:size(zeros2, 1) loop
        w := Complex.'abs'(zeros2[i]);
        w_min := min(w_min, w);
        w_max := max(w_max, w);
        w_found := true;
      end for;
    end frequencyRange;

## Modelica_LinearSystems2.ZerosAndPoles.Internal.frequencyRangeBode

Determine min. and max. frequencies for Bode plot

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
ZerosAndPoles | tf |  | ZerosAndPoles transfer function

### Outputs

Type | Name | Description
--- | --- | ---
AngularVelocity | w_min | [rad/s]
AngularVelocity | w_max | [rad/s]

### Modelica definition

    function frequencyRangeBode
      "Determine min. and max. frequencies for Bode plot"
      import Modelica;
      import Modelica_LinearSystems2.Math.Complex;
      import Modelica_LinearSystems2;
      import Modelica_LinearSystems2.ZerosAndPoles;

      input ZerosAndPoles tf "ZerosAndPoles transfer function";
      output Modelica.SIunits.AngularVelocity w_min;
      output Modelica.SIunits.AngularVelocity w_max;
    protected
      Real phi_min=Modelica.SIunits.Conversions.from_deg(3);
      Real real_min=1.0e-4;
      Real pi=Modelica.Constants.pi;
      Complex numZeros[:];
      Complex denZeros[:];
      Integer n_num;
      Integer n_den;
      Real w_min1;
      Real w_min2;
      Real w_max1;
      Real w_max2;
    algorithm
      // Compute zeros and poles
      (numZeros,denZeros) := ZerosAndPoles.Analysis.zerosAndPoles(tf);

      // Compute frequencies for numerator
      n_num := size(numZeros, 1);
      if n_num > 0 then
        (w_min1,w_max1) := Modelica_LinearSystems2.Internal.frequencyRangeZeros(
              numZeros,
              phi_min,
              real_min);
      end if;

      // Compute frequencies for denominator
      n_den := size(denZeros, 1);
      if n_den > 0 then
        (w_min2,w_max2) := Modelica_LinearSystems2.Internal.frequencyRangeZeros(
              denZeros,
              phi_min,
              real_min);
      end if;

      // Use largest range
      if n_num == 0 and n_den == 0 then
        w_min := 0.1;
        w_max := 10;
      elseif n_num == 0 then
        w_min := w_min2;
        w_max := w_max2;
      elseif n_den == 0 then
        w_min := w_min1;
        w_max := w_max1;
      else
        w_min := min(w_min1, w_min2);
        w_max := max(w_max1, w_max2);
      end if;
    end frequencyRangeBode;

## Modelica_LinearSystems2.ZerosAndPoles.Internal.filterToNormalized

Given [p^2,p] and [p] coefficients, transform to normalized form with highest power of p equal 1

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | c1[:] |  | [p] coefficients of polynomials (a*p + 1)
Real | c2[:, 2] |  | [p^2, p] coefficients of polynomials (b*p^2 + a*p + 1)

### Outputs

Type | Name | Description
--- | --- | ---
Real | n1[size(c1, 1)] | [p^0] coefficients of polynomials a*(p+1/a)
Real | n2[size(c2, 1), 2] | [p, p^0] coefficients of polynomials b*(p^2 + (a/b)*p + (1/b))
Real | k | Gain (product(1/a)*(1/b)

### Modelica definition

    function filterToNormalized
      "Given [p^2,p] and [p] coefficients, transform to normalized form with highest power of p equal 1"

      input Real c1[:] "[p] coefficients of polynomials (a*p + 1)";
      input Real c2[:,2] "[p^2, p] coefficients of polynomials (b*p^2 + a*p + 1)";
      output Real n1[size(c1, 1)] "[p^0] coefficients of polynomials a*(p+1/a)";
      output Real n2[size(c2, 1),2]
        "[p, p^0] coefficients of polynomials b*(p^2 + (a/b)*p + (1/b))";
      output Real k "Gain (product(1/a)*(1/b)";
    algorithm
      k := 1.0;
      for i in 1:size(c1, 1) loop
        k := k*c1[i];
        n1[i] := 1/c1[i];
      end for;

      for i in 1:size(c2, 1) loop
        k := k*c2[i, 1];
        n2[i, 1] := c2[i, 2]/c2[i, 1];
        n2[i, 2] := 1/c2[i, 1];
      end for;
    end filterToNormalized;

## Modelica_LinearSystems2.ZerosAndPoles.Internal.numberOfRealZeros

Calculate number of real zeros

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | poly1[:] |  |
Real | poly2[:, 2] |  |

### Outputs

Type | Name | Description
--- | --- | ---
Integer | result | Number of real zeros

### Modelica definition

    encapsulated function numberOfRealZeros
      "Calculate number of real zeros"
      import Modelica;
      import Modelica_LinearSystems2.ZerosAndPoles;

      input Real poly1[:];
      input Real poly2[:,2];
      output Integer result "Number of real zeros";
    protected
      Real D;
    algorithm
      result := size(poly1, 1);
      for i in 1:size(poly2, 1) loop
        D := (poly2[i, 1]/2)^2 - poly2[i, 2];
        if D >= 0 then
            // two real zeros
          result := result + 2;
        end if;
      end for;
    end numberOfRealZeros;

## Modelica_LinearSystems2.ZerosAndPoles.Internal.numberOfRealZeros2

Calculate number of real zeros

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
TransferFunction | tf |  | TransferFunction

### Outputs

Type | Name | Description
--- | --- | ---
Integer | result |

### Modelica definition

    function numberOfRealZeros2 "Calculate number of real zeros"
      import Modelica;
      import Modelica_LinearSystems2.Internal;
      import Modelica_LinearSystems2.TransferFunction;
      import Modelica_LinearSystems2.Math.Polynomial;

      input TransferFunction tf "TransferFunction";
      output Integer result=Internal.numberOfRealZeros(Polynomial.roots(Polynomial(tf.n)));
    algorithm
    end numberOfRealZeros2;

## Modelica_LinearSystems2.ZerosAndPoles.Internal.numberOfRealPoles

Calculate number of real poles

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
TransferFunction | tf |  | TransferFunction

### Outputs

Type | Name | Description
--- | --- | ---
Integer | result |

### Modelica definition

    function numberOfRealPoles "Calculate number of real poles"
      import Modelica;
      import Modelica_LinearSystems2.Internal;
      import Modelica_LinearSystems2.TransferFunction;
      import Modelica_LinearSystems2.Math.Polynomial;

      input TransferFunction tf "TransferFunction";
    output Integer result=Internal.numberOfRealZeros(Polynomial.roots(Polynomial(tf.d)));
    algorithm
    end numberOfRealPoles;

## Modelica_LinearSystems2.ZerosAndPoles.Internal.isRoot

check if frequency is an element of the complex vector zeros

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Complex | zeros[:] |  |
Complex | p |  |
Real | eps | 0 |

### Outputs

Type | Name | Description
--- | --- | ---
Boolean | result |

### Modelica definition

    function isRoot
      "check if frequency is an element of the complex vector zeros"

      import Modelica_LinearSystems2.Math.Complex;
      //import Modelica_LinearSystems2;

      input Complex zeros[:];
      input Complex p;
    //Never used
      input Real eps(min=0) = 0;
      output Boolean result;

    protected
      Integer sz=size(zeros, 1);
      Integer i;

    algorithm
      i := 1;
      result := false;
      while i <= sz and not result loop

    //The implementation of Complex.'==' does not take a third parameter.
        result := p == zeros[i];
        i := i + 1;
      end while;

    end isRoot;

## Modelica_LinearSystems2.ZerosAndPoles.Internal.numberOfRealZerosAndPoles_zp

Get the number of first oder polynomials (n1, d1) and second order polynomials (n2, d2) of zeros and poles from zeros and poles written in a MAT-file

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
String | fileName | DataDir + "/zp.mat" | Name of the zeros and poles data file
Real | eps | Modelica.Constants.eps |

### Outputs

Type | Name | Description
--- | --- | ---
Integer | n1n2d1d2[4] |

### Modelica definition

    function numberOfRealZerosAndPoles_zp
      "Get the number of first oder polynomials (n1, d1) and second order polynomials (n2, d2) of zeros and poles from zeros and poles written in a MAT-file"

      import Modelica_LinearSystems2.DataDir;
      import Modelica_LinearSystems2.Internal;

      input String fileName=DataDir + "/zp.mat"
        "Name of the zeros and poles data file";
      input Real eps=Modelica.Constants.eps;

      output Integer n1n2d1d2[4];

    protected
      Integer n1;
      Integer d1;

      Integer zSize[2]=readMatrixSize(fileName, "z");
      Integer pSize[2]=readMatrixSize(fileName, "p");

      Real zerosMatrix[zSize[1],zSize[2]]=readMatrix(
              fileName,
              "z",
              zSize[1],
              zSize[2]) "zeros in rows of real parts and imaginary parts";
      Real polesMatrix[pSize[1],pSize[2]]=readMatrix(
              fileName,
              "p",
              pSize[1],
              pSize[2]) "poles in rows of real parts and imaginary parts";

    algorithm
      n1 := zSize[1];
      d1 := pSize[1];
      for i in 1:zSize[1] loop
        if abs(zerosMatrix[i, 2]) >= eps then
          n1 := n1 - 1;
        end if;
      end for;
      for i in 1:pSize[1] loop
        if abs(polesMatrix[i, 2]) >= eps then
          d1 := d1 - 1;
        end if;
      end for;
      n1n2d1d2[1] := n1;
      n1n2d1d2[2] := div((zSize[1] - n1), 2);
      n1n2d1d2[3] := d1;
      n1n2d1d2[4] := div((pSize[1] - d1), 2);

    end numberOfRealZerosAndPoles_zp;

## Modelica_LinearSystems2.ZerosAndPoles.Internal.numberOfRealZerosAndPoles_pc

Generate a zeros and poles data record by reading the polynomial coefficients from a file (default file name is zp.mat)

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
String | fileName | "pc.mat" | Name of the zeros and poles data file

### Outputs

Type | Name | Description
--- | --- | ---
Integer | n1n2d1d2[4] |

### Modelica definition

    function numberOfRealZerosAndPoles_pc
      "Generate a zeros and poles data record by reading the polynomial coefficients from a file (default file name is zp.mat)"
      import Modelica_LinearSystems2.ZerosAndPoles;
      import Modelica_LinearSystems2.Math.Complex;
      import Modelica;

      input String fileName="pc.mat" "Name of the zeros and poles data file";
      output Integer n1n2d1d2[4];

    protected
      Integer n1Size[2]=readMatrixSize(fileName, "n1");
      Integer n2Size[2]=readMatrixSize(fileName, "n2");
      Integer d1Size[2]=readMatrixSize(fileName, "d1");
      Integer d2Size[2]=readMatrixSize(fileName, "d2");

    algorithm
      n1n2d1d2[1] := n1Size[1];
      n1n2d1d2[2] := n2Size[1];
      n1n2d1d2[3] := d1Size[1];
      n1n2d1d2[4] := d2Size[1];

    end numberOfRealZerosAndPoles_pc;

## Modelica_LinearSystems2.ZerosAndPoles.Internal.isControllableAndObservableSISO

To check whether a SISO system is controllable and observable

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
ZerosAndPoles | zp |  |

### Outputs

Type | Name | Description
--- | --- | ---
Boolean | controllableAndObservable |

### Modelica definition

    encapsulated function isControllableAndObservableSISO
      "To check whether a SISO system is controllable and observable"

      import Modelica_LinearSystems2;
      import Modelica_LinearSystems2.StateSpace;
      import Modelica_LinearSystems2.ZerosAndPoles;

    input ZerosAndPoles zp;

      output Boolean controllableAndObservable;
    protected
      StateSpace ss=StateSpace(zp);

    algorithm
      controllableAndObservable := StateSpace.Internal.isControllableAndObservableSISO(ss=ss);

    end isControllableAndObservableSISO;