updating syntaxe to julia 1.0

README.md

@@ -90,7 +90,7 @@ plot(time_result)
 ![Plot real part of acoustic field](example/intro/plot_time_result.png)
 Or for a Gaussian impulse in time:
 ```julia
-t_vec = linspace(0.,700.,400)
+t_vec = LinRange(0.,700.,400)
 time_result = frequency_to_time(result; t_vec = t_vec, impulse = GaussianImpulse(max_ω))
 plot(time_result)
 ```

REQUIRE

@@ -1,6 +1,10 @@
-julia 0.5
+julia 1.0
 RecipesBase
 ProgressMeter
 SpecialFunctions
+Printf
+LinearAlgebra
+Random
+Statistics
 StaticArrays
 OffsetArrays

docs/src/example/box_size/box_size.jl

@@ -2,7 +2,7 @@ using MultipleScattering
 using Plots; pyplot()
 
 function box_size_convergence(m=4, volfrac = 0.05,
-    radius = 1.0, times = [20.0,30.0,40.0,50.0,60.0,70.0], k_arr=collect(linspace(0.01,1.0,100)) )
+    radius = 1.0, times = [20.0,30.0,40.0,50.0,60.0,70.0], k_arr=collect(LinRange(0.01,1.0,100)) )
 
     listener_position = [-10.0,0.0]
     bigshape = TimeOfFlight(listener_position,maximum(times))
@@ -17,12 +17,12 @@ function box_size_convergence(m=4, volfrac = 0.05,
         return FrequencySimulation(particles, k_arr; seed=seed, hankel_order=m)
     end
 
-    return map(m->TimeSimulation(m;time_arr=linspace(0.0,100.0,201)),simulations)
+    return map(m->TimeSimulation(m;time_arr=LinRange(0.0,100.0,201)),simulations)
 end
 
 function plot_box_size_convergence(simulations;times = [20,30,40,50,60,70])
 
-    colors = linspace(RGB(0.6,1,0.6),RGB(0,0.4,0),length(times))
+    colors = LinRange(RGB(0.6,1,0.6),RGB(0,0.4,0),length(times))
 
     plot(xlab="Time (t)",ylab="Response")
     for s in eachindex(simulations)

docs/src/example/gaussian_impulse/gaussian_impulse.jl

@@ -3,7 +3,7 @@ using Plots
 pyplot()
 
 # A very sharp sinc function (dicrete delta)
-w_arr = collect(linspace(0.0,1.0,1000))
+w_arr = collect(LinRange(0.0,1.0,1000))
 x0 = 10.;
 t_arr = 4.:0.01:22;
 fs = reshape([exp(im*w*x0) for w in w_arr],length(w_arr),1)

docs/src/example/hankel_convergence/convergence.jl

@@ -2,15 +2,15 @@ using MultipleScattering
 using Plots; pyplot()
 
 function hankel_order_convergence(m=[0,1,2,3,4,5,6,7,8,9,10], volfrac = 0.1,
-    radius = 1.0, maxtime = 40.0, k_arr=collect(linspace(0.01,1.0,100)) )
+    radius = 1.0, maxtime = 40.0, k_arr=collect(LinRange(0.01,1.0,100)) )
 
     listener_position = [-10.0,0.0]
     shape = TimeOfFlight(listener_position,maxtime)
 
     seed = MersenneTwister(1).seed
     particles = random_particles(volfrac, radius, shape; seed = seed)
 
-    simulations = Vector{FrequencySimulation{Float64}}(length(m))
+    simulations = Vector{FrequencySimulation{Float64}}(undef,length(m))
 
     for i = eachindex(m)
         simulations[i] = FrequencySimulation(particles, k_arr; seed=seed, hankel_order=m[i])
