Add smooth option for sources in Blocks

- For `Sine`, `Cosine`, `ExpSine`, `Step`, `Ramp`

Add Square and Triangular sources

Add a non-zero start_time to smooth Sine, Ramp, ExpSine tests & non-zero `offset` to smooth ExpSine

Add "almost" exact versions of Triangular and Square waves
The formulae from  https://en.wikipedia.org/wiki/Square_wave & https://en.wikipedia.org/wiki/Triangle_wave
 are modified to accomodate any amplitude, frequency, offset and start_time

translate t on x-axis with start_time for square and triangular waves

trigger workflow

Merge branch 'vk/smooth_blocks' of https://github.com/ven-k/ModelingToolkitStandardLibrary.jl into vk/smooth_blocks

thermal tests are no longer broken; so `@test_broken` -> `@test`

Author: ven-k

src/Blocks/Blocks.jl

@@ -16,7 +16,7 @@ export Abs, Sign, Sqrt, Sin, Cos, Tan, Asin, Acos, Atan, Atan2, Sinh, Cosh, Tanh
 export Log, Log10
 include("math.jl")
 
-export Constant, Sine, Cosine, ContinuousClock, Ramp, Step, ExpSine
+export Constant, Sine, Cosine, ContinuousClock, Ramp, Step, ExpSine, Square, Triangular
 include("sources.jl")
 
 export Limiter, DeadZone, SlewRateLimiter

src/Blocks/sources.jl

@@ -1,3 +1,57 @@
+# Define and register smooth functions
+# These are "smooth" aka differentiable and avoid Gibbs effect
+# These follow: `offset` + `smooth_wave` * `smooth_step` with zero output for `t < start_time`
+function smooth_cos(x, δ, f, amplitude, ϕ, offset, start_time)
+    offset + amplitude * cos(2*π*f*(x - start_time) + ϕ) * smooth_step(x, δ, one(x), zero(x), start_time)
+end
+
+function smooth_damped_sin(x, δ, f, amplitude, damping, ϕ, offset, start_time)
+    offset + exp((start_time - x)*damping)*amplitude*sin(2*π*f*(x - start_time) + ϕ) * smooth_step(x, δ, one(x), zero(x), start_time)
+end
+
+function smooth_ramp(x, δ, height, duration, offset, start_time)
+    offset + height/(duration) * (smooth_xH(x, δ, start_time) - smooth_xH(x, δ, start_time+duration))
+end
+
+function smooth_sin(x, δ, f, amplitude, ϕ, offset, start_time)
+    offset + amplitude * sin(2*pi*f*(x - start_time) + ϕ) * smooth_step(x, δ, one(x), zero(x), start_time)
+end
+
+function smooth_square(x, δ, f, amplitude, offset, start_time)
+    offset + amplitude*2atan(sin(2π*(x - start_time)*f)/δ)/π * smooth_step(x, δ, one(x), zero(x), start_time)
+end
+
+function smooth_step(x, δ, height, offset, start_time)
+    offset + height*(atan((x - start_time)/δ)/π + 0.5)
+end
+
+function smooth_triangular(x, δ, f, amplitude, offset, start_time)
+    offset + amplitude * (1-2acos((1 - δ)sin(2π*(x - start_time)*f))/π) * smooth_step(x, δ, one(x), zero(x), start_time)
+end
+
+function smooth_xH(x, δ, tₒ)
+    0.5*(x-tₒ) * (1+((x-tₒ)/sqrt((x-tₒ)^2+δ^2)))
+end
+
+function square(x, f, amplitude, offset, start_time)
+    offset + (x > start_time) * (amplitude * (4*floor(f*(x - start_time)) - 2*floor(2*(x - start_time)*f) + 1))
+end
+
+function triangular(x, f, amplitude, offset, start_time)
+    p = 1/f # period
+    offset + (x > start_time) * (4 * amplitude * f * abs(abs((x - p/4 - start_time) % p) - p/2) - amplitude)
+end
+
+@register_symbolic smooth_cos(x, δ, f, amplitude, ϕ, offset, start_time)
+@register_symbolic smooth_damped_sin(x, δ, f, amplitude, damping, ϕ, offset, start_time)
+@register_symbolic smooth_ramp(x, δ, height, duration, offset, start_time)
+@register_symbolic smooth_sin(x, δ, f, amplitude, ϕ, offset, start_time)
+@register_symbolic smooth_square(x, δ, f, amplitude, offset, start_time)
+@register_symbolic smooth_step(x, δ, height, offset, start_time)
+@register_symbolic smooth_triangular(x, δ, f, amplitude, offset, start_time)
+@register_symbolic triangular(x, f, amplitude, offset, start_time)
+@register_symbolic square(x, f, amplitude, offset, start_time)
+
 """
 Generate constant signal.
 
@@ -22,25 +76,36 @@ Generate sine signal.
 # Parameters:
 - `frequency`: [Hz] Frequency of sine wave
 - `amplitude`: Amplitude of sine wave
-- `phase`: [rad] Phase of sine wave 
+- `phase`: [rad] Phase of sine wave
 - `offset`: Offset of output signal
 - `start_time`: [s] Output `y = offset` for `t < start_time`
+- `smooth`:  If `true`, returns a smooth wave. Defaults to `false`
+             It uses a smoothing factor of `δ=1e-5`
 
 # Connectors:
 - `output`
 """
