diff --git a/lib/ControlSystemsBase/src/analysis.jl b/lib/ControlSystemsBase/src/analysis.jl
index 022018e32..c54a5801b 100644
--- a/lib/ControlSystemsBase/src/analysis.jl
+++ b/lib/ControlSystemsBase/src/analysis.jl
@@ -460,8 +460,8 @@ function sisomargin(sys::LTISystem, w::AbstractVector{<:Real}; full=false, allMa
     if !allMargins #Only output the smallest margins
         gm, idx = findmin([gm;Inf])
         wgm = [wgm;NaN][idx]
-        fi = [fi;NaN][idx]
         pm, idx = findmin([abs.(pm);Inf])
+        fi = [fi;NaN][idx]
         wpm = [wpm;NaN][idx]
         if full
             if !isnan(fi) #fi may be NaN, fullPhase is a scalar
@@ -509,19 +509,19 @@ function _findCrossings(w, n, res)
     tcross = Vector{eltype(w)}()
     for i in 1:(length(w)-1)
         if res[i] == 0
-            wcross = [wcross; w[i]]
-            tcross = [tcross; i]
+            push!(wcross, w[i])
+            push!(tcross, i)
         elseif n[i] != n[i+1]
             #Interpolate to approximate crossing
             t = res[i]/(res[i]-res[i+1])
-            tcross = [tcross; i+t]
+            push!(tcross, i+t)
             wt = w[i]+t*(w[i+1]-w[i])
-            wcross = [wcross; wt]
+            push!(wcross, wt)
         end
     end
     if res[end] == 0 #Special case if multiple points
-        wcross = [wcross; w[end]]
-        tcross = [tcross; length(w)]
+        push!(wcross, w[end])
+        push!(tcross, length(w))
     end
     wcross, tcross
 end
diff --git a/lib/ControlSystemsBase/test/test_analysis.jl b/lib/ControlSystemsBase/test/test_analysis.jl
index d29668e0b..444c46a6e 100644
--- a/lib/ControlSystemsBase/test/test_analysis.jl
+++ b/lib/ControlSystemsBase/test/test_analysis.jl
@@ -229,6 +229,17 @@ dm = delaymargin(P)[]
 
 @test delaymargin(tf(0.1, [1, 1])) == Inf
 
+# https://github.com/JuliaControl/ControlSystems.jl/issues/941
+C = 0.6 + 30 * tf([1, 0], [1])
+G = tf([0.04, 0.0001, 1.1], [1, 0.03, 254.9])
+H = tf([0.25], [1, 1, 0.25])
+L = C * G * H * delay(1)
+m = margin(L)
+Lw = freqresp(L, m[1][])[]
+@test imag(Lw) ≈ 0 atol = 1e-6 # Test definition of gain margin
+@test inv(-real(Lw)) ≈ m[2][] atol = 1e-6 # Test definition of gain margin
+
+
 # RGA
 a = 10
 P = ss([0 a; -a 0], I(2), [1 a; -a 1], 0)