diff --git a/README.md b/README.md
index 5a6f7d4..4f9c831 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-# ParticleScattering.jl
+# ParticleScattering
 
 [![Travis](https://travis-ci.org/bblankrot/ParticleScattering.jl.svg?branch=master)](https://travis-ci.org/bblankrot/ParticleScattering.jl)
 [![AppVeyor](https://ci.appveyor.com/api/projects/status/p0p636vtrx95ch8m/branch/master?svg=true)](https://ci.appveyor.com/project/bblankrot/particlescattering-jl/branch/master)
@@ -11,11 +11,16 @@ those containing a large number of penetrable smooth
 particles. Provides the ability to optimize over the
 particle parameters for various design problems.
 
-# Installation
+### Installation
 
-Currently, only Julia 0.6 is supported, and the package must be installed
-manually. See directions at [the official Julia documentation](https://docs.julialang.org/en/stable/manual/packages/#Installing-Unregistered-Packages-1).
+ParticleScattering can be installed using `Pkg.add`. Currently, only Julia 0.6 is supported.
+
+```julia
+Pkg.add("ParticleScattering")
+using ParticleScattering
+```
 
 ### Community
 
-The easiest way to contribute is by opening issues! Of course, we'd be more than happy if you implement any fixes and send a PR. If you have any relevant scattering problems that would make good examples for the docs, feel free to open an issue for that as well.
+The easiest way to contribute is by opening issues! Of course, we'd be more than happy if you implement any fixes and send a PR.
+If you have any relevant scattering problems that would make good examples for the docs, feel free to open an issue for that as well.
diff --git a/REQUIRE b/REQUIRE
index b85a102..3ea19cb 100644
--- a/REQUIRE
+++ b/REQUIRE
@@ -6,4 +6,4 @@ Optim 0.14.0
 PyPlot
 DataFrames
 CSV
-PGFPlotsX
+PGFPlotsX 0.2.1
diff --git a/docs/make.jl b/docs/make.jl
index d79ebce..22f41f6 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -13,6 +13,7 @@ makedocs(format = :html,
             "tutorial_optim_radius.md"
         ],
         "Choosing Minimal N and P" => "minimalNP.md",
+        "Adding New Shapes" => "new_shapes.md",
         "API" => "api.md"
     ]
 )
diff --git a/docs/src/index.md b/docs/src/index.md
index a335e15..6bd3ab7 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -7,9 +7,10 @@ particles. Provides the ability to optimize over particle parameters for various
 
 ## Installation
 
-Currently, only Julia 0.6 is supported. Once Julia is set up, ParticleScattering can be installed by running
+ParticleScattering can be installed using `Pkg.add`. Currently, only Julia 0.6 is supported.
+
 ```julia
-Pkg.clone(https://github.com/bblankrot/ParticleScattering.jl.git)
+Pkg.add("ParticleScattering")
 using ParticleScattering
 ```
 which also installs the following dependencies:
diff --git a/docs/src/new_shapes.md b/docs/src/new_shapes.md
new file mode 100644
index 0000000..89aee42
--- /dev/null
+++ b/docs/src/new_shapes.md
@@ -0,0 +1,17 @@
+# Adding New Shapes
+
+ParticleScattering includes functions for drawing squircles, rounded stars, and ellipses.
+New shape functions can  be added, provided they have the following structure:
+```julia
+function my_shape(args, N)
+    t = Float64[π*j/N for j = 0:(2*N-1)] # or t = 0:π/N:π*(2-1/N)
+    ft =  [x    y]
+    dft = [dx/dt    dy/dt]
+    ShapeParams(t, ft, dft)
+end
+```
+
+Where `t` is the parametrization variable, `ft[i,:] = [x(t[i]) y(t[i])]` contains the coordinates, and dft contains the derivative of `ft` with respect to `t`.
+In particular, the quadrature used by ParticleScattering assumes `t`
+are equidistantly distributed in ``[0, 2\pi)``, and that none of the points `ft`
+lie on the origin.
diff --git a/examples/sim0_min_NP.jl b/examples/sim0_min_NP.jl
index 16f964d..c6a8ad5 100644
--- a/examples/sim0_min_NP.jl
+++ b/examples/sim0_min_NP.jl
@@ -121,20 +121,16 @@ JLD.@load dirname(@__FILE__) * "/mindata_5_new.jld"
 JLD.@load dirname(@__FILE__) * "/mindata_s_new.jld"
 
 pgf.@pgf begin
-    N_plot5 = pgf.Plot(pgf.Coordinates(errN5,N5),
-        {blue, dashdotdotted, no_markers, thick},
-		label = "\$N_{\\mathrm{min}} \\, \\mathrm{(star)}\$")
-    P_plot5 = pgf.Plot(pgf.Coordinates(errN5,P5),
-        {"green!50!black", no_markers, thick},
-		label = "\$P_{\\mathrm{min}} \\, \\mathrm{(star)}\$")
-    N_plot = pgf.Plot(pgf.Coordinates(errNs,Ns),
-        {red, dotted, no_markers, thick},
-		label = "\$N_{\\mathrm{min}} \\, \\mathrm{(squircle)}\$")
-    P_plot = pgf.Plot(pgf.Coordinates(errNs,Ps),
-        {black, dashed, no_markers, thick},
-		label = "\$P_{\\mathrm{min}} \\, \\mathrm{(squircle)}\$")
-    ax = pgf.Axis([N_plot5, N_plot, P_plot5, P_plot],
-        {
+    N_plot5 = pgf.Plot({blue, dashdotdotted, no_markers, thick},
+						pgf.Coordinates(errN5,N5))
+	N_plot = pgf.Plot({red, dotted, no_markers, thick},
+						pgf.Coordinates(errNs,Ns))
+	P_plot5 = pgf.Plot({"green!50!black", no_markers, thick},
+						pgf.Coordinates(errN5,P5))
+    P_plot = pgf.Plot({black, dashed, no_markers, thick},
+						pgf.Coordinates(errNs,Ps))
+    ax = pgf.Axis(
+    	{
             xmin = 5e-10,
             xmax = 1e-1,
 			xlabel = "\$\\Delta u\$",
@@ -144,7 +140,13 @@ pgf.@pgf begin
             legend_pos = "north east",
             legend_style = "font = \\footnotesize",
 			legend_cell_align = "left"
-        })
+        },
+		[N_plot5, N_plot, P_plot5, P_plot],
+		pgf.Legend(["\$N_{\\mathrm{min}} \\, \\mathrm{(star)}\$",
+					"\$N_{\\mathrm{min}} \\, \\mathrm{(squircle)}\$",
+					"\$P_{\\mathrm{min}} \\, \\mathrm{(star)}\$",
+					"\$P_{\\mathrm{min}} \\, \\mathrm{(squircle)}\$"]))
+
 end
 pgf.save(dirname(@__FILE__) * "/tikz/minNP.tex", ax ,include_preamble = false)
 
@@ -152,9 +154,8 @@ shape = myshapefun5(200)
 pgf_shape = pgf.Coordinates([shape.ft[:,1];shape.ft[1,1]],
                             [shape.ft[:,2];shape.ft[1,2]])
 pgf.@pgf begin
-    ax2 = pgf.Axis(pgf.Plot(pgf_shape, {thick, black, no_markers, fill = "black!20"}),
-        {   axis_equal,
-            ticks = "none"})
+    ax2 = pgf.Axis({axis_equal, ticks = "none"}, pgf.Plot({thick, black,
+			no_markers, fill = "black!20"}, pgf_shape))
     Rd = shape.R
     push!(ax2, "\\addplot [black, dashed, thick, domain=0:2*pi,samples=100]({$(shape.R)*cos(deg(x))},{$(shape.R)*sin(deg(x))});")
     push!(ax2, "\\node at (axis cs: -0.28,-0.15) {\$D\$};")
@@ -169,10 +170,9 @@ shape = myshapefun_squircle(200)
 pgf_shape = pgf.Coordinates([shape.ft[:,1];shape.ft[1,1]],
                             [shape.ft[:,2];shape.ft[1,2]])
 pgf.@pgf begin
-    ax3 = pgf.Axis(pgf.Plot(pgf_shape, {thick, black, no_markers, fill = "black!20"}),
-        {   axis_equal,
-            ticks = "none"})
-    Rd = shape.R
+	ax3 = pgf.Axis({axis_equal, ticks = "none"}, pgf.Plot({thick, black,
+					no_markers, fill = "black!20", label="bla"}, pgf_shape))
+	Rd = shape.R
     push!(ax3, "\\addplot [black, dashed, thick, domain=0:2*pi,samples=100]({$(shape.R)*cos(deg(x))},{$(shape.R)*sin(deg(x))});")
     push!(ax3, "\\node at (axis cs: -0.28,-0.15) {\$D\$};")
     push!(ax3, "\\node at (axis cs: 0.25,-0.25) {\$k_0\$};")
diff --git a/examples/sim1_complexity.jl b/examples/sim1_complexity.jl
index e7758a5..9c17a83 100644
--- a/examples/sim1_complexity.jl
+++ b/examples/sim1_complexity.jl
@@ -39,7 +39,6 @@ for i=1:simlen-1
 end
 dt0 /= simlen
 
-
 for is = 1:simlen, it = 1:trials
     #compute shape variables
     begin #setup
@@ -123,27 +122,27 @@ pgf.@pgf begin
             ylabel = "\$ \\mathrm{Run time} \\ [\\mathrm{s}]\$",
             xmode = "linear",
             ymode = "log",
-            width = "\\figurewidth",
+            width = "10cm",
             legend_pos = "north west",
-            legend_style = "font = \\footnotesize"})
-    push!(ax, pgf.Plot(pgf.Coordinates(M_vec, res_vec),
-            {blue, "only marks", mark = "*"};
-            label = "Elapsed time (solution)"))
-            temp2 = floor(log10(10^a_total))
-            temp1 = 10^a_total/10^temp2
-    push!(ax, pgf.Plot(pgf.Coordinates(M_vec, res_ana),
-            {blue, thick, no_markers};
-            label = @sprintf("\$%.1f \\cdot 10^{%d} \\cdot M^{%.2f}\$",
-                temp1, temp2, b_total)))
-    push!(ax, pgf.Plot(pgf.Coordinates(M_vec, mvp_vec),
-            {red, only_marks, mark = "triangle*", mark_options = {fill = "red"}};
-            label= "Elapsed time (MVP)"))
-            temp2 = floor(log10(10^a_mvp))
-            temp1 = 10^a_mvp/10^temp2
-    push!(ax, pgf.Plot(pgf.Coordinates(M_vec, mvp_ana),
-            {red, thick, dashed, no_markers};
-            label = @sprintf("\$%.1f \\cdot 10^{%d} \\cdot M^{%.2f}\$",
-                temp1, temp2, b_mvp)))
+            legend_style = "font = \\footnotesize",
+            legend_cell_align = "left"})
+    push!(ax, pgf.Plot({blue, "only marks", mark = "*"},
+                        pgf.Coordinates(M_vec, res_vec)))
+    tmp_total = floor(log10(10^a_total))
+    push!(ax, pgf.Plot({blue, thick, no_markers},
+                        pgf.Coordinates(M_vec, res_ana)))
+    push!(ax, pgf.Plot({red, only_marks, mark = "triangle*",
+                        mark_options = {fill = "red"}},
+                        pgf.Coordinates(M_vec, mvp_vec)))
+    tmp_mvp = floor(log10(10^a_mvp))
+    push!(ax, pgf.Plot({red, thick, dashed, no_markers},
+                        pgf.Coordinates(M_vec, mvp_ana)))
+    push!(ax, pgf.Legend([
+                    "Elapsed time (solution)",
+                    @sprintf("\$%.1f \\cdot 10^{%d} \\cdot M^{%.2f}\$",
+                            10^a_total/10^tmp_total, tmp_total, b_total),
+                    "Elapsed time (MVP)",
+                    @sprintf("\$%.1f \\cdot 10^{%d} \\cdot M^{%.2f}\$",
+                            10^a_mvp/10^tmp_mvp, tmp_mvp, b_mvp)]))
 end
