bug fix causing incorrect eclipse calculation

Project.toml

@@ -1,7 +1,7 @@
 name = "EclipsingBinaryStars"
 uuid = "298f00ba-5fcf-566c-9100-f8b551e1d946"
 authors = ["m-wells"]
-version = "0.4.3"
+version = "0.4.4"
 
 [deps]
 DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"

src/binary.jl

@@ -15,16 +15,16 @@ Binary(M1, R1, M2, R2, x; kwargs...) = Binary(Star(M1, R1), Star(M2, R2), x; kwa
 get_star1(b::Binary) = b.pri
 get_star2(b::Binary) = b.sec
 get_orbit(b::Binary) = b.orb
-get_m1(b) = get_m(get_star1(b))
-get_m2(b) = get_m(get_star2(b))
-get_r1(b) = get_r(get_star1(b))
-get_r2(b) = get_r(get_star2(b))
-get_a(b) = get_a(get_orbit(b))
-get_P(b) = get_P(get_orbit(b))
-get_e(b) = get_e(get_orbit(b))
-get_i(b) = get_i(get_orbit(b))
-get_ω(b) = get_ω(get_orbit(b))
-get_Ω(b) = get_Ω(get_orbit(b))
+get_m1(b::Binary) = get_m(get_star1(b))
+get_m2(b::Binary) = get_m(get_star2(b))
+get_r1(b::Binary) = get_r(get_star1(b))
+get_r2(b::Binary) = get_r(get_star2(b))
+get_a(b::Binary) = get_a(get_orbit(b))
+get_P(b::Binary) = get_P(get_orbit(b))
+get_e(b::Binary) = get_e(get_orbit(b))
+get_i(b::Binary) = get_i(get_orbit(b))
+get_ω(b::Binary) = get_ω(get_orbit(b))
+get_Ω(b::Binary) = get_Ω(get_orbit(b))
 
 Base.show(io::IO, b::Binary) = printfields(io, b)
 Base.show(io::IO, ::MIME"text/plain", b::Binary) = print(io, typeof(b), b)

src/ebs_widget.jl

@@ -109,7 +109,6 @@ function ebs_widget(;
  pcolor = "tab:orange",
  scolor = "tab:blue",
 )
- ux, uy, uz = unit_sphere()
 
  r1_rsun = ustrip(u"Rsun", r1)
  r2_rsun = ustrip(u"Rsun", r2)
@@ -217,7 +216,14 @@ function ebs_widget(;
  xmin, xmax = get_lims(xmin, xmax)
  ax_pos.set_xlim(xmin, xmax)
  ax_pos.set_ylim(get_lims(ymin, ymax)...)
- ax_pos.axvspan(xmin, xmax; facecolor="black", alpha=0.2, zorder=0)
+ ax_pos.add_patch(
+ patches.Rectangle((0,0), 1, 1;
+ transform = ax_pos.transAxes,
+ zorder = 0,
+ color = "black",
+ alpha = 0.2
+ )
+ )
 
  xyz1, xyz2 = pos(ν, b)
  x,y,z = ustrip.(u"Rsun", xyz1)

src/eclipse.jl

@@ -13,49 +13,72 @@ end
 Eclipse(ν::Angle) = Eclipse(uconvert(u"°", ν))
 Eclipse(ν::Real) = Eclipse(rad2deg(ν)*u"°")
 
-get_eclipse_ν(eclip::Eclipse) = eclip.ν
-has_eclipse(eclip::Eclipse) = !isnan(get_eclipse_ν(eclip))
+abstract type Solver end
+
+struct Fast <: Solver
+end
+export Fast
 
-# assuming r1, r2, and a are given in same units
-#function get_eclipses(r1::Real, r2::Real, a::Real, e::Real, i::Real, ω::Real; atol=sqrt(eps()))
-function get_eclipses(b; atol=sqrt(eps()))
- a = ustrip(u"Rsun", get_a(b))
- e = get_e(b)
- i = ustrip(u"rad", get_i(b))
- ω = ustrip(u"rad", get_ω(b))
- r1 = ustrip(u"Rsun", get_r1(b))
- r2 = ustrip(u"Rsun", get_r2(b))
+struct Accurate <: Solver
+end
+export Accurate
+
+function get_eclipses(r1, r2, a, e, i, ω, ::Accurate; atol=sqrt(eps()), νpad=π/12)
+ r1 = ustrip(u"Rsun", r1)
+ r2 = ustrip(u"Rsun", r2)
  r1r2² = (r1 + r2)^2
- if (Δ²(0, a, e, i, π/2) ≥ r1r2²) && (Δ²(0, a, e, i, 3π/2) ≥ r1r2²)
+ a = ustrip(u"Rsun", a)
+ i = ustrip(u"rad", i)
+ ω = ustrip(u"rad", ω)
+
+ r1r2² = (r1 + r2)^2
+
+ if Δ²(0, a, e, i, π/2) ≥ r1r2² 
  (return Eclipse(NaN), Eclipse(NaN))
  end
- #Δ²(0, a, e, i, π/2) ≥ (r1 + r2)^2 && (return Eclipse(NaN), Eclipse(NaN))
- f(x) = Δ²(x, a, e, i, ω) - r1r2²
- f(x::AbstractArray) = f(first(x))
- g(x) = dΔ²_dν(x, a, e, i, ω)
- function g!(xout, xin::AbstractArray)
- xout[1] = g(first(xin))
- return nothing
- end
-
- x₁, x₂, x₃ = [-ω], [π - ω], [2π - ω]
- ν₁, ν₂ = [sconj(ω)], [iconj(ω)]
-
- res = optimize(f, g!, x₁, x₂, ν₁)
- ν₁ = mod2pi(first(Optim.minimizer(res)))
- ν₁ = (f(ν₁) < 0) && isapprox(g(ν₁), 0; atol=atol) ? ν₁ : NaN
+
+ f(x) = Δ²(x, a, e, i, ω)
 
- res = optimize(f, g!, x₂, x₃, ν₂)
- ν₂ = mod2pi(first(Optim.minimizer(res)))
- ν₂ = (f(ν₂) < 0) && isapprox(g(ν₂), 0; atol=atol) ? ν₂ : NaN
+ x₂ = sconj(ω)
+ x₁ = x₂ - νpad
+
+ res = optimize(f, x₁, x₂)
+ ν₁ = mod2pi(Optim.minimizer(res))
+ ν₁ = f(ν₁) < r1r2² ? ν₁ : NaN
+
+ x₁ = iconj(ω)
+ x₂ = x₁ + νpad
+ res = optimize(f, x₁, x₂)
+ ν₂ = mod2pi(Optim.minimizer(res))
+ ν₂ = f(ν₂) < r1r2² ? ν₂ : NaN
+
+ return Eclipse(ν₁), Eclipse(ν₂)
+end
+
+function get_eclipses(r1, r2, a, e, i, ω, ::Fast)
+ r1 = ustrip(u"Rsun", r1)
+ r2 = ustrip(u"Rsun", r2)
+ a = ustrip(u"Rsun", a)
+ i = ustrip(u"rad", i)
+ ω = ustrip(u"rad", ω)
+
+ cosi = cos(i)
+ sumradi = r1 + r2
+ ν₁ = sconj(ω)
+ ν₁ = abs(r(ν₁, a, e)*cosi) < sumradi ? ν₁ : NaN
+ ν₂ = iconj(ω)
+ ν₂ = abs(r(ν₂, a, e)*cosi) < sumradi ? ν₂ : NaN
  return Eclipse(ν₁), Eclipse(ν₂)
 end
 
-#get_eclipses(b) = get_eclipses(
-# ustrip.(u"Rsun", (get_r1(b), get_r2(b), get_a(b)))...,
-# get_e(b),
-# ustrip.(u"rad", (get_i(b), get_ω(b)))...
-#)
+# default method for get_eclipses is Accurate
+get_eclipses(r1, r2, a, e, i, ω; kwargs...) = get_eclipses(
+ r1, r2, a, e, i, ω, Accurate(); kwargs...
+)
+
+get_eclipses(b, args...; kwargs...) = get_eclipses(
+ get_r1(b), get_r2(b), get_a(b), get_e(b), get_i(b), get_ω(b), args...; kwargs...
+)
 
 struct EclipsingBinary{T}
  bin::Binary{T}
@@ -68,8 +91,10 @@ function EclipsingBinary(b::Binary{T1}, p::Eclipse{T2}, s::Eclipse{T3}) where {T
  return EclipsingBinary(convert(Binary{T}, b), convert.(Eclipse{T}, (p, s))...)
 end
 
-EclipsingBinary(b::Binary) = EclipsingBinary(b, get_eclipses(b)...)
-EclipsingBinary(p::Star, s::Star, x; kwargs...) = EclipsingBinary(Binary(p, s, x; kwargs...))
+EclipsingBinary(b::Binary, args...; kwargs...) = EclipsingBinary(
+ b, get_eclipses(b, args...; kwargs...)...
+)
+EclipsingBinary(args...; kwargs...) = EclipsingBinary(Binary(args...; kwargs...))
 
 Base.show(io::IO, b::EclipsingBinary) = printfields(io, b)
 Base.show(io::IO, ::MIME"text/plain", b::EclipsingBinary) = print(io, typeof(b), b)
@@ -78,12 +103,17 @@ get_binary(b::EclipsingBinary) = b.bin
 get_star1(b::EclipsingBinary) = get_star1(get_binary(b))
 get_star2(b::EclipsingBinary) = get_star2(get_binary(b))
 get_orbit(b::EclipsingBinary) = get_orbit(get_binary(b))
+
 get_eclipse1(b::EclipsingBinary) = b.pecl
 get_eclipse2(b::EclipsingBinary) = b.secl
 get_eclipses(b::EclipsingBinary) = get_eclipse1(b), get_eclipse2(b)
+
+get_eclipse_ν(eclip::Eclipse) = eclip.ν
 get_eclipse1_ν(b) = get_eclipse_ν(get_eclipse1(b))
 get_eclipse2_ν(b) = get_eclipse_ν(get_eclipse2(b))
-get_eclipses_ν(b) = get_eclipse_ν.(get_eclipses(b))
+get_eclipses_ν(b, args...; kwargs...) = get_eclipse_ν.(get_eclipses(b, args...; kwargs...))
+
+has_eclipse(eclip::Eclipse) = !isnan(get_eclipse_ν(eclip))
 has_eclipse1(b) = has_eclipse(get_eclipse1(b))
 has_eclipse2(b) = has_eclipse(get_eclipse2(b))
 has_eclipse(b) = has_eclipse1(b) || has_eclipse2(b)

src/geometry.jl

@@ -60,6 +60,9 @@ dΔ²_dν(ν, b) = dΔ²_dν(ν, get_a(b), get_e(b), get_i(b), get_ω(b))
 sconj(ω) = 0.5π - ω
 iconj(ω) = 1.5π - ω
 
+sconj(b::Binary) = sconj(get_ω(b))
+iconj(b::Binary) = iconj(get_ω(b))
+
 periastron(a, e) = r(0, a, e)
 periastron(b) = r(0, b)
 

test/runtests.jl

@@ -87,3 +87,54 @@ end
  @test !has_eclipse2(eb)
  @test !has_eclipse(eb)
 end
+
+function _get(xmin, xmax)
+ xmin > xmax && error("xmin > xmax\n\txmin = ", xmin, "\n\txmax = ", xmax)
+ xran = xmax - xmin
+ return rand()*xran + xmin
+end
+
+@testset "accurate versus fast" begin
+ n = 10000
+
+ rmin = 0.1
+ rmax = 10.0
+ amax = 10000.0
+ mmin = 0.1
+ mmax = 10.0
+
+ eclips_accurate = falses(n)
+ eclips_fast = falses(n)
+
+ for j in 1:n
+ m1 = _get(mmin, mmax)u"Msun"
+ m2 = _get(mmin, mmax)u"Msun"
+ r1 = _get(rmin, rmax)u"Rsun"
+ r2 = _get(rmin, rmax)u"Rsun"
+ amin = Inf
+ e = NaN
+ check_e = true
+ while check_e
+ e = _get(0, 1)
+ amin = ustrip(u"Rsun", 1.5*(r1 + r2)/(1-e))
+ check_e = amin > amax
+ end
+ a = _get(amin, amax)u"Rsun"
+ i = _get(0, π)u"rad"
+ ω = _get(0, 2π)u"rad"
+
+ b = Binary(m1, r1, m2, r2, a; e=e, i=i, ω=ω)
+ eb = EclipsingBinary(b, Fast())
+ eclips_fast[j] = has_eclipse(eb)
+
+ eb = EclipsingBinary(b)
+ eclips_accurate[j] = has_eclipse(eb)
+
+ end
+ @test sum(eclips_accurate) ≥ sum(eclips_fast)
+ mask = xor.(eclips_accurate, eclips_fast)
+ inds = findall(mask)
+
+ @test all(eclips_accurate[inds])
+ @test !any(eclips_fast[inds])
+end