Wrote plot(Simulation, linetype=:contour) #20

Am struggling to get the right dispatch behaviour...

src/plot/plot.jl

@@ -1,6 +1,42 @@
 include("plot_domain.jl")
 # include("plot_moments.jl")
 
+"Plot just the particles and source"
+@recipe function plot(simres::SimulationResult, linetype = :line, field_apply = real)
+
+    if linetype == :line
+        @series begin
+            labs = [ "real x=$x" for x in simres.x];
+            labs = [labs ; [ "imag x=$x" for x in simres.x]]
+            ωt = getfield(simres, fieldnames(simres)[end])
+
+            labels --> reshape(labs,1,length(labs))
+            (transpose(ωt), transpose([real.(field(simres)); imag.(field(simres))]))
+        end
+
+    elseif linetype == :contour || linetype == :field
+        @series begin
+            ind_ωt = 1 # choose which time or frequency
+
+            x_pixels = union([x[1] for x in field_sim.x])
+            y_pixels = union([x[2] for x in field_sim.x])
+
+            # Make the vector field(field_sim)[:,ind_ωt] into a matrix for heatmap
+            field_mat = transpose(
+                reshape(field(simres)[:,ind_ωt], (length(x_pixels), length(y_pixels)))
+            )
+            linetype --> :contour
+            fill --> true
+            aspect_ratio := 1.0
+            fillcolor --> :pu_or
+            # title --> "Field at ω=$ω"
+
+            (x_pixels, y_pixels, field_apply.(field_mat))
+        end
+    else error("Unknown linestyle = $linestyle for $res")
+    end
+end
+
 "Plot just the particles and source"
 @recipe function plot(sim::FrequencySimulation; bounds = :auto,
                          drawparticles=true, drawsource=true)
@@ -18,17 +54,16 @@ end
 
 "Plot the field for a particular wavenumber"
 @recipe function plot(sim::FrequencySimulation, ω::Number; res=10, xres=res, yres=res,
-                         field_apply=real, build_field=true, bounds = :auto,
+                         field_apply=real, bounds = :auto, build_field = true,
                          drawparticles=false)
 
     @series begin
         # find a box which covers everything
         if bounds == :auto bounds = bounding_rectangle(sim.particles) end
 
         if build_field
-          field_sim = run(sim, bounds, [ω]; xres=xres, yres=yres)
-        else
-          field_sim = sim
+            field_sim = run(sim, bounds, [ω]; xres=xres, yres=yres)
+        else field_sim = sim
         end
 
         xy_mat = reshape(field_sim.x, (xres+1, yres+1))

src/result.jl

