Merge branch 'dev' into master

docs/make.jl

@@ -13,6 +13,7 @@ makedocs(format = :html,
             "tutorial_optim_radius.md"
         ],
         "Choosing Minimal N and P" => "minimalNP.md",
+        "Incident Field Types" => "incident_fields.md",
         "Adding New Shapes" => "new_shapes.md",
         "API" => "api.md"
     ]

docs/src/incident_fields.md

@@ -0,0 +1,54 @@
+# Incident Field Types
+
+ParticleScattering solves the multiple-scattering equation in the cylindrical
+wave domain, and thus a transformation is necessary in order to solve for arbitrary incident fields. The following field types are currently supported:
+
+## Plane Wave
+
+A TM plane wave with angle ``\theta_i`` relative to the x-axis is represented by
+```math
+E_z^{\mathrm{inc}} (\mathbf{r}) = e^{ik(\cos \theta_i, \, \sin \theta_i) \cdot \mathbf{r}},
+```
+and constructed by calling `PlaneWave(θ_i)`.
+
+## Line Current
+
+```math
+E_z^{\mathrm{inc}} (\mathbf{r}) = \frac{i}{4} H_0^{(1)} (k|\mathbf{r} - (x_0, y_0)|)
+```
+
+Following standard notation, this is a line current with total current ``I = i/k \eta``. This source cannot be placed inside a particle, and is constructed with `LineSource(x0, y0)`.
+
+## Current Source
+
+A straight line containing an arbitrary single potential density `\sigma`, producing the following incident field:
+
+```math
+E_z^{\mathrm{inc}} (\mathbf{r}) = \frac{i}{4} \int \sigma(\mathbf{r}') H_0^{(1)} (k|\mathbf{r} - \mathbf{r}'|) \, \mathrm{d}\mathbf{r}'
+```
+
+The current density of this source is ``J_z = i \sigma /k \eta ``. Current
+sources must lie completely outside all particles, and are contructed with `CurrentSource(x1, y1, x2, y2, σ)`, where `(x1,y1)` and `(x2,y2)` denote the start and end points, and `σ` is an array containing the single layer potential density at equidistant points along the source. For example,
+```julia
+using ParticleScattering, PyPlot
+λ = 1
+yc = linspace(-0.5λ, 0.5λ, 100)
+ui = CurrentSource(0, -0.5λ, 0, 0.5λ, cos.(π*yc/λ))
+#now plot
+x_points = 100; y_points = 100
+x = linspace(-3λ, 3λ, x_points + 1)
+y = linspace(-3λ, 3λ, y_points + 1)
+xgrid = repmat(x', y_points + 1, 1)
+ygrid = repmat(y, 1, x_points + 1)
+points = cat(2, vec(xgrid[1:y_points, 1:x_points]) + 3λ/x_points,
+            vec(ygrid[1:y_points, 1:x_points]) + 3λ/y_points)
+u = uinc(2π/λ, points, ui)
+pcolormesh(xgrid, ygrid, abs.(reshape(u, y_points, x_points)))
+colorbar()
+```
+yields the figure
+```@raw html
+<div style="text-align:center">
+<img alt=halfwave src="./assets/halfwave.png" style="width:80%; height:auto; margin:1%; max-width: 400px">
+</div>
+```

docs/src/tutorial1.md

@@ -24,7 +24,7 @@ relationship between these parameters and the various resulting errors, see
 k0 = 2π/λ0
 kin = 3k0
 θ_i = 0.0 #incident wave is left->right
-
+pw = PlaneWave(θ_i)
 N = 260
 P = 10
 shapes = [rounded_star(0.1λ0, 0.05λ0, 5, N)]
@@ -38,7 +38,7 @@ Now that the scattering problem is set up, we solve for the cylindrical harmonic
 coefficients and potential densities, respectively, by using
 
 ```julia
-beta,inner = solve_particle_scattering(k0, kin, P, sp, θ_i)
+beta,inner = solve_particle_scattering(k0, kin, P, sp, pw)
 ```
 
 These can be used to calculate the scattered field at any point in space using
