Merge pull request #38 from JuliaAstro/ml/orbitsjl

.JuliaFormatter.toml

@@ -0,0 +1 @@
+style = "blue"

.gitignore

@@ -5,4 +5,5 @@ Manifest.toml
 /docs/build/
 
 dev
-test/output
+test/output
+test.jl

Project.toml

@@ -12,9 +12,9 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 KeywordCalls = "4d827475-d3e4-43d6-abe3-9688362ede9f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Orbits = "88f3d70c-3a4c-11ed-1389-4902f2e49de8"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -30,11 +30,11 @@ ConcreteStructs = "0.2"
 Distributions = "0.22, 0.23, 0.24, 0.25"
 FastGaussQuadrature = "0.4"
 KeywordCalls = "0.2"
-Reexport = "0.2, 1"
+Orbits = "0.1"
 Rotations = "1"
 SpecialFunctions = "0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 1, 2"
 StaticArrays = "0.12, 1"
-StatsFuns = "0.9"
+StatsFuns = "0.9, 1"
 Unitful = "1.9"
 UnitfulAstro = "1"
 julia = "1.5"

bench/speed.jl

@@ -7,13 +7,12 @@ using Statistics
 using StatsPlots, ColorSchemes
 using ProgressLogging
 
-
 N_pts = [100, 316, 1000, 3160, 10000, 31600, 100000]
 N_us = [1 2 3 5 8 13 21 34 55 90 144]
 
 function benchmark(N_u, N_pts; r=0.1)
     u = ones(N_u) ./ N_u
-    b = [sqrt(((i - N_pts/2) * 2/N_pts * (1 + 2 * r))^2) for i in 1:N_pts]
+    b = [sqrt(((i - N_pts / 2) * 2 / N_pts * (1 + 2 * r))^2) for i in 1:N_pts]
 
     ld = PolynomialLimbDark(u)
     output = similar(b)