-
 pgf.save(dirname(@__FILE__) * "/sim1.tex", ax ,include_preamble = false)
diff --git a/examples/sim2_LuneburgOptim.jl b/examples/sim2_LuneburgOptim.jl
index e5c33ee..6ca26c4 100644
--- a/examples/sim2_LuneburgOptim.jl
+++ b/examples/sim2_LuneburgOptim.jl
@@ -45,21 +45,21 @@ 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,
         opt = fmm_options, border = border)
-plotNearField_pgf(filename1, k0, kin, P, sp1, θ_i; opt = fmm_options,
+plot_near_field_pgf(filename1, k0, kin, P, sp1, θ_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,
             opt = fmm_options, border = border)
-plotNearField_pgf(filename2, k0, kin, P, sp2, θ_i; opt = fmm_options,
+plot_near_field_pgf(filename2, k0, kin, P, sp2, θ_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,
         opt = fmm_options, border = border)
-plotNearField_pgf(filename3, k0, kin, P, sp3, θ_i; opt = fmm_options,
+plot_near_field_pgf(filename3, k0, kin, P, sp3, θ_i; opt = fmm_options,
     x_points = 201, y_points = 201, border = border)
 
 #plot convergence
@@ -81,23 +81,26 @@ JLD.@save output_dir * "/luneburg_optim.jld"
 
 import PGFPlotsX; const pgf = PGFPlotsX
 pgf.@pgf begin
-    fobj_plot = pgf.Plot(pgf.Coordinates(0:inner_iters-1, fobj),
-                {blue, thick, no_markers},
-                label = "\$f_{\\mathrm{obj}}\$")
-    gobj_plot = pgf.Plot(pgf.Coordinates(0:inner_iters-1, gobj),
-                {red, thick, dashed, no_markers},
-                label = "\$\\|\\mathbf{g}_{\\mathrm{obj}}\\|_{\\infty}\$")
-    fobj_outer = pgf.Plot(pgf.Coordinates((0:inner_iters-1)[rng], fobj[rng]),
-                {blue, only_marks, mark = "*", mark_options = {fill = "blue"}})
-    gobj_outer = pgf.Plot(pgf.Coordinates((0:inner_iters-1)[rng], gobj[rng]),
-                {red, only_marks, mark = "triangle*", mark_options = {fill = "red"}})
-    ax = pgf.Axis([fobj_plot;gobj_plot;fobj_outer;gobj_outer],
-        {
-            width = "\\figurewidth",
-            xlabel = "Iterations",
-            legend_pos = "north east",
-            legend_style = "font = \\footnotesize"
-        })
+    fobj_plot = pgf.Plot({blue, thick, no_markers},
+                        pgf.Coordinates(0:inner_iters-1, fobj))
+    gobj_plot = pgf.Plot({red, thick, dashed, no_markers},
+                        pgf.Coordinates(0:inner_iters-1, gobj))
+    fobj_outer = pgf.Plot({blue, only_marks, mark = "*",
+                            mark_options = {fill = "blue"}},
+                        pgf.Coordinates((0:inner_iters-1)[rng], fobj[rng]))
+    gobj_outer = pgf.Plot({red, only_marks, mark = "triangle*",
+                            mark_options = {fill = "red"}},
+                        pgf.Coordinates((0:inner_iters-1)[rng], gobj[rng]))
+    ax = pgf.Axis({ xmin = 0,
+                    width = "6cm",
+                    xlabel = "Iterations",
+                    legend_pos = "north east",
+                    legend_style = "font = \\footnotesize",
+                    legend_cell_align = "left"},
+                fobj_plot, gobj_plot, fobj_outer, gobj_outer)
+    push!(ax, pgf.Legend(["\$f_{\\mathrm{obj}}\$";
+                        "\$\\|\\mathbf{g}_{\\mathrm{obj}}\\|_{\\infty}\$"]))
+    ax
 end
 pgf.save(output_dir * "/opt_r_conv.tex", ax ,include_preamble = false)
 
diff --git a/examples/sim3_AngleOpt.jl b/examples/sim3_AngleOpt.jl
index ba24e07..0422b72 100644
--- a/examples/sim3_AngleOpt.jl
+++ b/examples/sim3_AngleOpt.jl
@@ -74,18 +74,18 @@ plot_near_field(k0, kin, P, sp_before, θ_i,
                 x_points = 600, y_points = 200, border = plot_border);
 colorbar()
 clim([0;5])
-plotNearField_pgf(output_dir * "/opt_phi_before.tex", k0, kin, P,
-    sp_before, θ_i; opt = fmm_options, x_points = 600, y_points = 200,
-    border = plot_border, downsample = 4)
+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,
+    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,
                 x_points = 600, y_points = 200, border = plot_border)
 colorbar()
 clim([0;5])
-plotNearField_pgf(output_dir * "/opt_phi_after.tex", k0, kin, P,
-    sp_after, θ_i; opt = fmm_options, x_points = 600, y_points = 200,
-    border = plot_border, downsample = 4)
+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,
+    border = plot_border, downsample = 10, include_preamble = true)
 
 inner_iters = length(test_max.trace)
 fobj = -[test_max.trace[i].value for i=1:inner_iters]
@@ -95,20 +95,17 @@ figure()
 plot(0:inner_iters-1, fobj)
 plot(0:inner_iters-1, gobj)
 
-import PGFPlotsX; const pgf = PGFPlotsX
 pgf.@pgf begin
-    fobj_plot = pgf.Plot(pgf.Coordinates(0:inner_iters-1, fobj),
-                {blue, thick, no_markers},
-                label = "\$f_{\\mathrm{obj}}\$")
-    gobj_plot = pgf.Plot(pgf.Coordinates(0:inner_iters-1, gobj),
-                {red, thick, dashed, no_markers},
-                label = "\$\\|\\mathbf{g}_{\\mathrm{obj}}\\|\$")
-    ax = pgf.Axis([fobj_plot;gobj_plot],
-        {
-            width = "\\figurewidth",
-            xlabel = "Iterations",
-            legend_pos = "north east",
-            legend_style = "font = \\footnotesize"
-        })
+    fobj_plot = pgf.Plot({blue, thick, no_markers},
+                        pgf.Coordinates(0:inner_iters-1, fobj))
+    gobj_plot = pgf.Plot({red, thick, dashed, no_markers},
+                        pgf.Coordinates(0:inner_iters-1, gobj))
+    ax = pgf.Axis({ width = "6cm",
+                    xlabel = "Iterations",
+                    legend_pos = "north east",
+                    legend_style = "font = \\footnotesize"},
+                fobj_plot, gobj_plot)
+    push!(ax, pgf.Legend(["\$f_{\\mathrm{obj}}\$";
+                        "\$\\|\\mathbf{g}_{\\mathrm{obj}}\\|\$"]))
 end
 pgf.save(output_dir * "/opt_phi_conv.tex", ax ,include_preamble = false)
diff --git a/paper/fig.png b/paper/fig.png
new file mode 100644
index 0000000..2e4a0ab
Binary files /dev/null and b/paper/fig.png differ
diff --git a/paper/paper.bib b/paper/paper.bib
new file mode 100644
index 0000000..a3ae6be
--- /dev/null
+++ b/paper/paper.bib
@@ -0,0 +1,32 @@
+@article{ar:blankrot2018,
+	author = {Boaz Blankrot and Clemens Heitzinger},
+	title = {Efficient computational design and optimization of dielectric metamaterial devices},
+	journal = {submitted for publication},
+	year = {2018}
+}
+
+@article{ar:optimjl,
+	doi = {10.21105/joss.00615},
+	url = {https://doi.org/10.21105/joss.00615},
+	year  = {2018},
+	month = {apr},
+	publisher = {The Open Journal},
+	volume = {3},
+	number = {24},
+	pages = {615},
+	author = {Patrick K Mogensen and Asbj{\o}rn N Riseth},
+	title = {Optim: A mathematical optimization package for Julia},
+	journal = {Journal of Open Source Software}
+}
+
+@article{ar:alldielectric,
+	title={All-dielectric metamaterials},
+	author={Jahani, Saman and Jacob, Zubin},
+	journal={Nature Nanotechnology},
+	volume={11},
+	number={1},
+	pages={23--36},
+	year={2016},
+	publisher={Nature Research},
+	doi={10.1038/nnano.2015.304}
+}
diff --git a/paper/paper.md b/paper/paper.md
new file mode 100644
index 0000000..779c925
--- /dev/null
+++ b/paper/paper.md
@@ -0,0 +1,56 @@
+---
+title: 'ParticleScattering: Solving and optimizing multiple-scattering problems in Julia'
+tags:
+- Computational electromagnetics
+- Scattering
+- Julia
+authors:
+- name: Boaz Blankrot
+  orcid: 0000-0003-3364-9298
+  affiliation: 1 # (Multiple affiliations must be quoted: "1, 2")
+- name: Clemens Heitzinger
+  orcid: 0000-0003-1613-5164
+  affiliation: 1
+affiliations:
+- name: Vienna University of Technology, A-1040 Vienna, Austria
+  index: 1
+date: 8 April 2018
+bibliography: paper.bib
+---
+
+# Summary
+
+[ParticleScattering](https://github.com/bblankrot/ParticleScattering.jl) is a Julia package for computing the electromagnetic fields scattered by a large number of two-dimensional particles, as well as optimizing particle parameters for various applications.
+Such problems naturally arise in the design and analysis of metamaterials, including photonic crystals [@ar:alldielectric].
+Unlike most solvers for these problems, ours does not require a periodic structure and is scalable to a large number of particles.
+In particular, this software is designed for scattering problems involving TM plane waves impinging on a collection of homogeneous dielectric particles with arbitrary smooth shapes.
+Our code performs especially well when the number of particles is substantially larger than the number of distinct shapes, where particles are considered indistinct if they are identical up to rotation.
+
+## Solver overview
+Given a scattering problem consisting of a collection of penetrable particles in a homogeneous medium, the software performs the following steps to calculate the total electric field:
+
+- For each distinct non-circular shape, a single- and double-layer potential formulation is constructed.
+- The potential formulations are transformed to a multipole basis of Hankel functions, reducing the degrees of freedom by at least an order of magnitude.
+- Analytical multipole basis is computed for circular particles.
+- A multiple-scattering system of equations is constructed, and then solved with the Fast Multipole Method.
+- Electric field is computed at any point of interest.
+
+In addition, ParticleScattering can plot near- and far-field results using the popular framework [PyPlot](https://github.com/JuliaPy/PyPlot.jl), create publication-level plots with [PGFPlots](https://github.com/KristofferC/PGFPlotsX.jl) integration, and compute minimum parameters for a desired error level.
+
+## Optimization
+
+ParticleScattering is especially targeted at users who wish to design metamaterials belonging to the class described above.
+While the large number of variables such metamaterials contain allows for a variety of devices that meet different objectives, it also creates a large search space for choosing them. Therefore, a fast and automated approach can be beneficial for both inventing new designs and improving existing ones.
+As the results of many ParticleScattering computations can be recycled between optimization iterations, a large number of parameters can be optimized simultaneously in reasonable time.
+ParticleScattering performs gradient-based optimization of rotation angle for arbitrarily-shaped particles, and radius of circular particles, in conjunction with the Optim optimization package [@ar:optimjl], where the objective is to minimize or maximize the electric field intensity at chosen points.
+Figure 1 shows an example of angle optimization of 20 particles, where the objective is the electric field intensity at the origin. From left to right, we see the electric field before optimization, after minimization, and after maximization. The field intensity at the origin is clearly different in both optimization results, with minimization decreasing the intensity by 95%, and maximization increasing it by over 700%. The total run time for both optimizations and all necessary precomputations was 35 seconds.
+
+![Scattering problem before optimization, after minimization, and after maximization.](fig.png)
+
+For a detailed description of our approach, including several numerical examples generated by ParticleScattering, see our recent publication [@ar:blankrot2018].
+
+# Acknowledgments
+
+This work was supported by the Austrian Science Fund (FWF) through the START Project Y 660 *PDE Models for Nanotechnology*.
+
+# References
diff --git a/sim1_complexity.jld b/sim1_complexity.jld
new file mode 100644
index 0000000..ca996d0
Binary files /dev/null and b/sim1_complexity.jld differ
diff --git a/src/ParticleScattering.jl b/src/ParticleScattering.jl
index 7dfa316..224f00e 100644
--- a/src/ParticleScattering.jl
+++ b/src/ParticleScattering.jl
@@ -27,8 +27,8 @@ include("utilities.jl")
 include("minimum_N_P.jl")
 
 #methods, shapes.jl
-export rounded_star, squircle, ellipse, square_grid, rect_grid, randpoints,
-    luneburg_grid, verify_min_distance
+export rounded_star, squircle, ellipse, square_grid, rect_grid, hex_grid,
+    randpoints, luneburg_grid, verify_min_distance
 #methods, multipole.jl
 export solve_particle_scattering, scattered_field_multipole
 #methods, minimum_N_P.jl
@@ -48,13 +48,11 @@ export find_border, uniqueind
 export optimize_φ
 #methods, optimize_rs.jl
 export optimize_radius
-
-
 # For advanced plotting with pgfplots
 import DataFrames, CSV, PGFPlotsX; const pgf = PGFPlotsX
 include("visualization_pgf.jl")
 #methods, visualization_pgf.jl
-export plotNearField_pgf, drawShapes_pgf
+export plot_near_field_pgf
 
 #temp
 export divideSpace, FMMtruncation, particleExpansion, FMMbuildMatrices
diff --git a/src/scattering.jl b/src/scattering.jl
index 65c7acf..a08dcfb 100644
--- a/src/scattering.jl
+++ b/src/scattering.jl
@@ -17,11 +17,10 @@ function get_potential(kout, kin, P, t, ft, dft)
 
     #assuming the wave is sampled on the shape
     nz = sqrt.(sum(abs2,ft,2))
+    θ = atan2.(ft[:,2], ft[:,1])
     ndz = sqrt.(sum(abs2,dft,2))
     nzndz = nz.*ndz
-    #precompute? or negative = (-1)^n positive?
-    #no precompute
-	wro = dft[:,2].*ft[:,1] - dft[:,1].*ft[:,2]
+    wro = dft[:,2].*ft[:,1] - dft[:,1].*ft[:,2]
 	zz = dft[:,1].*ft[:,1] + dft[:,2].*ft[:,2]
 
     bessp = besselj.(-P-1, kout*nz)
@@ -31,13 +30,46 @@ function get_potential(kout, kin, P, t, ft, dft)
     for p = -P:P
         bess[:] = besselj.(p, kout*nz)
 		du[:] = kout*bessp.*wro - (p*bess./nz).*(wro + 1im*zz)
-        rhs[:] = -[bess.*exp.(1.0im*p*t);
-               (du./nzndz).*exp.(1.0im*p*t)]
+        rhs[:] = -[bess.*exp.(1.0im*p*θ);
+               (du./nzndz).*exp.(1.0im*p*θ)]
         sigma_mu[:,p + P + 1] = LU\rhs
         copy!(bessp, bess)
     end
     return sigma_mu
 end
+#
+# function get_potential_deprecated(kout, kin, P, t, ft, dft)
+# #problematic assumption here that t=θ in the Jacobi-Anger expansion (for rhs).
+# #this means that we assume all shapes are given in polar form (r(θ)cosθ,r(θ)sinθ)),
+# #which is not the case for the simple ellipse (r1cosθ,r2sinθ).
+#     N = length(t) #N here is 2N elesewhere.
+#
+#     A = SDNTpotentialsdiff(kout, kin, t, ft, dft)
+#     LU = lufact(A)
+#
+#     sigma_mu = Array{Complex{Float64}}(2*N, 2*P+1)
+#
+#     #assuming the wave is sampled on the shape
+#     nz = sqrt.(sum(abs2,ft,2))
+#     ndz = sqrt.(sum(abs2,dft,2))
+#     nzndz = nz.*ndz
+#     wro = dft[:,2].*ft[:,1] - dft[:,1].*ft[:,2]
+# 	zz = dft[:,1].*ft[:,1] + dft[:,2].*ft[:,2]
+#
+#     bessp = besselj.(-P-1, kout*nz)
+#     bess = similar(bessp)
+#     du = Array{Complex{Float64}}(length(bessp))
+#     rhs = Array{Complex{Float64}}(2*length(bessp))
+#     for p = -P:P
+#         bess[:] = besselj.(p, kout*nz)
+# 		du[:] = kout*bessp.*wro - (p*bess./nz).*(wro + 1im*zz)
+#         rhs[:] = -[bess.*exp.(1.0im*p*t);
+#                (du./nzndz).*exp.(1.0im*p*t)]
+#         sigma_mu[:,p + P + 1] = LU\rhs
+#         copy!(bessp, bess)
+#     end
+#     return sigma_mu
+# end
 
 """
     get_potential(kout, kin, P, t, ft, dft) -> sigma_mu
@@ -201,10 +233,12 @@ function scatteredfield(sigma_mu, k, t, ft, dft, p)
     u_s = SDout*sigma_mu
 end
 
+
 function shapeMultipoleExpansion(k, t, ft, dft, P)
     #unlike others (so far), this does *not* assume t_j=pi*j/N
     N = div(length(t),2)
     nz = vec(sqrt.(sum(abs2,ft,2)))
+    θ = atan2.(ft[:,2], ft[:,1])
     ndz = vec(sqrt.(sum(abs2,dft,2)))
     AB = Array{Complex{Float64}}(2*P + 1, 4*N)
     bessp = besselj.(-P-1,k*nz)
@@ -212,20 +246,46 @@ function shapeMultipoleExpansion(k, t, ft, dft, P)
     for l = -P:0
         bess[:] = besselj.(l,k*nz)
         for j = 1:2*N
-            AB[l+P+1,j] = 0.25im*(π/N)*bess[j]*exp(-1.0im*l*t[j])*ndz[j]
-            l!=0 && (AB[-l+P+1,j] = 0.25im*((-1.0)^l*π/N)*bess[j]*exp(1.0im*l*t[j])*ndz[j])
+            AB[l+P+1,j] = 0.25im*(π/N)*bess[j]*exp(-1.0im*l*θ[j])*ndz[j]
+            l != 0 && (AB[-l+P+1,j] = 0.25im*((-1.0)^l*π/N)*bess[j]*exp(1.0im*l*θ[j])*ndz[j])
             wro = ft[j,1]*dft[j,2] - ft[j,2]*dft[j,1]
             zdz = -1.0im*(ft[j,1]*dft[j,1] + ft[j,2]*dft[j,2])
             b1 = (-l*bess[j]/nz[j])*(zdz + wro)
             b1_ = (-l*bess[j]/nz[j])*(zdz - wro)
             b2 = k*bessp[j]*wro
-            AB[l+P+1,j+2*N] = 0.25im*(π/N)*(exp(-1.0im*l*t[j])/nz[j])*(b1 + b2)
-            l!=0 && (AB[-l+P+1,j+2*N] = 0.25im*((-1.0)^l*π/N)*(exp(1.0im*l*t[j])/nz[j])*(-b1_ + b2))
+            AB[l+P+1,j+2*N] = 0.25im*(π/N)*(exp(-1.0im*l*θ[j])/nz[j])*(b1 + b2)
+            l != 0 && (AB[-l+P+1,j+2*N] = 0.25im*((-1.0)^l*π/N)*(exp(1.0im*l*θ[j])/nz[j])*(-b1_ + b2))
         end
         copy!(bessp,bess)
     end
     return AB
 end
+#
+# function shapeMultipoleExpansion_deprecated(k, t, ft, dft, P)
+#     #unlike others (so far), this does *not* assume t_j=pi*j/N
+#     N = div(length(t),2)
+#     nz = vec(sqrt.(sum(abs2,ft,2)))
+#     ndz = vec(sqrt.(sum(abs2,dft,2)))
+#     AB = Array{Complex{Float64}}(2*P + 1, 4*N)
+#     bessp = besselj.(-P-1,k*nz)
+#     bess = similar(bessp)
+#     for l = -P:0
+#         bess[:] = besselj.(l,k*nz)
+#         for j = 1:2*N
+#             AB[l+P+1,j] = 0.25im*(π/N)*bess[j]*exp(-1.0im*l*t[j])*ndz[j]
+#             l!=0 && (AB[-l+P+1,j] = 0.25im*((-1.0)^l*π/N)*bess[j]*exp(1.0im*l*t[j])*ndz[j])
+#             wro = ft[j,1]*dft[j,2] - ft[j,2]*dft[j,1]
+#             zdz = -1.0im*(ft[j,1]*dft[j,1] + ft[j,2]*dft[j,2])
+#             b1 = (-l*bess[j]/nz[j])*(zdz + wro)
+#             b1_ = (-l*bess[j]/nz[j])*(zdz - wro)
+#             b2 = k*bessp[j]*wro
+#             AB[l+P+1,j+2*N] = 0.25im*(π/N)*(exp(-1.0im*l*t[j])/nz[j])*(b1 + b2)
+#             l!=0 && (AB[-l+P+1,j+2*N] = 0.25im*((-1.0)^l*π/N)*(exp(1.0im*l*t[j])/nz[j])*(-b1_ + b2))
+#         end
+#         copy!(bessp,bess)
+#     end
+#     return AB
+# end
 
 function solvePotential_forError(kin, kout, shape, ls_pos, ls_amp, θ_i)
     #plane wave outside, line sources inside
diff --git a/src/shapes.jl b/src/shapes.jl
index a2ebd8b..0bf59b2 100644
--- a/src/shapes.jl
+++ b/src/shapes.jl
@@ -40,8 +40,7 @@ Return a `ShapeParams` object containing the shape parametrized by
 `(x/r1)^2 + (y/r2)^2 = 1` with 2`N` nodes.
 """
 function ellipse(r1, r2, N)
-	warn("ellipse currently exhibiting abnormal behavior in cylindrical harmonics (minimumP does not converge)")
-    t = Float64[pi*j/N for j=0:(2*N-1)]
+	t = Float64[pi*j/N for j=0:(2*N-1)]
     ft = [r1*cos.(t) r2*sin.(t)]
     dft = [(-r1)*sin.(t) r2*cos.(t)]
     return ShapeParams(t, ft, dft)
@@ -75,6 +74,42 @@ function rect_grid(a::Integer, b::Integer, dx, dy)
     centers = [repeat(xpoints, outer=[b]) 	repeat(ypoints, inner=[a])]
 end
 
+"""
+	hex_grid(a::Integer, b::Integer, d)
+
+Return `centers`, an `(M,2)` array  containing the points on a hexagonal lattice
+with horizontal rows, with `a` points in each row and `rows` rows, distanced `d`.
+If `minus1` is true, the last point in every odd row is omitted.
+"""
+function hex_grid(a::Integer, rows::Integer, d; minus1 = false)
+	h = d*sqrt(0.75) #row height
+	M = minus1 ? a*rows - div(rows,2) : a*rows
+	centers = Array{Float64}(M,2)
+	ind = 1
+	for r = 0:rows-1
+		if mod(r,2) == 0
+			for i = 1:a
+				centers[ind,1] = i*d
+				centers[ind,2] = (r - 1)*h
+				ind += 1
+			end
+		else
+			for i = 1:a-1
+				centers[ind,1] = (i + 0.5)*d
+				centers[ind,2] = (r - 1)*h
+				ind += 1
+			end
+			if !minus1
+				centers[ind,1] = (a + 0.5)*d
+				centers[ind,2] = (r - 1)*h
+				ind += 1
+			end
+		end
+	end
+	offset = mean(centers, 1)
+    centers .-= offset
+end
+
 """
 	randpoints(M, dmin, width, height; failures = 100)
 
@@ -129,8 +164,8 @@ function randpoints(M, dmin, width, height, points; failures = 100)
         dist2 <= dmin2 && error("randpoints: given points have distance <= dmin")
     end
 
-    x_res = [rand(Float64)*width]
-    y_res = [rand(Float64)*height]
+    x_res = Array{Float64}(0)
+    y_res = Array{Float64}(0)
     fail = 0
     while fail < failures && length(x_res) < M
         accepted = true
diff --git a/src/visualization_pgf.jl b/src/visualization_pgf.jl
index 4cef640..d04555e 100644
--- a/src/visualization_pgf.jl
+++ b/src/visualization_pgf.jl
@@ -1,4 +1,27 @@
-function plotNearField_pgf(filename, k0, kin, P, sp::ScatteringProblem, θ_i = 0;
+"""
+    plot_near_field_pgf(filename, k0, kin, P, sp::ScatteringProblem, θ_i = 0;
+                        opt::FMMoptions = FMMoptions(), use_multipole = true,
+                        x_points = 201, y_points = 201, border = find_border(sp),
+                        downsample = 1, include_preamble = false, normalize = 1.0)
+
+Plots the total electric field as a result of a plane wave with incident
+angle `θ_i` scattering from the ScatteringProblem `sp`, using pgfplots's
+`surf`. Can accept number of sampling points in each direction, and either
+a given `border` or calculate it automatically. The plots of the shapes (but not
+the field) can be downsampled by setting an integer `downsample`, since pgfplots
+slows down dramatically when drawing many shapes with many nodes.
+
+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
+harmonics space. Normalizes all distances and sizes in plot by `normalize`.
+
+Saves the generated pgfplots file to `filename`, with just a surrounding `tikzpicture`
+environment if `include_preamble=false`, and a compilable tandalone document
+otherwise.
+"""
+function plot_near_field_pgf(filename, k0, kin, P, sp::ScatteringProblem, θ_i = 0;
                             opt::FMMoptions = FMMoptions(), use_multipole = true,
                             x_points = 201, y_points = 201, border = find_border(sp),
                             downsample = 1, include_preamble = false, normalize = 1.0)
@@ -21,31 +44,30 @@ function plotNearField_pgf(filename, k0, kin, P, sp::ScatteringProblem, θ_i = 0
     dt = DataFrames.DataFrame(x = points[:,1]/normalize, y = points[:,2]/normalize, z = abs.(Ez))
     CSV.write(filename * ".dat", dt, delim = '\t')
     pgf.@pgf begin
-        p = pgf.Plot3(pgf.Table(basename(filename) * ".dat"),
-            {   surf,
-                no_markers,
-                shader = "interp",
-                "mesh/rows" = y_points}, incremental = false)
-        ax = pgf.Axis(p,
-         {  xmin = x_min/normalize,
-            xmax = x_max/normalize,
-            ymin = y_min/normalize,
-            ymax = y_max/normalize,
-            xlabel = "\$x\$",
-            ylabel = "\$y\$",
-            scale_only_axis,
-            width = "\\linewidth",
-            height = "$aspect*\\linewidth",
-            view = "{0}{90}",
-            "mesh/ordering" = "y varies",
-            colorbar,
-            point_meta_max = maximum(abs.(Ez))})
+        ax = pgf.Axis({ xmin = x_min/normalize,
+                        xmax = x_max/normalize,
+                        ymin = y_min/normalize,
+                        ymax = y_max/normalize,
+                        xlabel = "\$x/$normalize\$",
+                        ylabel = "\$y/$normalize\$",
+                        scale_only_axis,
+                        width = "\\linewidth",
+                        height = "$aspect*\\linewidth",
+                        view = "{0}{90}",
+                        "mesh/ordering" = "y varies",
+                        colorbar,
+                        point_meta_max = maximum(abs.(Ez))})
+        #new release of pgfplotsx isn't printing table correctly (too many newlines for pgfplots)
+        push!(ax, "\\addplot3[surf, no markers, shader={interp},
+                    mesh/rows={$(y_points)}] table{$(basename(filename)).dat};")
+
     end
-    drawShapes_pgf(shapes, centers, ids, φs, ax, 1.1*maximum(abs.(Ez)), downsample, normalize = normalize)
+    draw_shapes_pgf(shapes, centers, ids, φs, ax, 1.1*maximum(abs.(Ez)), downsample, normalize = normalize)
     pgf.save(filename, ax ,include_preamble = include_preamble)
 end
 
-function drawShapes_pgf(shapes, centers, ids, φs, ax, floating, downsample; normalize = 1.0)
+function draw_shapes_pgf(shapes, centers, ids, φs, ax, floating, downsample; normalize = 1.0)
+    #draw in 3d so it is "above" surf plot
     for ic = 1:size(centers,1)
         if typeof(shapes[ids[ic]]) == ShapeParams
             if φs[ic] == 0.0
@@ -54,9 +76,10 @@ function drawShapes_pgf(shapes, centers, ids, φs, ax, floating, downsample; nor
                 Rot = [cos(φs[ic]) sin(φs[ic]);-sin(φs[ic]) cos(φs[ic])]
                 ft_rot = shapes[ids[ic]].ft[1:downsample:end,:]*Rot .+ centers[ic,:]'
             end
-            push!(ax, pgf.Plot3(pgf.Coordinates([ft_rot[:,1];ft_rot[1,1]]/normalize,
-                [ft_rot[:,2];ft_rot[1,2]]/normalize, floating*ones(Float64,size(ft_rot,1)+1)),
-                "black", "no markers", "thick", incremental = false))
+            co = pgf.Coordinates([ft_rot[:,1];ft_rot[1,1]]/normalize,
+                                 [ft_rot[:,2];ft_rot[1,2]]/normalize,
+                                 floating*ones(size(ft_rot,1) + 1))
+            pgf.@pgf push!(ax, pgf.Plot3({black, no_markers, thick}, co))
         else
             x = centers[ic,1]/normalize
             y = centers[ic,2]/normalize
diff --git a/test/scatteredfield_test.jl b/test/scatteredfield_test.jl
index 612b4bc..2698790 100644
--- a/test/scatteredfield_test.jl
+++ b/test/scatteredfield_test.jl
@@ -6,7 +6,8 @@
     N = 1000
     shapes = [squircle(1, N);
             rounded_star(0.8, 0.2, 3, N);
-            rounded_star(0.75, 0.25, 5, N)]
+            rounded_star(0.75, 0.25, 5, N);
+            ellipse(0.2, 0.9, N)]
     for s in shapes
         sigma_mu = get_potentialPW(k0, kin, s, θ_i)
         dx = 2*norm(s.ft[1,:] - s.ft[2,:])