diff --git a/README.md b/README.md
index 9668d24f..f43bd91f 100644
--- a/README.md
+++ b/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)
 ```
diff --git a/REQUIRE b/REQUIRE
index 5f9fafa3..8cc0c8e3 100644
--- a/REQUIRE
+++ b/REQUIRE
@@ -1,6 +1,10 @@
-julia 0.5
+julia 1.0
 RecipesBase
 ProgressMeter
 SpecialFunctions
+Printf
+LinearAlgebra
+Random
+Statistics
 StaticArrays
 OffsetArrays
diff --git a/docs/src/example/box_size/box_size.jl b/docs/src/example/box_size/box_size.jl
index 656811be..8fcf0d5c 100644
--- a/docs/src/example/box_size/box_size.jl
+++ b/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)
diff --git a/docs/src/example/gaussian_impulse/gaussian_impulse.jl b/docs/src/example/gaussian_impulse/gaussian_impulse.jl
index d42006c6..71afc6bc 100644
--- a/docs/src/example/gaussian_impulse/gaussian_impulse.jl
+++ b/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)
diff --git a/docs/src/example/hankel_convergence/convergence.jl b/docs/src/example/hankel_convergence/convergence.jl
index cc74149a..cb0b56e6 100644
--- a/docs/src/example/hankel_convergence/convergence.jl
+++ b/docs/src/example/hankel_convergence/convergence.jl
@@ -2,7 +2,7 @@ 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)
@@ -10,7 +10,7 @@ function hankel_order_convergence(m=[0,1,2,3,4,5,6,7,8,9,10], volfrac = 0.1,
     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(
diff --git a/docs/src/example/lens/lens.jl b/docs/src/example/lens/lens.jl
index 0bd9ab8b..4f8fc217 100644
--- a/docs/src/example/lens/lens.jl
+++ b/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)
diff --git a/docs/src/example/moments/README.md b/docs/src/example/moments/README.md
index 77358459..88ec9754 100644
--- a/docs/src/example/moments/README.md
+++ b/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)
 ```
