From 252b0f5e414fb6ec685b18aae9004ec31274c64a Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Mon, 16 Sep 2024 16:07:51 +0200
Subject: [PATCH] add Position source (#330)

* add Position source

* add docstring
---
 .../TranslationalModelica.jl                  |  2 +-
 .../TranslationalModelica/sources.jl          | 43 +++++++++++++++++++
 test/Mechanical/translational_modelica.jl     | 24 ++++++++++-
 3 files changed, 67 insertions(+), 2 deletions(-)

diff --git a/src/Mechanical/TranslationalModelica/TranslationalModelica.jl b/src/Mechanical/TranslationalModelica/TranslationalModelica.jl
index ee16578d1..fd3ad41a1 100644
--- a/src/Mechanical/TranslationalModelica/TranslationalModelica.jl
+++ b/src/Mechanical/TranslationalModelica/TranslationalModelica.jl
@@ -13,7 +13,7 @@ include("utils.jl")
 export Fixed, Mass, Spring, Damper, SpringDamper
 include("components.jl")
 
-export Force
+export Force, Position
 include("sources.jl")
 
 end
diff --git a/src/Mechanical/TranslationalModelica/sources.jl b/src/Mechanical/TranslationalModelica/sources.jl
index 4657e226f..f3ce5645c 100644
--- a/src/Mechanical/TranslationalModelica/sources.jl
+++ b/src/Mechanical/TranslationalModelica/sources.jl
@@ -13,3 +13,46 @@ Input signal acting as external force on a flange
         flange.f ~ -f.u
     end
 end
+
+"""
+    Position(; name, exact = false, f_crit = 50)
+
+Forced movement of a flange according to a reference position
+
+The input signal `s_ref` defines the reference position in [m]. Flange flange is forced to move relative to the support connector according to this reference motion. According to parameter `exact`, this is done in the following way:
+
+- `exact=true`: The reference position is treated exactly. This is only possible, if the input signal is defined by an analytical function which can be differentiated at least twice. If this prerequisite is fulfilled, the Modelica translator will differentiate the input signal twice in order to compute the reference acceleration of the flange.
+- `exact=false`: The reference position is filtered and the second derivative of the filtered curve is used to compute the reference acceleration of the flange. This second derivative is not computed by numerical differentiation but by an appropriate realization of the filter. For filtering, a second order Bessel filter is used. The critical frequency (also called cut-off frequency) of the filter is defined via parameter `f_crit` in [Hz]. This value should be selected in such a way that it is higher as the essential low frequencies in the signal.
+
+The input signal can be provided from one of the signal generator blocks of the block library `Blocks.Sources`.
+"""
+@mtkmodel Position begin
+    @extend (s,) = ptf = PartialElementaryOneFlangeAndSupport2()
+    @structural_parameters begin
+        exact = false
+    end
+    @parameters begin
+        f_crit = 50
+    end
+    @variables begin
+        v(t)
+        a(t)
+    end
+    @components begin
+        s_ref = RealInput()
+    end
+    begin
+        w_crit = 2π * f_crit
+        af = 1.3617
+        bf = 0.6180
+    end
+    @equations begin
+        if exact
+            s ~ s_ref.u
+        else
+            a ~ ((s_ref.u - s) * w_crit - af * v) * (w_crit / bf)
+        end
+        v ~ D(s)
+        a ~ D(v)
+    end
+end
diff --git a/test/Mechanical/translational_modelica.jl b/test/Mechanical/translational_modelica.jl
index 9155a9a08..f22894c55 100644
--- a/test/Mechanical/translational_modelica.jl
+++ b/test/Mechanical/translational_modelica.jl
@@ -4,7 +4,8 @@ using ModelingToolkit: t_nounits as t, D_nounits as D
 using ModelingToolkitStandardLibrary.Blocks: Sine
 using ModelingToolkitStandardLibrary.Mechanical.TranslationalModelica: Damper, Spring, Mass,
                                                                        Fixed, Force,
-                                                                       SpringDamper
+                                                                       SpringDamper,
+                                                                       Position
 
 @testset "spring damper mass fixed" begin
     @mtkmodel SpringDamperMassFixed begin
@@ -85,3 +86,24 @@ end
     @test -lb≈ub atol=1e-2
     @test -0.11 < lb < -0.1
 end
+
+@testset "Position source" begin
+    @mtkmodel TestPositionSource begin
+        @components begin
+            p1 = Position(exact = true)
+            source = Sine(frequency = 3, amplitude = 2)
+            mass = Mass(m = 1, v = 1, s = 0)
+        end
+
+        @equations begin
+            connect(source.output, p1.s_ref)
+            connect(p1.flange, mass.flange_a)
+        end
+    end
+
+    @mtkbuild sys = TestPositionSource()
+    prob = ODEProblem(sys, [], (0, 2pi))
+    sol = solve(prob, Rodas4())
+    tv = 0:0.1:(2pi)
+    @test sol(tv, idxs = sys.mass.s)≈@.(2sin(2pi * tv * 3)) atol=1e-2
+end