Skip to content

Commit 013597b

Browse files
committed
using works
1 parent 4029231 commit 013597b

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

62 files changed

+10826
-554
lines changed

src/ApproxFunOrthogonalPolynomials.jl

Lines changed: 24 additions & 554 deletions
Large diffs are not rendered by default.

src/Domains/Arc.jl

Lines changed: 109 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,109 @@
1+
export Arc
2+
3+
mobius(a,b,c,d,z) = (a*z+b)/(c*z+d)
4+
mobius(a,b,c,d,z::Number) = isinf(z) ? a/c : (a*z+b)/(c*z+d)
5+
mobiusinv(a,b,c,d,z) = mobius(d,-b,-c,a,z)
6+
7+
8+
mobiusD(a,b,c,d,z) = (a*d-b*c)/(d+c*z)^2
9+
mobiusinvD(a,b,c,d,z) = mobiusD(d,-b,-c,a,z)
10+
"""
11+
Arc(c,r,(θ₁,θ₂))
12+
13+
represents the arc centred at `c` with radius `r` from angle `θ₁` to `θ₂`.
14+
"""
15+
struct Arc{T,V<:Real,TT} <: SegmentDomain{TT}
16+
center::T
17+
radius::V
18+
angles::Tuple{V,V}
19+
Arc{T,V,TT}(c,r,a) where {T,V,TT} = new{T,V,TT}(convert(T,c),convert(V,r),Tuple{V,V}(a))
20+
end
21+
22+
23+
Arc(c::T,r::V,t::Tuple{V1,V2}) where {T<:Number,V<:Real,V1<:Real,V2<:Real} =
24+
Arc{promote_type(T,V,V1,V2),
25+
promote_type(real(T),V,V1,V2),
26+
Complex{promote_type(real(T),V,V1,V2)}}(c,r,t)
27+
Arc(c::Vec{2,T},r::V,t::Tuple{V,V}) where {T<:Number,V<:Real} =
28+
Arc{Vec{2,promote_type(T,V)},
29+
promote_type(real(T),V),
30+
Vec{2,promote_type(real(T),V)}}(c,r,t)
31+
Arc(c::Vec{2,T},r::V,t::Tuple{V1,V2}) where {T<:Number,V<:Real,V1<:Real,V2<:Real} =
32+
Arc{Vec{2,promote_type(T,V,V1,V2)},
33+
promote_type(real(T),V,V1,V2),
34+
Vec{2,promote_type(T,V,V1,V2)}}(c,r,t)
35+
36+
Arc(c::Tuple,r,t) = Arc(Vec(c...),r,t)
37+
Arc(c,r,t0,t1) = Arc(c,r,(t0,t1))
38+
39+
40+
complex(a::Arc{V}) where {V<:Vec} = Arc(complex(a.center...),a.radius,a.angles)
41+
42+
isambiguous(d::Arc) =
43+
isnan(d.center) && isnan(d.radius) && isnan(d.angles[1]) && isnan(d.angles[2])
44+
convert(::Type{Arc{T,V}},::AnyDomain) where {T<:Number,V<:Real} =
45+
Arc{T,V}(NaN,NaN,(NaN,NaN))
46+
convert(::Type{IT},::AnyDomain) where {IT<:Arc} =
47+
Arc(NaN,NaN,(NaN,NaN))
48+
49+
isempty(d::Arc) = false
50+
51+
reverseorientation(a::Arc) = Arc(a.center,a.radius,reverse(a.angles))
52+
53+
arclength(d::Arc) = d.radius*(d.angles[2]-d.angles[1])
54+
55+
56+
function mobiuspars(z0,r,t0,t1)
57+
c=exp(im*t0/2)+exp(im*t1/2)
58+
f=exp(im*t0/2)-exp(im*t1/2)
59+
d=exp(im*t0/2+im*t1/2)*r
60+
-c,c*(z0+d),f,f*(d-z0)
61+
end
62+
63+
mobiuspars(a::Arc) = mobiuspars(a.center,a.radius,a.angles...)
64+
65+
for OP in (:mobius,:mobiusinv,:mobiusD,:mobiusinvD)
66+
@eval $OP(a::Arc,z) = $OP(mobiuspars(a)...,z)
67+
end
68+
69+
70+
tocanonical(a::Arc{T},x) where {T<:Number} = real(mobius(a,x))
71+
tocanonicalD(a::Arc{T},x) where {T<:Number} = mobiusD(a,x)
72+
fromcanonical(a::Arc{T,V,TT},x) where {T<:Number,V<:Real,TT<:Complex} =
73+
mobiusinv(a,x)
74+
fromcanonicalD(a::Arc{T},x) where {T<:Number} = mobiusinvD(a,x)
75+
76+
tocanonical(a::Arc{V},x::Vec) where {V<:Vec} =
77+
tocanonical(complex(a),complex(x...))
78+
fromcanonical(a::Arc{V},x::Number) where {V<:Vec} =
79+
Vec(reim(fromcanonical(complex(a),x))...)
80+
fromcanonicalD(a::Arc{V},x::Number) where {V<:Vec} =
81+
Vec(reim(fromcanonicalD(complex(a),x))...)
82+
83+
## information
84+
85+
issubset(a::Arc,b::Arc) =
86+
a.center==b.center && a.radius==b.radius &&
87+
(b.angles[1]a.angles[1]a.angles[2]b.angles[2] ||
88+
b.angles[1]a.angles[1]a.angles[2]b.angles[2])
89+
90+
91+
# Algebra
92+
93+
*(c::Real,d::Arc) =
94+
Arc(c*d.center,c*d.radius,(sign(c)*d.angles[1],sign(c)*d.angles[2]))
95+
*(d::Arc,c::Real) = *(c,d)
96+
97+
for op in (:+,:-)
98+
@eval begin
99+
$op(c::Number,d::Arc) = Arc($op(c,d.center),d.radius,d.angles)
100+
$op(d::Arc,c::Number) = Arc($op(d.center,c),d.radius,d.angles)
101+
end
102+
end
103+
104+
105+
# allow exp(im*Segment(0,1)) for constructing arc
106+
function exp(d::IntervalOrSegment{<:Complex})
107+
@assert isapprox(real(leftendpoint(d)),0) && isapprox(real(rightendpoint(d)),0)
108+
Arc(0,1,(imag(leftendpoint(d)),imag(rightendpoint(d))))
109+
end

