Fix bug in calculation of complex transfer function

HansOlsson · Mar 17, 2021 · ab5b5fb · ab5b5fb
1 parent 3bcfa2b
commit ab5b5fb
Show file tree

Hide file tree

Showing 2 changed files with 20 additions and 14 deletions.
diff --git a/Modelica/ComplexBlocks/ComplexMath/TransferFunction.mo b/Modelica/ComplexBlocks/ComplexMath/TransferFunction.mo
@@ -3,9 +3,9 @@ block TransferFunction "Complex Transfer Function"
   extends Modelica.ComplexBlocks.Interfaces.ComplexSISO;
   import Modelica.ComplexMath.j;
   parameter Real b[:]={1}
-    "Numerator coefficients of transfer function (e.g., 2*s+3 is specified as {2,3})";
+    "Numerator coefficients of transfer function (e.g., 2*(jw)+3 is specified as {2,3})";
   parameter Real a[:]={1}
-    "Denominator coefficients of transfer function (e.g., 5*s+6 is specified as {5,6})";
+    "Denominator coefficients of transfer function (e.g., 5*(jw)+6 is specified as {5,6})";
   Modelica.Blocks.Interfaces.RealInput w(unit="rad/s") "Angular frequency input" annotation (Placement(transformation(
         extent={{-20,-20},{20,20}},
         rotation=90,
@@ -19,15 +19,15 @@ protected
     input Integer k;
     output Complex x;
   protected
-     Integer m=mod(k,4);
+    Integer m=mod(k,4);
   algorithm
    x:=if m==0 then Complex(1) elseif m==1 then j elseif m==2 then Complex(-1) else -j;
    annotation(Inline=true);
   end powerOfJ;
 equation
   // Avoid computing power of Complex numbers - since it fails for w==0
-  bw = {b[i]*(w)^(i-1)*powerOfJ(i-1) for i in 1:size(b,1)};
-  aw = {a[i]*(w)^(i-1)*powerOfJ(i-1) for i in 1:size(a,1)};
+  bw = {b[i]*(w)^(size(b,1)-i)*powerOfJ(size(b,1)-i) for i in 1:size(b,1)};
+  aw = {a[i]*(w)^(size(a,1)-i)*powerOfJ(size(a,1)-i) for i in 1:size(a,1)};
   bSum = Complex(sum(bw.re), sum(bw.im));
   aSum = Complex(sum(aw.re), sum(aw.im));
   y = u*bSum/aSum;
@@ -50,5 +50,9 @@ The complex input u is multiplied by the complex transfer function (depending on
 y(jw) = ------------------------------------------------- * u(jw)
         a[1]*(jw)^[na-1] + a[2]*(jw)^[na-2] + ... + a[na]
 </pre></blockquote>
+</html>", revisions="<html>
+<ul>
+<li>2021: <em>Important bug fix</em> of the order of coefficients which has been interpreted wrongly, see <a href=\"https://github.com/modelica/ModelicaStandardLibrary/issues/3651\">#3651</a></li>
+</ul>
 </html>"));
 end TransferFunction;
diff --git a/Modelica/ComplexBlocks/Examples/ShowTransferFunction.mo b/Modelica/ComplexBlocks/Examples/ShowTransferFunction.mo
@@ -1,9 +1,11 @@
 within Modelica.ComplexBlocks.Examples;
-model ShowTransferFunction "Test Complex Transfer Function Block"
+model ShowTransferFunction "Test complex transfer function block"
   extends Modelica.Icons.Example;
-  parameter Real d=1/sqrt(2) "Damping coefficient";
-  parameter Real b[:]={1} "Numerator polynomial coefficients of the transfer function";
-  parameter Real a[:]={1,2*d,1} "Denominator polynomial coefficients of the transfer function";
+  parameter Real d=0.01 "Damping coefficient in kg/s (not the damping ratio)";
+  parameter Real m=0.2 "Mass in kg";
+  parameter Real c=0.1 "Stiffness in N/m";
+  parameter Real b[:]={-m} "Numerator polynomial coefficients {-m} of the transfer function";
+  parameter Real a[:]={m,d,c} "Denominator polynomial coefficients {m,d,c} of the transfer function";
   parameter Real wMin=0.01 "Lower bound for frequency sweep";
   parameter Real wMax=100 "Upper bound for frequency sweep";
   Real lg_w=log10(logFrequencySweep.y) "Logarithm of frequency";
@@ -28,13 +30,13 @@ equation
   connect(transferFunction.y, complexToPolar.u)
     annotation (Line(points={{-19,0},{-2,0}}, color={85,170,255}));
   annotation (
-    experiment(StopTime=1, Interval=0.001), Documentation(info=
-               "<html>
-<p>This example shows the response of a PT2 defined by its transfer function</p>
+    experiment(StopTime=1, Interval=0.001), Documentation(info="<html>
+<p>This example shows the response of a PT2 (mechanical spring-mass-damper system with
+an acceleration acting on the mass) defined by its transfer function</p>
 <blockquote><pre>
-            1
+              -m
 H(jw)=-------------------
-      1 + 2 d jw + (jw)^2
+      m*(jw)^2 + d*jw + c
 </pre></blockquote>
 <p>Frequency performs a logarithmic ramp from 0.01 to 100 s^-1.</p>
 <p>