@@ -35,7 +34,7 @@ labels = reduce(hcat, "uₙ: $N" for N in N_us)
 
 plot(
     N_pts,
-    vals,
+    vals;
     scale=:log10,
     ls=:solid,
     m=:o,
@@ -44,13 +43,12 @@ plot(
     lab=labels,
     xlabel="number of points",
     ylabel="timing [s]",
-    leg=:topleft
+    leg=:topleft,
 )
 
-
 function benchmark_limbdark(N_u, N_pts; r=0.1)
     u = ones(N_u) ./ N_u
-    b = [sqrt(((i - N_pts/2) * 2/N_pts * (1 + 2 * r))^2) for i in 1:N_pts]
+    b = [sqrt(((i - N_pts / 2) * 2 / N_pts * (1 + 2 * r))^2) for i in 1:N_pts]
 
     trans = Limbdark.transit_init(r, b[1], u, true)
     output = similar(b)
@@ -71,21 +69,21 @@ end
 
 plot(
     N_pts,
-    vals_limbdark ./ vals,
+    vals_limbdark ./ vals;
     xscale=:log10,
     ls=:solid,
     m=:o,
     msw=0,
     palette=palette(:plasma, length(N_us) + 2),
     lab=labels,
     xlabel="number of points",
-    ylabel="relative speedup (times faster) - Limbdark.jl"
+    ylabel="relative speedup (times faster) - Limbdark.jl",
 )
 
-tmed = median(vals ./ vals[:, 1], dims=1)
+tmed = median(vals ./ vals[:, 1]; dims=1)
 plot(
-    N_us |> vec,
-    tmed |> vec,
+    vec(N_us),
+    vec(tmed);
     scale=:log10,
     ls=:solid,
     m=:o,
@@ -94,17 +92,16 @@ plot(
     xlabel="number of limb-darkening coefficients",
     ylabel="relative timing",
     ylims=(0.8, 10),
-    leg=:topleft
+    leg=:topleft,
 )
 plot!(
-    N_us |> vec, 
-    median(vals_limbdark ./ vals_limbdark[:, 1], dims=1) |> vec,
+    vec(N_us),
+    vec(median(vals_limbdark ./ vals_limbdark[:, 1]; dims=1));
     lab="Limbdark.jl",
     c=2,
     ls=:solid,
     m=:o,
-    msw=0
+    msw=0,
 )
-plot!(N_us |> vec, n -> n^(0.2), ls=:dash, c=:gray, lab=L"N^{0.2}")
-plot!(N_us |> vec, n -> tmed[end] * n / N_us[1, end], ls=:dash, c=:gray, lab=L"N^1")
-
+plot!(vec(N_us), n -> n^(0.2); ls=:dash, c=:gray, lab=L"N^{0.2}")
+plot!(vec(N_us), n -> tmed[end] * n / N_us[1, end]; ls=:dash, c=:gray, lab=L"N^1")

docs/Project.toml

@@ -1,10 +1,11 @@
 [deps]
-ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2"
+ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+Orbits = "88f3d70c-3a4c-11ed-1389-4902f2e49de8"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Transits = "2e59a628-7bac-4d38-8059-3a73ba0928ab"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
 Unitful = "1"
-Documenter = "0.26, 0.27"

docs/make.jl

@@ -1,14 +1,17 @@
-using Transits
 using Documenter
+using Orbits
+using Transits
 
 setup = quote
     using Transits
 end
 
-DocMeta.setdocmeta!(Transits, :DocTestSetup, setup; recursive = true)
+DocMeta.setdocmeta!(Transits, :DocTestSetup, setup; recursive=true)
+
 include("pages.jl")
+
 makedocs(;
-    modules=[Transits],
+    modules=[Transits, Orbits],
     authors="Miles Lucas <[email protected]> and contributors",
     repo="https://github.com/juliaastro/Transits.jl/blob/{commit}{path}#L{line}",
     sitename="Transits.jl",
@@ -17,11 +20,7 @@ makedocs(;
         canonical="https://juliaastro.github.io/Transits.jl",
         assets=String[],
     ),
-    pages = pages
+    pages=pages,
 )
 
-deploydocs(;
-    repo="github.com/JuliaAstro/Transits.jl",
-    push_preview=true,
-    devbranch="main"
-)
+deploydocs(; repo="github.com/JuliaAstro/Transits.jl", push_preview=true, devbranch="main")

docs/pages.jl

@@ -1,8 +1,8 @@
 
-pages=[
+pages = [
     "Home" => "index.md",
     "Introduction" => "introduction.md",
     "Getting Started" => "gettingstarted.md",
     "Benchmarks" => "bench.md",
-    "API/Reference" => "api.md"
-]
+    "API/Reference" => "api.md",
+]

docs/src/api.md

@@ -41,7 +41,7 @@ grad([0.1, 0.1, 0.4, 0.26])
 To help demonstrate the logic behind these chain rules, here we derive a simple gradient function manually.
 
 ```jldoctest grads
-using ChainRules
+using ChainRulesCore
 
 u_n = [0.4, 0.26]
 μ = 0.1
@@ -70,14 +70,6 @@ gradr([0.1, 0.1, 0.4, 0.26])
 
 For the most granular support for gradients and jacobians, peer into the depths of `polynomial/poly-grad.jl` and `polynomial/quad-grad.jl`. These functions are not part of the public API and are not guaranteed any stability according to [semantic versioning](https://semver.org/).
 
-### Orbits
-
-```@docs
-SimpleOrbit
-KeplerianOrbit
-Orbits.relative_position
-```
-
 ### Distributions
 
 ```@docs

docs/src/gettingstarted.md

@@ -4,6 +4,7 @@
 ## Usage
 
 ```julia
+using Orbits
 using Transits
 
 orbit = SimpleOrbit(period=3, duration=1)

src/Transits.jl

@@ -1,22 +1,19 @@
 module Transits
 
 using ChainRulesCore
-using Reexport
+using Orbits
+using Orbits: AbstractOrbit
 using StaticArrays
+using Unitful
 
 export AbstractLimbDark,
-       PolynomialLimbDark,
-       QuadLimbDark,
-       IntegratedLimbDark,
-       SecondaryLimbDark,
-       compute,
-       Orbits,
-       # distributions
-       Kipping13
-
-include("orbits/Orbits.jl")
-@reexport using .Orbits
-using .Orbits: AbstractOrbit
+    PolynomialLimbDark,
+    QuadLimbDark,
+    IntegratedLimbDark,
+    SecondaryLimbDark,
+    compute,
+    # distributions
+    Kipping13
 
 """
     AbstractLimbDark
@@ -59,6 +56,8 @@ Compute the relative flux by calculating the impact parameter at time `t` from t
 
 # Examples
 ```jldoctest orb
+julia> using Orbits
+
 julia> ld = PolynomialLimbDark([0.4, 0.26]);
 
 julia> orbit = SimpleOrbit(period=3, duration=1);
@@ -82,12 +81,14 @@ julia> ld(orbit, 0u"d", 0.1)
 """
 function compute(ld::AbstractLimbDark, orbit::AbstractOrbit, t, r)
     coords = Orbits.relative_position(orbit, t)
+    # strip units, if necessary (e.g., KeplerianOrbit)
+    x, y, z = ustrip.(coords)
     # make sure los is in front of star
-    if coords[3] > 0
-        μ = sqrt(coords[1]^2 + coords[2]^2)
+    if z > 0
+        μ = sqrt(x^2 + y^2)
         return compute(ld, μ, r)
     else
-        return one(eltype(coords))
+        return one(Base.promote_typeof(x, y, z))
     end
 end
 

src/distributions.jl

@@ -13,17 +13,17 @@ A non-informative prior for two-parameter limb-darkening coefficients using *tri
 # Examples
 
 ```jldoctest
-julia> using Random; rng = Random.seed!(10);
+julia> using StableRNGs; rng = StableRNG(10);
 
 julia> rand(rng, Kipping13())
 2-element Vector{Float64}:
- 0.24716310305467298
- 0.08836997249298882
+  0.3361047299132651
+ -0.025681638815114587
 
 julia> rand(rng, Kipping13(), 5)
 2×5 Matrix{Float64}:
- 0.0664907  0.124817   1.00732    0.806902  0.74165
- 0.520411   0.222718  -0.389412  -0.314755  0.0768429
+ 0.0621057   0.992689   1.77965    0.784055  0.186386
+ 0.0659477  -0.236613  -0.795884  -0.187791  0.592194
 ```
 
 # References
@@ -36,7 +36,7 @@ struct Kipping13 <: MultivariateDistribution{Continuous} end
 
 Base.length(::Kipping13) = 2
 
-function _rand!(rng::AbstractRNG, ::Kipping13, x::AbstractVector{T}) where T
+function _rand!(rng::AbstractRNG, ::Kipping13, x::AbstractVector{T}) where {T}
     q1, q2 = rand(rng, T, 2)
     sqrtq1 = sqrt(q1)
     twoq2 = 2 * q2

src/integrated.jl

@@ -1,7 +1,7 @@
 
 using FastGaussQuadrature
 
-struct IntegratedLimbDark{LD <: AbstractLimbDark,WT,NT} <: AbstractLimbDark
+struct IntegratedLimbDark{LD<:AbstractLimbDark,WT,NT} <: AbstractLimbDark
     driver::LD
     order::Int
     nodes::NT
@@ -58,8 +58,9 @@ function IntegratedLimbDark(u::AbstractVector; kwargs...)
     return IntegratedLimbDark(driver; kwargs...)
 end
 
-compute(ld::IntegratedLimbDark, orbit::AbstractOrbit, t, r; texp=nothing) =
-    compute(ld, orbit, t, r, texp)
+function compute(ld::IntegratedLimbDark, orbit::AbstractOrbit, t, r; texp=nothing)
+    return compute(ld, orbit, t, r, texp)
+end
 
 function compute(ld::IntegratedLimbDark, orbit::AbstractOrbit, t, r, texp)
     # perform change of interval
@@ -72,17 +73,18 @@ function compute(ld::IntegratedLimbDark, orbit::AbstractOrbit, t, r, texp)
     return 0.5 * flux
 end
 
-compute(ld::IntegratedLimbDark, orbit::AbstractOrbit, t, r, ::Nothing) =
-    compute(ld.driver, orbit, t, r)
+function compute(ld::IntegratedLimbDark, orbit::AbstractOrbit, t, r, ::Nothing)
+    return compute(ld.driver, orbit, t, r)
+end
 
-function Base.show(io::IO, ld::IntegratedLimbDark{L1}) where L1
+function Base.show(io::IO, ld::IntegratedLimbDark{L1}) where {L1}
     N = length(ld.nodes)
     np = L1.name.name
-    print(io, "IntegratedLimbDark($np, $N)")
+    return print(io, "IntegratedLimbDark($np, $N)")
 end
 
 function Base.show(io::IO, ::MIME"text/plain", ld::IntegratedLimbDark)
     p = ld.driver
     N = length(ld.nodes)
-    print(io, "IntegratedLimbDark\n driver: $p\n N: $N")
-end
+    return print(io, "IntegratedLimbDark\n driver: $p\n N: $N")
+end