Merge pull request #15 from bblankrot/dev

Weekly PR #4

docs/src/minimalNP.md

@@ -25,7 +25,7 @@ Specifically, this error is measured on a circle of radius
 R_multipole*maximum(hypot.(s.ft[:,1],s.ft[:,2]))
 ```
 
-as this is the scattering disk on which the translation to cylindrical harmonics
+as this is the scattering disc on which the translation to cylindrical harmonics
 (Bessel & Hankel functions) are performed (and beyond which any gain or loss of
 accuracy due to `N` is mostly irrelevant).
 

src/ParticleScattering.jl

@@ -11,7 +11,7 @@ using IterativeSolvers, LinearMaps, Optim
 import LineSearches
 #For plotting with PyPlot
 using PyPlot, PyCall
-@pyimport matplotlib.patches as patch #circles,polygons
+@pyimport matplotlib.patches as patch #circles, polygons
 
 include("PS_types.jl")
 include("shapes.jl")
@@ -44,9 +44,6 @@ export ScatteringProblem, OptimBuffer, FMMoptions, R_multipole,
     ShapeParams, CircleParams, AbstractShapeParams
 #methods, utilities.jl
 export find_border, uniqueind
-
-### documented till here####
-
 #methods, optimize_phis.jl
 export optimize_φ
 #methods, optimize_rs.jl
@@ -60,10 +57,7 @@ include("visualization_pgf.jl")
 export plotNearField_pgf, drawShapes_pgf
 
 #temp
-include("optimize_rs_old.jl")
-
 export divideSpace, FMMtruncation, particleExpansion, FMMbuildMatrices
-
 export FMMbuffer, FMMmatrices
 
 end

src/minimum_N_P.jl

@@ -10,7 +10,7 @@ Error is calculated on `N_points` points on the scattering disk (`s.R`), by
 assuming a fictitious line source and comparing its field to that produced
 by the resulting potential densities.
 
-Since for moderate wavelengths and errors, \$\varepsilon \\propto N^3\$, we
+Since the error scales with \$N^{-3}\$ for moderate wavelengths and errors, we
 estimate `N` using the error of `N_start`, then binary search based on that
 guess between `N_min` and `N_max`.
 """

src/multipole.jl

@@ -185,7 +185,8 @@ function innerFieldCircle(kin, gamma, center::Array{Float64,1}, points::Array{Fl
 	bess = [besselj(p,kin*rs_moved[ii]) for ii=1:len, p=0:P]
 	E = gamma[P + 1]*bess[:,1]
 	for p = 1:P
-		E += bess[:,p+1].*(gamma[p + P + 1]*exp.(1.0im*p*ts_moved) + (-1)^p*gamma[-p + P + 1]*exp.(-1.0im*p*ts_moved))
+		E += bess[:,p+1].*(gamma[p + P + 1]*exp.(1.0im*p*ts_moved) + 
+							(-1)^p*gamma[-p + P + 1]*exp.(-1.0im*p*ts_moved))
 	end
 	return E
 end

src/optimize_phis.jl

@@ -3,8 +3,8 @@
         optimopts::Optim.Options, minimize = true)
 
 Optimize the rotation angles of a particle collection for minimization or
-maximization of the field intensity at `points`, depending on `minimize`.
-`optimopts` and `method` define the optimization emthod, convergence criteria,
+maximization (depending on `minimize`) of the field intensity at `points`.
+`optimopts` and `method` define the optimization method, convergence criteria,
 and other optimization parameters.
 Returns an object of type `Optim.MultivariateOptimizationResults`.
 """

src/shapes.jl

@@ -126,7 +126,7 @@ function randpoints(M, dmin, width, height, points; failures = 100)
     dmin2 = dmin^2
     for i = 2:size(points,1), j = 1:i-1
         dist2 = sum(x -> x^2, points[i,:] - points[j,:])
-        dist2 <= dmin2 && error("randpoints: given points have distance > dmin")
+        dist2 <= dmin2 && error("randpoints: given points have distance <= dmin")
     end
 
     x_res = [rand(Float64)*width]

src/visualization.jl

@@ -9,7 +9,7 @@ angle `θ_i` scattering from the ScatteringProblem `sp`, using matplotlib's
 `pcolormesh`. Can accept number of sampling points in each direction plus
 bounding box or calculate automatically.
 
-Uses the FMM options given by `opt` (default behavious is disabled FMM);
+Uses the FMM options given by `opt` (FMM is disabled by default);
 `use_multipole` dictates whether electric field is calculated using the
 multipole/cylindrical harmonics (true) or falls back on potential densities
 (false). Either way, the multiple-scattering system is solved in the cylindrical

test/fmm_test.jl

@@ -14,10 +14,10 @@
     φs = 2π*rand(M^2) #random rotation angles
     dist = 2*maximum(shapes[i].R for i=1:2)
     try
-    centers = randpoints(M^2, dist, 5λ0, 5λ0, [0.0 0.0],
+    centers = randpoints(M^2, dist, 5λ0, 5λ0, [0.0 0.0; 0.01+dist 0.01],
                 failures = 1_000)
     sp = ScatteringProblem(shapes, ids, centers, φs)
-    @test verify_min_distance(sp, [0.0 0.0])
+    @test verify_min_distance(sp, [0.0 0.0; 0.01+dist 0.01])
     catch
     warn("Could not find random points (1)")
     end

test/multipole_test.jl

@@ -0,0 +1,56 @@
+@testset "multipole" begin
+    #let's check that rotation works properly
+    k0 = 0.1
+    kin = 100k0
+    λ0 = 2π/k0
+    θ_i = 0.0
+
+    N = 1000
+    P = 5
+    shapes = [squircle(2π/λ0, N)]
+    centers = [-1.5λ0 0; 1.5λ0 0]
+    Ns = size(centers,1)
+    φs = [0.0; 0.23]
+    ids = ones(Int, Ns)
+    sp = ScatteringProblem(shapes, ids, centers, φs)
+    points = λ0*[5 0; -2 3]
+
+    verify_min_distance(shapes, centers, ids, points)
+    β1, σ1 = solve_particle_scattering(k0, kin, P, sp, θ_i;
+        get_inner = true, verbose = false)
+    Ez1 = scattered_field_multipole(k0, β1, centers, points)
+
+    #now rotate everything by some angle and compare
+    θ_r = 1.23#2π*rand()
+    #notation is transposed due to structure of centers
+    shapes2 = [squircle(2π/λ0, N); CircleParams(15)] #for extra code coverage
+    centers2 = centers*([cos(θ_r) -sin(θ_r); sin(θ_r) cos(θ_r)].')
+    φs2 = θ_r + φs
+    sp2 = ScatteringProblem(shapes2, ids, centers2, φs2)
+    points2 = points*([cos(θ_r) -sin(θ_r); sin(θ_r) cos(θ_r)].')
+
+    β2, σ2 = solve_particle_scattering(k0, kin, P, sp2, θ_i + θ_r;
+        get_inner = true, verbose = true)
+    Ez2 = scattered_field_multipole(k0, β2, centers2, points2)
+
+    @test Ez1 ≈ Ez2
+    β2_r = copy(β2)
+    for ic = 1:Ns, p = -P:P
+        β2_r[(ic-1)*(2P+1) + p + P + 1] *= exp(1.0im*θ_r*p)
+    end
+    @test β1 ≈ β2_r
+    @test σ1 ≈ σ2
+end
+
+@testset "circle" begin
+    k0 = 1
+    kin = 0.1
+    R = 0.2
+    P = 10
+    S, gamma = ParticleScattering.circleScatteringMatrix(k0, kin, R, P; gamma = true)
+    center = [0 2R]
+    points = center .+ [0 0.5R; 0.1R -0.2R]
+    E1 = ParticleScattering.innerFieldCircle(kin, gamma, center[1,:], points)
+    E2 = ParticleScattering.innerFieldCircle(kin, gamma, center[1,:], points[1,:])
+    @test E1[1] == E2
+end

test/runtests.jl

@@ -2,7 +2,9 @@ using ParticleScattering
 using Base.Test
 
 include("scatteredfield_test.jl")
+include("multipole_test.jl")
+include("fmm_test.jl")
 include("optimize_radius_test.jl")
 include("optimize_angle_test.jl")
 include("minimum_N_P_test.jl")
-include("fmm_test.jl")
+include("utilities_test.jl")

test/utilities_test.jl

@@ -0,0 +1,13 @@
+@testset "utilities" begin
+    s = rounded_star(1, 0.1, 7, 100)
+    @test ParticleScattering.pInPolygon([0.0 0.0], s.ft)
+    @test !ParticleScattering.pInPolygon([0.84 0.35], s.ft)
+    @test ParticleScattering.pInPolygon([0.84 0.35], s.ft+0.1)
+
+    sp = ScatteringProblem([s],[1],[0.1 0.1], [15π])
+    border = find_border(sp)
+    border2 = find_border(sp,[0.1 0.1])
+    @test border == border2
+    border3 = find_border(sp,[0.1 0.9999])
+    @test border != border3
+end