k = 50 # Set wavenumber and fundamental solution for Helmholtz equation
g1 = (x,y) -> -besselj0(k*abs(y-x))/2
g2 = (x,y) -> x == y ? -(log(k/2)+γ)/2/π + im/4 : im/4*hankelh1(0,k*abs(y-x)) - g1(x,y).*logabs(y-x)/π

ui = (x,y) -> exp(im*k*(x-y)/sqrt(2))    # Incident plane wave at 45°

dom = ChebyshevInterval()                     # Set the domain
sp = Space(dom)                      # Canonical space on the domain
⨍ = DefiniteLineIntegral(dom)        # Line integration functional
uiΓ = Fun(t->ui(real(t),imag(t)),sp) # Incident wave on Γ

# Instantiate the fundamental solution
G = GreensFun(g1,CauchyWeight(sp⊗sp,0)) + GreensFun(g2,sp⊗sp)

# Solve for the density
∂u∂n = ⨍[G]\uiΓ

# What is the forward error?
norm(⨍[G]*∂u∂n-uiΓ)

# Represent the scattered field
us = (x,y) -> -logkernel(g1,∂u∂n,complex(x,y))-linesum(g2,∂u∂n,complex(x,y))

ui = (x,y) -> logabs(complex(x,y)-2)     # Single source at (2,0) of strength 2π

N,r = 10,1e-1
cr = exp.(im*2π*(0:N-1)/N)
crl,crr = (1-2im*r)cr,(1+2im*r)cr
dom = ∪(Segment.(crl,crr))            # Set the shielding domain

sp = Space(dom)                      # Canonical space on the domain
⨍ = DefiniteLineIntegral(dom)        # Line integration functional
uiΓ = Fun(t->ui(real(t),imag(t)),sp) # Action of source on shields

# Instantiate the fundamental solution
G = GreensFun((x,y)->0.5,CauchyWeight(sp⊗sp,0))

# The first column augments the system for global unknown constant charge φ0
# The first row ensure constant charge φ0 on all plates
φ0,∂u∂n=[0 ⨍;1 ⨍[G]]\[0.;uiΓ]

# Represent the scattered field
us = (x,y) -> -logkernel(∂u∂n,complex(x,y))/2

G=Fun(z->[-1 -3; -3 -1]/z +
         [ 2  2;  1 -3] +
         [ 2 -1;  1  2]*z,Circle())

C  = Cauchy(-1)
V  = (I+(I-G)*C)\(G-I)

L  = ToeplitzOperator(inv(I+C*V))
U  = ToeplitzOperator(I+V+C*V)

α = 2.5 # ∈ [0, d+2), where d is the number of dimensions
        # α is the strength of the singularity of the algebraic kernel
δ = 0.1 # the horizon of the nonlocal integral operator
L = NonlocalLaplacian(Fourier(), α, δ)