SimpleOrbit updates

src/orbits/simple.jl

@@ -9,7 +9,7 @@ struct SimpleOrbit{T,V,S} <: AbstractOrbit
 end
 
 """
-    SimpleOrbit(; period, duration, t0=0, b=0)
+    SimpleOrbit(; period, duration, t0=0, b=0.0)
 
 Circular orbit parameterized by the basic observables of a transiting system.
 
@@ -19,50 +19,66 @@ Circular orbit parameterized by the basic observables of a transiting system.
 * `t0` - The midpoint time of the reference transit, similar units as `period`
 * `b` - The impact parameter of the orbit, unitless
 """
-function SimpleOrbit(;period, duration, t0=zero(period), b=0)
+function SimpleOrbit(;period, duration, t0=zero(period), b=0.0)
     half_period = 0.5 * period
     duration > half_period && error("duration cannot be longer than half the period")
-    speed = 2 * sqrt(1 - b^2) / duration
+    speed = 2.0 * sqrt(1.0 - b^2) / duration
     ref_time =  t0 - half_period
     # normalize time types
     period, t0, duration, half_period, ref_time = promote(period, t0, duration, half_period, ref_time)
     SimpleOrbit(period, t0, b, duration, speed, half_period, ref_time)
 end
 
-
 period(orbit::SimpleOrbit) = orbit.period
 duration(orbit::SimpleOrbit) = orbit.duration
 
-relative_time(orbit::SimpleOrbit, t) = (t - orbit.ref_time) % period(orbit) - orbit.half_period
+relative_time(orbit::SimpleOrbit, t) = mod(t - orbit.ref_time, period(orbit)) - orbit.half_period
 
 function relative_position(orbit::SimpleOrbit, t)
     Δt = relative_time(orbit, t)
     x = orbit.speed * Δt
     y = orbit.b
-    z = abs(Δt) < 0.5 * orbit.duration ? one(x) : -one(x)
+    z = abs(Δt) < 0.5 * duration(orbit) ? one(x) : -one(x)
     return SA[x, y, z]
 end
 
+# TODO: if texp, ϵ += 0.5 * texp
+#function in_transit(orbit::SimpleOrbit, t)
+#    Δt = relative_time.(orbit, t)
+#    ϵ = 0.5 * duration(orbit)
+#    return findall(x -> abs(x) < ϵ, Δt)
+#end
+#function in_transit(orbit::SimpleOrbit, t, r)
+#    Δt = relative_time.(orbit, t)
+#    ϵ = √((1.0 + r)^2.0 - orbit.b^2.0) / orbit.speed
+#    return findall(x -> abs(x) < ϵ, Δt)
+#end
+
 function flip(orbit::SimpleOrbit, ror)
     t0 = orbit.t0 + orbit.half_period
     b = orbit.b / ror
     speed = orbit.speed / ror
     ref_time = orbit.t0
-    return SimpleOrbit(orbit.period, t0, b, orbit.duration, speed, orbit.half_period, ref_time)
+    return SimpleOrbit(period(orbit), t0, b, duration(orbit), speed, orbit.half_period, ref_time)
 end
 
 function Base.show(io::IO, orbit::SimpleOrbit)
-    T = orbit.duration
-    P = orbit.period
+    T = duration(orbit)
+    P = period(orbit)
     b = orbit.b
     t0 = orbit.t0
-    print(io, "SimpleOrbit(P=", P, ", T=", T, ", t0=", t0, ", b=", b, ")")
+    print(io, "SimpleOrbit(P=$P, T=$T, t0=$t0, b=$b)")
 end
 
 function Base.show(io::IO, ::MIME"text/plain", orbit::SimpleOrbit)
-    T = orbit.duration
-    P = orbit.period
+    T = duration(orbit)
+    P = period(orbit)
     b = orbit.b
     t0 = orbit.t0
-    print(io, "SimpleOrbit\n period: ", P, "\n duration: ", T, "\n t0: ", t0, "\n b: ", b)
+    print(io,
+        "SimpleOrbit\n period: ", P,
+        "\n duration: ", T,
+        "\n t0: ", t0,
+        "\n b: ", b
+    )
 end

test/orbits/simple.jl

@@ -0,0 +1,47 @@
+using Transits.Orbits: relative_position
+
+@testset "simple orbit" begin
+    period = 3.456
+    t0 = 1.45
+    b = 0.5
+    duration = 0.12
+    t = t0 .+ range(-2.0*period, 2.0*period, length=5000)
+    m0 = @. abs(mod(t - t0 + 0.5*period, period) - 0.5*period) < 0.5 * duration
+
+    orbit = SimpleOrbit(;period, t0, b, duration)
+    pos = @. relative_position(orbit, t)
+    bs = map(v -> sqrt(v[1]^2.0 + v[2]^2.0), pos)
+    zs = map(v->v[3], pos)
+    m = @. (bs ≤ 1.0) & (zs > 0.0)
+    @test orbit.ref_time == -0.278
+    @test m ≈ m0
+
+    #idxs = in_transit(orbit, t)
+    #@test all(bs[idxs] .≤ 1.0)
+    #@test all(zs[idxs] .> 0.0)
+
+    # TODO: Do we want this?
+    #pos  = @. star_position(orbit, t)
+    #for vec in pos
+    #    @test all(vec .≈ 0.0)
+    #end
+end
+
+@testset "simple limbdark" begin
+    period = 3.456
+    t0 = 1.45
+    b = 0.5
+    duration = 0.12
+
+    t = t0 .+ range(-2.0*period, 2.0*period, length=5000)
+    m0 = @. abs(mod(t - t0 + 0.5*period, period) - 0.5*period) < 0.5 * duration
+    orbit = SimpleOrbit(;period, t0, b, duration)
+    ld = PolynomialLimbDark([0.2, 0.3])
+
+    lc1 = @. ld.(orbit, t, 0.01)
+    @test all(lc1[.!m0] .== 1.0)
+
+    ldt = IntegratedLimbDark(ld)
+    lc2 = @. ldt.(orbit, t, 0.01, 0.01)
+    @assert lc1 ≉ lc2
+end

test/runtests.jl

@@ -17,4 +17,5 @@ rng = StableRNG(2752)
     include("show.jl")
     include("distributions.jl")
     include("grads.jl")
+    include("orbits/simple.jl")
 end

test/show.jl

@@ -55,12 +55,12 @@ end
 @testset "SimpleOrbit" begin
     orbit = SimpleOrbit(duration=1, period=3)
 
-    @test sprint(show, orbit) == "SimpleOrbit(P=3.0, T=1.0, t0=0.0, b=0)"
+    @test sprint(show, orbit) == "SimpleOrbit(P=3.0, T=1.0, t0=0.0, b=0.0)"
 
     @test sprint(show, "text/plain", orbit) == """
     SimpleOrbit
      period: 3.0
      duration: 1.0
      t0: 0.0
-     b: 0"""
+     b: 0.0"""
 end