@@ -20,10 +20,10 @@ function hankel_order_convergence(m=[0,1,2,3,4,5,6,7,8,9,10], volfrac = 0.1,
 end
 
 function plot_hankel_order_convergence(simulations)
-    responses = Vector{Vector{Complex{Float64}}}(length(simulations))
-    m = Vector{Int64}(length(simulations))
+    responses = Vector{Vector{Complex{Float64}}}(undef,length(simulations))
+    m = Vector{Int64}(undef,length(simulations))
 
-    labels = Matrix{String}(1,0)
+    # labels = Matrix{String}(undef,1,0)
     for i = eachindex(simulations)
         responses[i] = reshape(simulations[i].response, size(simulations[i].response,1))
         m[i] = simulations[i].hankel_order
@@ -33,9 +33,9 @@ function plot_hankel_order_convergence(simulations)
     error = [r .- responses[end] for r in responses[1:end-1]]
     integrated_error = norm.(error).*map(m->((m.k_arr[end]-m.k_arr[1])/length(m.k_arr)),simulations[1:end-1])
 
-    colors = reshape(linspace(RGB(0.6,1,0.6),RGB(0,0.4,0),length(m)),1,length(m))
-    realcolors = reshape(linspace(RGB(0.6,0.6,1),RGB(0,0,0.4),length(m)),1,length(m))
-    imagcolors = reshape(linspace(RGB(1,0.6,0.6),RGB(0.4,0,0),length(m)),1,length(m))
+    colors = reshape(LinRange(RGB(0.6,1,0.6),RGB(0,0.4,0),length(m)),1,length(m))
+    realcolors = reshape(LinRange(RGB(0.6,0.6,1),RGB(0,0,0.4),length(m)),1,length(m))
+    imagcolors = reshape(LinRange(RGB(1,0.6,0.6),RGB(0.4,0,0),length(m)),1,length(m))
 
     absvec(v) = abs.(v)
     plot(

docs/src/example/lens/lens.jl

@@ -9,7 +9,7 @@ function run_lens(;
         innertime=34.0
     )
 
-    k_arr=collect(linspace(0.01,2.0,100))
+    k_arr=collect(LinRange(0.01,2.0,100))
 
     # Generate particles which are at most outertime away from our listener
     outershape = TimeOfFlight(listener_position,outertime)

docs/src/example/moments/README.md

@@ -27,7 +27,7 @@ plot_shape = annotate!([(listener_position[1], listener_position[2] -2., "Receiv
 ## Calculate the moments of the scattered wave
 The code below chooses a random (uniform distribution) configuration of particles inside `shape` and calculates the received signal at `listener` for wavenumbers `k_arr`,
 ```julia
-k_arr = collect(linspace(0.01,1.0,100))
+k_arr = collect(LinRange(0.01,1.0,100))
 simulation = FrequencySimulation(volfrac,radius,k_arr; shape=shape, listener_positions = listener, seed = 1)
 plot(simulation)
 ```

docs/src/example/moments/moments.jl

@@ -8,12 +8,12 @@ function moments_example()
     volfrac = 0.01
     radius = 1.0
     num_particles = 10
-    k_arr = collect(linspace(0.01,1.0,100))
+    k_arr = collect(LinRange(0.01,1.0,100))
 
     # region to place the particles
     shape = Rectangle(volfrac, radius, num_particles)
     # Holder for our simulations
-    simulations = Vector{FrequencySimulation{Float64}}(10)
+    simulations = Vector{FrequencySimulation{Float64}}(undef,10)
     for i=1:10
         simulations[i] = FrequencySimulation(volfrac,radius,k_arr;seed=[0x7f5def91, 0x82255da3, 0xc5e461c7, UInt32(i)])
     end

docs/src/example/near_surface_backscattering/backscattering.jl

@@ -48,7 +48,7 @@ widths = [round(s.shape.topright[1]) for s in simulations]
 listener_position = simulations[1].listener_positions[:]
 times = 2*(widths .- listener_position[1]) # time if takes for an incident plane wave to reach the furthest particles and then return to the receiver
 
-# [Int.(round.(linspace(1,length(num_particles)-1,5))); length(num_particles)]
+# [Int.(round.(LinRange(1,length(num_particles)-1,5))); length(num_particles)]
 plot()
 for i in [1,3,6,9,12,13]
     plot!(time_simulations[i],label="$(num_particles[i]) particles"

docs/src/example/random_particles/README.md

@@ -15,7 +15,7 @@ which wavenumbers (k) to evaluate at
 ```julia
 volfrac = 0.01
 radius = 1.0
-k_arr = collect(linspace(0.01,1.0,100))
+k_arr = collect(LinRange(0.01,1.0,100))
 simulation = FrequencySimulation(volfrac,radius,k_arr)
 ```
 

docs/src/example/random_particles/random_particles.jl

@@ -10,7 +10,7 @@ using MultipleScattering
 # which wavenumbers (k) to evaluate at
 volfrac = 0.01
 radius = 1.0
-k_arr = collect(linspace(0.01,1.0,100))
+k_arr = collect(LinRange(0.01,1.0,100))
 simulation = FrequencySimulation(volfrac,radius,k_arr)
 
 # We use the `Plots` package to plot both the response at the listener position

docs/src/example/random_particles_in_circle/README.md

@@ -5,7 +5,7 @@ The code [particles_in_circle.jl](particles_in_circle.jl) compares the scattered
 ```julia
 using MultipleScattering
 
-k_arr = collect(linspace(0.1,1.0,10))
+k_arr = collect(LinRange(0.1,1.0,10))
 
 # You can also pick your own shape, an generate random particles inside it
 # with a certain radius ands volume fraction

docs/src/example/random_particles_in_circle/particles_in_circle.jl

@@ -1,6 +1,6 @@
 using MultipleScattering
 
-k_arr = collect(linspace(0.1,1.0,10))
+k_arr = collect(LinRange(0.1,1.0,10))
 
 # You can also pick your own shape, an generate random particles inside it
 # with a certain radius ands volume fraction

docs/src/example/time_response_single_particle/time_response_single_particle.jl

@@ -2,20 +2,18 @@ using MultipleScattering
 using Plots
 
 function run_time_response_single_particle(;
-        k_arr = collect(linspace(0.001,50.0,2000)),
+        k_arr = collect(LinRange(0.001,50.0,2000)),
         particle_x = 100.0
     )
 
-    # Vector of one particle
-    particles = Vector{Particle{Float64}}(1)
     # Radius of particle
     radius = 1.0
     # Phase speed inside the particle
     c = 0.5 + 0.0*im
     # Density of particle
     ρ = 10.0
-    # Define positions and radii
-    particles[1] = Particle{Float64}([particle_x,0.0],radius,c,ρ)
+    # Define positions and radii of particle
+    particles = [Particle{Float64}([particle_x,0.0],radius,c,ρ)]
     # Simulate a single particle in frequency space
     freq_simulation = FrequencySimulation(particles,k_arr; hankel_order = 10)
     # Convert the frequency simulation into a time simulation

docs/src/example/transducers/transducer_scatterer.jl

@@ -8,7 +8,7 @@ a2_host = Acoustic(1.0,1.0 + 0.0im,2)
 # Create a finite transducer by adding point sources
 n = 100
 amp = 10.0/n
-transducer = sum(point_source(a2_host, [-2.,y], amp) for y in linspace(-2.,2.,n))
+transducer = sum(point_source(a2_host, [-2.,y], amp) for y in LinRange(-2.,2.,n))
 sim = FrequencySimulation(a2_host, transducer)
 
 using Plots; gr()
@@ -22,7 +22,7 @@ timres = frequency_to_time(simres);
 
 plot(timres, 10., seriestype=:contour)
 
-ts = filter(t -> t<16, transpose(timres.t))
+ts = filter(t -> t<16, timres.t)
 anim = @animate for t in ts
     plot(timres,t,seriestype=:contour, clim=(0.,1.25), c=:balance)
 end

example/plot/README.md

@@ -27,7 +27,7 @@ bottomleft = [-25.,-max_width]
 bounds = Rectangle(bottomleft,topright)
 result = run(simulation, bounds, [ω]; res=100)
 
-ts = linspace(0.,2pi/ω,30)
+ts = LinRange(0.,2pi/ω,30)
 
 maxc = round(10*maximum(real.(field(result))))/10
 minc = round(10*minimum(real.(field(result))))/10

example/random_particles/README.md

@@ -33,7 +33,7 @@ x = [-10.,0.]
 host_medium = Acoustic(1.0, 1.0, 2)
 source =  plane_source(host_medium; position = x, direction = [1.0,0.])
 
-ωs = linspace(0.01,1.0,100)
+ωs = LinRange(0.01,1.0,100)
 
 simulation = FrequencySimulation(host_medium, particles, source)
 result = run(simulation, x, ωs)

example/random_particles/random_particles.jl

@@ -25,7 +25,7 @@ x = [-10.,0.]
 host_medium = Acoustic(2; ρ=1.0, c=1.0)
 source =  plane_source(host_medium; position = x, direction = [1.0,0.])
 
-ωs = linspace(0.01,1.0,100)
+ωs = LinRange(0.01,1.0,100)
 
 simulation = FrequencySimulation(host_medium, particles, source)
 result = run(simulation, x, ωs)

example/random_particles_in_circle/README.md

@@ -50,7 +50,7 @@ Resulting in the figures:
 If we compare the response measured at the listener `[-10., 0.]`, they should be very similar:
 ```julia
 # define angular frequency range
-ωs = collect(linspace(0.1,1.0,10))
+ωs = collect(LinRange(0.1,1.0,10))
 result = run(simulation, x, ωs)
 big_result = run(big_particle_simulation, x, ωs)
 

example/random_particles_in_circle/particles_in_circle.jl

@@ -23,7 +23,7 @@ source =  plane_source(host_medium; position = x, direction = [1.0,0.])
 simulation = FrequencySimulation(host_medium, particles, source)
 
 # define angular frequency range
-ωs = collect(linspace(0.1,1.0,10))
+ωs = collect(LinRange(0.1,1.0,10))
 result = run(simulation, x, ωs)
 
 big_particle = Particle(particle_medium, circle)
@@ -52,7 +52,7 @@ savefig("plot_field.png")
 
 pyplot(leg=false, size=(1.8*height,height))
 
-ωs = collect(linspace(0.1,1.0,10))
+ωs = collect(LinRange(0.1,1.0,10))
 result = run(simulation, x, ωs)
 big_result = run(big_particle_simulation, x, ωs)
 

src/MultipleScattering.jl

@@ -25,12 +25,13 @@ export t_matrix, get_t_matrices
 export scattering_matrix
 
 import SpecialFunctions: besselj, hankelh1
+import Printf: @printf
 
 import StaticArrays: SVector
 
 import OffsetArrays: OffsetArray
 
-using RecipesBase
+using Random, LinearAlgebra, RecipesBase, Statistics
 using ProgressMeter
 
 # Generic machinery common to all physical models

src/physics/acoustics/acoustics.jl

@@ -37,7 +37,7 @@ include("boundary_data.jl")
 function basis_function(medium::Acoustic{T,2}, ω::T) where {T}
     return function acoustic_basis_function(m::Integer, x::SVector{2,T})
         r = norm(x)
-        θ = atan2(x[2],x[1])
+        θ = atan(x[2],x[1])
         k = ω/medium.c
         hankelh1(m,k*r)*exp(im*θ*m)
     end
@@ -47,7 +47,7 @@ end
 function basis_function(p::Particle{T,2,Acoustic{T,2}}, ω::T) where {T}
     return function acoustic_basis_function(m::Integer, x::SVector{2,T})
         r = norm(x)
-        θ = atan2(x[2],x[1])
+        θ = atan(x[2],x[1])
         k = ω/p.medium.c
         besselj(m,k*r)*exp(im*θ*m)
     end

src/physics/acoustics/boundary_data.jl

@@ -21,15 +21,15 @@ function boundary_data(shape::Shape{T,2}, inside_medium::Acoustic{T,2}, outside_
     in_pressure = run(sim, inside_points, ωs; kws...)
 
     fields = (in2_results.field - in1_results.field)/(dr * in_m.ρ)
-    in_displace = FrequencySimulationResult(fields, inside_points, RowVector(ωs))
+    in_displace = FrequencySimulationResult(fields, inside_points, Vector(ωs))
 
     # results from just outside particles
     out1_results = run(sim, outside1_points, ωs; kws...)
     out2_results = run(sim, outside2_points, ωs; kws...)
     out_pressure  = run(sim, outside_points, ωs; kws...)
 
     fields = (out2_results.field - out1_results.field)/(dr * out_m.ρ)
-    out_displace = FrequencySimulationResult(fields, outside_points, RowVector(ωs))
+    out_displace = FrequencySimulationResult(fields, outside_points, Vector(ωs))
 
     return ([in_pressure, out_pressure], [in_displace, out_displace])
 end

src/physics/acoustics/source.jl

@@ -5,7 +5,7 @@ function besselj_field(source::Source{Acoustic{T,2},T}, medium::Acoustic{T,2}, c
     centre = SVector{2,T}(centre)
 
     return (x,ω) -> sum(
-        source.coef(n,centre,ω)*besselj(n,ω/medium.c*norm(x - centre))*exp(im*n*atan2(x[2] - centre[2],x[1] - centre[1]))
+        source.coef(n,centre,ω)*besselj(n,ω/medium.c*norm(x - centre))*exp(im*n*atan(x[2] - centre[2],x[1] - centre[1]))
     for n = -basis_order:basis_order)
 
 end
@@ -15,7 +15,7 @@ end
 
 Create 2D [`Acoustic`](@ref) point [`Source`](@ref) (zeroth Hankel function of first type)
 """
-function point_source{T}(medium::Acoustic{T,2}, source_position, amplitude::Union{T,Complex{T}} = one(T))::Source{Acoustic{T,2},T}
+function point_source(medium::Acoustic{T,2}, source_position, amplitude::Union{T,Complex{T}} = one(T))::Source{Acoustic{T,2},T} where T <: AbstractFloat
 
     # Convert to SVector for efficiency and consistency
     source_position = SVector{2,T}(source_position)
@@ -25,7 +25,7 @@ function point_source{T}(medium::Acoustic{T,2}, source_position, amplitude::Unio
     function source_coef(n,centre,ω)
         k = ω/medium.c
         r = norm(centre - source_position)
-        θ = atan2(centre[2]-source_position[2], centre[1]-source_position[1])
+        θ = atan(centre[2]-source_position[2], centre[1]-source_position[1])
         # using Graf's addition theorem
         return (amplitude*im)/4 * hankelh1(-n,k*r) * exp(-im*n*θ)
     end
@@ -61,7 +61,7 @@ function plane_source(medium::Acoustic{T,2}, position, direction = SVector(one(T
 
     function source_coef(n,centre,ω)
         # Jacobi-Anger expansion
-        θ = atan2(direction[2],direction[1])
+        θ = atan(direction[2],direction[1])
         source_field(centre,ω) * exp(im * n *(T(pi)/2 -  θ))
     end
 

src/plot/plot.jl

@@ -15,7 +15,7 @@ include("plot_field.jl")
         @series begin
             label --> "$apply x=$(simres.x[x_ind])"
             xlabel --> xlab
-            (transpose(getfield(simres, 3)[ω_indices]), apply_field)
+            (getfield(simres, 3)[ω_indices], apply_field)
         end
     end
 end

src/plot/plot_field.jl

@@ -1,7 +1,7 @@
 # Plot the result in space (across all x) for a specific angular frequency
 @recipe function plot(simres::FrequencySimulationResult, ω::AbstractFloat;
         time = 0.0, # does not except "t" as name of variable..
-        x_indices = indices(simres.x,1),
+        x_indices = axes(simres.x,1),
         ω_index = findmin(abs.(getfield(simres, 3) .- ω))[2],
         field_apply = real, seriestype = :surface)
 
@@ -35,7 +35,7 @@ end
 
 # Plot the result in space (across all x) for a specific angular frequency
 @recipe function plot(timres::TimeSimulationResult, t::AbstractFloat;
-        x_indices = indices(timres.x,1),
+        x_indices = axes(timres.x,1),
         t_index = findmin(abs.(getfield(timres, 3) .- t))[2],
         field_apply = real, seriestype = :surface)