-function Sine(;name, 
-    frequency, 
+function Sine(;name,
+    frequency,
     amplitude=1,
     phase=0,
     offset=0,
-    start_time=0)
+    start_time=0,
+    smooth=false)
 
     @named output = RealOutput()
     pars = @parameters offset=offset start_time=start_time amplitude=amplitude frequency=frequency phase=phase
+    equation = if smooth == false
+        offset + ifelse(t < start_time, 0, amplitude* sin(2*pi*frequency*(t - start_time) + phase))
+    else
+        δ = 1e-5
+        smooth_sin(t, δ, frequency, amplitude, phase, offset, start_time)
+    end
+
     eqs = [
-        output.u ~ offset + ifelse(t < start_time, 0, amplitude* sin(2*pi*frequency*(t - start_time) + phase))
+        output.u ~ equation
     ]
+
     compose(ODESystem(eqs, t, [], pars; name=name), [output])
 end
 
@@ -50,32 +115,43 @@ Generate cosine signal.
 # Parameters:
 - `frequency`: [Hz] Frequency of sine wave
 - `amplitude`: Amplitude of sine wave
-- `phase`: [rad] Phase of sine wave 
+- `phase`: [rad] Phase of sine wave
 - `offset`: Offset of output signal
 - `start_time`: [s] Output `y = offset` for `t < start_time`
+- `smooth`:  If `true`, returns a smooth wave. Defaults to `false`
+             It uses a smoothing factor of `δ=1e-5`
 
 # Connectors:
 - `output`
 """
-function Cosine(;name, 
-    frequency, 
+
+function Cosine(;name,
+    frequency,
     amplitude=1,
     phase=0,
     offset=0,
-    start_time=0)
+    start_time=0,
+    smooth=false)
 
     @named output = RealOutput()
     pars = @parameters offset=offset start_time=start_time amplitude=amplitude frequency=frequency phase=phase
+    equation = if smooth == false
+        offset + ifelse(t < start_time, zero(t), amplitude* cos(2*pi*frequency*(t - start_time) + phase))
+    else
+        δ = 1e-5
+        smooth_cos(t, δ, frequency, amplitude, phase, offset, start_time)
+    end
     eqs = [
-        output.u ~ offset + ifelse(t < start_time, 0, amplitude* cos(2*pi*frequency*(t - start_time) + phase))
+        output.u ~ equation
     ]
+
     compose(ODESystem(eqs, t, [], pars; name=name), [output])
 end
 
 """
 Generate current time signal.
 
-# Parameters: 
+# Parameters:
 - `offset`: Offset of output signal
 - `start_time`: [s] Output `y = offset` for `t < start_time`
 
@@ -86,8 +162,9 @@ function ContinuousClock(;name, offset=0, start_time=0)
     @named output = RealOutput()
     pars = @parameters offset=offset start_time=start_time
     eqs = [
-        output.u ~ offset + ifelse(t < start_time, 0, t - start_time)
+        output.u ~ offset + ifelse(t < start_time, zero(t), t - start_time)
     ]
+
     compose(ODESystem(eqs, t, [], pars; name=name), [output])
 end
 
@@ -99,22 +176,73 @@ Generate ramp signal.
 - `duration`: [s] Duration of ramp (= 0.0 gives a Step)
 - `offset`: Offset of output signal
 - `start_time`: [s] Output `y = offset` for `t < start_time`