@@ -51,7 +51,7 @@ appropriate method:
 ```julia
 #calculate field on the x-axis passing through the particle
 points = [linspace(-0.5λ0, 0.5λ0, 200)  zeros(200)]
-u = calc_near_field(k0, kin, P, sp, points, θ_i)
+u = calc_near_field(k0, kin, P, sp, points, pw)
 ```
 
 We use `PyPlot` to display the result:
@@ -64,7 +64,7 @@ plot(points[:,1]/λ0, abs.(u))
 
 Similarly, a 2D plot can be drawn of the total field around the scatterer:
 ```julia
-plot_near_field(k0, kin, P, sp, θ_i;
+plot_near_field(k0, kin, P, sp, pw;
     x_points = 201, y_points = 201, border = 0.5λ0*[-1;1;-1;1])
 ```
 
@@ -107,7 +107,7 @@ lead to unnecessary computations.
 
 Plotting the near field with the code
 ```julia
-data = plot_near_field(k0, kin, P, sp, θ_i)
+data = plot_near_field(k0, kin, P, sp, PlaneWave(θ_i))
 colorbar()
 ```
 yields the following near-field plot:

docs/src/tutorial2.md

@@ -69,7 +69,7 @@ Calculating and plotting the near or far fields with FMM is just as in the
 `FMMoptions` object:
 
 ```julia
-plot_near_field(k0, kin, P, sp, θ_i, opt = fmm_options,
+plot_near_field(k0, kin, P, sp, PlaneWave(θ_i), opt = fmm_options,
                 border = [-12;12;-10;10], x_points = 480, y_points = 400)
 colorbar()
 ```
@@ -111,9 +111,9 @@ inners = Array{Vector}(Pmax)
 inners_FMM = Array{Vector}(Pmax)
 fmmopts = ParticleScattering.FMMoptions(true, nx = 1, acc = 9)
 for P = 1:Pmax
