diff --git a/Project.toml b/Project.toml
index 15dbeaa1..dc85223f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -22,6 +22,7 @@ julia = "1"
 [extras]
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 
 [targets]
-test = ["Test", "SparseArrays"]
+test = ["Test", "FFTW", "SparseArrays"]
diff --git a/src/elements.jl b/src/elements.jl
index d1cab768..ca701a58 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -439,21 +439,41 @@ Pins: `gate`, `source`, `drain`
         end)
 end
 
-"""
-    opamp()
+@doc doc"""
+    opamp(;maxgain=Inf, gain_bw_prod=Inf)
 
-Creates an ideal operational amplifier. It enforces the voltage between the
-input pins to be zero without sourcing any current while sourcing arbitrary
-current on the output pins wihtout restricting their voltage.
+Creates a linear operational amplifier as a voltage-controlled voltage source.
+The input current is zero while the input voltage is mapped to the output
+voltage according to the transfer function
 
-Note that the opamp has two output pins, one of which will typically be
-connected to a ground node and has to provide the current sourced on the other
-output pin.
+$H(f) = \frac{A_\text{max}}{\sqrt{A_\text{max}^2-1} i \frac{f}{f_\text{UG}} + 1}$
+
+where $f$ is the signal frequency, $A_\text{max}$ (`maxgain`) is the maximum
+open loop gain and $f_\text{UG}$ (`gain_bw_prod`) is the gain/bandwidth
+product (unity gain bandwidth). For `gain_bw_prod=Inf` (the default), this
+corresponds to a frequency-independent gain of `maxgain`. For `maxgain=Inf`
+(the default), the amplifier behaves as a perfect integrator.
+
+For both `maxgain=Inf` and `gain_bw_prod=Inf`, i.e. just `opamp()`, an ideal
+operational amplifier is obtained that enforces the voltage between the input
+pins to be zero while sourcing arbitrary current on the output pins without
+restricting their voltage.
+
+Note that the opamp has two output pins, where the negative one will typically
+be connected to a ground node and has to provide the current sourced on the
+positive one.
 
 Pins: `in+` and `in-` for input, `out+` and `out-` for output
-"""
-opamp() = Element(mv=[0 0; 1 0], mi=[1 0; 0 0],
-                  ports=["in+" => "in-", "out+" => "out-"])
+""" ->
+opamp(;maxgain=Inf, gain_bw_prod=Inf) =
+    if gain_bw_prod==Inf # special case to avoid unnecessary state
+        Element(mv=[0 0; 1 -1/maxgain], mi=[1 0; 0 0],
+                ports=["in+" => "in-", "out+" => "out-"])
+    else
+        Element(mv=[0 0; -1/sqrt(1-1/maxgain^2) 0; 0 -1], mi=[1 0; 0 0; 0 0],
+                mx=[0; 1/sqrt(maxgain^2-1); 1], mxd=[0; 1/(2π*gain_bw_prod); 0],
+                ports=["in+" => "in-", "out+" => "out-"])
+    end
 
 @doc doc"""
     opamp(Val{:macak}, gain, vomin, vomax)
diff --git a/test/runtests.jl b/test/runtests.jl
index 565bb2b7..1c269ad4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -5,6 +5,7 @@ include("checklic.jl")
 
 using ACME
 using Test: @test, @test_broken, @test_logs, @test_throws, @testset
+using FFTW: rfft
 using ProgressMeter
 using SparseArrays: sparse, spzeros
 
@@ -519,6 +520,33 @@ end
     end
 end
 
+@testset "op amp" begin
+    for Amax in (10, Inf), GBP in (50e3, Inf)
+        # test circuit: non-inverting amplifier in high shelving configuration
+        circ = @circuit begin
+            input = voltagesource(), [-] ⟷ gnd
+            op = opamp(maxgain=Amax, gain_bw_prod=GBP), ["in+"] ⟷ input[+], ["out-"] ⟷ gnd
+            r1 = resistor(109e3), [1] ⟷ op["out+"], [2] ⟷ op["in-"]
+            r2 = resistor(1e3), [1] ⟷ op["in-"]
+            c = capacitor(22e-9), [1] ⟷ r2[2], [2] ⟷ gnd
+            output = voltageprobe(), [+] ⟷ op["out+"], [-] ⟷ gnd
+        end
+        model = DiscreteModel(circ, 1/44100)
+        # obtain impulse response / transfer function
+        u = [1; zeros(4095)]'
+        y = run!(model, u)[1,:]
+        Y = rfft(y)
+        # inverse of op amp transfer function
+        G⁻¹(s) = sqrt(1-1/Amax^2)*s/(2π*GBP) + 1/Amax
+        #  feedback transfer function
+        H(s) = (1e3*22e-9*s + 1) / ((109e3+1e3)*22e-9*s + 1)
+        # overall transfer function evaluated taking frequency warping of
+        # bilinear transform into account
+        Yref = [let ω=2*44100*tan(π*k/length(y)); 1/(G⁻¹(im*ω) + H(im*ω)); end for k in eachindex(Y).-1]
+        @test Y ≈ Yref
+    end
+end
+
 function checksteady!(model)
     x_steady = steadystate!(model)
     for s in model.solvers