+- `smooth`:  If `true`, returns a smooth wave. Defaults to `false`
+             It uses a smoothing factor of `δ=1e-5`
 
 # Connectors:
 - `output`
 """
-function Ramp(;name, 
-    offset=0,
+function Ramp(;name,
     height=1,
-    duration=1, 
-    start_time=0)
+    duration=1,
+    offset=0,
+    start_time=0,
+    smooth=false)
 
     @named output = RealOutput()
     pars = @parameters offset=offset start_time=start_time height=height duration=duration
-    eqs = [
-        output.u ~ offset + ifelse(t < start_time, 0, 
+    equation = if smooth == false
+        offset + ifelse(t < start_time, 0,
             ifelse(t < (start_time + duration), (t - start_time) * height / duration, height))
+    else
+        δ = 1e-5
+        smooth_ramp(t, δ, height, duration, offset, start_time)
+    end
+
+    eqs = [
+        output.u ~ equation
+    ]
+
+    compose(ODESystem(eqs, t, [], pars; name=name), [output])
+end
+
+"""
+Generate smooth square signal.
+
+# Parameters:
+- `frequency`: [Hz] Frequency of square wave
+- `amplitude`: Amplitude of square wave
+- `offset`: Offset of output signal
+- `start_time`: [s] Output `y = offset` for `t < start_time`
+- `smooth`:  If `true`, returns a smooth wave. Defaults to `false`
+             It uses a smoothing factor of `δ=1e-5`
+
+# Connectors:
+- `output`
+"""
+function Square(; name, frequency=1.0, amplitude=1.0,
+    offset=0.0, start_time=0.0, smooth=false)
+    δ = 1e-5
+
+    @named output = RealOutput()
+    pars = @parameters begin
+        frequency=frequency
+        amplitude=amplitude
+        offset=offset
+        start_time=start_time
+    end
+
+    equation = if smooth == false
+        square(t, frequency, amplitude, offset, start_time)
+    else
+        δ = 1e-5
+        smooth_square(t, δ, frequency, amplitude, offset, start_time)
+    end
+
+    eqs = [
+        output.u ~ equation
     ]
+
     compose(ODESystem(eqs, t, [], pars; name=name), [output])
 end
 
@@ -125,16 +253,26 @@ Generate step signal.
 - `height`: Height of step
 - `offset`: Offset of output signal
 - `start_time`: [s] Output `y = offset` for `t < start_time`
+- `smooth`:  If `true`, returns a smooth wave. Defaults to `false`
+             It uses a smoothing factor of `δ=1e-5`
 
 # Connectors:
 - `output`
 """
-function Step(;name, offset=0, height=1, start_time=0)
+function Step(;name, height=1, offset=0, start_time=0, smooth=true)
     @named output = RealOutput()
     pars = @parameters offset=offset start_time=start_time height=height
+    equation = if smooth == false
+        offset + ifelse(t < start_time, zero(t), height)
+    else
+        δ = 1e-5
+        smooth_step(t, δ, height, offset, start_time)
+    end
+
     eqs = [
-        output.u ~ offset + ifelse(t < start_time, 0, height)
+        output.u ~ equation
     ]
+
     compose(ODESystem(eqs, t, [], pars; name=name), [output])
 end
 
@@ -145,26 +283,81 @@ Generate exponentially damped sine signal.
 - `frequency`: [Hz] Frequency of sine wave
 - `amplitude`: Amplitude of sine wave
 - `damping`: [1/s] Damping coefficient of sine wave
-- `phase`: [rad] Phase of sine wave 
+- `phase`: [rad] Phase of sine wave
 - `offset`: Offset of output signal
 - `start_time`: [s] Output `y = offset` for `t < start_time`
+- `smooth`:  If `true`, returns a smooth wave. Defaults to `false`
+             It uses a smoothing factor of `δ=1e-5`
 
 # Connectors:
 - `output`
 """
-function ExpSine(;name, 
-    frequency, 
+function ExpSine(; name,
+    frequency,
     amplitude=1,
     damping=0.1,
     phase=0,
     offset=0,
-    start_time=0)
+    start_time=0,
+    smooth=false)
 
     @named output = RealOutput()
     pars = @parameters offset=offset start_time=start_time amplitude=amplitude frequency=frequency phase=phase damping=damping
+
+    equation = if smooth == false
+        offset + ifelse(t < start_time, 0, amplitude * exp(-damping * (t - start_time)) * sin(2*pi*frequency*(t - start_time) + phase))
+    else
+        δ = 1e-5
+        smooth_damped_sin(t, δ, frequency, amplitude, damping, phase, offset, start_time)
+    end
+
+    eqs = [
+        output.u ~ equation
+    ]
+
+    compose(ODESystem(eqs, t, [], pars; name=name), [output])
+end
+
+"""
+Generate smooth triangular signal for frequencies less than or equal to 25 Hz
+
+# Parameters:
+- `frequency`: [Hz] Frequency of square wave
+- `amplitude`: Amplitude of square wave
+- `offset`: Offset of output signal.
+- `start_time`: [s] Output `y = offset` for `t < start_time`
+- `smooth`:  If `true`, returns a smooth wave. Defaults to `false`
+             It uses a smoothing factor of `δ=1e-5`
+
+# Connectors:
+- `output`
+"""
+function Triangular(; name, amplitude=1.0, frequency=1.0,
+    offset=0.0, start_time=0.0, smooth=false)
+
+    if smooth
+        frequency > 25 && @warn "`frequency > 25` can lead to non-triangular wave pattern"
+    end
+
+    @named output = RealOutput()
+    pars = @parameters begin
+        amplitude=amplitude
+        frequency=frequency
+        offset=offset
+        start_time=start_time
+    end
+
+    equation = if smooth == false
+        triangular(t, frequency, amplitude, offset, start_time)
+    else
+        δ = 1e-5
+        smooth_triangular(t, δ, frequency, amplitude, offset, start_time)
+    end
+
     eqs = [
-        output.u ~ offset + ifelse(t < start_time, 0, amplitude * exp(-damping * (t - start_time)) * sin(2*pi*frequency*(t - start_time) + phase))
+        output.u ~ equation
     ]
+
     compose(ODESystem(eqs, t, [], pars; name=name), [output])
 end