src/Domains/Domains.jl

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
include("Ray.jl")
2+
include("Arc.jl")
3+
include("Line.jl")
4+
include("IntervalCurve.jl")
5+
6+
# sort
7+
8+
isless(d1::IntervalOrSegment{T1},d2::Ray{false,T2}) where {T1<:Real,T2<:Real} = d1 d2.center
9+
isless(d2::Ray{true,T2},d1::IntervalOrSegment{T1}) where {T1<:Real,T2<:Real} = d2.center d1

src/Domains/IntervalCurve.jl

Lines changed: 49 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,49 @@
1+
2+
3+
4+
5+
struct IntervalCurve{S<:Space,T,VT} <: SegmentDomain{T}
6+
curve::Fun{S,T,VT}
7+
end
8+
9+
==(a::IntervalCurve, b::IntervalCurve) = a.curve == b.curve
10+
11+
isempty(::IntervalCurve) = false
12+
13+
points(c::IntervalCurve, n::Integer) = c.curve.(points(domain(c.curve),n))
14+
15+
checkpoints(d::IntervalCurve) = fromcanonical.(Ref(d),checkpoints(domain(d.curve)))
16+
17+
for op in (:(leftendpoint),:(rightendpoint),:(rand))
18+
@eval $op(c::IntervalCurve) = c.curve($op(domain(c.curve)))
19+
end
20+
21+
22+
23+
canonicaldomain(c::IntervalCurve) = domain(c.curve)
24+
25+
fromcanonical(c::IntervalCurve{S,T},x) where {S<:Space,T<:Number} = c.curve(x)
26+
function tocanonical(c::IntervalCurve,x)
27+
rts=roots(c.curve-x)
28+
@assert length(rts)==1
29+
first(rts)
30+
end
31+
32+
33+
fromcanonicalD(c::IntervalCurve,x)=differentiate(c.curve)(x)
34+
35+
function indomain(x,c::IntervalCurve)
36+
rts=roots(c.curve-x)
37+
if length(rts)  1
38+
false
39+
else
40+
in(first(rts),canonicaldomain(c))
41+
end
42+
end
43+
44+
isambiguous(d::IntervalCurve) = ncoefficients(d.curve)==0 && isambiguous(domain(d.curve))
45+
46+
reverseorientation(d::IntervalCurve) = IntervalCurve(reverseorientation(d.curve))
47+
convert(::Type{IntervalCurve{S,T}},::AnyDomain) where {S,T}=Fun(S(AnyDomain()),[NaN])
48+
49+
arclength(d::IntervalCurve) = linesum(ones(d))

src/Domains/Line.jl

