diff --git a/src/ParticleInCell.jl b/src/ParticleInCell.jl
index b7d3040..814de14 100644
--- a/src/ParticleInCell.jl
+++ b/src/ParticleInCell.jl
@@ -44,7 +44,7 @@ end
 export step!
 
 include("field_utils.jl")
-export FiniteDifferenceToEdges, AverageEdgesToNodes
+export FiniteDifferenceToEdges, AverageEdgesToNodes, ZeroField, MultiplyField
 
 include("poisson.jl")
 export PoissonSolveFFT
diff --git a/src/field_utils.jl b/src/field_utils.jl
index 2a81e28..d747bdf 100644
--- a/src/field_utils.jl
+++ b/src/field_utils.jl
@@ -68,3 +68,10 @@ struct ZeroField{F} <: AbstractSimulationStep
 end
 
 step!(step::ZeroField) = step.field.values .= 0
+
+struct MultiplyField{F,T} <: AbstractSimulationStep
+    field::F
+    x::T
+end
+
+step!(step::MultiplyField) = step.field.values .*= step.x
diff --git a/src/interpolation.jl b/src/interpolation.jl
index 174c92b..865c30d 100644
--- a/src/interpolation.jl
+++ b/src/interpolation.jl
@@ -163,6 +163,16 @@ struct BSplineChargeInterpolation{S,F,IF} <: AbstractSimulationStep
             interp_func,
         )
     end
+
+    function BSplineChargeInterpolation(
+        species::S,
+        rho::Field{T,N,NodeOffset},
+        interp_width,
+        interp_func::F,
+    ) where {T,N,S,F}
+
+        new{S,typeof(rho),F}(species, rho, -1, interp_width, interp_func)
+    end
 end
 
 function step!(step::BSplineChargeInterpolation)
diff --git a/src/particle_updaters/electrostatic.jl b/src/particle_updaters/electrostatic.jl
index 29bb943..9390798 100644
--- a/src/particle_updaters/electrostatic.jl
+++ b/src/particle_updaters/electrostatic.jl
@@ -56,7 +56,7 @@ function step!(step::ElectrostaticParticlePush)
         )
 
         for I in Is
-            grid_cell_coord = cell_index_to_cell_coords(step.E, I)
+            grid_cell_coord = cell_index_to_cell_coords(step.E, I, 1)
             dist = Tuple(particle_cell_coord .- grid_cell_coord)
             interp_weights = step.interpolation_function.(dist)
             interp_weight = prod(interp_weights)
diff --git a/src/poisson.jl b/src/poisson.jl
index cdb265b..e2ce48d 100644
--- a/src/poisson.jl
+++ b/src/poisson.jl
@@ -1,3 +1,7 @@
+field_solve_three_pt(k, dx) = k^2 * sinc(k * dx / 2 / pi)^2
+field_solve_lagrange(k, dx) = k^2 * sinc(k * dx / 2 / pi)^2 * (2 + cos(k * dx)) / 3
+field_solve_five_pt(k, dx) = -1 * k^2 * sinc(k * dx / 2 / pi)^2 * (-7 + cos(k * dx)) / 6
+
 struct PoissonSolveFFT{T,D,G,P,F<:AbstractField} <: AbstractSimulationStep
     rho::F
     phi::F
@@ -7,7 +11,11 @@ struct PoissonSolveFFT{T,D,G,P,F<:AbstractField} <: AbstractSimulationStep
 
     fft_plan::P
 
-    function PoissonSolveFFT(rho::F, phi::F) where {T,N,O,D,G,F<:AbstractField{T,N,O,D,G}}
+    function PoissonSolveFFT(
+        rho::F,
+        phi::F,
+        k2s = field_solve_three_pt,
+    ) where {T,D,G,F<:AbstractField{T,D,G}}
         # This restriction could possibly be relaxed to just require compatible grids...
         @assert rho.grid === phi.grid
         # Currently only supports periodic boundary conditions...
@@ -31,9 +39,12 @@ struct PoissonSolveFFT{T,D,G,P,F<:AbstractField} <: AbstractSimulationStep
                 continue
             end
 
-            ks = 2π .* (It .- 1) ./ sim_lengths
-            grid_angles = ks .* cell_lengths ./ 2
-            inv_Ksqs = (cell_lengths ./ (2 .* sin.(grid_angles))) .^ 2 ./ epsilon_0
+            # ks = 2π .* (It .- 1) ./ sim_lengths
+            ks =
+                2π .*
+                ifelse.(It .<= size(Ksq_inv) ./ 2, It .- 1, It .- size(Ksq_inv) .- 1) ./
+                sim_lengths
+            inv_Ksqs = (k2s.(ks, cell_lengths)) .^ (-1) ./ epsilon_0
 
             Ksq_inv[I] = prod(inv_Ksqs)
         end