diff --git a/doc/capsule-boundary-conditions.wl b/doc/capsule-boundary-conditions.wl
index 1aed7371..2e62c86c 100644
--- a/doc/capsule-boundary-conditions.wl
+++ b/doc/capsule-boundary-conditions.wl
@@ -30,7 +30,6 @@
 (*F[f_,g_][x_,y_]:= f[x] g[y]-f[y] g[x]*)
 (*Fd[f_,g_][x_,y_]:= f[x] g'[y]-f'[y] g[x]*)
 (**)
-(**)
 (*Yd[x_,y_]=  Fd[Subscript[H, n],Subscript[J, n]] [x,y];*)
 (*Y[x_,y_]=  F[Subscript[H, n],Subscript[J, n]] [x,y];*)
 (*Ydd[x_,y_]=  F[Subscript[H, n]',Subscript[J, n]'] [x,y];*)
@@ -214,9 +213,32 @@
 (*eqs = {\[Psi][1]  == \[Psi][0]/.r-> a[0],D[\[Psi][1],r]/\[Rho][1]  == D[\[Psi][0],r]/\[Rho][0]/.r-> a[0], \[Psi][s] + \[Psi]o  == \[Psi][1]/.r-> a[1]};*)
 (*subsol = Solve[eqs,{f[0,n],f[1,n],A[1,n]}] /.{\[Rho]o -> qo ko,\[Rho][1] -> q[1] k[1],\[Rho][0] -> q[0] k[0]};*)
 (**)
