add BO equation

.gitignore

@@ -1,4 +1,4 @@
-examples/.ipynb_checkpoints/Manipulate SIEs-checkpoint.ipynb
+examples/.ipynb_checkpoints/*.ipynb
 src/CUDA/**/*.ptx
 deps/**/*.o
 deps/**/*.dll

examples/Scatteraux.jl

@@ -22,6 +22,12 @@ if !isdefined(:scatteraux_loaded)
     end
     @vectorize_2arg Number g3neumann
 
+    function logheatmap(A::AbstractMatrix,ϵ;kwds...)
+        ltϵ = log10(ϵ)
+        LGABSAt = [isinf(log10(abs(A[j,i]))) ? ltϵ : max(log10(abs(A[j,i])),ltϵ) for i=1:size(A,2),j=1:size(A,1)]
+        heatmap(flipdim(LGABSAt,1);kwds...)
+    end
+
     function makegif(x,y,u,L;plotfunction=Main.Plots.PyPlot.contourf,seconds=1,cmap="seismic",vert=1)
         tm=string(time_ns())
         dr = pwd()*"/"*tm*"mov"

src/Operators/Hilbert.jl

@@ -41,17 +41,21 @@ for Op in (:PseudoHilbert,:Hilbert,:SingularIntegral)
         end
 
         bandinds{s,DD}(::$ConcOp{Hardy{s,DD}})=0,0
-        domainspace{s,DD}(H::$ConcOp{Hardy{s,DD}})=H.space
         rangespace{s,DD}(H::$ConcOp{Hardy{s,DD}})=H.space
 
         bandinds{DD}(::$ConcOp{Laurent{DD}})=0,0
         domainspace{DD}(H::$ConcOp{Laurent{DD}})=H.space
         rangespace{DD}(H::$ConcOp{Laurent{DD}})=H.space
 
         bandinds{DD}(H::$ConcOp{Fourier{DD}})=-H.order,H.order
-        domainspace{DD}(H::$ConcOp{Fourier{DD}})=H.space
         rangespace{DD}(H::$ConcOp{Fourier{DD}})=H.space
 
+        bandinds{DD<:PeriodicInterval}(H::$ConcOp{CosSpace{DD}})=(0,1)
+        bandinds{DD<:PeriodicInterval}(H::$ConcOp{SinSpace{DD}})=(-1,0)
+        rangespace{DD<:PeriodicInterval}(H::$ConcOp{CosSpace{DD}})=SinSpace(domain(H))
+        rangespace{DD<:PeriodicInterval}(H::$ConcOp{SinSpace{DD}})=CosSpace(domain(H))
+
+
         function rangespace{DD}(H::$ConcOp{JacobiWeight{Chebyshev{DD},DD}})
             @assert domainspace(H).α==domainspace(H).β==-0.5
             H.order==0 ? Chebyshev(domain(H)) : Ultraspherical(H.order,domain(H))
@@ -302,6 +306,71 @@ function getindex{DD<:Circle,OT,T}(H::ConcreteSingularIntegral{Fourier{DD},OT,T}
     end
 end
 
+## PeriodicInterval
+
+function getindex{DD<:PeriodicInterval,OT,T}(H::ConcreteHilbert{Taylor{DD},OT,T},k::Integer,j::Integer)
+    d = domain(H)
+    a,b = d.a,d.b
+    if H.order == 1
+        if k>1
+            if k == j
+                T(im)
+            else
+                zero(T)
+            end
+        else
+            zero(T)
+        end
+    else
+        error("Hilbert order $(H.order) not implemented for Taylor")
+    end
+end
+
+function getindex{DD<:PeriodicInterval,OT,T}(H::ConcreteHilbert{Hardy{false,DD},OT,T},k::Integer,j::Integer)
+    d = domain(H)
+    a,b = d.a,d.b
+    if H.order == 1
+        if k>1
+            if k == j
+                T(-im)
+            else
+                zero(T)
+            end
+        else
+            zero(T)
+        end
+    else
+        error("Hilbert order $(H.order) not implemented for Hardy{false}")
+    end
+end
+
+function getindex{DD<:PeriodicInterval,OT,T}(H::ConcreteHilbert{CosSpace{DD},OT,T},k::Integer,j::Integer)
+    d = domain(H)
+    a,b = d.a,d.b
+    if H.order == 1
+        if k == j-1
+            -one(T)
+        else
+            zero(T)
+        end
+    else
+        error("Hilbert order $(H.order) not implemented for CosSpace")
+    end
+end
+
+function getindex{DD<:PeriodicInterval,OT,T}(H::ConcreteHilbert{SinSpace{DD},OT,T},k::Integer,j::Integer)
+    d = domain(H)
+    a,b = d.a,d.b
+    if H.order == 1
+        if k == j+1
+            one(T)
+        else
+            zero(T)
+        end
+    else
+        error("Hilbert order $(H.order) not implemented for CosSpace")
+    end
+end
 
 ## JacobiWeight
 

src/periodicline.jl

@@ -11,7 +11,7 @@ end
 
 function hilbert{L<:PeriodicLine,SS}(S::Space{SS,L},f,z::Number)
     S2=setdomain(S,Circle())
-    hilbert(S2,f,mappoint(domain(f),Circle(),z))-hilbert(S2,f,-1)
+    hilbert(S2,f,mappoint(domain(S),Circle(),z))-hilbert(S2,f,-1)
 end
 
 

src/stieltjesmoment.jl

@@ -28,7 +28,7 @@ stieltjesjacobimoment(α::Real,β::Real,z) = stieltjesjacobimoment(α,β,0,z)
 function logjacobimoment(α::Real,β::Real,n::Int,z)
     x = 2./(1-z)
     if n == 0
-        2normalization(0,α,β)*(log.(z-1)-dualpart(_₂F₁(dual(zero(α)+eps(α+β),one(β)),α+1,α+β+2,x)))
+        2normalization(0,α,β)*(log.(z-1)-dualpart.(_₂F₁(dual(zero(α)+eps(α+β),one(β)),α+1,α+β+2,x)))
         # For testing purposes only, should be equivalent to above within radius of convergence
         #2normalization(0,α,β)*(log(z-1)-(α+1)/(α+β+2)*x.*_₃F₂(α+2,α+β+3,x))
     else

test/runtests.jl

@@ -31,6 +31,17 @@ for s in (true,false)
 end
 
 
+#
+# Hilbert transform on the real line for periodic functions on [a,b).
+# It is involutory, so H^2*f = - f, except it zeros the constant.
+#
+
+H = Hilbert()
+f = Fun(ones(20),Laurent(PeriodicInterval()))
+@test norm(H^2*f+f-1) == 0.0
+f = Fun(ones(20),Fourier(PeriodicInterval()))
+@test norm(H^2*f+f-1) == 0.0
+
 x=Fun()
 @test sum(logabs(x-2.0)) ≈ logabslegendremoment(2.0)
 @test sum(logabs(x-(2.0+im))) ≈ logabslegendremoment(2.0+im)