Add newtoniansquare routine

Author: dlfivefifty

README.md

@@ -1,2 +1,18 @@
 # MultivariateSingularIntegrals.jl
  Julia package for computing multivariate singular integrals
+
+
+The function `newtoniansquare([x,y], n)` computes a matrix of  Newtonian potentials of Legendre polynomials on the unit square $Ω := [-1,1]^2$. That is it computes
+```math
+\begin{pmatrix}
+L_{00}(𝐱) & \cdots & L_{0,p-1}(𝐱) & L_{0p}(𝐱) \\
+L_{10}(𝐱) & \cdots & L_{1,p-1}(𝐱) \\
+\vdots & \iddots \\
+L_{p0}(𝐱)
+\end{pmatrix}
+```
+where $𝐱 = [x,y]$ is any point in $ℝ²$ (including on or near the unit square $Ω$) and 
+```math
+L_{k,j}(𝐱) := ∬_Ω P_k(s) P_j(t) \log(\| [s,t] - [x,y] \|) \, ds \, dt
+```
+where $P_k$ and $P_j$ are the Legendre polynomials of degree $k$ and $j$, respectively.

src/MultivariateSingularIntegrals.jl

@@ -2,14 +2,17 @@ module MultivariateSingularIntegrals
 using LinearAlgebra, ClassicalOrthogonalPolynomials, SingularIntegrals, FillArrays
 import Base: size
 
-export logkernelsquare, stieltjessquare, logkernelsquare!, stieltjessquare!
+export newtoniansquare, logkernelsquare, stieltjessquare, logkernelsquare!, stieltjessquare!
 
 
 function imlogkernel_vec(n, z)
     T = float(real(typeof(z)))
     transpose(complexlogkernel(Legendre{T}(), -im*float(z)))[1:n] + m_const_vec(n, float(z))
 end
 
+newtoniansquare(z::AbstractVector, n) = real(logkernelsquare(complex(z...), n))
+newtoniansquare(z::Number, n) = real(logkernelsquare(z, n))
+
 include("stieltjessquare.jl")
 include("logkernelsquare.jl")
 

test/runtests.jl

@@ -57,6 +57,10 @@ end
     end
 end
 
+@testset "newtonian" begin
+    @test newtoniansquare(0.1+0.2im, 10) == newtoniansquare([0.1,0.2], 10) == real(logkernelsquare(0.1+0.2im, 10))
+end
+
 # n = 4000
 # z = 2.0
 # @profview logkernelsquare(z, n)