-	betas[P], inners[P] = solve_particle_scattering(k0, kin, P, sp, 0.0;
+	betas[P], inners[P] = solve_particle_scattering(k0, kin, P, sp, PlaneWave();
                             verbose = false)
-	res, inners_FMM[P] = solve_particle_scattering_FMM(k0, kin, P, sp, 0.0,
+	res, inners_FMM[P] = solve_particle_scattering_FMM(k0, kin, P, sp, PlaneWave(),
                             fmmopts; verbose = false)
 	betas_FMM[P] = res[1]
 end

docs/src/tutorial_optim_angle.md

@@ -24,6 +24,7 @@ First we set up our scattering problem:
 k0 = 2π/λ0
 kin = 0.5k0
 θ_i = 0 #incident wave e^{i k_0 (1/sqrt{2},1/sqrt{2}) \cdot \mathbf{r}}
+pw = PlaneWave(θ_i)
 M = 20
 shapes = [rounded_star(0.35λ0, 0.1λ0, 4, 202)]
 P = 12
@@ -57,9 +58,9 @@ optim_method = Optim.BFGS(;linesearch = LineSearches.BackTracking())
 We now run both minimization and maximization:
 
 ```julia
-res_min = optimize_φ(φs0, points, P, θ_i, k0, kin, shapes, centers, ids,
+res_min = optimize_φ(φs0, points, P, pw, k0, kin, shapes, centers, ids,
             fmm_options, optim_options, optim_method; minimize = true)
-res_max = optimize_φ(φs0, points, P, θ_i, k0, kin, shapes, centers, ids,
+res_max = optimize_φ(φs0, points, P, pw, k0, kin, shapes, centers, ids,
             fmm_options, optim_options, optim_method; minimize = false)
 sp_min = ScatteringProblem(shapes, ids, centers, res_min.minimizer)
 sp_max = ScatteringProblem(shapes, ids, centers, res_max.minimizer)
@@ -71,19 +72,19 @@ following PyPlot code:
 
 ```julia
 plts = Array{Any}(3)
-plts[1] = plot_near_field(k0, kin, P, sp, θ_i, x_points = 100, y_points = 300,
+plts[1] = plot_near_field(k0, kin, P, sp, pw, x_points = 100, y_points = 300,
         opt = fmm_options, border = find_border(sp, points))
-plts[2] = plot_near_field(k0, kin, P, sp_min, θ_i, x_points = 100, y_points = 300,
+plts[2] = plot_near_field(k0, kin, P, sp_min, pw, x_points = 100, y_points = 300,
         opt = fmm_options, border = find_border(sp, points))
-plts[3] = plot_near_field(k0, kin, P, sp_max, θ_i, x_points = 100, y_points = 300,
+plts[3] = plot_near_field(k0, kin, P, sp_max, pw, x_points = 100, y_points = 300,
         opt = fmm_options, border = find_border(sp, points))
 close("all")
 
 fig, axs = subplots(ncols=3); msh = 0
 for (i, spi) in enumerate([sp;sp_min;sp_max])
     msh = axs[i][:pcolormesh](plts[i][2][1], plts[i][2][2], abs.(plts[i][2][3]),
                         vmin = 0, vmax = 3.4, cmap="viridis")
-    draw_shapes(spi.shapes, spi.centers, spi.ids, spi.φs, axs[i])
+    draw_shapes(spi, ax = axs[i])
     axs[i][:set_aspect]("equal", adjustable = "box")
     axs[i][:set_xlim]([border[1];border[2]])
     axs[i][:set_ylim]([border[3];border[4]])

docs/src/tutorial_optim_radius.md

@@ -16,6 +16,7 @@ k0 = 2π
 kin = sqrt(er)*k0
 a = 0.2*2π/k0     #wavelength/5
 θ_i = 0.0
+pw = PlaneWave(θ_i)
 P = 5
 centers = square_grid(5, a)
 φs = zeros(size(centers,1))
@@ -60,7 +61,7 @@ assert(verify_min_distance([CircleParams(r_max[i]) for i = 1:J],
 The optimization process is initiated by running:
 
 ```julia
-res = optimize_radius(rs0, r_min, r_max, points, ids, P, θ_i, k0, kin,
+res = optimize_radius(rs0, r_min, r_max, points, ids, P, pw, k0, kin,
                 centers, fmm_options, optim_options, minimize = true)
 rs = res.minimizer
 ```
@@ -70,14 +71,14 @@ initial and optimized radii:
 
 ```julia
 sp1 = ScatteringProblem([CircleParams(rs0[i]) for i = 1:J], ids, centers, φs)
-plot_near_field(k0, kin, P, sp1, θ_i, x_points = 150, y_points = 150,
+plot_near_field(k0, kin, P, sp1, pw, x_points = 150, y_points = 150,
         opt = fmm_options, border = 0.9*[-1;1;-1;1], normalize = a)
 colorbar()
 clim([0;2.5])
 xlabel("x/a")
 ylabel("y/a")
 sp2 = ScatteringProblem([CircleParams(rs[i]) for i = 1:J], ids, centers, φs)
-plot_near_field(k0, kin, P, sp2, θ_i, x_points = 150, y_points = 150,
+plot_near_field(k0, kin, P, sp2, pw, x_points = 150, y_points = 150,
         opt = fmm_options, border = 0.9*[-1;1;-1;1], normalize = a)
 colorbar()
 clim([0;2.5])

examples/sim2_LuneburgOptim.jl

@@ -17,7 +17,8 @@ P = 5
 
 fmm_options = FMMoptions(true, acc = 6, dx = 2*a_lens, method = "pre")
 
-optim_options =  Optim.Options(x_tol = 1e-6, iterations = 100,
+optim_options =  Optim.Options(x_tol = 1e-6, outer_x_tol = 1e-6,
+                        iterations = 100, outer_iterations = 100,
                         store_trace = true, extended_trace = true,
                         show_trace = false, allow_f_increases = true)
 
@@ -27,10 +28,10 @@ r_min = (a_lens*1e-3)*ones(size(centers,1))
 rs0 = (0.25*a_lens)*ones(size(centers,1))
 
 ids_max = collect(1:length(rs0))
-test_max = optimize_radius(rs0, r_min, r_max, points, ids_max, P, θ_i, k0, kin, #precompile
+test_max = optimize_radius(rs0, r_min, r_max, points, ids_max, P, PlaneWave(θ_i), k0, kin, #precompile
                 centers, fmm_options, optim_options, minimize = false)
 tic()
-test_max = optimize_radius(rs0, r_min, r_max, points, ids_max, P, θ_i, k0, kin,
+test_max = optimize_radius(rs0, r_min, r_max, points, ids_max, P, PlaneWave(θ_i), k0, kin,
                 centers, fmm_options, optim_options, minimize = false)
 optim_time = toq()
 rs_max = test_max.minimizer
@@ -43,23 +44,23 @@ border = (R_lens + a_lens)*[-1;1;-1;1]
 
 sp1 = ScatteringProblem([CircleParams(rs_lnbrg[i]) for i in eachindex(rs_lnbrg)],
         ids_lnbrg, centers, φs)
-Ez1 = plot_near_field(k0, kin, P, sp1, θ_i, x_points = 150, y_points = 150,
+Ez1 = plot_near_field(k0, kin, P, sp1, PlaneWave(θ_i), x_points = 150, y_points = 150,
         opt = fmm_options, border = border)
-plot_near_field_pgf(filename1, k0, kin, P, sp1, θ_i; opt = fmm_options,
+plot_near_field_pgf(filename1, k0, kin, P, sp1, PlaneWave(θ_i); opt = fmm_options,
     x_points = 201, y_points = 201, border = border)
 
 sp2 = ScatteringProblem([CircleParams(rs_max[i]) for i in eachindex(rs_max)],
         ids_max, centers, φs)
-Ez2 = plot_near_field(k0, kin, P, sp2, θ_i, x_points = 150, y_points = 150,
+Ez2 = plot_near_field(k0, kin, P, sp2, PlaneWave(θ_i), x_points = 150, y_points = 150,
             opt = fmm_options, border = border)
-plot_near_field_pgf(filename2, k0, kin, P, sp2, θ_i; opt = fmm_options,
+plot_near_field_pgf(filename2, k0, kin, P, sp2, PlaneWave(θ_i); opt = fmm_options,
     x_points = 201, y_points = 201, border = border)
 
 sp3 = ScatteringProblem([CircleParams(rs0[i]) for i in eachindex(rs0)],
         collect(1:length(rs0)), centers, φs)
-Ez3 = plot_near_field(k0, kin, P, sp3, θ_i, x_points = 150, y_points = 150,
+Ez3 = plot_near_field(k0, kin, P, sp3, PlaneWave(θ_i), x_points = 150, y_points = 150,
         opt = fmm_options, border = border)
-plot_near_field_pgf(filename3, k0, kin, P, sp3, θ_i; opt = fmm_options,
+plot_near_field_pgf(filename3, k0, kin, P, sp3, PlaneWave(θ_i); opt = fmm_options,
     x_points = 201, y_points = 201, border = border)
 
 #plot convergence
@@ -114,26 +115,26 @@ r_min = (a_lens*1e-3)*ones(J)
 rs0 = (0.25*a_lens)*ones(J)
 
 tic()
-test_max_sym = optimize_radius(rs0, r_min, r_max, points, ids_sym, P, θ_i, k0, kin,
+test_max_sym = optimize_radius(rs0, r_min, r_max, points, ids_sym, P, PlaneWave(θ_i), k0, kin,
                 centers, fmm_options, optim_options, minimize = false)
 sym_time = toq()
 rs_sym = test_max_sym.minimizer
 JLD.@save output_dir * "/luneburg_optim_sym.jld" test_max_sym sym_time
 
 sp4 = ScatteringProblem([CircleParams(rs_sym[i]) for i in eachindex(rs_sym)],
         ids_sym, centers, φs)
-Ez4 = plot_near_field(k0, kin, P, sp4, θ_i, x_points = 150, y_points = 150,
+Ez4 = plot_near_field(k0, kin, P, sp4, PlaneWave(θ_i), x_points = 150, y_points = 150,
         opt = fmm_options, border = border)
-
-u1 = calc_near_field(k0, kin, 7, sp1, points, θ_i; opt = fmm_options)
-u2 = calc_near_field(k0, kin, 7, sp2, points, θ_i; opt = fmm_options)
-u3 = calc_near_field(k0, kin, 7, sp3, points, θ_i; opt = fmm_options)
-u4 = calc_near_field(k0, kin, 7, sp4, points, θ_i; opt = fmm_options)
+pw = PlaneWave(θ_i)
+u1 = calc_near_field(k0, kin, 7, sp1, points, pw; opt = fmm_options)
+u2 = calc_near_field(k0, kin, 7, sp2, points, pw; opt = fmm_options)
+u3 = calc_near_field(k0, kin, 7, sp3, points, pw; opt = fmm_options)
+u4 = calc_near_field(k0, kin, 7, sp4, points, pw; opt = fmm_options)
 abs.([u1[1];u2[1];u3[1];u4[1]])
 
 #####################################
 # selfconsistent err P calculation
-Ez_4 = calc_near_field(k0, kin, 4, sp1, points, θ_i; opt = fmm_options)
-Ez_5 = calc_near_field(k0, kin, 5, sp1, points, θ_i; opt = fmm_options)
-Ez_6 = calc_near_field(k0, kin, 6, sp1, points, θ_i; opt = fmm_options)
-Ez_7 = calc_near_field(k0, kin, 7, sp1, points, θ_i; opt = fmm_options)
+Ez_4 = calc_near_field(k0, kin, 4, sp1, points, pw; opt = fmm_options)
+Ez_5 = calc_near_field(k0, kin, 5, sp1, points, pw; opt = fmm_options)
+Ez_6 = calc_near_field(k0, kin, 6, sp1, points, pw; opt = fmm_options)
+Ez_7 = calc_near_field(k0, kin, 7, sp1, points, pw; opt = fmm_options)

examples/sim3_AngleOpt.jl

@@ -43,7 +43,7 @@ begin
     if draw_fig
         figure()
         #draw shapes and points
-        draw_shapes(shapes, centers, ids, φs)
+        draw_shapes(shapes, ids, centers, φs)
         plot(points[:,1], points[:,2], "r*")
         tight_layout()
         ax = gca()
@@ -63,28 +63,28 @@ optim_options =  Optim.Options(f_tol = 1e-6,
 optim_method = Optim.BFGS(;linesearch = LineSearches.BackTracking())
 
 tic()
-test_max = optimize_φ(φs, points, P, θ_i, k0, kin, shapes,
+test_max = optimize_φ(φs, points, P, PlaneWave(θ_i), k0, kin, shapes,
             centers, ids, fmm_options, optim_options, optim_method; minimize = false)
 optim_time = toq()
 
 # %%
 
 sp_before = ScatteringProblem(shapes, ids, centers, φs)
-plot_near_field(k0, kin, P, sp_before, θ_i,
+plot_near_field(k0, kin, P, sp_before, PlaneWave(θ_i),
                 x_points = 600, y_points = 200, border = plot_border);
 colorbar()
 clim([0;5])
 plot_near_field_pgf(output_dir * "/opt_phi_before.tex", k0, kin, P,
-    sp_before, θ_i; opt = fmm_options, x_points = 150, y_points = 50,
+    sp_before, PlaneWave(θ_i); opt = fmm_options, x_points = 150, y_points = 50,
     border = plot_border, downsample = 10, include_preamble = true)
 
 sp_after = ScatteringProblem(shapes, ids, centers, test_max.minimizer)
-plot_near_field(k0, kin, P, sp_after, θ_i,
+plot_near_field(k0, kin, P, sp_after, PlaneWave(θ_i),
                 x_points = 600, y_points = 200, border = plot_border)
 colorbar()
 clim([0;5])
 plot_near_field_pgf(output_dir * "/opt_phi_after.tex", k0, kin, P,
-    sp_after, θ_i; opt = fmm_options, x_points = 150, y_points = 50,
+    sp_after, PlaneWave(θ_i); opt = fmm_options, x_points = 150, y_points = 50,
     border = plot_border, downsample = 10, include_preamble = true)
 
 inner_iters = length(test_max.trace)