Lines changed: 226 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,226 @@
1+
2+
3+
export Line
4+
5+
6+
7+
8+
## Standard interval
9+
10+
11+
12+
# angle is π*a where a is (false==0) and (true==1)
13+
# or ranges from (-1,1]. We use 1 as 1==true.
14+
15+
"""
16+
Line{a}(c)
17+
18+
represents the line at angle `a` in the complex plane, centred at `c`.
19+
"""
20+
struct Line{angle,T<:Number} <: SegmentDomain{T}
21+
center::T
22+
α::Float64
23+
β::Float64
24+
25+
#TODO get this inner constructor working again.
26+
Line{angle,T}(c,α,β) where {angle,T} = new{angle,T}(c,α,β)
27+
Line{angle,T}(c) where {angle,T} = new{angle,T}(c,-1.,-1.)
28+
Line{angle,T}() where {angle,T} = new{angle,T}(zero(T),-1.,-1.)
29+
end
30+
31+
const RealLine{T} = Union{Line{false,T},Line{true,T}}
32+
33+
Line{a}(c,α,β) where {a} = Line{a,typeof(c)}(c,α,β)
34+
Line{a}(c::Number) where {a} = Line{a,typeof(c)}(c)
35+
Line{a}() where {a} = Line{a,Float64}()
36+
37+
angle(d::Line{a}) where {a} = a*π
38+
39+
reverseorientation(d::Line{true}) = Line{false}(d.center,d.β,d.α)
40+
reverseorientation(d::Line{false}) = Line{true}(d.center,d.β,d.α)
41+
reverseorientation(d::Line{a}) where {a} = Line{a-1}(d.center,d.β,d.α)
42+
43+
# ensure the angle is always in (-1,1]
44+
Line(c,a,α,β) = Line{mod(a/π-1,-2)+1,typeof(c)}(c,α,β)
45+
Line(c,a) = Line(c,a,-1.,-1.)
46+
# true is negative orientation, false is positive orientation
47+
# this is because false==0 and we take angle=0
48+
Line(b::Bool) = Line{b}()
49+
Line() = Line(false)
50+
51+
52+
isambiguous(d::Line)=isnan(d.center)
53+
convert(::Type{Domain{T}},d::Line{a}) where {a,T<:Number} = Line{a,T}(d.center, d.α, d.β)
54+
convert(::Type{Line{a,T}},::AnyDomain) where {a,T<:Number} = Line{a,T}(NaN)
55+
convert(::Type{IT},::AnyDomain) where {IT<:Line}=Line(NaN,NaN)
56+
57+
## Map interval
58+
59+
60+
##TODO non-1 alpha,beta
61+
62+
isempty(::Line) = false
63+
64+
function line_tocanonical(α,β,x)
65+
@assert α==β==-1. || α==β==-.5
66+
67+
if α==β==-1.
68+
2x/(1+sqrt(1+4x^2))
69+
elseif α==β==-.5
70+
x/sqrt(1 + x^2)
71+
end
72+
end
73+
74+
function line_tocanonicalD(α,β,x)
75+
@assert α==β==-1. || α==β==-.5
76+
77+
if α==β==-1.
78+
2/(1+4x^2+sqrt(1+4x^2))
79+
elseif α==β==-0.5
80+
(1 + x^2)^(-3/2)
81+
end
82+
end
83+
function line_fromcanonical(α,β,x)
84+
#TODO: why is this consistent?
85+
if α==β==-1.
86+
x/(1-x^2)
87+
else
88+
x*(1 + x)^α*(1 - x)^β
89+
end
90+
end
91+
function line_fromcanonicalD(α,β,x)
92+
if α==β==-1.
93+
(1+x^2)/(1-x^2)^2
94+
else
95+
(1 --α)x -+α+1)x^2)*(1+x)^-1)*(1-x)^-1)
96+
end
97+
end
98+
99+
function line_invfromcanonicalD(α,β,x)
100+
if α==β==-1.
101+
(1-x^2)^2/(1+x^2)
102+
else
103+
1/(1 --α)x -+α+1)x^2)*(1+x)^(1-α)*(1-x)^(1-β)
104+
end
105+
end
106+
107+
108+
tocanonical(d::Line,x) = line_tocanonical(d.α,d.β,cis(-angle(d)).*(x-d.center))
109+
tocanonical(d::Line{false},x) = line_tocanonical(d.α,d.β,x-d.center)
110+
tocanonical(d::Line{true},x) = line_tocanonical(d.α,d.β,d.center-x)
111+
112+
tocanonicalD(d::Line,x) = cis(-angle(d)).*line_tocanonicalD(d.α,d.β,cis(-angle(d)).*(x-d.center))
113+
tocanonicalD(d::Line{false},x) = line_tocanonicalD(d.α,d.β,x-d.center)
114+
tocanonicalD(d::Line{true},x) = -line_tocanonicalD(d.α,d.β,d.center-x)
115+
116+
fromcanonical(d::Line,x) = cis(angle(d))*line_fromcanonical(d.α,d.β,x)+d.center
117+
fromcanonical(d::Line{false},x) = line_fromcanonical(d.α,d.β,x)+d.center
118+
fromcanonical(d::Line{true},x) = -line_fromcanonical(d.α,d.β,x)+d.center
119+
120+
fromcanonicalD(d::Line,x) = cis(angle(d))*line_fromcanonicalD(d.α,d.β,x)
121+
fromcanonicalD(d::Line{false},x) = line_fromcanonicalD(d.α,d.β,x)
122+
fromcanonicalD(d::Line{true},x) = -line_fromcanonicalD(d.α,d.β,x)
123+
124+
invfromcanonicalD(d::Line,x) = cis(-angle(d))*line_invfromcanonicalD(d.α,d.β,x)
125+
invfromcanonicalD(d::Line{false},x) = line_invfromcanonicalD(d.α,d.β,x)
126+
invfromcanonicalD(d::Line{true},x) = -line_invfromcanonicalD(d.α,d.β,x)
127+
128+
129+
130+
131+
132+
133+
==(d::Line{a},m::Line{a}) where {a} = d.center == m.center && d.β == m.β &&d.α == m.α
134+
135+
136+
137+
# algebra
138+
*(c::Real,d::Line{false}) = Line{sign(c)>0 ? false : true}(isapprox(d.center,0) ? d.center : c*d.center,d.α,d.β)
139+
*(c::Real,d::Line{true}) = Line{sign(c)>0 ? true : false}(isapprox(d.center,0) ? d.center : c*d.center,d.α,d.β)
140+
*(c::Number,d::Line) = Line(isapprox(d.center,0) ? d.center : c*d.center,angle(d)+angle(c),d.α,d.β)
141+
*(d::Line,c::Number) = c*d
142+
for OP in (:+,:-)
143+
@eval begin
144+
$OP(c::Number,d::Line{a}) where {a} = Line{a}($OP(c,d.center),d.α,d.β)
145+
$OP(d::Line{a},c::Number) where {a} = Line{a}($OP(d.center,c),d.α,d.β)
146+
end
147+
end
148+
149+
150+
151+
152+
153+
154+
# algebra
155+
arclength(d::Line) = Inf
156+
leftendpoint(d::Line) = -Inf
157+
rightendpoint(d::Line) = Inf
158+
complexlength(d::Line) =Inf
159+
160+
## vectorized
161+
162+
163+
function convert(::Type{Line},d::ClosedInterval)
164+
a,b=d.left,d.right
165+
@assert abs(a) == abs(b) == Inf
166+
167+
if isa(a,Real) && isa(b,Real)
168+
if a==Inf
169+
@assert b==-Inf
170+
Line(true)
171+
else
172+
@assert a==-Inf&&b==Inf
173+
Line(false)
174+
end
175+
elseif abs(real(a)) < Inf
176+
@assert real(a)==real(b)
177+
@assert sign(imag(a))==-sign(imag(b))
178+
179+
Line(real(b),angle(b))
180+
elseif isnan(real(a)) && isnan(real(b)) # hack for -im*Inf
181+
Line([imag(a)*im,imag(b)*im])
182+
elseif isnan(real(a)) # hack for -im*Inf
183+
Line([real(b)+imag(a)*im,b])
184+
elseif isnan(real(b)) # hack for -im*Inf
185+
Line([a,real(a)+imag(b)*im])
186+
elseif abs(imag(a)) < Inf
187+
@assert imag(a)==imag(b)
188+
@assert sign(real(a))==-sign(real(b))
189+
190+
Line(imag(b),angle(b))
191+
else
192+
@assert angle(b) == -angle(a)
193+
194+
Line(0.,angle(b))
195+
end
196+
end
197+
198+
199+
200+
issubset(a::IntervalOrSegment, b::Line) = leftendpoint(a)b && rightendpoint(a)b
201+
issubset(a::Ray{angle}, b::Line{angle}) where angle = leftendpoint(a) b
202+
issubset(a::Ray{true}, b::Line{false}) = true
203+
issubset(a::Ray{false}, b::Line{true}) = true
204+
205+
function intersect(a::Union{Interval,Segment,Ray},b::Line)
206+
@assert a b
207+
a
208+
end
209+
210+
function union(a::Union{Interval,Segment,Ray},b::Line)
211+
@assert a b
212+
b
213+
end
214+
215+
intersect(b::Line,a::Union{Interval,Segment,Ray}) = intersect(a,b)
216+
union(b::Line,a::Union{Interval,Segment,Ray}) = union(a,b)
217+
218+
219+
function setdiff(b::Line,a::Segment)
220+
@assert a b
221+
if leftendpoint(a)>rightendpoint(a)
222+
b\reverseorientation(a)
223+
else
224+
Ray([leftendpoint(b),leftendpoint(a)]) Ray([rightendpoint(a),rightendpoint(b)])
225+
end
226+
end

0 commit comments

Comments
 (0)