diff --git a/src/BallAndBeam.jl b/src/BallAndBeam.jl
index 55692e7daf7fa47193e96e4fa7e5ed557560aa0e..ab87053654472b0bd9eafe1d50d205e140ba8905 100644
--- a/src/BallAndBeam.jl
+++ b/src/BallAndBeam.jl
@@ -26,9 +26,9 @@ using LabProcesses, Plots, Polynomials, ControlSystems, ProgressMeter
 	run_experiment(P::AbstractBeam, ω, duration, settling_time, amplitude, bias)
 Perform fra-experiemnt For a single frequency `ω`. Called from inside `fra`
 """
-function run_experiment(P::AbstractBeam, ω, duration, settling_time, amplitude)
-	h = sampletime(P)
-	data = zeros(0:h:duration)
+function run_experiment(P::Beam, ω, duration, settling_time, amplitude)
+	h       = sampletime(P)
+	data    = zeros(0:h:duration)
 	LabProcesses.initialize(P)
 	for t = 0:h:settling_time
 		@periodically h begin
@@ -47,6 +47,23 @@ function run_experiment(P::AbstractBeam, ω, duration, settling_time, amplitude)
 	data
 end
 
+function run_experiment(P::BeamSimulator, ω, duration, settling_time, amplitude)
+	h       = sampletime(P)
+	data    = zeros(0:h:duration)
+	LabProcesses.initialize(P)
+	for t = 0:h:settling_time
+		u = amplitude*sin(ω*t)
+		control(P, u + bias(P))
+	end
+	for (i,t) = enumerate(0:h:duration)
+		data[i] = measure(P)
+		u       = amplitude*sin(ω*t)
+		control(P, u + bias(P))
+	end
+	LabProcesses.finalize(P)
+	data
+end
+
 """
 	∫ = integrate(fun::Function, data::Vector, ω::Real,h)
 Calculate sin or cos channel, call like this `integrate(sin, data, ω,h)` for sin-channel
@@ -91,12 +108,14 @@ function sortfqs(G)
 end
 
 function bopl!(G; kwargs...)
-	plot!(G[:,1].*(180/π), abs.(G[:,2]); ylabel="Magnitude", title="Bode plot", xscale=:log10, yscale=:log10, subplot=1, kwargs...)
-	plot!(G[:,1].*(180/π), angle.(G[:,2]).*(180/π); xlabel="Frequency [1/s]", ylabel="Phase", xscale=:log10, subplot=2, kwargs...)
+	plot!(G[:,1], abs.(G[:,2]); ylabel="Magnitude", title="Bode plot", xscale=:log10, yscale=:log10, subplot=1, kwargs...)
+	phase = angle.(G[:,2]).*(180/π)
+	phase[phase .> 0] .-= 360
+	plot!(G[:,1], phase; xlabel="Frequency [rad/s]", ylabel="Phase", xscale=:log10, subplot=2, kwargs...)
     plot!(link=:x)
 end
 
-nypl!(G; kwargs...) = plot!(real.(G[:,2]), imag.(G[:,2]),
+nypl!(G; kwargs...) = plot!(real.(G[:,2]), imag.(G[:,2]);
 	xlabel="\$Re(G)\$", ylabel="\$Im(G)\$", title="Nyquist plot", kwargs...)
 
 """
diff --git a/test/runtests.jl b/test/runtests.jl
index 22258e3260d9a9adc9d1179ceb0304faf6e114ec..b4f2b660f65094ca00ea4518a7c2dbb6c4b69f49 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -47,3 +47,24 @@ gain     = 1
 sysFBc,L,T,C = fbdesign(G12, polevect, zerovect, gain)
 sysFFc,YR,FF = ffdesign(T, polevect, zerovect, gain)
 sysFB,sysFF  = c2d(sysFBc,h),c2d(sysFFc,h)
+
+
+
+
+# Test BeamSimulator
+bias = 0.01                     # Change this if your process drifts over time
+P    = LabProcesses.BeamSimulator(0.001)  # Replace for BeamSimulator to run simulations
+h    = sampletime(P)
+
+settling_time  = 2
+nbr_of_periods = 2
+
+w1_100 = logspace(-1,log10(300),500)
+G1     = fra(P, w1_100, amplitude=1, nbr_of_periods=nbr_of_periods, settling_time=settling_time)
+
+true_resp = freqresp(P.sys, w1_100)
+@test sum(abs, log.(abs.(G1[:,2])) - log.(abs.(true_resp[1][:]))) < 3.2 # Some numerical errors expected
+phase_id = angle.(G1[:,2]) |> ControlSystems.unwrap
+phase_true = angle.(true_resp[1][:]) |> ControlSystems.unwrap
+@test sum(abs, (phase_id-phase_true)[1:400]) < 0.065
+@test sum(abs, (phase_id-phase_true)[400:end]) < 5.3 # Some numerical errors for high freqs