Sphere in a point charge field.

In general the LaplaceBIE allows to find field on the object if the potential in absence of object(s) is known. To demonstrate that let's consider a field on the dielectric sphere (radius $r_1$) due to a point charge at distance $\zeta$ away. The solution inside the sphere is given as series [Straton page 221]:

\[\psi = \sum_{n=0}^{\infty} \frac{2n + 1}{\epsilon n + n + 1} \frac{r_1^n}{\zeta^{n+1}} L_n(\cos \theta)\]

where $L_n$ is Legendre polynomial and $\theta$ angle from line connecting point charge and spehre.

using Jacobi
function ψt(cosθ,ϵ,ζ,r1)
    s = 0
    for n in 0:25
        s += (2*n + 1)/(n*ϵ + n + 1)*r1^n/ζ^(n+1)*legendre(cosθ,n)
    end
    return s
end

ψt (generic function with 1 method)

The normal derivative can also be given for the sphere which coincides with derivative in the radial direction

\[\frac{\partial \psi}{\partial n} = \sum_{n=0}^{\infty} \frac{n(2n + 1)}{\epsilon n + n + 1} \frac{r_1^{n-1}}{\zeta^{n+1}} L_n(\cos \theta)\]

function ∇ψn(cosθ,ϵ,ζ,r1)
    s = 0
    for n in 1:25
        s += n*(2*n + 1)/(n*ϵ + n + 1)*r1^(n-1)/ζ^(n+1) * legendre(cosθ,n)
    end
    return s
end

∇ψn (generic function with 1 method)

First we can define the exterior potential in absence of objects

using LinearAlgebra

ζ = 1.2
r1 = 1
y = [0,0,ζ]

ψ0(x) = 1/norm(x.-y)
∇ψ0(x) = -(x.-y)/norm(x.-y)^3

∇ψ0 (generic function with 1 method)

Let's now place a sphere in the field with relative permitivity $\epsilon=10$.

ϵ = 10

include("sphere.jl")
msh = unitsphere(2)
vertices, faces = msh.vertices, msh.faces
n = vertices

162-element Array{GeometryTypes.Point{3,Float64},1}:
 [-0.5257311121191336, 0.85065080835204, 0.0]
 [0.5257311121191336, 0.85065080835204, 0.0]
 [-0.5257311121191336, -0.85065080835204, 0.0]
 [0.5257311121191336, -0.85065080835204, 0.0]
 [0.0, -0.5257311121191336, 0.85065080835204]
 [0.0, 0.5257311121191336, 0.85065080835204]
 [0.0, -0.5257311121191336, -0.85065080835204]
 [0.0, 0.5257311121191336, -0.85065080835204]
 [0.85065080835204, 0.0, -0.5257311121191336]
 [0.85065080835204, 0.0, 0.5257311121191336]
 ⋮
 [0.25989191300775444, 0.4338885645526948, -0.8626684804161862]
 [0.26286555605956685, 0.16245984811645317, -0.9510565162951536]
 [0.26286555605956685, -0.16245984811645317, -0.9510565162951536]
 [0.9619383577839176, 0.0, 0.2732665289126717]
 [0.9510565162951536, 0.26286555605956685, 0.16245984811645317]
 [0.8626684804161862, 0.25989191300775444, -0.4338885645526948]
 [0.9510565162951536, 0.26286555605956685, -0.16245984811645317]
 [0.6937804775604491, 0.7020464447761631, 0.16062203564002314]
 [0.8506508083520399, 0.5257311121191337, 0.0]

With LaplaceBIE to calculate the surface field we execute theese three lines:

using LaplaceBIE

ψ = surfacepotential(vertices,n,faces,ϵ,ψ0);
P∇ψ = tangentderivatives(vertices,n,faces,ψ);
n∇ψ = normalderivatives(vertices,n,faces,P∇ψ,ϵ,∇ψ0);

162-element Array{Float64,1}:
 -0.06205534122187392
 -0.062055341221873825
 -0.06205534122187398
 -0.06205534122187367
  0.09154529388081109
  0.09154529388081074
 -0.06365584711404197
 -0.0636558471140412
 -0.06367958215679821
 -0.04611694351733614
  ⋮
 -0.06373122123157465
 -0.06363392295695912
 -0.06363392295695924
 -0.05852127628110055
 -0.06092044651491243
 -0.06376068401697482
 -0.06308197836685187
 -0.06102685443985177
 -0.06187692442427096

Now to compare with analytics we use azimuthal symetry which in numerics was choosem as x axis.

using Winston

cosθ = range(-1,1,length=100)
cosθs = [x[3]/1 for x in vertices]

162-element Array{Float64,1}:
  0.0
  0.0
  0.0
  0.0
  0.85065080835204
  0.85065080835204
 -0.85065080835204
 -0.85065080835204
 -0.5257311121191336
  0.5257311121191336
  ⋮
 -0.8626684804161862
 -0.9510565162951536
 -0.9510565162951536
  0.2732665289126717
  0.16245984811645317
 -0.4338885645526948
 -0.16245984811645317
  0.16062203564002314
  0.0

Potential on the surface vs analytics:

scatter(cosθs,ψ)
oplot(cosθ,ψt.(cosθ,ϵ,ζ,r1) )
savefig("potential.svg")

Normal derivatives on the surface vs analytics:

scatter(cosθs,n∇ψ)
oplot(cosθ,∇ψn.(cosθ,ϵ,ζ,r1) )
savefig("nderivatives.svg")