-(*f[0,n]/.subsol [[1]]//Simplify*)
-(*f[1,n]/.subsol [[1]]//Simplify*)
-(*A[1,n]/.subsol [[1]]//Simplify*)
+(**)
+(**)
+(**)
+
+
+(* ::InheritFromParent:: *)
+(**)
+
+
+(* ::InheritFromParent:: *)
+(*force = Ao[n] Subscript[H, n][ko a[1]]+fo[n] Subscript[J, n][ko a[1]];*)
+(**)
+(*numer = q[0]/q[1]  Yd[a[0] k[1],a[0] k[1]];*)
+(*denom = q[0]/q[1]  Subscript[J, n][a[0] k[0]] Yd[a[1] k[1],a[0] k[1]] -Derivative[1][Subscript[J, n]][a[0] k[0]] Y[a[1] k[1],a[0] k[1]] ;*)
+(*force numer/denom  == f[0,n]/.subsol [[1]]//Simplify*)
+(**)
+(*numer =  q[0]/q[1] Subscript[J, n][a[0] k[0]] Derivative[1][Subscript[H, n]][a[0] k[1]]- Subscript[H, n][a[0] k[1]] Derivative[1][Subscript[J, n]][a[0] k[0]];*)
+(*denom = Y[a[1] k[1],a[0] k[1]] Derivative[1][Subscript[J, n]][a[0] k[0]]-q[0]/q[1] Subscript[J, n][a[0] k[0]] Yd[a[1] k[1],a[0] k[1]];*)
+(*force numer/denom  ==f[1,n]/.subsol [[1]]//Simplify*)
+(**)
+(*numer =  Subscript[J, n][a[0] k[1]] Derivative[1][Subscript[J, n]][a[0] k[0]]-q[0]/q[1] Subscript[J, n][a[0] k[0]] Derivative[1][Subscript[J, n]][a[0] k[1]];*)
+(*denom = Y[a[1] k[1],a[0] k[1]] Derivative[1][Subscript[J, n]][a[0] k[0]]-q[0]/q[1] Subscript[J, n][a[0] k[0]] Yd[a[1] k[1],a[0] k[1]];*)
+(*force  numer/denom== A[1,n]/.subsol [[1]]//Simplify*)
+(**)
+(**)
+(* Y[a[1] k[1],a[0] k[1]]*)
 
 
 (* ::Input:: *)
@@ -354,6 +376,8 @@ the background material.
 
 
 
+
+
 (* ::Section:: *)
 (*Addition Theorems*)
 
diff --git a/src/physics/acoustics/concentric_capsule.jl b/src/physics/acoustics/concentric_capsule.jl
index ed9f298a..978abbd7 100644
--- a/src/physics/acoustics/concentric_capsule.jl
+++ b/src/physics/acoustics/concentric_capsule.jl
@@ -3,20 +3,57 @@
 function internal_field(x::SVector{2,T}, p::CapsuleParticle{T,2,Acoustic{T,2},Circle{T}}, sim::FrequencySimulation{T,2,Acoustic{T,2}}, ω::T, scattering_coefficients::AbstractVector;
         basis_order::Int=5) where T
 
-    medium = sim.medium
     Nh = basis_order
 
-    if iszero(p.medium.c) || isinf(abs(p.medium.c))
+    if iszero(p.outer.medium.c) || isinf(abs((p.outer.medium.c)))
         return zero(Complex{T})
-    else
-        r = outer_radius(p)
-        k = ω/medium.c
-        kp = ω/p.medium.c
-        Z = - t_matrix(p.shape, p.medium, medium, ω, basis_order)
+    elseif norm(x - origin(p)) > outer_radius(p)
+        warn("point $x not insider particle $p. Returning zero.")
+        return zero(Complex{T})
+    end
+
+    k = ω / medium.c
+    k0 = ω / p.inner.medium.c
+    k1 = ω / p.outer.medium.c
+    a0 = outer_radius(p.inner)
+    a1 = outer_radius(p.outer)
+    q0 = (p.inner.medium.ρ*p.inner.medium.c)/(p.outer.medium.ρ*p.outer.medium.c)
+
+    Yn(n::Integer) = hankelh1(n,k1*a1)*besselj(n,k1*a0) - hankelh1(n,k1*a0)*besselj(n,k1*a1)
+    Ydn(n::Integer,x::Complex{T},y::Complex{T}) = hankelh1(n,x)*diffbesselj(n,y) - diffhankelh1(n,y)*besselj(n,x)
+
+    force(n) = scattering_coefficients[m+Nh+1] * hankelh1(n, k*a1) +
+        sim.source.coef(n,origin(p),ω) * besselj(n, k*a1);
+
+    if norm(x - origin(p)) <= outer_radius(p.inner)
+        function coef(n::Int)
+            numer = q0*Ydn(n, a0*k1, a0*k1)
+            denom = q0*besselj(n, a0*k0)*Ydn(n, a1*k1, a0*k1) - diffbesselj(n, a0*k0)*Yn(n)
+            return force(n) * numer / denom
+        end
+        basis = basis_function(p.inner, ω)
         return map(-Nh:Nh) do m
-             scattering_coefficients[m+Nh+1] / (Z[m+Nh+1,m+Nh+1]*besselj(m,kp*r)) * (Z[m+Nh+1,m+Nh+1]*hankelh1(m,k*r) - besselj(m,k*r))
+             basis(m,x) * coef(n)
         end
     end
+    # else calculate field in mid region
+    function J_coef(n::Int)
+        numer = q0*besselj(n, a0*k0)*diffhankelh1(n, a0*k1) - hankelh1(n, a0*k1)*diffbesselj(n, a0*k0)
+        denom = Yn(n)*diffbesselj(n, a0*k0) - q0*besselj(n, a0*k0)*Ydn(n, a1*k1, a0*k1)
+        return force(n) * numer / denom
+    end
+    function H_coef(n::Int)
+        numer =  besselj(n, a0*k1)*diffbesselj(n,a0*k0) - q0*besselj(n, a0*k0)*diffbesselj(n,a0*k1)
+        denom = Yn(n)*diffbesselj(n,a0*k0) - q0*besselj(n, a0*k0)*Ydn(n, a1*k1, a0*k1)
+        return force(n) * numer / denom
+    end
+
+    J_basis = basis_function(p.outer, ω)
+    H_basis = basis_function(p.outer.medium, ω)
+    # return map(-Nh:Nh) do m
+    #      basis(m,x) * coef(n)
+    # end
+
 end
 
 function t_matrix(cap::CapsuleParticle{T,2,Acoustic{T,2},Circle{T}}, medium::Acoustic{T,2}, ω::T, M::Integer)::Diagonal{Complex{T}} where T <: AbstractFloat