diff --git a/docs/src/example/moments/moments.jl b/docs/src/example/moments/moments.jl
index 39d9ca1a..ae17caac 100644
--- a/docs/src/example/moments/moments.jl
+++ b/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
diff --git a/docs/src/example/near_surface_backscattering/backscattering.jl b/docs/src/example/near_surface_backscattering/backscattering.jl
index 7bfd506d..41219207 100644
--- a/docs/src/example/near_surface_backscattering/backscattering.jl
+++ b/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"
diff --git a/docs/src/example/random_particles/README.md b/docs/src/example/random_particles/README.md
index e63aa3c2..89d40e86 100644
--- a/docs/src/example/random_particles/README.md
+++ b/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)
 ```
 
diff --git a/docs/src/example/random_particles/random_particles.jl b/docs/src/example/random_particles/random_particles.jl
index 0006c44e..3452c658 100644
--- a/docs/src/example/random_particles/random_particles.jl
+++ b/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
diff --git a/docs/src/example/random_particles_in_circle/README.md b/docs/src/example/random_particles_in_circle/README.md
index bee1dab8..7ea79f93 100644
--- a/docs/src/example/random_particles_in_circle/README.md
+++ b/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
diff --git a/docs/src/example/random_particles_in_circle/particles_in_circle.jl b/docs/src/example/random_particles_in_circle/particles_in_circle.jl
index c13a895b..ef522188 100644
--- a/docs/src/example/random_particles_in_circle/particles_in_circle.jl
+++ b/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
diff --git a/docs/src/example/time_response_single_particle/time_response_single_particle.jl b/docs/src/example/time_response_single_particle/time_response_single_particle.jl
index 0694764e..4178c0eb 100644
--- a/docs/src/example/time_response_single_particle/time_response_single_particle.jl
+++ b/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
diff --git a/docs/src/example/transducers/transducer_scatterer.jl b/docs/src/example/transducers/transducer_scatterer.jl
index f83c1545..22808110 100644
--- a/docs/src/example/transducers/transducer_scatterer.jl
+++ b/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
diff --git a/example/plot/README.md b/example/plot/README.md
index ffa9efe3..5dfffb9f 100644
--- a/example/plot/README.md
+++ b/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
diff --git a/example/random_particles/README.md b/example/random_particles/README.md
index aaf0def4..33734f87 100644
--- a/example/random_particles/README.md
+++ b/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)
diff --git a/example/random_particles/random_particles.jl b/example/random_particles/random_particles.jl
index d709dabd..31d9464f 100644
--- a/example/random_particles/random_particles.jl
+++ b/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)
diff --git a/example/random_particles_in_circle/README.md b/example/random_particles_in_circle/README.md
index 6c8e1ab3..19589353 100644
--- a/example/random_particles_in_circle/README.md
+++ b/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)
 
diff --git a/example/random_particles_in_circle/particles_in_circle.jl b/example/random_particles_in_circle/particles_in_circle.jl
index 3ec01370..356bcf8f 100644
--- a/example/random_particles_in_circle/particles_in_circle.jl
+++ b/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)
 
diff --git a/src/MultipleScattering.jl b/src/MultipleScattering.jl
index 22e8e2a6..31eea59a 100644
--- a/src/MultipleScattering.jl
+++ b/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
diff --git a/src/physics/acoustics/acoustics.jl b/src/physics/acoustics/acoustics.jl
index 54aa5a16..24e6ab1c 100644
--- a/src/physics/acoustics/acoustics.jl
+++ b/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
diff --git a/src/physics/acoustics/boundary_data.jl b/src/physics/acoustics/boundary_data.jl
index be5fb5ec..334cddd3 100644
--- a/src/physics/acoustics/boundary_data.jl
+++ b/src/physics/acoustics/boundary_data.jl
@@ -21,7 +21,7 @@ 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...)
@@ -29,7 +29,7 @@ function boundary_data(shape::Shape{T,2}, inside_medium::Acoustic{T,2}, outside_
     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
diff --git a/src/physics/acoustics/source.jl b/src/physics/acoustics/source.jl
index 5646f926..d061b0dc 100644
--- a/src/physics/acoustics/source.jl
+++ b/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
 
diff --git a/src/plot/plot.jl b/src/plot/plot.jl
index 7a3f951c..f19a232d 100644
--- a/src/plot/plot.jl
+++ b/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
diff --git a/src/plot/plot_field.jl b/src/plot/plot_field.jl
index 204889e5..944ce546 100644
--- a/src/plot/plot_field.jl
+++ b/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)
 
diff --git a/src/plot/plot_time.jl b/src/plot/plot_time.jl
index 59ed8979..ac895c6a 100644
--- a/src/plot/plot_time.jl
+++ b/src/plot/plot_time.jl
@@ -10,7 +10,7 @@ end
 
 # "Plot the field for an array of time"
 # @recipe function plot{T}(TimeSimulation::TimeSimulation{T}; res=res)
-#     map(indices(TimeSimulation.time_arr)) do i
+#     map(axes(TimeSimulation.time_arr,1)) do i
 #         simulation = deepcopy(TimeSimulation)
 #         simulation.response = reshape(TimeSimulation.response[i,:],1,:)
 #         plot(simulation,TimeSimulation.time_arr[i]; res=res, build_field=false)
@@ -46,8 +46,8 @@ end
         end
 
         # For this we sample at the centre of each pixel
-        x_pixels = linspace(bounds.bottomleft[1], bounds.topright[1], xres+1)
-        y_pixels = linspace(bounds.bottomleft[2], bounds.topright[2], yres+1)
+        x_pixels = LinRange(bounds.bottomleft[1], bounds.topright[1], xres+1)
+        y_pixels = LinRange(bounds.bottomleft[2], bounds.topright[2], yres+1)
 
         # NOTE only plots the first time plot for now...
         response_mat = transpose(reshape(field_TimeSimulation.response[1,:], (xres+1, yres+1)))
diff --git a/src/random/moments.jl b/src/random/moments.jl
index f9139930..8f367fc2 100644
--- a/src/random/moments.jl
+++ b/src/random/moments.jl
@@ -11,7 +11,7 @@ function calculate_moments(results::AbstractVector{SimRes}, num_moments::Int; ap
     # Number of realisations
     R = length(results)
 
-    moments = [Matrix{T}(X,K) for i=1:num_moments]
+    moments = [Matrix{T}(undef,X,K) for i=1:num_moments]
 
     # For each wavenumber or timestep, calculate each moment up to num_moments
     for i = 1:X
diff --git a/src/random/random_particles.jl b/src/random/random_particles.jl
index 3991c2e7..8077f4cb 100644
--- a/src/random/random_particles.jl
+++ b/src/random/random_particles.jl
@@ -16,7 +16,7 @@ function random_particles(particle_medium::PhysicalProperties{T,Dim}, particle_s
 end
 
 function random_particles(particle_medium::P, particle_shape::S,
-    box_shape::Shape{T,Dim}, volfrac::AbstractFloat; seed=Base.Random.make_seed()) where {T,Dim,P<:PhysicalProperties{T,Dim},S<:Shape{T,Dim}}
+    box_shape::Shape{T,Dim}, volfrac::AbstractFloat; seed = Random.make_seed()) where {T,Dim,P<:PhysicalProperties{T,Dim},S<:Shape{T,Dim}}
 
     N = Int(round(volume(box_shape) / volume(particle_shape) * volfrac))
     return random_particles(particle_medium, particle_shape, box_shape, N; seed=seed)
@@ -24,9 +24,9 @@ end
 
 """
     random_particles(particle_medium, particle_shape, box_shape, volume_fraction::Number;
-        seed=Base.Random.make_seed())
+        seed=Random.make_seed())
     random_particles(particle_medium, particle_shape, box_shape, N::Integer;
-        seed=Base.Random.make_seed())
+        seed=Random.make_seed())
 
 Generate `N` random particles that fit inside `box_shape` (or fill with `volume_fraction`)
 
@@ -35,7 +35,7 @@ the bounding rectangle of `box_shape` and discards particle if it overlaps (base
 or does not lies completely in box box.
 """
 function random_particles(particle_medium::P, particle_shape::S,
-    box_shape::Shape{T,Dim}, N::Integer; seed=Base.Random.make_seed()) where {T,Dim,P<:PhysicalProperties{T,Dim},S<:Shape{T,Dim}}
+    box_shape::Shape{T,Dim}, N::Integer; seed=Random.make_seed()) where {T,Dim,P<:PhysicalProperties{T,Dim},S<:Shape{T,Dim}}
 
     # Check volume fraction is not impossible
     volfrac = N * volume(particle_shape) / volume(box_shape)
@@ -64,7 +64,7 @@ function random_particles(particle_medium::P, particle_shape::S,
     )
 
     # Allocate memory for particles
-    particles = Vector{Particle{T,Dim,P,S}}(N)
+    particles = Vector{Particle{T,Dim,P,S}}(undef,N)
 
     for n = 1:N
 
@@ -74,7 +74,7 @@ function random_particles(particle_medium::P, particle_shape::S,
 
             outside_box = true
             while outside_box
-                x = bounding_rect_size .* (1 .- 2.*rand(randgen,T,2)) .+ origin(bounding_rect)
+                x = bounding_rect_size .* (1 .- 2 .* rand(randgen,T,2)) .+ origin(bounding_rect)
                 particles[n] = Particle(particle_medium, congruent(particle_shape, x))
                 outside_box = !(particles[n] ⊆ box_shape)
             end
diff --git a/src/result.jl b/src/result.jl
index 6c57b384..016a0bd0 100644
--- a/src/result.jl
+++ b/src/result.jl
@@ -11,14 +11,14 @@ struct FrequencySimulationResult{T<:AbstractFloat,Dim,FieldDim} <: SimulationRes
     "Positions"
     x::Vector{SVector{Dim,T}}
     "Angular frequencies"
-    ω::RowVector{T}
+    ω::Vector{T}
 end
 
 function FrequencySimulationResult(field::Union{Matrix{Complex{T}}, Matrix{AbstractVector{Complex{T}}}}, x::AbstractVector{SVector{Dim,T}}, ω::AbstractVector{T}) where {Dim,T}
 
     field = [SVector(d...) for d in field]
     FieldDim = size(field[1],1)
-    FrequencySimulationResult{T,Dim,FieldDim}(field, Vector(x), RowVector(ω))
+    FrequencySimulationResult{T,Dim,FieldDim}(field, Vector(x), Vector(ω))
 end
 
 """
@@ -30,14 +30,14 @@ struct TimeSimulationResult{T<:AbstractFloat,Dim,FieldDim} <: SimulationResult{T
     "Positions"
     x::Vector{SVector{Dim,T}}
     "Times"
-    t::RowVector{T}
+    t::Vector{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))
+    TimeSimulationResult{T,Dim,FieldDim}(time_field, Vector(x), Vector(t))
 end
 
 """
diff --git a/src/scattering_matrix.jl b/src/scattering_matrix.jl
index f095747a..996b4a89 100644
--- a/src/scattering_matrix.jl
+++ b/src/scattering_matrix.jl
@@ -8,7 +8,7 @@ function scattering_matrix(medium::PhysicalProperties, particles::AbstractPartic
     # No particles means no scattering
     if P == 0
         # warn("You have computed the scattering matrix with no particles, are you sure something hasn't gone wrong?")
-        return Matrix{Complex{T}}(0,0)
+        return Matrix{Complex{T}}(undef,0,0)
     end
 
     # Number of hankel basis function at each particle
@@ -19,7 +19,7 @@ function scattering_matrix(medium::PhysicalProperties, particles::AbstractPartic
     # Faire: this could potentially return an MMatrix
     function S_block(j,l)
         if j == l
-            return eye(Complex{T}, H, H)
+            return Matrix{Complex{T}}(I, H, H)
         else
             x_lj = origin(particles[j]) .- origin(particles[l])
             # Faire: basis functions could be more efficient if it returned a vector
@@ -33,7 +33,7 @@ function scattering_matrix(medium::PhysicalProperties, particles::AbstractPartic
     S_blocks = [S_block(j,l) for j in 1:P, l in 1:P]
 
     # Reshape S_blocks into big matrix
-    S = Matrix{Complex{T}}(P*H, P*H)
+    S = Matrix{Complex{T}}(undef, P*H, P*H)
     for i in 1:P
         for j in 1:P
             S[((i-1)*H+1):(i*H), ((j-1)*H+1):(j*H)] .= S_blocks[i,j]
diff --git a/src/shapes/shapes.jl b/src/shapes/shapes.jl
index 2d17f363..902e20cd 100644
--- a/src/shapes/shapes.jl
+++ b/src/shapes/shapes.jl
@@ -47,7 +47,7 @@ include("sphere.jl")
 function boundary_points(shape::Shape{T,Dim}, num_points::Int = 4; dr = zero(T)) where {Dim,T}
     x, y = boundary_functions(shape)
     v(τ) = SVector(x(τ),y(τ)) + dr * (SVector(x(τ),y(τ)) - origin(shape))
-    return [ v(τ) for τ in linspace(zero(T),one(T),num_points+1)[1:end-1] ]
+    return [ v(τ) for τ in LinRange(zero(T),one(T),num_points+1)[1:end-1] ]
 end
 
 "Returns rectangle which completely encloses the shapes"
diff --git a/src/simulation.jl b/src/simulation.jl
index 03d63b6d..fb1531ab 100644
--- a/src/simulation.jl
+++ b/src/simulation.jl
@@ -26,7 +26,7 @@ end
 
 # A simulation with just sources is perfectly reasonable
 function FrequencySimulation(medium::P, source::Source{P,T}) where {Dim,T,P<:PhysicalProperties{T,Dim}}
-    FrequencySimulation{T,Dim,P}(medium, Vector{AbstractParticle{T,Dim}}(0), source)
+    FrequencySimulation{T,Dim,P}(medium, Vector{AbstractParticle{T,Dim}}(undef,0), source)
 end
 
 import Base.run
@@ -44,7 +44,7 @@ function run(sim::FrequencySimulation{T,Dim,P}, x_vec::Union{Vector{Vector{T}},V
 
     # Construct results object
     field_vec = reshape(map(f->SVector{FieldDim,Complex{T}}(f), field_vec), :, 1)
-    return FrequencySimulationResult{T,Dim,FieldDim}(field_vec, x_vec, RowVector([ω]))
+    return FrequencySimulationResult{T,Dim,FieldDim}(field_vec, x_vec, Vector([ω]))
 
 end
 
@@ -62,7 +62,7 @@ function run(sim::FrequencySimulation{T,Dim,P}, x_vec::Union{Vector{Vector{T}},V
     if basis_order_vec == [-1]
         max_basis_order = max(basis_order,min_basis_order)
         basis_order_vec = Int.(round.(
-            linspace(min_basis_order, max_basis_order, length(ωs))
+            LinRange(min_basis_order, max_basis_order, length(ωs))
         ))
         basis_order_vec = basis_order_vec[sortperm(ωs)]
     end
@@ -89,12 +89,12 @@ function run(sim::FrequencySimulation{T,Dim,P}, x_vec::Union{Vector{Vector{T}},V
     end
 
     if !result_in_time
-        FrequencySimulationResult(fields,x_vec,RowVector(ωs))
+        FrequencySimulationResult(fields,x_vec,Vector(ωs))
     else
         if isempty(ts) ts = ω_to_t(ωs) end
 
         # better to use the defaults of TimeSimulationResult's Constructor.
-        frequency_to_time(FrequencySimulationResult(fields,x_vec,RowVector(ωs)); t_vec = reshape(ts,length(ts)), time_kws...)
+        frequency_to_time(FrequencySimulationResult(fields,x_vec,Vector(ωs)); t_vec = reshape(ts,length(ts)), time_kws...)
     end
 end
 
@@ -145,7 +145,7 @@ Create forcing vector from source, forms the right hand side of matrix equation
 """
 function forcing(source::Source{P,T}, particles::AbstractParticles, ω::T, Nh::Integer)::Vector{Complex{T}} where {P,T}
     mat = [source.coef(n,origin(p),ω) for n in -Nh:Nh, p in particles]
-    f = Vector{Complex{T}}(prod(size(mat)))
+    f = Vector{Complex{T}}(undef,prod(size(mat)))
     H = 2Nh + 1
     for i in eachindex(particles)
         f[((i-1)*H+1):(i*H)] .= mat[:,i]
@@ -174,7 +174,7 @@ function basis_coefficients(sim::FrequencySimulation{T,Dim,P}, ω::T; basis_orde
 
     # reshape and multiply by t-matrix to get the scattering coefficients
     a = reshape(a,2basis_order+1,length(sim.particles))
-    for i in indices(a,2)
+    for i in axes(a,2)
         a[:,i] = t_matrices[i] * a[:,i]
     end
     a
@@ -196,8 +196,8 @@ function field(sim::FrequencySimulation{T,Dim,P}, ω::T, x_vec::Vector{SVector{D
         end
     end
     map(x_vec) do x
-        j = findfirst(p->x∈p, sim.particles)
-        if iszero(j)
+        j = findfirst(p -> x∈p, sim.particles)
+        if typeof(j) === Nothing
             sim.source.field(x,ω) + (isempty(sim.particles) ? zero(Complex{T}) : sum_basis(x))
         else
             p = sim.particles[j]
diff --git a/src/t_matrix.jl b/src/t_matrix.jl
index a2ec490b..203c4be2 100644
--- a/src/t_matrix.jl
+++ b/src/t_matrix.jl
@@ -15,11 +15,11 @@ vector.
 """
 function get_t_matrices(medium::PhysicalProperties, particles::AbstractParticles, ω::AbstractFloat, Nh::Integer)::Vector
 
-    t_matrices = Vector{AbstractMatrix}(length(particles))
+    t_matrices = Vector{AbstractMatrix}(undef,length(particles))
 
     # Vector of particles unique up to congruence, and the respective T-matrices
-    unique_particles = Vector{AbstractParticle}(0)
-    unique_t_matrices = Vector{AbstractMatrix}(0)
+    unique_particles = Vector{AbstractParticle}(undef,0)
+    unique_t_matrices = Vector{AbstractMatrix}(undef,0)
 
     for p_i in eachindex(particles)
         p = particles[p_i]
diff --git a/src/time_simulation.jl b/src/time_simulation.jl
index 99de9378..cf859fda 100644
--- a/src/time_simulation.jl
+++ b/src/time_simulation.jl
@@ -6,11 +6,11 @@ Assumes only positive frequencies and a real time signal
 function frequency_to_time(simres::FrequencySimulationResult{T,Dim,FieldDim};
         t_vec::AbstractVector{T} = ω_to_t(simres.ω),
         impulse::ContinuousImpulse{T} = TimeDiracImpulse(zero(T)), #GaussianImpulse(maximum(simres.ω)),
-        discrete_impulse::DiscreteImpulse{T} = continuous_to_discrete_impulse(impulse, t_vec, transpose(simres.ω)),
+        discrete_impulse::DiscreteImpulse{T} = continuous_to_discrete_impulse(impulse, t_vec,simres.ω),
         method = :dft
     ) where {Dim,FieldDim,T}
 
-    time_field = frequency_to_time(transpose(field(simres)), transpose(simres.ω), t_vec;
+    time_field = frequency_to_time(transpose(field(simres)), simres.ω, t_vec;
         discrete_impulse = discrete_impulse, method = method)
 
     return TimeSimulationResult(transpose(time_field), simres.x, t_vec)
@@ -23,11 +23,11 @@ Assumes only positive frequencies and a real time signal
 function time_to_frequency(timres::TimeSimulationResult{T,Dim,FieldDim};
         ω_vec = t_to_ω(timres.t),
         impulse::ContinuousImpulse{T} = TimeDiracImpulse(zero(T)), #GaussianImpulse(maximum(ω_vec)),
-        discrete_impulse::DiscreteImpulse{T} = continuous_to_discrete_impulse(impulse, transpose(timres.t), ω_vec),
+        discrete_impulse::DiscreteImpulse{T} = continuous_to_discrete_impulse(impulse,timres.t, ω_vec),
         method =:dft
     ) where {Dim,FieldDim,T}
 
-    freq_field = time_to_frequency(transpose(field(timres)), transpose(timres.t), ω_vec;
+    freq_field = time_to_frequency(transpose(field(timres)), timres.t, ω_vec;
         discrete_impulse = discrete_impulse, method = method)
 
     return FrequencySimulationResult(transpose(freq_field), deepcopy(timres.x), ω_vec)
@@ -46,7 +46,7 @@ function ω_to_t(ω_arr::AbstractArray{T}) where T <: AbstractFloat
         error("expected only non-negative values for the frequencies")
     end
     dω = ω_arr[2] - ω_arr[1]
-    t_arr = linspace(zero(T),2π/dω,2N+2)[1:(2N+1)]
+    t_arr = LinRange(zero(T),2π/dω,2N+2)[1:(2N+1)]
     return t_arr
 end
 
@@ -55,7 +55,7 @@ function t_to_ω(t_arr::AbstractArray{T}) where T <: AbstractFloat
     N = Int(round((length(t_arr)-one(T))/T(2)))
     maxt = t_arr[2] - t_arr[1] + t_arr[end] - t_arr[1] # subtract t_arr[1] in case t_arr[1] != zero(T)
     maxω = N*2π/maxt
-    ω_vec = linspace(zero(T),maxω,N+1)
+    ω_vec = LinRange(zero(T),maxω,N+1)
     return ω_vec
 end
 
@@ -99,7 +99,7 @@ function frequency_to_time(field_mat::AbstractArray{Complex{T}}, ω_vec::Abstrac
         fs
     end
     inverse_fourier_integral = (t,j) -> numerical_integral(ω_vec, f(t,j), method)
-    u = [inverse_fourier_integral(t,j) for t in t_vec, j in indices(field_mat,2)]
+    u = [inverse_fourier_integral(t,j) for t in t_vec, j in axes(field_mat,2)]
 
     return real.(u)/pi # a constant 1/(2pi) appears due to our Fourier convention, but because we only use positive frequencies, and assume a real time signal, this becomes 1/pi.
 end
@@ -118,7 +118,7 @@ function time_to_frequency(field_mat::AbstractArray{T}, t_vec::AbstractVector{T}
     # to use an impulse below in time we would need to do a discrete convolution, which we decided against.
     f(ω::T, j::Int) = field_mat[:,j].*exp.((im*ω).*t_vec)
     fourier_integral = (ω,j) -> numerical_integral(t_vec, f(ω,j), method)
-    uhat = [discrete_impulse.in_freq[i]*fourier_integral(ω_vec[i],j) for i in eachindex(ω_vec), j in indices(field_mat,2)]
+    uhat = [discrete_impulse.in_freq[i]*fourier_integral(ω_vec[i],j) for i in eachindex(ω_vec), j in axes(field_mat,2)]
 
     return uhat
 end
diff --git a/test/acoustic_physics_tests.jl b/test/acoustic_physics_tests.jl
index 49e15de9..fb23b4bc 100644
--- a/test/acoustic_physics_tests.jl
+++ b/test/acoustic_physics_tests.jl
@@ -31,7 +31,7 @@ end
 
     # Test return type is satisfied for valid input
     t = t_matrix(Particle(a2,circle), a2_host, ω, N_basis)
-    @test typeof(t) == Diagonal{Complex{Float64}}
+    @test typeof(t) <: Diagonal{Complex{Float64}}
 
     p = Particle(Acoustic(Inf, 0.0im, 2),circle)
     @test_throws DomainError t_matrix(p, Acoustic(1.0, 1.0+0.0im, 2), ω, N_basis)
@@ -53,7 +53,7 @@ end
     source_position = SVector(0.0,1.0)
     amplitude = 1.0
     s1 = point_source(a2, source_position, amplitude)
-    s2 = point_source(a2, 2.*source_position, amplitude)
+    s2 = point_source(a2, 2.0*source_position, amplitude)
 
     # Create new souce as a linear combination of two other sources
     s3 = 2*s1 + s2
@@ -66,7 +66,7 @@ end
     ω = 0.8
     centre =  SVector(1.0,0.0)
     s3_besselj = besselj_field(s3, a2, centre; basis_order = 7)
-    xs = [centre + 0.1.*[cos(τ),sin(τ)] for τ = 0.0:0.3:1.5]
+    xs = [centre + 0.1 .* [cos(τ),sin(τ)] for τ = 0.0:0.3:1.5]
     @test norm([s3.field(x,ω) - s3_besselj(x,ω) for x in xs]) < 1e-7*norm([s3.field(x,ω) for x in xs])
 
     source = plane_source(a2_host, SVector(-10.0,0.0), SVector(1.0,0.0), 1.0)
diff --git a/test/boundary_conditions.jl b/test/boundary_conditions.jl
index 5a589992..4970ce9c 100644
--- a/test/boundary_conditions.jl
+++ b/test/boundary_conditions.jl
@@ -6,7 +6,7 @@ function boundary_conditions_test(numberofparticles::Int=4, seed = 1 )
     particles = [Particle([0.,0.], rand(0.1:0.1:2.0), rand(0.2:0.1:10)+0.0im, rand(0.2:0.1:10)) for i=1:numberofparticles]
     shape = Rectangle(0.1, mapreduce(p->p.r,max,particles), numberofparticles)
     random_particles!(particles, shape; seed=seed) # choose random positions
-    k_arr = collect(linspace(0.01,1.0,10));
+    k_arr = collect(LinRange(0.01,1.0,10));
     simulation = FrequencySimulation(particles,k_arr);
     boundary_data= map(4:6) do m
       simulation.hankel_order = m
@@ -103,5 +103,5 @@ end
 
 "returns points on along a radial axes from the centre of the particle."
 function points_on_radial_axes{T}(p::Particle{T}, R=T(2)*p.r, θ=zero(T); numberofpoints::Int = 20)
-    return [p.x + r*[cos(θ),sin(θ)] for r in linspace(zero(T),R,numberofpoints+1)[2:end]]
+    return [p.x + r*[cos(θ),sin(θ)] for r in LinRange(zero(T),R,numberofpoints+1)[2:end]]
 end
diff --git a/test/capsule_test.jl b/test/capsule_test.jl
index 7f3b335b..86521030 100644
--- a/test/capsule_test.jl
+++ b/test/capsule_test.jl
@@ -56,7 +56,7 @@ concen = CapsuleParticle(
 )
 
 n=12
-ps = [Particle(Acoustic(0.1,0.1,2),Circle((4*cos(θ),4*sin(θ)),0.5)) for θ in linspace(0.,2pi,n+1)[1:n] ]
+ps = [Particle(Acoustic(0.1,0.1,2),Circle((4*cos(θ),4*sin(θ)),0.5)) for θ in LinRange(0.,2pi,n+1)[1:n] ]
 ps = [concen; ps]
 plot(ps)
 
diff --git a/test/fourier_tests.jl b/test/fourier_tests.jl
index 4ad01424..546f36b1 100644
--- a/test/fourier_tests.jl
+++ b/test/fourier_tests.jl
@@ -1,11 +1,13 @@
-import Base.Test: @testset, @test, @test_throws
+import Test: @testset, @test, @test_throws
 import StaticArrays: SVector
+import Statistics: mean
 
 using MultipleScattering
+using LinearAlgebra
 
 @testset "Fourier transform" begin
 
-    ω_arr = collect(linspace(0.0,100.0,200))
+    ω_arr = collect(LinRange(0.0,100.0,200))
     t_arr = ω_to_t(ω_arr)
     N = length(ω_arr) - 1
     t0 = t_arr[N]
@@ -24,18 +26,18 @@ using MultipleScattering
     @test Fs ≈ time_to_frequency(fs,t_arr,ω_arr)
 
     # sinc to rectangle
-    rect(t) = (abs(t)<0.5) ? 1.0:0.
+    rect(t) = (abs(t)<0.5) ? 1.0 : 0.
     as = [0.7]; #only positive values
-    Fs = [ (ω == 0.0im) ? 1.0+0.0im:sin(a*pi*ω)*exp(im*t0*ω)/(a*pi*ω) for ω in ω_arr, a in as]
+    Fs = [ (ω == 0.0im) ? 1.0+0.0im : sin(a*pi*ω)*exp(im*t0*ω)/(a*pi*ω) for ω in ω_arr, a in as]
     fs = frequency_to_time(Fs,ω_arr,t_arr) # only correct for half of time_arr
     true_fs = (1/(2pi))*[(1/abs(a))*rect((t-t0)/(2*pi*a)) for t in t_arr, a in as]
     @test mean(abs.(fs-true_fs)) < 0.02*mean(abs.(true_fs))
 
     # trapezoidal integration
-    ω_arr = linspace(0.0,130.0,600)
+    ω_arr = LinRange(0.0,130.0,600)
     t_arr = ω_to_t(ω_arr)
 
-    fs = cos.(t_arr).*exp.(-(t_arr - t_arr[end]/2).^2/(t_arr[end]^2/25))
+    fs = cos.(t_arr).*exp.(-(t_arr .- t_arr[end]/2).^2/(t_arr[end]^2/25))
     fs = reshape(fs,length(fs),1)
 
     Fs_trap = time_to_frequency(fs, t_arr; method=:trapezoidal)
@@ -46,22 +48,24 @@ using MultipleScattering
     fs_trap = frequency_to_time(Fs_trap,ω_arr; method=:trapezoidal)
     @test mean(abs.(fs-fs_trap))< 1e-5*mean(abs.(fs))
 
-    # rectangle to sinc
+    ## rectangle to sinc
     # as = [1/(2*(maximum(ω_arr)+ω_arr[2]))]; #only positive values
     # # t0 = t_arr[end]/3
     # Fs = Complex{Float64}[ rect(a*ω)*exp(im*t0*ω) for ω in ω_arr, a in as]
-    # fs =  frequency_to_time(Fs,ω_arr,t_arr)
+    # fs =  frequency_to_time(Fs, ω_arr, t_arr)
     # true_fs = [ (t==t0) ? 1/(2pi*a) : sin((t-t0)/(2a))/(pi*(t-t0)) for t in t_arr, a in as]
-    # mean(abs(fs-true_fs)) < 0.04*mean(abs.(true_fs))
+    # @test maximum(abs.(fs-true_fs)) < 1e-3*maximum(abs.(true_fs))
 
     # pole to wavy heaviside
-    # ω_arr = collect(linspace(0.0,10.0,200)) # need fine mesh near ω=0
+    # ω_arr = collect(LinRange(0.0,10.0,4000)) # need fine mesh near ω=0
     # dω = (ω_arr[2]-ω_arr[1])
     # N = length(ω_arr)-1;
-    # time_arr = linspace(0,2π/dω,2N+2)[1:(2N+1)]
+    # time_arr = LinRange(0,2π/dω,2N+2)[1:(2N+1)]
     # as = [0.1]; #only positive values
     # Fs = Complex{Float64}[1/(a-im*ω) for ω in ω_arr, a in as]
     # fs =  frequency_to_time(Fs,ω_arr,time_arr);
     # true_fs = [exp(-a*t) for t in time_arr, a in as]
+    # mean(abs.(fs-true_fs))
+    #  1e-3*mean(abs.(true_fs))
 
 end
diff --git a/test/runoldtests.jl b/test/runoldtests.jl
index a3ee15de..4e708ca9 100644
--- a/test/runoldtests.jl
+++ b/test/runoldtests.jl
@@ -19,7 +19,7 @@ using Plots
         # Define everything as a Float32
         volfrac = 0.01f0
         radius = 1.0f0
-        k_arr = collect(linspace(0.01f0,1.0f0,100))
+        k_arr = collect(LinRange(0.01f0,1.0f0,100))
         simulation = FrequencySimulation(volfrac,radius,k_arr)
         @test typeof(simulation.response[1]) == Complex{Float32}
     end
@@ -40,7 +40,7 @@ using Plots
         # Just run it to see if we have any errors (yes thats a very low bar)
         volfrac = 0.01
         radius = 2.0
-        k_arr = collect(linspace(0.2,1.0,5))
+        k_arr = collect(LinRange(0.2,1.0,5))
         simulation = FrequencySimulation(volfrac,radius,k_arr)
 
         plot(simulation)
diff --git a/test/runtests.jl b/test/runtests.jl
index d6539ab0..e914f4a8 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,7 +1,8 @@
-import Base.Test: @testset, @test, @test_throws
+import Test: @testset, @test, @test_throws
 import StaticArrays: SVector
 
 using MultipleScattering
+using LinearAlgebra
 using Plots
 
 include("shapetests.jl")
diff --git a/test/single_scatter.jl b/test/single_scatter.jl
index 1ab008bf..0ce2b3b7 100644
--- a/test/single_scatter.jl
+++ b/test/single_scatter.jl
@@ -5,7 +5,7 @@ function single_scatter_test()
     xp = 10*rand(2);
     a1 = 10*rand();
     a2 = 10*rand();
-    k_arr = linspace(0.01,2.0,200)
+    k_arr = LinRange(0.01,2.0,200)
     M1 = Int(round(a1)); #largest hankelh1 order
     M2 = Int(round(a2)); #largest hankelh1 order
 
@@ -15,7 +15,7 @@ function single_scatter_test()
     # known solution for the scattering coefficients for a single scattering
     function exact_single_scatter(x, k, a, M)
         A(n,k) = -exp(im * k * (xp[1] - listener[1]) + im * n * pi / 2);
-        θ = atan2(x[2] - xp[2], x[1] - xp[1])
+        θ = atan(x[2] - xp[2], x[1] - xp[1])
         r = norm(x - xp)
         u = 0.0im
         for m = -M:M
diff --git a/test/time_tests.jl b/test/time_tests.jl
index 4246cb24..bb475509 100644
--- a/test/time_tests.jl
+++ b/test/time_tests.jl
@@ -1,8 +1,5 @@
-import Base.Test: @testset, @test, @test_throws
 import StaticArrays: SVector
 
-using MultipleScattering
-
 @testset "Time Result" begin
     sound_p = Acoustic(.1, 0.1 + 0.0im,2)
     particles = [Particle(sound_p,Circle([10.5,0.0], .5))]
@@ -12,7 +9,7 @@ using MultipleScattering
     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
+    @test LinRange(ω_vec) == t_to_ω(ω_to_t(ω_vec)) # only exact for length(ω_vec) = even number
 
     # invertability of dft
         x_vec = [ [0.0,0.0], [3.0,0.0]]
@@ -49,7 +46,7 @@ using MultipleScattering
         timres1 = frequency_to_time(simres; t_vec = t_vec, method=:trapezoidal, impulse = GaussianImpulse(maximum(simres.ω)))
         timres2 = frequency_to_time(simres; t_vec = t_vec, method=:dft, impulse = GaussianImpulse(maximum(simres.ω)))
         @test norm(field(timres1) - field(timres2))/norm(field(timres1)) < 2e-5
-        # plot(timres1.t', [field(timres1)[1,:]-field(timres2)[1,:]])
+        # plot(timres1.t, [field(timres1)[1,:]-field(timres2)[1,:]])
 end
 
 @testset "Time shift signal" begin