@@ -22,7 +22,7 @@ struct TimeSimulationResult{T<:AbstractFloat,Dim,FieldDim} <: SimulationResult{T
 end
 
 function TimeSimulationResult(time_field::Union{Matrix{T},Matrix{AbstractVector{T}}}, x::AbstractVector{SVector{Dim,T}}, t::AbstractVector{T}) where {Dim,T}
-    
+
     time_field = [SVector(d...) for d in time_field]
     FieldDim = size(time_field[1],1)
     TimeSimulationResult{T,Dim,FieldDim}(time_field, Vector(x), RowVector(t))

src/simulation.jl

@@ -35,21 +35,27 @@ function run(sim::FrequencySimulation{T,Dim,P}, x_vec::Vector{SVector{Dim,T}}, 
 
 end
 
-function run(sim::FrequencySimulation{T,Dim,P}, x_vec::Vector{SVector{Dim,T}}, ωs::AbstractVector{T} = T[]; ts::AbstractVector{T} = T[], result_in_time = !isempty(ts),
+function run(sim::FrequencySimulation{T,Dim,P},x_vec::Vector{SVector{Dim,T}},ωs::AbstractArray{T}=T[];
+        ts::AbstractArray{T} = T[], result_in_time = !isempty(ts),
         kws...)::(SimulationResult{T,Dim,FieldDim} where FieldDim)  where {Dim,P,T}
 
     if isempty(ωs) ωs = t_to_ω(ts) end
 
+    # ugly bit of code to seperate keywords for simulating frequencies
+    ks = [:basis_order]
+    freq_kws = filter(k -> contains(==,ks,k[1]), kws)
+    time_kws = setdiff(kws,freq_kws)
+
     # if user asks for ω = 0, then we provide
     if first(ωs) == zero(T)
         # Compute for each angular frequency, then join up all the results
-        fields = mapreduce(ω->run(sim,x_vec,ω; kws...).field, hcat, ωs[2:end])
+        fields = mapreduce(ω->run(sim,x_vec,ω; freq_kws...).field, hcat, ωs[2:end])
 
         # extrapolate field at ω = 0, which should be real when the time signal is real
         zeroresponse = (ωs[3].*fields[:,1] - ωs[2].*fields[:,2])./(ωs[3]-ωs[2])
         fields = reshape([real.(zeroresponse); fields[:]], length(zeroresponse), size(fields,2)+1)
     else
-        fields = mapreduce(ω->run(sim,x_vec,ω; kws...).field, hcat, ωs)
+        fields = mapreduce(ω->run(sim,x_vec,ω; freq_kws...).field, hcat, ωs)
     end
 
     if !result_in_time
@@ -58,7 +64,7 @@ function run(sim::FrequencySimulation{T,Dim,P}, x_vec::Vector{SVector{Dim,T}}, 
         if isempty(ts) ts = ω_to_t(ωs) end
 
         # better to use the defaults of TimeSimulationResult's Constructor.
-        TimeSimulationResult(FrequencySimulationResult(fields,x_vec,RowVector(ωs)); t_vec = ts, kws...)
+        TimeSimulationResult(FrequencySimulationResult(fields,x_vec,RowVector(ωs)); t_vec = reshape(ts,length(ts)), time_kws...)
     end
 end
 

src/time_simulation.jl

@@ -132,7 +132,7 @@ function frequency_to_time(field_mat::AbstractArray{Complex{T}}, ω_vec::Abstrac
         fs = impulse_vec[:].*field_mat[:,j].*exp.(-(im*t).*ω_vec)
         if method == :dft && ω_vec[1] == zero(T)
             fs[1] = fs[1]/T(2) # so as to not add ω=0 tωice in the integral of ω over [-inf,inf]
-        end    
+        end
         fs
     end
     inverse_fourier_integral = (t,j) -> numerical_integral(ω_vec, f(t,j); method=method)

test/time_tests.jl

@@ -14,21 +14,27 @@ sim = FrequencySimulation(sound_sim, particles, source)
 ω_vec = 0.0:0.01:5.01
 @test ω_vec == t_to_ω(ω_to_t(ω_vec)) # only exact for length(ω_vec) = even number
 
-x_vec = [SVector(0.0,0.0), SVector(3.0,0.0)]
 x_vec = [SVector(0.0,0.0)]
+x_vec = [SVector(0.0,0.0), SVector(3.0,0.0)]
+ω_vec = 0.0:0.1:1.01
 simres = run(sim, x_vec, ω_vec)
+timres = run(sim, x_vec, ω_vec; result_in_time=true, impulse = delta_freq_impulse, method=:trapezoidal)
+
+using Plots; pyplot()
+
+bounds = Rectangle([-2.,2.], [3.,4.])
+timres = run(sim, bounds, ω_vec; result_in_time=true)
 
-# timres = run(sim, x_vec; ts = ω_to_t(ω_vec))
+plot(timres, linetype=:contour)
 
 # timres = TimeSimulationResult(simres; t_vec = 0.0:0.2:50., method=:trapezoidal);
-timres = TimeSimulationResult(simres; impulse = delta_freq_impulse, method=:trapezoidal);
+timres2 = TimeSimulationResult(simres; impulse = delta_freq_impulse, method=:trapezoidal);
 d1 = transpose(field(timres));
 
 timres = TimeSimulationResult(simres; method=:dft, impulse = delta_freq_impulse);
 d2 = transpose(field(timres));
 norm(d1 - d2)/norm(d1)
 
-# using Plots; pyplot()
 # plot(timres.t', [d1 d2])
 
 simres2 = FrequencySimulationResult(timres; method=:dft, impulse = delta_freq_impulse)