Support Unitful in element construction (#93)

* Support `Unitful` in element contructors via package extension

* Briefly mention Unitful support in documentation

Author: martinholters

Project.toml

@@ -19,10 +19,17 @@ ProgressMeter = "0.6, 0.7, 0.8, 0.9, 1"
 StaticArrays = "0.8, 0.9, 0.10, 0.11, 0.12, 1.0"
 julia = "1.4"
 
+[extensions]
+UnitfulExt = "Unitful"
+
 [extras]
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [targets]
-test = ["Test", "FFTW", "SparseArrays"]
+test = ["Test", "FFTW", "SparseArrays", "Unitful"]
+
+[weakdeps]
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"

docs/src/ug.md

@@ -5,6 +5,20 @@
 All circuit elements are created by calling corresponding functions; see the
 [Element Reference](@ref) for details.
 
+### Unitful elements
+
+ACME provides a package extension for
+[Unitful](https://github.com/PainterQubits/Unitful.jl) to support quantities
+with units when constructing elements. E.g. `resistor(4.7e3)` and
+`resistor(4.7u"kΩ")` are equivalent after `Unitful` has been loaded. This can
+increase readability and help catch bugs (e.g. `resistor(5u"V")` will throw an
+error). The input and output signals of the curcuit models will still be
+unitless, however.
+
+!!! compat "Julia 1.9"
+    Package extensions require Julia 1.9 or later. Consequently, unitful
+    quantities are not supported on earlier Julia versions.
+
 ## Circuit Description
 
 Circuits are described using `Circuit` instances, which are most easily created

ext/UnitfulExt.jl

@@ -0,0 +1,117 @@
+# Copyright 2023 Martin Holters
+# See accompanying license file.
+
+module UnitfulExt
+
+import ACME
+using Unitful: Unitful, @u_str, NoUnits, Quantity, Units, uconvert
+
+remove_unit(unit::Units, q::Number) = NoUnits(uconvert(unit, q) / unit)
+
+ACME.resistor(r::Quantity) = ACME.resistor(remove_unit(u"Ω", r))
+
+ACME.potentiometer(r::Quantity, pos) = ACME.potentiometer(remove_unit(u"Ω", r), pos)
+ACME.potentiometer(r::Quantity) = ACME.potentiometer(remove_unit(u"Ω", r))
+
+ACME.capacitor(c::Quantity) = ACME.capacitor(remove_unit(u"F", c))
+
+ACME.inductor(l::Quantity) = ACME.inductor(remove_unit(u"H", l))
+
+function ACME.transformer(
+    l1::Quantity, l2::Quantity;
+    coupling_coefficient=1,
+    mutual_coupling::Quantity=coupling_coefficient*sqrt(l1*l2)
+)
+    return ACME.transformer(
+        remove_unit(u"H", l1), remove_unit(u"H", l2);
+        mutual_coupling = remove_unit(u"H", mutual_coupling)
+    )
+end
+
+function ACME._transformer_ja(ns, α, c, kwargs::NamedTuple{K, <:Tuple{Vararg{Union{Real,Quantity}}}}) where {K}
+    units = map(
+        (k) -> begin
+            if k === :D
+                return u"m"
+            elseif k === :A
+                return u"m^2"
+            elseif k === :a || k === :k || k === :Ms
+                return u"A/m"
+            else
+                throw(ArgumentError("transformer: got unsupported keyword argument \"$(k)\""))
+            end
+        end,
+        K
+    )
+    return ACME._transformer_ja(
+        ns, α, c,
+        NamedTuple{K}(map((unit, value) -> remove_unit(unit, value), units, values(kwargs))),
+    )
+end
+
+ACME.voltagesource(v::Quantity; rs::Quantity=0u"Ω") =
+    ACME.voltagesource(remove_unit(u"V", v); rs=remove_unit(u"Ω", rs))
+ACME._voltagesource(rs::Quantity) = ACME._voltagesource(remove_unit(u"Ω", rs))
+
+ACME.currentsource(i::Quantity; gp::Quantity=0u"S") =
+    ACME.currentsource(remove_unit(u"A", i); gp=remove_unit(u"S", gp))
+ACME._currentsource(gp::Quantity) = ACME._currentsource(remove_unit(u"S", gp))
+
+ACME._voltageprobe(gp::Quantity) = ACME._voltageprobe(remove_unit(u"S", gp))
+
+ACME._currentprobe(rs::Quantity) = ACME._currentprobe(remove_unit(u"Ω", rs))
+
+ACME._diode(is::Quantity, η::Real) = ACME._diode(remove_unit(u"A", is), η)
+
+function ACME._bjt(typ, kwargs::NamedTuple{K, <:Tuple{Vararg{Union{Real,Quantity}}}}) where {K}
+    units = map(
+        (k) -> begin
+            if k === :is || k === :isc || k === :ise || k === :ilc || k === :ile || k === :ikf || k === :ikr
+                return u"A"
+            elseif k === :vaf || k === :var
+                return u"V"
+            elseif k === :re || k === :rc || k === :rb
+                return u"Ω"
+            elseif k === :η || k === :ηc || k === :ηe || k === :βf || k === :βr || k === :ηcl || k === :ηel
+                return NoUnits
+            else
+                throw(ArgumentError("bjt: got unsupported keyword argument \"$(k)\""))
+            end
+        end,
+        K
+    )
+    return ACME._bjt(
+        typ,
+        NamedTuple{K}(map((unit, value) -> remove_unit(unit, value), units, values(kwargs))),
+    )
+end
+
+_mosfet_remove_units(unit::Units, q::Number) = remove_unit(unit, q)
+_mosfet_remove_units(unit::Units, q::Tuple{Vararg{Number,N}}) where {N} =
+    ntuple(n -> remove_unit(unit / u"V"^(n-1), q[n]), Val(N))
+
+function ACME._mosfet(typ, kwargs::NamedTuple{K, <:Tuple{Vararg{Union{Union{Real,Quantity},Tuple{Vararg{Union{Real,Quantity}}}}}}}) where {K}
+    units = map(
+        (k) -> begin
+            if k === :vt
+                return u"V"
+            elseif k === :α
+                return u"A/V^2"
+            elseif k === :λ
+                return u"V^-1"
+            else
+                throw(ArgumentError("bjt: got unsupported keyword argument \"$(k)\""))
+            end
+        end,
+        K
+    )
+    return ACME._mosfet(
+        typ,
+        NamedTuple{K}(map((unit, q) -> _mosfet_remove_units(unit, q), units, values(kwargs))),
+    )
+end
+
+ACME.opamp(::Type{Val{:macak}}, gain::Real, vomin::Quantity, vomax::Quantity) =
+    ACME.opamp(Val{:macak}, gain, remove_unit(u"V", vomin), remove_unit(u"V", vomax))
+
+end

src/elements.jl

@@ -1,4 +1,4 @@
-# Copyright 2015, 2016, 2017, 2018, 2019, 2020, 2021 Martin Holters
+# Copyright 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2023 Martin Holters
 # See accompanying license file.
 
 export resistor, potentiometer, capacitor, inductor, transformer,
@@ -97,8 +97,12 @@ Magnetization"](http://dafx.de/paper-archive/2016/dafxpapers/08-DAFx-16_paper_10
 Pins: `1` and `2` for primary winding, `3` and `4` for secondary winding, and so
 on
 """
-function transformer(::Type{Val{:JA}}; D=2.4e-2, A=4.54e-5, ns=[],
-                     a=14.1, α=5e-5, c=0.55, k=17.8, Ms=2.75e5)
+transformer(::Type{Val{:JA}}; ns=[], α=5e-5, c=0.55, kwargs...) =
+    _transformer_ja(ns, α, c, (; kwargs...))
+_transformer_ja(ns, α, c, kwargs::NamedTuple{<:Any, <:Tuple{Vararg{Real}}}) =
+    __transformer_ja(; ns=ns, α=α, c=c, kwargs...)
+function __transformer_ja(; D=2.4e-2, A=4.54e-5, ns=[],
+                        a=14.1, α=5e-5, c=0.55, k=17.8, Ms=2.75e5)
     μ0 = 1.2566370614e-6
     nonlinear_eq = @inline function (q)
         coth_q1 = coth(q[1])
@@ -175,7 +179,8 @@ Pins: `+` and `-` with `v` being measured from `+` to `-`
 """
 function voltagesource end
 voltagesource(v; rs=0) = Element(mv=1, mi=-rs, u0=v, ports=[:+ => :-])
-voltagesource(; rs=0) = Element(mv=1, mi=-rs, mu=1, ports=[:+ => :-])
+voltagesource(; rs=0) = _voltagesource(rs)
+_voltagesource(rs) = Element(mv=1, mi=-rs, mu=1, ports=[:+ => :-])
 
 """
     currentsource(; gp=0)
@@ -190,7 +195,8 @@ Pins: `+` and `-` where `i` measures the current leaving source at the `+` pin
 """
 function currentsource end
 currentsource(i; gp=0) = Element(mv=gp, mi=-1, u0=i, ports=[:+ => :-])
-currentsource(; gp=0) = Element(mv=gp, mi=-1, mu=1, ports=[:+ => :-])
+currentsource(; gp=0) = _currentsource(gp)
+_currentsource(gp) = Element(mv=gp, mi=-1, mu=1, ports=[:+ => :-])
 
 """
     voltageprobe()
@@ -201,7 +207,8 @@ defaults to zero.
 
 Pins: `+` and `-` with the output voltage being measured from `+` to `-`
 """
-voltageprobe(;gp=0) = Element(mv=-gp, mi=1, pv=1, ports=[:+ => :-])
+voltageprobe(;gp=0) = _voltageprobe(gp)
+_voltageprobe(gp) = Element(mv=-gp, mi=1, pv=1, ports=[:+ => :-])
 
 """
     currentprobe()
@@ -213,7 +220,8 @@ defaults to zero.
 Pins: `+` and `-` with the output current being the current entering the probe
 at `+`
 """
-currentprobe(;rs=0) = Element(mv=1, mi=-rs, pi=1, ports=[:+ => :-])
+currentprobe(;rs=0) = _currentprobe(rs)
+_currentprobe(rs=0) = Element(mv=1, mi=-rs, pi=1, ports=[:+ => :-])
 
 @doc raw"""
     diode(;is=1e-12, η = 1)
@@ -224,7 +232,8 @@ The reverse saturation current `is` has to be given in Ampere, the emission
 coefficient `η` is unitless.
 
 Pins: `+` (anode) and `-` (cathode)
-""" diode(;is::Real=1e-12, η::Real = 1) =
+""" diode(;is=1e-12, η = 1) = _diode(is, η)
+_diode(is::Real, η::Real) =
   Element(mv=[1;0], mi=[0;1], mq=[-1 0; 0 -1], ports=[:+ => :-], nonlinear_eq =
         @inline function(q)
             v, i = q
@@ -295,7 +304,9 @@ The parameters are set using named arguments:
 | `rb`  | Base terminal resistance
 
 Pins: `base`, `emitter`, `collector`
-""" function bjt(typ; is=1e-12, η=1, isc=is, ise=is, ηc=η, ηe=η, βf=1000, βr=10,
+""" bjt(typ; kwargs...) = _bjt(typ, (; kwargs...))
+_bjt(typ, kwargs::NamedTuple{<:Any, <:Tuple{Vararg{Real}}}) = __bjt(typ; kwargs...)
+function __bjt(typ; is=1e-12, η=1, isc=is, ise=is, ηc=η, ηe=η, βf=1000, βr=10,
              ile=0, ilc=0, ηcl=ηc, ηel=ηe, vaf=Inf, var=Inf, ikf=Inf, ikr=Inf,
              re=0, rc=0, rb=0)
     local polarity
@@ -419,7 +430,10 @@ respectively. E.g. with `vt=(0.7, 0.1, 0.02)`, the $v_{GS}$-dpendent threshold
 voltage $v_T = 0.7 + 0.1\cdot v_{GS} + 0.02\cdot v_{GS}^2$ will be used.
 
 Pins: `gate`, `source`, `drain`
-""" function mosfet(typ; vt=0.7, α=2e-5, λ=0)
+""" mosfet(typ; kwargs...) = _mosfet(typ, (; kwargs...))
+_mosfet(typ, kwargs::NamedTuple{<:Any, <:Tuple{Vararg{Union{Real,Tuple{Vararg{Real}}}}}}) =
+     __mosfet(typ; kwargs...)
+function __mosfet(typ; vt=0.7, α=2e-5, λ=0)
     if typ == :n
         polarity = 1
     elseif typ == :p
@@ -519,7 +533,7 @@ connected to a ground node and has to provide the current sourced on the other
 output pin.
 
 Pins: `in+` and `in-` for input, `out+` and `out-` for output
-""" function opamp(::Type{Val{:macak}}, gain, vomin, vomax)
+""" function opamp(::Type{Val{:macak}}, gain::Real, vomin::Real, vomax::Real)
     offset = 0.5 * (vomin + vomax)
     scale = 0.5 * (vomax - vomin)
     nonlinear_eq =

test/runtests.jl

@@ -794,3 +794,7 @@ end
         # TODO: further validate y
     end
 end
+
+if isdefined(Base, :get_extension) # Julia 1.10
+    include("unitful.jl")
+end

test/unitful.jl

@@ -0,0 +1,94 @@
+# Copyright 2023 Martin Holters
+# See accompanying license file.
+
+using Unitful: @u_str, DimensionError
+
+@testset "element constructors supporting Unitful" begin
+    @test resistor(10) == resistor(10u"Ω")
+    @test resistor(1000) == resistor(1u"kΩ")
+    @test_throws DimensionError resistor(1u"kV")
+
+    @test potentiometer(10) == potentiometer(10u"Ω")
+    @test potentiometer(1000) == potentiometer(1u"kΩ")
+    @test_throws DimensionError potentiometer(1u"kV")
+    @test potentiometer(10, 0.8) == potentiometer(10u"Ω", 0.8)
+    @test potentiometer(1000, 0.8) == potentiometer(1u"kΩ", 0.8)
+    @test_throws DimensionError potentiometer(1u"kV", 0.8)
+
+    @test capacitor(10) == capacitor(10u"F")
+    @test capacitor(1e-9) == capacitor(1.0u"nF")
+    @test_throws DimensionError capacitor(1u"kV")
+
+    @test inductor(10) == inductor(10u"H")
+    @test inductor(1e-9) == inductor(1.0u"nH")
+    @test_throws DimensionError inductor(1u"kV")
+
+    @test inductor(Val{:JA}; D=10e-3) == inductor(Val{:JA}; D=10.0u"mm")
+    @test inductor(Val{:JA}; D=10e-3, A=45e-6, n=200, a=14, α=4e-5, c=0.5, k=17, Ms=275e3) ==
+        inductor(Val{:JA}; D=10.0u"mm", A=45e-6u"m^2", n=200, a=14u"A/m", α=4e-5, c=0.5, k=17u"A/m", Ms=275.0u"kA/m")
+    @test_throws DimensionError inductor(Val{:JA}; D=10.0u"mm^2")
+
+    @test transformer(10, 5) == transformer(10u"H", 5u"H")
+    @test transformer(1e-9, 5e-9) == transformer(1.0u"nH", 5.0u"nH")
+    @test transformer(4, 9; coupling_coefficient = 0.8) == transformer(4u"H", 9u"H"; coupling_coefficient = 0.8)
+    @test transformer(1e-9, 5e-9; mutual_coupling = 2e-9) == transformer(1.0u"nH", 5.0u"nH"; mutual_coupling = 2.0u"nH")
+    @test_throws DimensionError transformer(1u"kV", 2u"kV")
+    @test_throws TypeError transformer(1u"nH", 2u"nH"; mutual_coupling = 2e-9)
+    @test_throws DimensionError transformer(1u"nH", 2u"nH"; mutual_coupling = 2u"kV")
+
+    @test transformer(Val{:JA}; D=10e-3) == transformer(Val{:JA}; D=10.0u"mm")
+    @test transformer(Val{:JA}; D=10e-3, A=45e-6, ns=[20], a=14, α=4e-5, c=0.5, k=17, Ms=275e3) ==
+        transformer(Val{:JA}; D=10.0u"mm", A=45e-6u"m^2", ns=[20], a=14u"A/m", α=4e-5, c=0.5, k=17u"A/m", Ms=275.0u"kA/m")
+    @test_throws DimensionError transformer(Val{:JA}; D=10.0u"mm^2")
+
+    @test voltagesource(10) == voltagesource(10u"V")
+    @test voltagesource(1e-3) == voltagesource(1.0u"mV")
+    @test voltagesource(5; rs=0.3) == voltagesource(5.0u"V"; rs=0.3u"Ω")
+    @test voltagesource(; rs=0.3) == voltagesource(; rs=0.3u"Ω")
+    @test_throws DimensionError voltagesource(1u"kA")
+    @test_throws TypeError voltagesource(1u"V"; rs=1)
+    @test_throws DimensionError voltagesource(1u"V"; rs=1u"V")
+    @test_throws DimensionError voltagesource(; rs=1u"V")
+
+    @test currentsource(10) == currentsource(10u"A")
+    @test currentsource(1e-3) == currentsource(1.0u"mA")
+    @test currentsource(5; gp=0.3) == currentsource(5.0u"A"; gp=0.3u"S")
+    @test currentsource(; gp=0.3) == currentsource(; gp=0.3u"S")
+    @test_throws DimensionError currentsource(1u"kV")
+    @test_throws TypeError currentsource(1u"A"; gp=1)
+    @test_throws DimensionError currentsource(1u"A"; gp=1u"A")
+    @test_throws DimensionError currentsource(; gp=1u"A")
+
+    @test voltageprobe(; gp=0.3) == voltageprobe(; gp=0.3u"S")
+    @test_throws DimensionError voltageprobe(; gp=1u"A")
+
+    @test currentprobe(; rs=0.3) == currentprobe(; rs=0.3u"Ω")
+    @test_throws DimensionError currentprobe(; rs=1u"V")
+
+    @test diode(; is=1e-15) == diode(; is=1.0u"fA")
+    @test_throws DimensionError diode(; is=1u"V")
+    @test_throws MethodError diode(; η=1u"V")
+
+    @test bjt(:npn; is=1e-15) == bjt(:npn; is=1.0u"fA")
+    @test bjt(:npn; isc=2e-12, ise=3e-12, ηc=1.1, ηe=1.2, βf=1000, βr=10, ile=4e-15,
+              ilc=5e-15, ηcl=1.15, ηel=1.18, vaf=40, var=30, ikf=5e-12, ikr=6e-12, re=1,
+              rc=2, rb=3) ==
+        bjt(:npn; isc=2e-12u"A", ise=3e-12u"A", ηc=1.1, ηe=1.2, βf=1000, βr=10,
+            ile=4e-15u"A", ilc=5e-15u"A", ηcl=1.15, ηel=1.18, vaf=40u"V", var=30u"V",
+            ikf=5e-12u"A", ikr=6e-12u"A", re=1u"Ω", rc=2u"Ω", rb=3u"Ω")
+    @test_throws DimensionError bjt(:npn; is=1u"V")
+    @test_throws DimensionError bjt(:npn; βr=1u"A")
+
+    @test mosfet(:n; vt=0.5) == mosfet(:n; vt=500.0u"mV")
+    @test mosfet(:n; vt=0.5, α=25e-6, λ=0.1) ==
+        mosfet(:n; vt=500.0u"mV", α=25.0e-6u"A/V^2", λ=0.1u"V^-1")
+    @test mosfet(:n; vt=(-1.2454, -0.199, -0.0483), α=(0.0205, -0.0017), λ=0.1) ==
+        mosfet(:n; vt=(-1.2454u"V", -0.199, -0.0483u"V^-1"), α=(0.0205u"A/V^2", -0.0017u"A/V^3"), λ=0.1u"V^-1")
+    @test_throws DimensionError mosfet(:n; vt=500.0u"mA")
+    @test_throws DimensionError mosfet(:n; vt=(500.0u"mV", 500.0u"mV"))
+
+    @test opamp(Val{:macak}, 10, -3, 5) == opamp(Val{:macak}, 10, -3u"V", 5u"V")
+    @test_throws DimensionError opamp(Val{:macak}, 10, -3u"V", 5u"A")
+    @test_throws MethodError opamp(Val{:macak}, 10, -3u"V", 5)
+    @test_throws MethodError opamp(Val{:macak}, 10u"V", -3u"V", 5u"V")
+end