From 07e138f725899369c4e8094e44ea97740daa086b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jorge=20P=C3=A9rez?= Date: Mon, 17 Jul 2023 21:46:33 +0200 Subject: [PATCH] Update Apophis script and some cleanup (#29) * HTTP.jl strikes back * Cleanup * Cleanup * Update src/observations/topocentric.jl * Revert changes in obs_pos_ECEF * Try larger timeouts * Update apophis.jl script * Add pha/Project.toml for script runs * Minor fix * More fixes in apophis.jl * Even more fixes * Fix typo * Fix in bwdfwdeph * Fix type signature in bwdfwdeph * Fix script and bwdfwdeph * Switch back to Downloads: low speed limit 300 seconds, 4 retries * Update apophis.jl * Update apophis.jl * Update apophis.jl * Revert "Switch back to Downloads: low speed limit 300 seconds, 4 retries" This reverts commit fa0796425bb793307774d41c188df8c19fe5bb6e. * Try setting connect_timeout and readtimeout to 180 in download_scratch * Update pha/Project.toml * Update README.md * Cleanup * Update docstrings * Update README.md * Update apophis.jl * Update Apophis notebook * Minor fixes * Update Project.toml * Update apophis.jl * Update apophis.jl --- Project.toml | 2 +- README.md | 24 +- notebooks/Apophis.ipynb | 1709 ++++--------------------- pha/Project.toml | 19 + pha/apophis.jl | 173 ++- src/NEOs.jl | 2 +- src/observations/catalogue_mpc.jl | 21 +- src/observations/jpl_eph.jl | 10 +- src/observations/topocentric.jl | 6 + src/orbit_determination/osculating.jl | 6 +- src/propagation/propagation.jl | 92 -- test/osculating.jl | 2 +- 12 files changed, 469 insertions(+), 1597 deletions(-) create mode 100644 pha/Project.toml diff --git a/Project.toml b/Project.toml index 601094ec..6ec448d8 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "NEOs" uuid = "b41c07a2-2abb-11e9-070a-c3c1b239e7df" authors = ["Jorge A. Pérez Hernández", "Luis Benet", "Luis Eduardo Ramírez Montoya"] -version = "0.7.0" +version = "0.7.1" [deps] Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" diff --git a/README.md b/README.md index ed6671b7..7b2a8fc4 100644 --- a/README.md +++ b/README.md @@ -5,8 +5,8 @@ [![Coverage Status](https://coveralls.io/repos/github/PerezHz/NEOs.jl/badge.svg?branch=main)](https://coveralls.io/github/PerezHz/NEOs.jl?branch=main) [![codecov](https://codecov.io/gh/PerezHz/NEOs.jl/branch/main/graph/badge.svg?token=F1IY79YP3J)](https://codecov.io/gh/PerezHz/NEOs.jl) -`NEOs.jl` is a Near-Earth Object orbital propagator and -fitter in Julia. `NEOs.jl` exploits jet transport techniques via +`NEOs.jl` is a Julia package for high-accuracy orbit determination and propagation of +Near-Earth Objects. `NEOs.jl` exploits jet transport techniques via [TaylorIntegration.jl](https://github.com/PerezHz/TaylorIntegration.jl). ## Authors @@ -27,8 +27,20 @@ The current version of this package may be installed in Julia pkg manager via: ## Usage -The `apophis.jl` file in the `pha` directory contains an example script. This -script may be called as: +The `pha` directory contains the `apophis.jl` script which performs an +orbit determination for asteroid (99942) Apophis from optical and radar astrometry. In order +to run this script, the environment corresponding to the `Project.toml` contained in the +`pha` directory has to be active and instantiated. This can be done, for example, by running +the following command from the `pha` directory: + +``julia -e `import Pkg; Pkg.activate(); Pkg.instantiate()` # run this from the `pha` directory `` + +Once the `pha` environment is active, this script may be called from the `pha` directory +with the default settings as: + +`julia --project apophis.jl` + +The `--help` option can be passed to see a list of the customizable settings `julia --project apophis.jl --help` @@ -43,6 +55,8 @@ resources provided by LANCAD-UNAM-DGTIC-284. ## References +- Pérez-Hernández, J.A., Benet, L. Non-zero Yarkovsky acceleration for near-Earth asteroid + (99942) Apophis. Commun Earth Environ 3, 10 (2022). https://doi.org/10.1038/s43247-021-00337-x - Pérez-Hernández, Jorge A., & Benet, Luis. (2023). [PerezHz/TaylorIntegration.jl](https://github.com/PerezHzTaylorIntegration.jl): - v0.13.0 (Version v0.13.0). Zenodo. http://doi.org/10.5281/zenodo.7953772 + v0.14.2 (Version v0.14.2). Zenodo. https://doi.org/10.5281/zenodo.8104080 diff --git a/notebooks/Apophis.ipynb b/notebooks/Apophis.ipynb index b8f68135..59f9c9a7 100644 --- a/notebooks/Apophis.ipynb +++ b/notebooks/Apophis.ipynb @@ -9,7 +9,8 @@ }, "outputs": [], "source": [ - "]st" + "import Pkg\n", + "Pkg.activate(\"../pha\")" ] }, { @@ -24,829 +25,22 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "using DataFrames, Query, TaylorSeries, PlanetaryEphemeris" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "using JLD" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "using NEOs" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "using Statistics" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "using LinearAlgebra" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "using Dates" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "using Gaston" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "using StatsBase" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "using DelimitedFiles" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "NEOs.loadjpleph()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "## 16 ast + radial nongrav (A1 / SRP)\n", - "#dq = set_variables(\"δx\", order=5, numvars=8)\n", - "#radec_bwd = JLD.load(\"radec_99942_2004_2020.dat.jld\", \"radec_table\")\n", - "#radec_fwd = JLD.load(\"radec_99942_2020_2021.dat.jld\", \"radec_table\")\n", - "#deldop_bwd = JLD.load(\"deldop_99942_RADAR_2005_2013.dat.jld\", \"deldop_table\")\n", - "#deldop_fwd_ = JLD.load(\"deldop_99942_RADAR_2021.dat.jld\", \"deldop_table\")\n", - "\n", - "# 32 ast (NO srp)\n", - "dq = set_variables(\"δx\", order=5, numvars=7)\n", - "radec_bwd = JLD.load(\"radec_99942_2004_2020_32AST.dat.jld\", \"radec_table\")\n", - "radec_fwd = JLD.load(\"radec_99942_2020_2021_32AST.dat.jld\", \"radec_table\")\n", - "deldop_bwd = JLD.load(\"deldop_99942_RADAR_2005_2013_32AST.dat.jld\", \"deldop_table\")\n", - "deldop_fwd_ = JLD.load(\"deldop_99942_RADAR_2021_32AST.dat.jld\", \"deldop_table\")\n", - "\n", - "## 16 ast (NO srp)\n", - "#dq = set_variables(\"δx\", order=5, numvars=7)\n", - "#radec_bwd = JLD.load(\"radec_99942_2004_2020_LAST_16ast.dat.jld\", \"radec_table\")\n", - "#radec_fwd = JLD.load(\"radec_99942_2020_2021_LAST_16ast.dat.jld\", \"radec_table\")\n", - "#deldop_bwd = JLD.load(\"deldop_99942_RADAR_2005_2013_LAST_16ast.dat.jld\", \"deldop_table\")\n", - "#deldop_fwd_ = JLD.load(\"deldop_99942_RADAR_2021_LAST_16ast.dat.jld\", \"deldop_table\")\n", - ";" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "radec_all_ = vcat(radec_bwd, radec_fwd)\n", - "radec_all = radec_all_ |> @filter( \n", - " (Date(_.dt_utc_obs) != Date(2021, 1, 28)) #&& (Date(_.dt_utc_obs) <= Date(2021, 3, 8))\n", - " ) |> DataFrame\n", - ";" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "#deldop_fwd = deldop_fwd_[1:end-1,:] # in paper results, only 3 out of 4 radar measurements are used\n", - "deldop_fwd = deldop_fwd_\n", - "\n", - "deldop_all = vcat(deldop_bwd, deldop_fwd)\n", - ";" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "### first and last dates (optical obs)\n", - "radec_all.dt_utc_obs[1], radec_all.dt_utc_obs[end]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "nrow(radec_all), nrow(deldop_all)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "del_bwd = deldop_bwd |> @filter(_.delay_index) |> DataFrame\n", - "dop_bwd = deldop_bwd |> @filter(_.doppler_index) |> DataFrame\n", - "\n", - "del_fwd = deldop_fwd |> @filter(_.delay_index) |> DataFrame\n", - "dop_fwd = deldop_fwd |> @filter(_.doppler_index) |> DataFrame\n", - "\n", - "del_all = deldop_all |> @filter(_.delay_index) |> DataFrame\n", - "dop_all = deldop_all |> @filter(_.doppler_index) |> DataFrame\n", - ";" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "res_τ = del_all |> @map( _.τ_obs - _.τ_comp ) |> collect\n", - "res_ν = dop_all |> @map( _.ν_obs - _.ν_comp ) |> collect\n", - "\n", - "σ_del_bwd = del_bwd.σ_τ\n", - "σ_dop_bwd = dop_bwd.σ_ν\n", - "\n", - "σ_del_fwd = del_fwd.σ_τ\n", - "σ_dop_fwd = dop_fwd.σ_ν\n", - "\n", - "σ_del_all = del_all.σ_τ\n", - "σ_dop_all = dop_all.σ_ν\n", - "\n", - "w_τ = del_all |> @map( 1/(_.σ_τ)^2 ) |> collect\n", - "w_ν = dop_all |> @map( 1/(_.σ_ν)^2 ) |> collect\n", - "\n", - "\n", - "# only Vokrouhlicky et al. (2015) data\n", - "res_τ_v15 = del_bwd |> @map( _.τ_obs - _.τ_comp ) |> collect\n", - "res_ν_v15 = dop_bwd |> @map( _.ν_obs - _.ν_comp ) |> collect\n", - "\n", - "w_τ_v15 = del_bwd |> @map( 1/(_.σ_τ)^2 ) |> collect\n", - "w_ν_v15 = dop_bwd |> @map( 1/(_.σ_ν)^2 ) |> collect\n", - ";" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# check if all RA/Dec observations are CCD:\n", - "\n", - "all(radec_all |> @map(begin x = _.j2000; x == \"C\" || x == \"X\" end))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# there are 494 RA/Dec obs without catalog info\n", - "nrow( radec_all |> @filter(_.catalog == \" \") |> DataFrame )" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "### WIP: outlier rejection\n", - "\n", - "### function chisqobs(row, Γ_β, x, outlier::Bool=false)\n", - "### A_α_TN = TaylorSeries.gradient(row.α_obs-row.α_corr-row.α_comp)\n", - "### A_δ_TN = TaylorSeries.gradient(row.δ_obs-row.δ_corr-row.δ_comp)\n", - "### A_i = hcat(A_α_TN(x), A_δ_TN(x))\n", - "### outlier_sign = outlier*2-1\n", - "### (A_i')*Γ_β*A_i\n", - "### end\n", - "\n", - "### chisqobs(res0, Γ_OR7, x_OR7)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "#radec_04XX = radec_0415\n", - "#radec_04XX = radec_0420\n", - "#radec_04XX = radec_0421\n", - "#radec_04XX = filter(x->( DateTime(x.yr, x.month, x.day) != DateTime(2021, 1, 28) ), radec_0421)\n", - "#radec_04XX = radec_0421\n", - "\n", - "#radec_04XX = filter(x->x.dt_utc_obs((x.α_obs-x.α_corr-x.α_comp())^2 + (x.δ_obs-x.δ_corr-x.δ_comp())^2)<χ2_rej, radec_04XX)\n", - "#radec_04XX = filter(x->((x.α_obs-x.α_corr-x.α_comp(x_OR7))^2 + (x.δ_obs-x.δ_corr-x.δ_comp(x_OR7))^2)/(x.σ^2) ≤ χ2_rej, radec_04XX)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "#χ2_rej = 4.0 #4.0 #4.25 #5.025 # Brozovic et al (2018) use Chesley et al (2010) outlier rejection, in turn from Carpino et al (2003)\n", - "radec_04XX = radec_all\n", - "#radec_04XX = filter(x->((x.α_obs-x.α_corr-x.α_comp(x_OR7))^2 + (x.δ_obs-x.δ_corr-x.δ_comp(x_OR7))^2)/(x.σ^2) ≤ χ2_rej, radec_04XX)\n", - "#;" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "#nrow(radec_all), nrow(radec_04XX)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "#radec_04XX.dt_utc_obs" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Vokrouhlicky et al. (2015) obs table\n", - "radec_vok15 = radec_table(\"../data/vokr15_mpc_formatted.dat\")\n", - "\n", - "# vector of RA values from Vokrouhlicky et al. (2015) observations (used for filtering)\n", - "vok15_rav = radec_vok15.α_obs\n", - "\n", - "vok15_ind_filt = findall(x->x ∈ vok15_rav, radec_all.α_obs)\n", - ";" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Vorkouhliky et al. (2015) RA/Dec table\n", - "radec_v15 = radec_all |> @filter(_.α_obs ∈ vok15_rav) |> DataFrame;" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "radec_v15_radec_corr = radec_vok15 |> @select( :α_corr, :δ_corr ) |> DataFrame;" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "res_α = radec_04XX |> @map(_.α_obs - _.α_corr - _.α_comp) |> collect\n", - "res_δ = radec_04XX |> @map(_.δ_obs - _.δ_corr - _.δ_comp) |> collect\n", - "\n", - "res_α_v15 = (radec_v15 |> @map(_.α_obs - _.α_comp)) .- radec_v15_radec_corr.α_corr\n", - "res_δ_v15 = (radec_v15 |> @map(_.δ_obs - _.δ_comp)) .- radec_v15_radec_corr.δ_corr\n", - ";" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "mean(res_α()), mean(res_δ())" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# read astrometric errors from Tholen et al. (2013)\n", - "tho13_errors = readdlm(\"../data/tholenetal2013_opterror.dat\", ',')\n", - "# compute weights\n", - "w_α_tho13 = 1 ./ (tho13_errors[:,1].^2 .+ tho13_errors[:,3].^2 .+ tho13_errors[:,5].^2)\n", - "w_δ_tho13 = 1 ./ (tho13_errors[:,2].^2 .+ tho13_errors[:,4].^2 .+ tho13_errors[:,6].^2)\n", - "\n", - "# Tholen et al. (2013) obs table\n", - "radec_tho13 = radec_table(\"../data/tholen13_mpc_formatted.dat\")\n", - "# vector of RA values from Tholen et al. (2013) observations (used for filtering)\n", - "tho13_rav = radec_tho13.α_obs\n", - "\n", - "tho13_ind_filt = findall(x->x ∈ tho13_rav, radec_all.α_obs)\n", - "tho13_ind_filt_v15 = findall(x->x ∈ tho13_rav, radec_v15.α_obs)\n", - "@show length(tho13_ind_filt), length(tho13_ind_filt_v15)\n", - ";" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "function w8sfarnocchia15(row::NamedTuple)\n", - " return w8sfarnocchia15(row.obscode, row.dt_utc_obs, row.catalog)\n", - "end\n", - "\n", - "function w8sfarnocchia15(df::DataFrame)\n", - " return w8sfarnocchia15.(df.obscode, df.dt_utc_obs, df.catalog)\n", - "end\n", - "\n", - "function w8sfarnocchia15(obscode::Union{Int,String}, dt_utc_obs::DateTime, catalog::String)\n", - " w = 1.0 # unit weight (arcseconds)\n", - " if obscode == \"H01\" # Magdalena Ridge\n", - " return 0.3w\n", - " elseif obscode == \"568\" # Mauna Kea (MPEC 2014-R71)\n", - " return 0.13w #0.2w #0.2\":Verevs17 weight;0.13\":Farnocchia15 weight\n", - " elseif obscode == \"F51\" # PanSTARRS-PS1\n", - " return 0.15w #0.2w #0.2\":Verevs17 weight;0.15\":Farnocchia15 weight\n", - " else\n", - " return w\n", - " end\n", - "end\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "# for each observation batch, count the number of observations made in the same night (?) by the same observatory\n", - "date_obscode_v = radec_04XX |> @map( (Date(_.dt_utc_obs),_.obscode) )\n", - "date_obscode_v_v15 = radec_v15 |> @map( (Date(_.dt_utc_obs),_.obscode) )\n", - "Nv = [count(x->x[1]==i[1] && x[2]==i[2], date_obscode_v) for i in date_obscode_v]\n", - "Nv_v15 = [count(x->x[1]==i[1] && x[2]==i[2], date_obscode_v_v15) for i in date_obscode_v_v15]\n", - "relax_factor = map(x->x>4.0 ? x/4.0 : 1.0, Nv)\n", - "relax_factor_v15 = map(x->x>4.0 ? x/5.0 : 1.0, Nv_v15) # map(x->x>2.0 ? x : 1.0, Nv_v15)\n", - "\n", - "# statistical weights according to Veres et al (2017)\n", - "σ_v17 = radec_04XX.σ\n", - "σ_v17_v15 = w8sfarnocchia15(radec_v15)\n", - "\n", - "w_α = 1 ./ (σ_v17.^2)\n", - "w_δ = 1 ./ (σ_v17.^2)\n", - "\n", - "w_α_v15 = 1 ./ (σ_v17_v15.^2)\n", - "w_δ_v15 = 1 ./ (σ_v17_v15.^2)\n", - "\n", - "# Tholen et al. (2013) astrometry uncertainties\n", - "w_α[tho13_ind_filt] .= w_α_tho13\n", - "w_δ[tho13_ind_filt] .= w_δ_tho13\n", - "w_α_v15[tho13_ind_filt_v15] .= w_α_tho13\n", - "w_δ_v15[tho13_ind_filt_v15] .= w_δ_tho13\n", - "\n", - "w_α = w_α./(relax_factor)\n", - "w_δ = w_δ./(relax_factor)\n", - "\n", - "w_α_v15 = w_α_v15./(relax_factor_v15)\n", - "w_δ_v15 = w_δ_v15./(relax_factor_v15)\n", - ";" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "length(w_τ), length(w_ν), length(w_α), length(w_δ)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "length(w_τ_v15), length(w_ν_v15), length(w_α_v15), length(w_δ_v15)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "nobs_R = nrow(del_all) + nrow(dop_all)\n", - "nobs_O = nrow(radec_04XX)\n", - "nobs = nobs_R + nobs_O\n", - "@show nobs_O, nobs_R, nobs\n", - ";" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "res_O = vcat(res_α, res_δ)\n", - "w_O = vcat(w_α, w_δ)\n", - "\n", - "res_R = vcat(res_τ, res_ν)\n", - "w_R = vcat(w_τ, w_ν)\n", - "\n", - "res_OR = vcat(res_R, res_O)\n", - "w_OR = vcat(w_R, w_O)\n", - ";" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# only Vokrouhliky et al. (2015) data\n", - "res_O_v15 = vcat(res_α_v15, res_δ_v15)\n", - "w_O_v15 = vcat(w_α_v15, w_δ_v15)\n", - "\n", - "res_R_v15 = vcat(res_τ_v15, res_ν_v15)\n", - "w_R_v15 = vcat(w_τ_v15, w_ν_v15)\n", - "\n", - "res_OR_v15 = vcat(res_R_v15, res_O_v15)\n", - "w_OR_v15 = vcat(w_R_v15, w_O_v15)\n", - ";" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# pre-2021, with SRP (\"OR8\")\n", - "\n", - "x_OR7_v15, Γ_OR7_v15 = newtonls(res_OR_v15, w_OR_v15, zeros(get_numvars()))\n", - "σ_OR7_v15 = sqrt.(diag(Γ_OR7_v15))\n", - "@show x_OR7_v15 σ_OR7_v15 x_OR7_v15[end]/σ_OR7_v15[end]\n", - "\n", - "x_OR6_v15, Γ_OR6_v15 = newtonls_6v(res_OR_v15, w_OR_v15, zeros(get_numvars()))\n", - "σ_OR6_v15 = sqrt.(diag(Γ_OR6_v15))\n", - "@show x_OR6_v15 σ_OR6_v15 x_OR6_v15[end]/σ_OR6_v15[end]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "@show x_OR7_v15;\n", - "@show Γ_OR7_v15;" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true, - "tags": [] - }, - "outputs": [], - "source": [ - "# post-2021, NO SRP, 32 ast\n", - "\n", - "x_OR7, Γ_OR7 = newtonls(res_OR, w_OR, zeros(get_numvars()))\n", - "σ_OR7 = sqrt.(diag(Γ_OR7))\n", - "@show x_OR7 σ_OR7 x_OR7[end]/σ_OR7[end]\n", - "\n", - "x_OR6, Γ_OR6 = newtonls_6v(res_OR, w_OR, zeros(get_numvars()))\n", - "σ_OR6 = sqrt.(diag(Γ_OR6))\n", - "@show x_OR6 σ_OR6 x_OR6[end]/σ_OR6[end]\n", - ";" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "(-2.898835108503474 - -2.8988788886209065), (-2.898835108503474 - -2.8988788886209065)/0.02474582241296382" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true, - "tags": [] - }, - "outputs": [], - "source": [ - "# post-2021, NO SRP, 16 ast\n", - "\n", - "x_OR7, Γ_OR7 = newtonls(res_OR, w_OR, zeros(get_numvars()))\n", - "σ_OR7 = sqrt.(diag(Γ_OR7))\n", - "@show x_OR7 σ_OR7 x_OR7[end]/σ_OR7[end]\n", - "\n", - "x_OR6, Γ_OR6 = newtonls_6v(res_OR, w_OR, zeros(get_numvars()))\n", - "σ_OR6 = sqrt.(diag(Γ_OR6))\n", - "@show x_OR6 σ_OR6 x_OR6[end]/σ_OR6[end]\n", - ";" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "@show x_OR7 Γ_OR7" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "#x_O7, Γ_O7 = newtonls(res_O, w_O, zeros(8))\n", - "#σ_O7 = sqrt.(diag(Γ_O7))\n", - "#@show x_O7 σ_O7 x_O7[end]/σ_O7[end]\n", - "#\n", - "#x_R7, Γ_R7 = newtonls(res_R, w_R, zeros(8))\n", - "#σ_R7 = sqrt.(diag(Γ_R7))\n", - "#@show x_R7 σ_R7 x_R7[end]/σ_R7[end]\n", - "#\n", - "#x_O6, Γ_O6 = newtonls(res_O, w_O, zeros(8))\n", - "#σ_O6 = sqrt.(diag(Γ_O6))\n", - "#@show x_O6 σ_O6 x_O6[end]/σ_O6[end]\n", - "#\n", - "#x_R6, Γ_R6 = newtonls(res_R, w_R, zeros(8))\n", - "#σ_R6 = sqrt.(diag(Γ_R6))\n", - "#@show x_R6 σ_R6 x_R6[end]/σ_R6[end]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "#@show x_OR7 Γ_OR7" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# 32 ast\n", - "#sqrt(Γ_OR7[7,7])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# 16 ast\n", - "sqrt(Γ_OR7[7,7])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "abs(-2.894053174513909 +2.8941071491385224)/sqrt(Γ_OR7[7,7])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "#@show x_OR7 Γ_OR7" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "@show x_OR7.*vcat(1e-8ones(6), 1e-14)\n", - ";" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "#@show x_OR7.*vcat(1e-8ones(6), 1e-14, 1e-13)\n", - "#;" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Fit quality statistics: $\\chi^2$, normalized RMS" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The gravity-only solution post-fit normalized RMS for the optical, radar, and combined optical and radar astrometry are, respectively, 0.448, 14.720 and 1.0259. Moreover, the gravity-only solution post-fit right-ascension weighted mean residual is 0.072$~$arcsec, while the declination weighted mean residual is 0.218$~$arcsec. The weighted time-delay and Doppler shift mean residual for the gravity-only solution are, respectively, -0.06 us and 0.043 Hz. On the other hand, the non-gravitational solution normalized RMS for the optical, radar, and combined optical and radar astrometry are, respectively, 0.313, 0.387 and 0.314. The right-ascension, declination, time-delay and Doppler shift weighted mean residuals, are, respectively, 0.008$~$arcsec, 0.003$~$arcsec, -0.0003$~$us and -0.31$~$Hz." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "mean(res_α(x_OR6), weights(w_α)), mean(res_δ(x_OR6), weights(w_δ))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "mean(res_α(x_OR7), weights(w_α)), mean(res_δ(x_OR7), weights(w_δ))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "mean(res_τ(x_OR6), weights(w_τ)), mean(res_ν(x_OR6), weights(w_ν))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "mean(res_τ(x_OR7), weights(w_τ)), mean(res_ν(x_OR7), weights(w_ν))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "nrms(res_O(x_OR6), w_O), nrms(res_R(x_OR6), w_R), nrms(res_OR(x_OR6), w_OR)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "nrms(res_O(x_OR7), w_O), nrms(res_R(x_OR7), w_R), nrms(res_OR(x_OR7), w_OR)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "chi2(res_OR(x_OR7), w_OR)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, + "metadata": { + "scrolled": true + }, "outputs": [], "source": [ - "chi2(res_OR(x_OR6), w_OR)" + "using NEOs\n", + "using DataFrames\n", + "using TaylorIntegration\n", + "using PlanetaryEphemeris\n", + "using JLD2\n", + "using Statistics\n", + "using LinearAlgebra\n", + "using Dates\n", + "using StatsBase\n", + "using Gaston\n", + "using DelimitedFiles" ] }, { @@ -855,25 +49,14 @@ "metadata": {}, "outputs": [], "source": [ - "# signal-to-noise ratio\n", - "abs(x_OR7[end]/σ_OR7[end])" + "dq = set_variables(\"δx\", order=5, numvars=7)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "# Semimajor axis drift" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "jd0 = datetime2julian(DateTime(2020,12,17)) #Julian date (TDB) of integration initial time\n", - "et0 = (jd0 - J2000)*daysec" + "## Load Sun, Earth and Apophis ephemeris" ] }, { @@ -882,7 +65,7 @@ "metadata": {}, "outputs": [], "source": [ - "q00 = [-0.18034747703273316, 0.9406910666200128, 0.3457360259054398, -0.016265942170279046, 4.392889725556651e-5, -0.00039519931615139716]" + "sol_bwd = JLD2.load(\"Apophis_bwd.jld2\", \"asteph\")" ] }, { @@ -891,11 +74,7 @@ "metadata": {}, "outputs": [], "source": [ - "if get_numvars() == 7\n", - " @show q00 .+ (vcat(1e-8*ones(6), 1e-14).*x_OR7)[1:6]\n", - "elseif get_numvars() == 8\n", - " q00 .+ (vcat(1e-8*ones(6), 1e-14, 1e-13).*x_OR7)[1:6]\n", - "end" + "sol_fwd = JLD2.load(\"Apophis_fwd.jld2\", \"asteph\")" ] }, { @@ -904,38 +83,61 @@ "metadata": {}, "outputs": [], "source": [ - "xas = q00 .- NEOs.kmsec2auday(NEOs.sun_pv(et0)) # .+ (vcat(1e-8*ones(6), 1e-14).*x_OR7)[1:6]" + "sseph::TaylorInterpolant{Float64,Float64,2} = loadpeeph(NEOs.sseph, sol_bwd.t0+sol_bwd.t[end], sol_fwd.t0+sol_fwd.t[end])\n", + "eph_su::TaylorInterpolant{Float64,Float64,2} = selecteph(sseph, su)\n", + "eph_ea::TaylorInterpolant{Float64,Float64,2} = selecteph(sseph, ea)\n", + "\n", + "# NEO\n", + "# Change t, x, v units, resp., from days, au, au/day to sec, km, km/sec\n", + "xva_bwd(et) = auday2kmsec(sol_bwd(et/daysec)[1:6])\n", + "xva_fwd(et) = auday2kmsec(sol_fwd(et/daysec)[1:6])\n", + "xva(et) = bwdfwdeph(et, sol_bwd, sol_fwd)\n", + "# Earth\n", + "# Change x, v units, resp., from au, au/day to km, km/sec\n", + "xve(et) = auday2kmsec(eph_ea(et/daysec))\n", + "# Sun\n", + "# Change x, v units, resp., from au, au/day to km, km/sec\n", + "xvs(et) = auday2kmsec(eph_su(et/daysec))" ] }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "scrolled": true + }, "outputs": [], "source": [ - "ϵ0_deg = 84381.448/3600;\n", - "m_eq2ecl = Rx(deg2rad(ϵ0_deg));\n", - "m_xv_eq2ecl = hcat(vcat(m_eq2ecl, zeros(3,3)), vcat(zeros(3,3), m_eq2ecl));" + "# get optical observations\n", + "radec_2004_2020 = read_radec_mpc(joinpath(pkgdir(NEOs), \"data\", \"99942_2004_2020.dat\"))\n", + "radec_2020_2021 = read_radec_mpc(joinpath(pkgdir(NEOs), \"data\", \"99942_2020_2021.dat\"))\n", + "radec = vcat(radec_2004_2020,radec_2020_2021)" ] }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "scrolled": true + }, "outputs": [], "source": [ - "kep0 = vcat([pv2kep( m_xv_eq2ecl*(xas .+ (vcat(1e-8*ones(6), 1e-14).*dq)[1:6]) , μ[1], jd0)...][1:6],1e-14(dq[7])); # x_OR6 +\n", - "#kep0 = vcat([pv2kep( m_xv_eq2ecl*(xas .+ (vcat(1e-8*ones(6), 1e-14, 1e-13).*dq)[1:6]) , μ[1], jd0)...][1:6],1e-14(dq[7]),1e-13(dq[8])); # x_OR6 +\n", - "kep0()" + "# get radar observations\n", + "deldop_2005_2013 = read_radar_jpl(joinpath(pkgdir(NEOs), \"data\", \"99942_RADAR_2005_2013.dat\"))\n", + "deldop_2021 = read_radar_jpl(joinpath(pkgdir(NEOs), \"data\", \"99942_RADAR_2021.dat\"))\n", + "deldop = vcat(deldop_2005_2013,deldop_2021)" ] }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "scrolled": true + }, "outputs": [], "source": [ - "kep_OR7 = kep0(x_OR7)" + "# Compute optical residuals\n", + "@time res_radec_all, w_radec_all = NEOs.residuals(radec; xvs, xve, xva);" ] }, { @@ -944,47 +146,81 @@ "metadata": {}, "outputs": [], "source": [ - "# WARNING: these are ecliptic elements!!! also, not sure about the order in the Ω ↔ ω angles\n", - "kep_JPL216 = [0.1914347568286823, 0.7460634895090179, 2459424.558511739019, 203.9662692226386, 126.5900238507279, 3.338874353538208, -2.901085583204654E-14, 5.E-13]\n", - "#σ_kep_JPL204 = [2.2197e-09, 2.2003e-09, 1.2886e-06, 1.4045e-05, 1.4205e-05, 2.4295e-07, 3.039E-16]\n", - ";" + "# Compute radar residuals\n", + "@time res_del, w_del, res_dop, w_dop = NEOs.residuals(deldop; xvs, xve, xva, niter=10, tord=10);" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "kep_JPL216" + "## Process optical data (filter, weight, debias)" ] }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "scrolled": true + }, "outputs": [], "source": [ - "orbel_OR7_TN = kep_OR7 .+ dq\n", - "#orbel_OR7_TN = kep_JPL216 .+ dq # length(kep_JPL216) == 8" + "# filter out biased observations from observatory 217 on 28-Jan-2021\n", + "df_radec = DataFrame(radec)\n", + "# add residuals and weights to optical astrometry DataFrame\n", + "df_radec[!, :res_α] .= res_radec_all[1:round(Int,length(res_radec_all)/2)]\n", + "df_radec[!, :res_δ] .= res_radec_all[1+round(Int,length(res_radec_all)/2):end]\n", + "df_radec[!, :w_α] .= w_radec_all[1:round(Int,length(res_radec_all)/2)]\n", + "df_radec[!, :w_δ] .= w_radec_all[1+round(Int,length(res_radec_all)/2):end]\n", + "filter!(\n", + " x->(Date(x.date) != Date(2021, 1, 28)),\n", + " df_radec\n", + ")\n", + "\n", + "# read astrometric errors from Tholen et al. (2013)\n", + "tho13_errors = readdlm(joinpath(pkgdir(NEOs), \"data\", \"tholenetal2013_opterror.dat\"), ',')\n", + "# compute weights\n", + "w_α_tho13 = 1 ./ (tho13_errors[:,1].^2 .+ tho13_errors[:,3].^2 .+ tho13_errors[:,5].^2)\n", + "w_δ_tho13 = 1 ./ (tho13_errors[:,2].^2 .+ tho13_errors[:,4].^2 .+ tho13_errors[:,6].^2)\n", + "# Tholen et al. (2013) obs table\n", + "radec_tho13 = DataFrame(read_radec_mpc(joinpath(pkgdir(NEOs), \"test\", \"data\", \"99942_Tholen_etal_2013.dat\")))\n", + "# vector of RA values from Tholen et al. (2013) observations (used for filtering)\n", + "tho13_α = radec_tho13[!,:α]\n", + "# set weights in Tholen et al. (2013) astrometry corresponding to associated uncertainties\n", + "df_radec[in.(df_radec.α, Ref(tho13_α)),:w_α] = w_α_tho13\n", + "df_radec[in.(df_radec.α, Ref(tho13_α)),:w_δ] = w_δ_tho13\n", + "\n", + "# Relaxation factor (account for correlations in optical astrometry data)\n", + "# for each observation batch, count the number of observations made in\n", + "# the same night by the same observatory\n", + "# Ref: Veres et al. (2017)\n", + "date_site_v = select(df_radec, :date => ByRow(Date), :observatory)\n", + "Nv = [count(x->x.date_Date==i.date_Date && x.observatory==i.observatory, eachrow(date_site_v)) for i in eachrow(date_site_v)]\n", + "relax_factor = map(x->x>4.0 ? x/4.0 : 1.0, Nv)\n", + "# inflate uncertainties (i.e., relax weights) by relaxation factor\n", + "df_radec[!, :w_α] .= (df_radec.w_α)./relax_factor\n", + "df_radec[!, :w_δ] .= (df_radec.w_δ)./relax_factor\n", + "\n", + "# update optical residuals and weights\n", + "res_radec = vcat(df_radec.res_α, df_radec.res_δ)\n", + "w_radec = vcat(df_radec.w_α, df_radec.w_δ)\n", + ";" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "var_names_str = prod(get_variable_names().*\" \")\n", - "vars_old = get_variables()\n", - "vars_new = set_variables(eltype(vars_old[1]), var_names_str[1:prevind(var_names_str, end)], order=vars_old[1].order)" + "# check if all RA/Dec observations are CCD:\n", + "\n", + "all(\n", + " select(\n", + " df_radec, \n", + " :obstech => ByRow(x->x == \"C\" || x == \"X\")\n", + " )[!,1]\n", + ")" ] }, { @@ -993,7 +229,9 @@ "metadata": {}, "outputs": [], "source": [ - "x_OR7[7], σ_OR7[7]" + "# there are 512 RA/Dec obs without catalog info\n", + "\n", + "length(df_radec[isunknown.(df_radec.catalogue),:catalogue])" ] }, { @@ -1002,8 +240,11 @@ "metadata": {}, "outputs": [], "source": [ - "# OR7, 16 ast\n", - "ady_kep_au_My = 1e6yr*yarkp2adot(orbel_OR7_TN[7], orbel_OR7_TN[2]/(1.0 - orbel_OR7_TN[1]), orbel_OR7_TN[1], μ[1]) # au/Myr" + "# Construct vector of residuals and weights\n", + "\n", + "res = vcat(res_radec, res_del, res_dop)\n", + "w = vcat(w_radec, w_del, w_dop)\n", + ";" ] }, { @@ -1012,8 +253,10 @@ "metadata": {}, "outputs": [], "source": [ - "# OR7, 16 ast\n", - "ady_kep_m_y = 1e3au*yr*yarkp2adot(orbel_OR7_TN[7], orbel_OR7_TN[2]/(1.0 - orbel_OR7_TN[1]), orbel_OR7_TN[1], μ[1]) #m/yr" + "### Perform 7-DOF orbital fit to optical and radar astrometry data\n", + "\n", + "@show success, δx_OR7, Γ_OR7 = newtonls(res, w, zeros(get_numvars()), 10)\n", + "x_OR7 = sol_fwd(sol_fwd.t0)(δx_OR7);" ] }, { @@ -1022,7 +265,10 @@ "metadata": {}, "outputs": [], "source": [ - "- 198.70516356300044 - (- 198.70886945315476)" + "### Perform 6-DOF orbital fit (i.e., without Yarkovsky effect) to optical and radar astrometry data\n", + "\n", + "δx_OR6, Γ_OR6 = newtonls_6v(res, w, zeros(get_numvars()), 10)\n", + "x_OR6 = sol_fwd(sol_fwd.t0)(δx_OR6);" ] }, { @@ -1031,26 +277,27 @@ "metadata": {}, "outputs": [], "source": [ - "t_car2kep(x) = TaylorSeries.jacobian(kep0, x)'\n", - "t_kep2ady(x, ady_kep) = TaylorSeries.gradient(ady_kep)(x)" + "### Print results\n", + "\n", + "# orbital fit\n", + "@show success\n", + "@show x_OR7 # cartesian state vector at solution reference epoch\n", + "@show sqrt.(diag(Γ_OR7)).*vcat(1e-8ones(6),1e-14) # uncertainties at solution reference epoch (1-sigma)\n", + ";" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "Γ_OR7_kep = (t_car2kep(x_OR7)')*Γ_OR7*(t_car2kep(x_OR7))" + "# Fit quality statistics: $\\chi^2$, normalized RMS" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "1e4sqrt( (t_kep2ady(kep_OR7, ady_kep_au_My)')*Γ_OR7_kep*t_kep2ady(kep_OR7, ady_kep_au_My) )" + "The results below correspond to Table 1 of Pérez-Hernández and Benet (2022)." ] }, { @@ -1059,16 +306,36 @@ "metadata": {}, "outputs": [], "source": [ - "sqrt( (t_kep2ady(kep_OR7, ady_kep_m_y)')*Γ_OR7_kep*t_kep2ady(kep_OR7, ady_kep_m_y) )" + "# post-fit statistics\n", + "\n", + "nradec = length(res_radec)\n", + "res_α = view(res_radec, 1:nradec÷2)\n", + "res_δ = view(res_radec, 1+nradec÷2:nradec)\n", + "res_τ = res_del\n", + "res_ν = res_dop\n", + "w_α = view(w_radec, 1:nradec÷2)\n", + "w_δ = view(w_radec, 1+nradec÷2:nradec)\n", + "res_α_OR6 = res_α(δx_OR6)\n", + "res_δ_OR6 = res_δ(δx_OR6)\n", + "res_α_OR7 = res_α(δx_OR7)\n", + "res_δ_OR7 = res_δ(δx_OR7)\n", + "\n", + "@show nrms_radec = nrms(res_radec(δx_OR7),w_radec)\n", + "@show nrms_radec = nrms(vcat(res_del,res_dop)(δx_OR7),vcat(w_del,w_dop))\n", + "@show nrms_optrad = nrms(res(δx_OR7),w)\n", + "@show mean_ra = mean(res_ra(δx_OR7), weights(w_ra))\n", + "@show mean_dec = mean(res_dec(δx_OR7), weights(w_dec))\n", + "@show mean_del = mean(res_del(δx_OR7), weights(w_del))\n", + "@show mean_dop = mean(res_dop(δx_OR7), weights(w_dop))\n", + "@show chi2_optrad = chi2(res(δx_OR7),w)\n", + ";" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "sqrt( (t_kep2ady(kep_OR7, ady_kep_m_y)')*Γ_OR7_kep*t_kep2ady(kep_OR7, ady_kep_m_y) )" + "# Semimajor axis drift" ] }, { @@ -1077,7 +344,8 @@ "metadata": {}, "outputs": [], "source": [ - "t_car2kep(x_OR7)" + "jd0 = J2000 + sol_fwd.t0 #Julian date (TDB) of integration initial time\n", + "et0 = (jd0 - J2000)*daysec" ] }, { @@ -1086,7 +354,8 @@ "metadata": {}, "outputs": [], "source": [ - "σ_OR7_kep = sqrt.(diag(Γ_OR7_kep))" + "# Apophis heliocentric position and velocity at solution reference epoch (plus autodiff perturbations)\n", + "xas = sol_fwd(sol_fwd.t0)[1:6] .- eph_su(jd0-J2000)" ] }, { @@ -1095,7 +364,8 @@ "metadata": {}, "outputs": [], "source": [ - "kep_OR7" + "# Orbital elements of new solution (ecliptic frame at J2000.0 epoch)\n", + "pv2kep(xas(δx_OR7), jd=jd0, frame=:ecliptic)" ] }, { @@ -1104,7 +374,17 @@ "metadata": {}, "outputs": [], "source": [ - "kep_JPL216" + "# JPL solution #216 for Apophis\n", + "kep_JPL216 = [\n", + " 0.1914347568286823, \n", + " 0.7460634895090179, \n", + " 2459424.558511739019, \n", + " 203.9662692226386, \n", + " 126.5900238507279, \n", + " 3.338874353538208, \n", + " -2.901085583204654E-14, \n", + " 5.E-13\n", + "]" ] }, { @@ -1113,17 +393,44 @@ "metadata": {}, "outputs": [], "source": [ - "σ_OR7_kep" + "kep_JT = pv2kep(xas, jd=jd0, frame=:ecliptic) # Orbital elements of JT solution at ref. epoch\n", + "kep_OR7 = pv2kep(xas(δx_OR7), jd=jd0, frame=:ecliptic) # Orbital elements of new solution at ref. epoch" ] }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "scrolled": true + }, "outputs": [], "source": [ - "Γ_OR7_kep = (t_car2kep(x_OR7)')*Γ_OR7*(t_car2kep(x_OR7));\n", - "Γ_OR6_kep = (t_car2kep(x_OR6)[1:6,1:6]')*Γ_OR6*t_car2kep(x_OR6)[1:6,1:6];" + "# Jacobian of transformation: cartesian state -> Keplerian elements\n", + "t_car2kep(x) = hcat(\n", + " TaylorSeries.gradient(kep_JT.e), \n", + " TaylorSeries.gradient(kep_JT.q),\n", + " TaylorSeries.gradient(kep_JT.tp),\n", + " TaylorSeries.gradient(kep_JT.Ω),\n", + " TaylorSeries.gradient(kep_JT.ω),\n", + " TaylorSeries.gradient(kep_JT.i),\n", + " TaylorSeries.gradient(kep_JT.M),\n", + " TaylorSeries.gradient(kep_JT.a),\n", + " )(x)\n", + "\n", + "# semimajor axis drift \\dot{a} due to Yarkovsky effect in units of au per day, as a function of initial condition\n", + "ady_kep = yarkp2adot(sol_fwd(sol_fwd.t0)[7], kep_JT.a, kep_JT.e)\n", + "\n", + "# semimajor axis drift, au per million years\n", + "ady_kep_au_My = 1e6yr*ady_kep\n", + "\n", + "# semimajor axis drift, meters per year\n", + "ady_kep_m_y = 1e3au*yr*ady_kep\n", + "\n", + "# Jacobian of transformation: cartesian state -> semimajor axis drift in au per million years\n", + "t_kep2ady_au_My(x) = TaylorSeries.gradient(ady_kep_au_My)(x)\n", + "\n", + "# Jacobian of transformation: cartesian state -> semimajor axis drift in meters per year\n", + "t_kep2ady_m_y(x) = TaylorSeries.gradient(ady_kep_m_y)(x)" ] }, { @@ -1132,7 +439,8 @@ "metadata": {}, "outputs": [], "source": [ - "jd0" + "# mean semimajor axis drift, au per million years\n", + "ady_kep_au_My_OR7 = ady_kep_au_My(δx_OR7) # au/Myr" ] }, { @@ -1141,35 +449,8 @@ "metadata": {}, "outputs": [], "source": [ - "for (x, Γ_kep, sol_str) in [(x_OR6, Γ_OR6_kep, \"OR6\"), (x_OR7, Γ_OR7_kep, \"OR7\")]\n", - " println(\"* ----===OOO===--- *\", sol_str)\n", - " @show kep0(x) sqrt.(diag(Γ_kep))\n", - "end" - ] - }, - { - "cell_type": "raw", - "metadata": {}, - "source": [ - "#OR6\n", - "\n", - "0.19150875(472) \\pm 1.28\\times 10^{-9}\n", - "0.74585316(035) \\pm 1.23\\times 10^{-9}\n", - "2459101.040820(608) \\pm 7.49\\times 10^{-7}\n", - "204.04254(376) \\pm 7.44\\times 10^{-6}\n", - "126.65331(533) \\pm 7.58\\times 10^{-6}\n", - "3.336760(462) \\pm 1.36\\times 10^{-7}\n", - " 0.0 (fixed) \n", - "\n", - "# OR7\n", - "\n", - "0.19150886(716) \\pm 1.60\\times 10^{-9}\n", - "0.74585305(033) \\pm 1.54\\times 10^{-9}\n", - "2459101.04092(537) \\pm 1.17\\times 10^{-6}\n", - "204.04199(116) \\pm 8.81\\times 10^{-6}\n", - "126.65396(094) \\pm 9.37\\times 10^{-6}\n", - "3.336773(201) \\pm 1.74\\times 10^{-7}\n", - "-2.8(988) \\pm 2.48\\times 10^{-2}" + "# mean semimajor axis drift, meters per year\n", + "ady_kep_m_y_OR7 = 1e3au*yr*ady_kep(δx_OR7) # m/yr" ] }, { @@ -1178,58 +459,18 @@ "metadata": {}, "outputs": [], "source": [ - "for (x, Γ, sol_str) in [(x_OR6, Γ_OR6, \"OR6\"), (x_OR7, Γ_OR7, \"OR7\")]\n", - " println(\"* ----===OOO===--- *\", sol_str)\n", - " σ_sol_ = sqrt.(diag(Γ))\n", - " σ_sol = length(σ_sol_)==6 ? 1e-8σ_sol_[1:6] : vcat(1e-8σ_sol_[1:6], 1e-14σ_sol_[7])\n", - " @show vcat(q00, 0).+vcat(1e-8x[1:6], 1e-14x[7]) σ_sol\n", - "end" - ] - }, - { - "cell_type": "raw", - "metadata": {}, - "source": [ - "#OR6\n", - "\n", - "-0.18034816(472) \\pm 7.05\\times 10^{-9} \n", - " 0.94069116(764) \\pm 1.71\\times 10^{-9} \n", - " 0.34573641(385) \\pm 4.03\\times 10^{-9} \n", - "-0.0162659357(312) \\pm 3.31\\times 10^{-11}\n", - " 4.39163(792)\\times 10^{-5} \\pm 6.68\\times 10^{-11}\n", - "-0.000395210(949) \\pm 1.24\\times 10^{-10}\n", - " 0.0 (fixed) \n", - "\n", - "#OR7\n", - "\n", - "-0.18034828(526) \\pm 7.12\\times 10^{-9} \n", - " 0.94069105(951) \\pm 1.94\\times 10^{-9} \n", - " 0.34573599(029) \\pm 5.41\\times 10^{-9} \n", - "-0.0162659397(886) \\pm 4.79\\times 10^{-11}\n", - " 4.39154(800)\\times 10^{-5} \\pm 6.72\\times 10^{-11}\n", - "-0.000395204(013) \\pm 1.38\\times 10^{-10}\n", - "-2.8(988)\\times 10^{-14} \\pm 2.48\\times 10^{-16}\n" + "# uncertainty in semimajor axis drift, au per million years\n", + "sqrt( t_kep2ady_m_y(δx_OR7)'*Γ*t_kep2ady_m_y(δx_OR7) )" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "res_α_OR6 = res_α(x_OR6)\n", - "res_δ_OR6 = res_δ(x_OR6)\n", - "\n", - "res_α_OR7 = res_α(x_OR7)\n", - "res_δ_OR7 = res_δ(x_OR7)\n", - ";" + "# uncertainty in semimajor axis drift, au per million years\n", + "sqrt( t_kep2ady_au_My(δx_OR7)'*Γ*t_kep2ady_au_My(δx_OR7) )" ] }, { @@ -1238,7 +479,8 @@ "metadata": {}, "outputs": [], "source": [ - "mean(res_α_OR6), std(res_α_OR6), mean(res_δ_OR6), std(res_δ_OR6)" + "# covariance matrix of Keplerian elements\n", + "Γ_kep = t_car2kep(δx_OR7)'*Γ*t_car2kep(δx_OR7)" ] }, { @@ -1247,7 +489,9 @@ "metadata": {}, "outputs": [], "source": [ - "mean(res_α_OR7), std(res_α_OR7), mean(res_δ_OR7), std(res_δ_OR7)" + "#formal uncertainties in Keplerian elements (1-sigma)\n", + "\n", + "sqrt.(diag(Γ_kep))" ] }, { @@ -1257,15 +501,6 @@ "# Plots: optical astrometry post-fit residuals" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "methods(fit)" - ] - }, { "cell_type": "code", "execution_count": null, @@ -1352,43 +587,6 @@ "mean(res_α_OR7), mean(res_δ_OR7)" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# use findall instead of filter, since we're looking for the indices of ZTF observations\n", - "#ztf_obs_indx = findall(x->x==\"I41\", radec_all.obscode)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "#res_α_OR7[ztf_obs_indx]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "#res_δ_OR7[ztf_obs_indx]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "#sqrt( mean(res_δ_OR7[ztf_obs_indx].^2 .+ res_α_OR7[ztf_obs_indx].^2) )" - ] - }, { "cell_type": "code", "execution_count": null, @@ -1486,176 +684,39 @@ " #linetype = :tab10,\n", " axesconf=\"set lmargin at screen 0.85+$gap; set rmargin at screen 0.99; set tmargin at screen 0.8; set bmargin at screen 0.1\",\n", " yrange = (ymin_plt, ymax_plt),\n", - " format = \"''\",\n", - " label=\"1 '(c)' at graph 0.15,0.925\"\n", - " ),\n", - " handle = 4\n", - ")\n", - "Gaston.plot!(\n", - " h_d.weights./2,\n", - " bc_d,\n", - " supp=[h_d.weights./2 bs_d./2],\n", - " curveconf=\"w boxxy notit fs transparent solid 0.4\"\n", - ")\n", - "\n", - "Gaston.plot([p1 nothing ; p2 p4])\n", - "\n", - "#Gaston.save(term = \"cairolatex\", output = \"radec_hist_NEW_SANSSERIF.tex\", saveopts = \"pdf standalone color dashed transparent size 7in,5in font ',12'\")\n", - "#run(`lualatex radec_hist_NEW_SANSSERIF.tex`)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true, - "tags": [] - }, - "outputs": [], - "source": [ - "?writedlm" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "open(\"optical_residuals_OR6.txt\", \"w\") do io\n", - " writedlm(io, [res_α_OR6 res_δ_OR6])\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "length(res_α)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "mean(w_α.*res_α(x_OR7)), std(w_α.*res_α(x_OR7))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "mean(w_δ.*res_δ(x_OR7)), std(w_δ.*res_δ(x_OR7))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "mean(res_τ(x_OR7)), std(res_τ(x_OR7))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "mean(res_ν(x_OR7)), std(res_ν(x_OR7))" + " format = \"''\",\n", + " label=\"1 '(c)' at graph 0.15,0.925\"\n", + " ),\n", + " handle = 4\n", + ")\n", + "Gaston.plot!(\n", + " h_d.weights./2,\n", + " bc_d,\n", + " supp=[h_d.weights./2 bs_d./2],\n", + " curveconf=\"w boxxy notit fs transparent solid 0.4\"\n", + ")\n", + "\n", + "Gaston.plot([p1 nothing ; p2 p4])\n", + "\n", + "#Gaston.save(term = \"cairolatex\", output = \"radec_hist.tex\", saveopts = \"pdf standalone color dashed transparent size 7in,5in font ',12'\")\n", + "#run(`lualatex radec_hist.tex`)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## OR6" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "mean(res_τ(x_OR6)), std(res_τ(x_OR6))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "mean(res_ν(x_OR6)), std(res_ν(x_OR6))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ind_del_0506 = 1:2\n", - "ind_del_1213 = 3:17\n", - "ind_del_2121 = 18:20 #18:19\n", - "\n", - "ind_dop_0506 = 1:5\n", - "ind_dop_1213 = 6:29\n", - "ind_dop_2121 = 30:30\n", - ";" + "# Plots: radar astrometry post-fit residuals" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "del_index = findall(x->x, deldop_all.delay_index)\n", - "dop_index = findall(x->x, deldop_all.doppler_index)\n", + "del_index = findall(x->x.Δτ_units==\"us\", deldop)\n", + "dop_index = findall(x->x.Δν_units==\"Hz\", deldop)\n", "length(del_index), length(dop_index)" ] }, @@ -1665,7 +726,7 @@ "metadata": {}, "outputs": [], "source": [ - "nrow(deldop_all)" + "length(deldop)" ] }, { @@ -1674,20 +735,15 @@ "metadata": {}, "outputs": [], "source": [ - "del_dates_plot_ = (deldop_all |> @filter(_.delay_index) |> DataFrame).dt_utc_obs\n", - "dop_dates_plot_ = (deldop_all |> @filter(_.doppler_index) |> DataFrame).dt_utc_obs\n", + "del_dates_plot_ = date.(deldop[del_index])\n", + "dop_dates_plot_ = date.(deldop[dop_index])\n", "\n", "del_dates_plot = Dates.format.(del_dates_plot_, \"dd/mm/yy HHMM\")\n", "dop_dates_plot = Dates.format.(dop_dates_plot_, \"dd/mm/yy HHMM\")\n", "\n", - "length(del_dates_plot), length(dop_dates_plot)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Plots: radar astrometry post-fit residuals" + "σ_del_all = delay_sigma.(deldop[del_index])\n", + "σ_dop_all = doppler_sigma.(deldop[dop_index])\n", + ";" ] }, { @@ -1698,7 +754,7 @@ "source": [ "p1 = Gaston.plot(\n", " del_index,\n", - " res_τ(x_OR6),\n", + " res_τ(δx_OR6),\n", " supp = σ_del_all,\n", " key=\"box spacing 1.2 height 0.4\", #at screen 0.835, screen 0.78 \n", " curveconf = \"u 1:2:3 w errorbars pt 6 lw 4 tit 'gravity-only'\",\n", @@ -1719,7 +775,7 @@ ")\n", "Gaston.plot!(\n", " del_index,\n", - " res_τ(x_OR7),\n", + " res_τ(δx_OR7),\n", " supp = σ_del_all,\n", " curveconf = \"u 1:2:3 w errorbars pt 4 lw 4 tit 'non-gravitational'\",\n", " lc = 2\n", @@ -1727,7 +783,7 @@ "\n", "p2 = Gaston.plot(\n", " dop_index,\n", - " res_ν(x_OR6),\n", + " res_ν(δx_OR6),\n", " supp = σ_dop_all,\n", " curveconf = \"u 1:2:3 w errorbars pt 6 lw 4 tit 'gravity-only'\",\n", " key=\"box spacing 1.2 height 0.4\",\n", @@ -1748,7 +804,7 @@ ")\n", "Gaston.plot!(\n", " dop_index,\n", - " res_ν(x_OR7),\n", + " res_ν(δx_OR7),\n", " supp = σ_dop_all,\n", " curveconf = \"u 1:2:3 w errorbars pt 4 lw 4 tit 'non-gravitational'\",\n", " lc = 2\n", @@ -1793,7 +849,7 @@ "\n", "p1 = Gaston.plot(\n", " del_dates_plot,\n", - " res_τ(x_OR6),\n", + " res_τ(δx_OR6),\n", " supp = σ_del_all,\n", " curveconf = \"u 1:3:4 w errorbars pt 6 lw $lw_ths_plt\",\n", " Axes(\n", @@ -1823,7 +879,7 @@ ")\n", "Gaston.plot!(\n", " del_dates_plot,\n", - " res_τ(x_OR7),\n", + " res_τ(δx_OR7),\n", " supp = σ_del_all,\n", " curveconf = \"u 1:3:4 w errorbars pt 4 lw $lw_ths_plt\",\n", " lc = 2\n", @@ -1832,7 +888,7 @@ "\n", "p2 = Gaston.plot(\n", " del_dates_plot,\n", - " res_τ(x_OR6),\n", + " res_τ(δx_OR6),\n", " supp = σ_del_all,\n", " curveconf = \"u 1:3:4 w errorbars pt 6 lw $lw_ths_plt\",\n", " Axes(\n", @@ -1861,7 +917,7 @@ ")\n", "Gaston.plot!(\n", " del_dates_plot,\n", - " res_τ(x_OR7),\n", + " res_τ(δx_OR7),\n", " supp = σ_del_all,\n", " curveconf = \"u 1:3:4 w errorbars pt 4 lw $lw_ths_plt\",\n", " lc = 2\n", @@ -1870,7 +926,7 @@ "\n", "p3 = Gaston.plot(\n", " del_dates_plot,\n", - " res_τ(x_OR6),\n", + " res_τ(δx_OR6),\n", " supp = σ_del_all,\n", " curveconf = \"u 1:3:4 w errorbars pt 6 lw $lw_ths_plt\",\n", " Axes(\n", @@ -1897,7 +953,7 @@ ")\n", "Gaston.plot!(\n", " del_dates_plot,\n", - " res_τ(x_OR7),\n", + " res_τ(δx_OR7),\n", " supp = σ_del_all,\n", " curveconf = \"u 1:3:4 w errorbars pt 4 lw $lw_ths_plt\",\n", " lc = 2\n", @@ -1913,17 +969,11 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "@show res_ν(x_OR7)\n", - ";" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, + "metadata": { + "jupyter": { + "source_hidden": true + } + }, "outputs": [], "source": [ "tmargin_ths_plt = 0.5 - tmargin_bmargin_gap\n", @@ -1932,7 +982,7 @@ "Gaston.set(preamble=\"unset grid\")\n", "p4 = Gaston.plot(\n", " dop_dates_plot,\n", - " res_ν(x_OR6),\n", + " res_ν(δx_OR6),\n", " supp = σ_dop_all,\n", " curveconf = \"u 1:3:4 w errorbars pt 6 lw $lw_ths_plt\",\n", " Axes(\n", @@ -1959,7 +1009,7 @@ ")\n", "Gaston.plot!(\n", " dop_dates_plot,\n", - " res_ν(x_OR7),\n", + " res_ν(δx_OR7),\n", " supp = σ_dop_all,\n", " curveconf = \"u 1:3:4 w errorbars pt 4 lw $lw_ths_plt\",\n", " lc = 2\n", @@ -1968,7 +1018,7 @@ "\n", "p5 = Gaston.plot(\n", " dop_dates_plot,\n", - " res_ν(x_OR6),\n", + " res_ν(δx_OR6),\n", " supp = σ_dop_all,\n", " curveconf = \"u 1:3:4 w errorbars pt 6 lw $lw_ths_plt\",\n", " Axes(\n", @@ -1995,7 +1045,7 @@ ")\n", "Gaston.plot!(\n", " dop_dates_plot,\n", - " res_ν(x_OR7),\n", + " res_ν(δx_OR7),\n", " supp = σ_dop_all,\n", " curveconf = \"u 1:3:4 w errorbars pt 4 lw $lw_ths_plt\",\n", " lc = 2\n", @@ -2004,7 +1054,7 @@ "\n", "p6 = Gaston.plot(\n", " dop_dates_plot,\n", - " res_ν(x_OR6),\n", + " res_ν(δx_OR6),\n", " supp = σ_dop_all,\n", " curveconf = \"u 1:3:4 w errorbars pt 6 lw $lw_ths_plt\",\n", " Axes(\n", @@ -2030,7 +1080,7 @@ ")\n", "Gaston.plot!(\n", " dop_dates_plot,\n", - " res_ν(x_OR7),\n", + " res_ν(δx_OR7),\n", " supp = σ_dop_all,\n", " curveconf = \"u 1:3:4 w errorbars pt 4 lw $lw_ths_plt\",\n", " lc = 2\n", @@ -2055,223 +1105,6 @@ "#run(`lualatex deldop_residuals_dates.tex`)" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# $\\Gamma$ matrix eigenvalues/vectors" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "Hermitian(Γ_OR7) - Γ_OR7" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "norm(Hermitian(Γ_OR7) - Γ_OR7)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "Γ_OR7_cholf = cholesky(Hermitian(Γ_OR7))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "Γ_OR7_cU = Γ_OR7_cholf.U\n", - "Γ_OR7_cL = Γ_OR7_cholf.L" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Compute nominal initial condition in au, au/day, au/day^2\n", - "q0_OR7 = [\n", - " -0.18034747703273316,\n", - " 0.9406910666200128,\n", - " 0.3457360259054398,\n", - " -0.016265942170279046,\n", - " 4.392889725556651e-5,\n", - " -0.00039519931615139716,\n", - " 0.0, 0.0\n", - "] .+ vcat(1e-8ones(6), 1e-14, 1e-13).*x_OR7\n", - "@show q0_OR7\n", - ";" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# relative scalings analogous to Milani (2005)\n", - "scaling_milani_2005 = vcat( (1/norm(q0_OR7[1:3]))*ones(3), (1/norm(q0_OR7[4:6]))*ones(3), 1/norm(q0_OR7[7]), 1/norm(q0_OR7[8]) )\n", - "\n", - "# jet transport propagation scaling\n", - "jt_scaling = vcat(1e-8ones(6), 1e-14, 1e-13)\n", - ";" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# compute cov matrix with relative scalings\n", - "Γ_OR7_scm05 = diagm(scaling_milani_2005) * ( diagm(jt_scaling)*Γ_OR7*diagm(jt_scaling) ) * diagm(scaling_milani_2005)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# compute normal matrix with relative scalings\n", - "C_OR7_scm05 = inv(Γ_OR7_scm05)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# compute eigenvalues, eigenvectors of normal matrix with relative scalings\n", - "C_eig_OR7_scm05 = eigen(C_OR7_scm05)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# k₁ * V₁ = 1/√λ₁ * V₁\n", - "k1_V1 = C_eig_OR7_scm05.vectors[:,1]/sqrt(C_eig_OR7_scm05.values[1])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# revert eigenvector of \"weak\" direction to JT scaling\n", - "( k1_V1 ./ scaling_milani_2005 ) ./ jt_scaling" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "Taylor1(10) * ( k1_V1 ./ scaling_milani_2005 )" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "@show ( ( k1_V1 ./ scaling_milani_2005 ) ) #./ jt_scaling" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "σ_OR7" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "@show x_OR7 .* jt_scaling" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "#@show Γ_OR7" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, { "cell_type": "code", "execution_count": null, @@ -2282,15 +1115,15 @@ ], "metadata": { "kernelspec": { - "display_name": "Julia 1.6.2", + "display_name": "Julia-4-threads 1.9.1", "language": "julia", - "name": "julia-1.6" + "name": "julia-4-threads-1.9" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.6.2" + "version": "1.9.2" } }, "nbformat": 4, diff --git a/pha/Project.toml b/pha/Project.toml new file mode 100644 index 00000000..ad9456c0 --- /dev/null +++ b/pha/Project.toml @@ -0,0 +1,19 @@ +[deps] + +ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" +DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" +Gaston = "4b11ee91-296f-5714-9832-002c20994614" +JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +NEOs = "b41c07a2-2abb-11e9-070a-c3c1b239e7df" +Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +PlanetaryEphemeris = "d83715d0-7e5f-11e9-1a59-4137b20d8363" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" +TaylorIntegration = "92b13dbe-c966-51a2-8445-caca9f8a7d42" + +[compat] + +NEOs = "0.7" \ No newline at end of file diff --git a/pha/apophis.jl b/pha/apophis.jl index 6faae846..fdca5417 100644 --- a/pha/apophis.jl +++ b/pha/apophis.jl @@ -1,4 +1,32 @@ -using ArgParse, NEOs, PlanetaryEphemeris, Dates, TaylorIntegration, JLD2 +# This file is part of the NEOs.jl package; MIT licensed + +### This script can be run either as a standalone script or via ArgParse. In any case, +### this folder's Project.toml environment has to be already active and instantiated +### before running this script. Uncomment the following three lines to activate and +### instantiate this folder's environment: + +# import Pkg +# Pkg.activate(".") +# Pkg.instantiate() + +### If ran as a standalone and the environment defined in this folder's Project.toml is +### already active and instantiated, then this script can be run from this folder with the +### default settings simply as: +### $ julia -t --project=. apophis.jl +### Finally, this script can be run via the ArgParse mechanism. Help can be displayed doing: +### $ julia --project=. apophis.jl --help + +using ArgParse +using NEOs +using Dates +using TaylorIntegration +using JLD2 +using PlanetaryEphemeris +using DataFrames +using DelimitedFiles +using LinearAlgebra: diag +using Statistics +using StatsBase # Load JPL ephemeris loadjpleph() @@ -47,13 +75,13 @@ function parse_commandline() end s.epilog = """ - examples:\n + Examples (run from the `pha` folder):\n \n # Multi-threaded\n - julia -t 4 --project apophis.jl --maxsteps 100 --nyears_bwd -0.02 --nyears_fwd 0.02 --parse_eqs true\n + julia -t 4 --project=. apophis.jl --maxsteps 100 --nyears_bwd -0.02 --nyears_fwd 0.02 --parse_eqs true\n \n # Single-threaded\n - julia --project apophis.jl --maxsteps 100 --nyears_bwd -0.02 --nyears_fwd 0.02 --parse_eqs true\n + julia --project=. apophis.jl --maxsteps 100 --nyears_bwd -0.02 --nyears_fwd 0.02 --parse_eqs true\n \n """ @@ -79,19 +107,20 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T, q00 = kmsec2auday(apophisposvel197(datetime2et(jd0_datetime))) # Perturbation to nominal initial condition (Taylor1 jet transport) - # vcat(fill(1e-8, 6), 1e-14, 1e-13) are the scaling factors for jet transport perturbation, + # vcat(fill(1e-8, 6), 1e-14, 1e-15) are the scaling factors for jet transport perturbation, # these are needed to ensure expansion coefficients remain small. # The magnitudes correspond to the typical order of magnitude of errors in - # position/velocity (1e-8), Yarkovsky (1e-13) and radiation pressure (1e-14) + # position/velocity (1e-8), Yarkovsky (1e-14) and radiation pressure (1e-15) + scalings = vcat(fill(1e-8, 6), 1e-14, 1e-15) if varorder == 0 dq = zeros(8) else - dq = NEOs.scaled_variables("δx", vcat(fill(1e-8, 6), 1e-14), order = varorder) + dq = NEOs.scaled_variables("δx", scalings, order = varorder) end - q0 = vcat(q00, 0.0, 0.0) .+ vcat(dq, zero(dq[1])) + q0 = vcat(q00, 0.0, 0.0) .+ dq - # Initial date (in julian days) + # Initial date (in Julian days) jd0 = datetime2julian(jd0_datetime) print_header("Integrator warmup", 2) @@ -105,19 +134,25 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T, println("• Final time of integration: ", julian2datetime(jd0 + tmax)) sol_bwd = NEOs.propagate(dynamics, maxsteps, jd0, nyears_bwd, q0, Val(true); - order = order, abstol = abstol, parse_eqs = parse_eqs) - PE.save2jld2andcheck("Apophis_bwd.jld2", (asteph = sol_bwd,)) + order, abstol, parse_eqs) + jldsave("Apophis_bwd.jld2"; sol_bwd) + # sol_bwd = JLD2.load("Apophis_bwd.jld2", "asteph") tmax = nyears_fwd*yr println("• Initial time of integration: ", string(jd0_datetime)) println("• Final time of integration: ", julian2datetime(jd0 + tmax)) sol_fwd = NEOs.propagate(dynamics, maxsteps, jd0, nyears_fwd, q0, Val(true); - order = order, abstol = abstol, parse_eqs = parse_eqs) - PE.save2jld2andcheck("Apophis_fwd.jld2", (asteph = sol_fwd,)) - + order, abstol, parse_eqs) + jldsave("Apophis_fwd.jld2"; sol_fwd) + # sol_fwd = JLD2.load("Apophis_fwd.jld2", "asteph") println() + # load Solar System ephemeris + sseph::TaylorInterpolant{Float64,Float64,2} = loadpeeph(NEOs.sseph, sol_bwd.t0+sol_bwd.t[end], sol_fwd.t0+sol_fwd.t[end]) + eph_su::TaylorInterpolant{Float64,Float64,2} = selecteph(sseph, su) + eph_ea::TaylorInterpolant{Float64,Float64,2} = selecteph(sseph, ea) + # NEO # Change t, x, v units, resp., from days, au, au/day to sec, km, km/sec xva_bwd(et) = auday2kmsec(sol_bwd(et/daysec)[1:6]) @@ -125,39 +160,107 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T, xva(et) = bwdfwdeph(et, sol_bwd, sol_fwd) # Earth # Change x, v units, resp., from au, au/day to km, km/sec - eph_ea = selecteph(NEOs.sseph, ea) xve(et) = auday2kmsec(eph_ea(et/daysec)) # Sun # Change x, v units, resp., from au, au/day to km, km/sec - eph_su = selecteph(NEOs.sseph, su) xvs(et) = auday2kmsec(eph_su(et/daysec)) - radec_2004_2020 = read_radec_mpc("/Users/Jorge/projects/NEOs/data/99942_2004_2020.dat") - radec_2020_2021 = read_radec_mpc("/Users/Jorge/projects/NEOs/data/99942_2020_2021.dat") + radec_2004_2020 = read_radec_mpc(joinpath(pkgdir(NEOs), "data", "99942_2004_2020.dat")) + radec_2020_2021 = read_radec_mpc(joinpath(pkgdir(NEOs), "data", "99942_2020_2021.dat")) radec = vcat(radec_2004_2020,radec_2020_2021) - deldop_2005_2013 = read_radar_jpl("/Users/Jorge/projects/NEOs/data/99942_RADAR_2005_2013.dat") - deldop_2021 = read_radar_jpl("/Users/Jorge/projects/NEOs/data/99942_RADAR_2021.dat") - # TODO: make radar astrometry residuals work with functions!!! - # deldop = vcat(deldop_2005_2013,deldop_2021) + deldop_2005_2013 = read_radar_jpl(joinpath(pkgdir(NEOs), "data", "99942_RADAR_2005_2013.dat")) + deldop_2021 = read_radar_jpl(joinpath(pkgdir(NEOs), "data", "99942_RADAR_2021.dat")) + deldop = vcat(deldop_2005_2013,deldop_2021) # Compute optical residuals - print_header("Compute residuals", 2) - res_radec, w_radec = residuals(radec, xvs = xvs, xve = xve, xva = xva) + res_radec_all, w_radec_all = NEOs.residuals(radec; xvs, xve, xva) + jldsave("Apophis_res_w_radec.jld2"; res_radec_all, w_radec_all) + # JLD2.@load "Apophis_res_w_radec.jld2" # Compute radar residuals - #### TODO: make radar astrometry residuals work with functions!!! - #### res_deldop, w_deldop = residuals(deldop, xvs = xvs, xve = xve, xva = xva_bwd) - res_del_bwd, w_del_bwd, res_dop_bwd, w_dop_bwd = residuals(deldop_2005_2013, xvs = xvs, xve = xve, xva = xva_bwd, niter=5) - res_del_fwd, w_del_fwd, res_dop_fwd, w_dop_fwd = residuals(deldop_2021, xvs = xvs, xve = xve, xva = xva_fwd, niter=5) - - PE.save2jld2andcheck("resw_radec.jld2", (;res_radec,w_radec,res_del_bwd,w_del_bwd,res_dop_bwd,w_dop_bwd,res_del_fwd,w_del_fwd,res_dop_fwd,w_dop_fwd)) - - @show NEOs.cte(res) NEOs.cte(w) - - nothing - + res_del, w_del, res_dop, w_dop = NEOs.residuals(deldop; xvs, xve, xva, niter=10, tord=10) + jldsave("Apophis_res_w_radec.jld2"; res_del, w_del, res_dop, w_dop) + # JLD2.@load "Apophis_res_w_deldop.jld2" + + ### Process optical astrometry (filter, weight, debias) + + # filter out biased observations from observatory 217 on 28-Jan-2021 + df_radec = DataFrame(radec) + # add residuals and weights to optical astrometry DataFrame + df_radec[!, :res_α] .= res_radec_all[1:round(Int,length(res_radec_all)/2)] + df_radec[!, :res_δ] .= res_radec_all[1+round(Int,length(res_radec_all)/2):end] + df_radec[!, :w_α] .= w_radec_all[1:round(Int,length(res_radec_all)/2)] + df_radec[!, :w_δ] .= w_radec_all[1+round(Int,length(res_radec_all)/2):end] + filter!( + x->(Date(x.date) != Date(2021, 1, 28)), + df_radec + ) + + # read astrometric errors from Tholen et al. (2013) + tho13_errors = readdlm(joinpath(pkgdir(NEOs), "data", "tholenetal2013_opterror.dat"), ',') + # compute weights + w_α_tho13 = 1 ./ (tho13_errors[:,1].^2 .+ tho13_errors[:,3].^2 .+ tho13_errors[:,5].^2) + w_δ_tho13 = 1 ./ (tho13_errors[:,2].^2 .+ tho13_errors[:,4].^2 .+ tho13_errors[:,6].^2) + # Tholen et al. (2013) obs table + radec_tho13 = DataFrame(read_radec_mpc(joinpath(pkgdir(NEOs), "test", "data", "99942_Tholen_etal_2013.dat"))) + # vector of RA values from Tholen et al. (2013) observations (used for filtering) + tho13_α = radec_tho13[!,:α] + # set weights in Tholen et al. (2013) astrometry corresponding to associated uncertainties + df_radec[in.(df_radec.α, Ref(tho13_α)),:w_α] = w_α_tho13 + df_radec[in.(df_radec.α, Ref(tho13_α)),:w_δ] = w_δ_tho13 + + # Relaxation factor (account for correlations in optical astrometry data) + # for each observation batch, count the number of observations made in + # the same night by the same observatory + # Ref: Veres et al. (2017) + date_site_v = select(df_radec, :date => ByRow(Date), :observatory) + Nv = [count(x->x.date_Date==i.date_Date && x.observatory==i.observatory, eachrow(date_site_v)) for i in eachrow(date_site_v)] + relax_factor = map(x->x>4.0 ? x/4.0 : 1.0, Nv) + # inflate uncertainties (i.e., relax weights) by relaxation factor + df_radec[!, :w_α] .= (df_radec.w_α)./relax_factor + df_radec[!, :w_δ] .= (df_radec.w_δ)./relax_factor + + # update optical residuals and weights + res_radec = vcat(df_radec.res_α, df_radec.res_δ) + w_radec = vcat(df_radec.w_α, df_radec.w_δ) + + ### Perform orbital fit to optical and radar astrometry data + + res = vcat(res_radec, res_del, res_dop) + w = vcat(w_radec, w_del, w_dop) + + success, δx_OR8, Γ_OR8 = newtonls(res, w, zeros(get_numvars()), 10) + x_OR8 = sol_fwd(sol_fwd.t0)(δx_OR8) + σ_OR8 = sqrt.(diag(Γ_OR8)).*scalings + + nradec = length(res_radec) + res_ra = view(res_radec, 1:nradec÷2) + res_dec = view(res_radec, 1+nradec÷2:nradec) + w_ra = view(w_radec, 1:nradec÷2) + w_dec = view(w_radec, 1+nradec÷2:nradec) + + ### Print results + + print_header("Orbital fit (8-DOF) and post-fit statistics", 2) + + # orbital fit + println("Success flag : ", success, "\n") + println("Nominal solution [au,au,au,au/d,au/d,au/d,au/d²,au/d²]: ", x_OR8, "\n") + println("1-sigma formal uncertainties [au,au,au,au/d,au/d,au/d,au/d²,au/d²]: ", σ_OR8, "\n") + + # post-fit statistics + println("Normalized RMS (optical-only) [adimensional] : ", nrms(res_radec(δx_OR8),w_radec)) + println("Normalized RMS (radar-only) [adimensional] : ", nrms(vcat(res_del,res_dop)(δx_OR8),vcat(w_del,w_dop))) + println("Normalized RMS (combined optical and radar) [adimensional] : ", nrms(res(δx_OR8),w), "\n") + println("Mean weighted right-ascension residual [arcseconds] : ", mean(res_ra(δx_OR8), weights(w_ra))) + println("Mean weighted declination residual [arcseconds] : ", mean(res_dec(δx_OR8), weights(w_dec))) + println("Mean weighted time-delay residual [micro-seconds]: ", mean(res_del(δx_OR8), weights(w_del))) + println("Mean weighted Doppler-shift residual [Hz] : ", mean(res_dop(δx_OR8), weights(w_dop)), "\n") + println("Chi-squared statistic (χ²): [adimensional] : ", chi2(res(δx_OR8),w)) + + return sol_bwd, sol_fwd, res_radec, res_del, res_dop, w_radec, w_del, w_dop end function main() diff --git a/src/NEOs.jl b/src/NEOs.jl index e907d3f0..bf969524 100644 --- a/src/NEOs.jl +++ b/src/NEOs.jl @@ -61,7 +61,7 @@ export gauss_method # Asteroid dynamical models export RNp1BP_pN_A_J23E_J2S_ng_eph_threads!, RNp1BP_pN_A_J23E_J2S_eph_threads! # Propagate -export propagate, propagate_lyap, propagate_root, save2jldandcheck +export propagate, propagate_lyap, propagate_root export valsecchi_circle, nrms, chi2, newtonls, newtonls_6v, diffcorr, newtonls_Q, bopik diff --git a/src/observations/catalogue_mpc.jl b/src/observations/catalogue_mpc.jl index 339e8685..6a240432 100644 --- a/src/observations/catalogue_mpc.jl +++ b/src/observations/catalogue_mpc.jl @@ -155,29 +155,18 @@ function write_catalogues_mpc(cats::Vector{CatalogueMPC}, filename::String) end end -const downloader = Downloads.Downloader() -downloader.easy_hook = (easy, info) -> Downloads.Curl.setopt(easy, Downloads.Curl.CURLOPT_LOW_SPEED_TIME, 60) - @doc raw""" - download_scratch(url::String, filename::String) + download_scratch(url::String, filename::String; connect_timeout=180, readtimeout=180) Download `url` and save the output to NEOs scratch space as `filename`. Return the local path and the contents of the file as a `String`. """ -function download_scratch(url::String, filename::String) +function download_scratch(url::String, filename::String; connect_timeout=180, readtimeout=180) # Local file path = joinpath(scratch_path[], filename) - # Download function - f = retry(delays = fill(0.1, 3)) do - try - download(url, path; downloader) - catch - rethrow() - end - end - # Download source_file - f() + # Get raw html (HTTP.get retries four times by default) + resp = get(url; connect_timeout, readtimeout) # Read local file - txt = read(path, String) + txt = String(resp.body) return path, txt end diff --git a/src/observations/jpl_eph.jl b/src/observations/jpl_eph.jl index 4d4f9920..bba0c18f 100644 --- a/src/observations/jpl_eph.jl +++ b/src/observations/jpl_eph.jl @@ -125,11 +125,11 @@ const ttmtdb::TaylorInterpolant{Float64, Float64, 1} = TaylorInterpolant(sseph.t @doc raw""" loadpeeph(eph::TaylorInterpolant{Float64, Float64, 2} = sseph, t_0::Real = 0.0, t_f::Real = 36525.0) -Load ephemeris produced by `PlanetaryEphemeris.jl` in timerange `[t_0, t_f] ⊆ [0.0, 36525.0]` where `t` must have units of TDB +Load ephemeris produced by `PlanetaryEphemeris.jl` in timerange `[t_0, t_f] ⊆ [0.0, 36525.0]` where `t` must have units of TDB days since J2000. The available options for `eph` are: -- `NEOs.sseph`: Solar system ephemeris. -- `NEOs.acceph`: accelerations ephemeris. +- `NEOs.sseph`: Solar system ephemeris. +- `NEOs.acceph`: accelerations ephemeris. - `NEOs.poteph`: newtonian potentials ephemeris. !!! warning @@ -143,12 +143,12 @@ function loadpeeph(eph::TaylorInterpolant{Float64, Float64, 2} = sseph, t_0::Rea end @doc raw""" - bwdfwdeph(et::Union{T, TaylorN{T}}, bwd::TaylorInterpolant{T, U, 2}, + bwdfwdeph(et::Union{T,Taylor1{T},TaylorN{T},Taylor1{TaylorN{T}}}, bwd::TaylorInterpolant{T, U, 2}, fwd::TaylorInterpolant{T, U, 2}) where {T <: AbstractFloat, U <: Union{T, TaylorN{T}}} Paste a backward and a forward integration, evaluate at `et` and convert from [au, au/day] -> [km, km/sec]. """ -function bwdfwdeph(et::Union{T,TaylorN{T}}, +function bwdfwdeph(et::Union{T,Taylor1{T},TaylorN{T},Taylor1{TaylorN{T}}}, bwd::TaylorInterpolant{T,U,2}, fwd::TaylorInterpolant{T,U,2} ) where {T<:AbstractFloat, U<:Union{T,TaylorN{T}}} diff --git a/src/observations/topocentric.jl b/src/observations/topocentric.jl index 66d94a01..b9d5bb3e 100644 --- a/src/observations/topocentric.jl +++ b/src/observations/topocentric.jl @@ -39,6 +39,9 @@ end obs_pos_ECEF(x::RadecMPC{T}) where {T <: AbstractFloat} = obs_pos_ECEF(x.observatory) obs_pos_ECEF(x::RadarJPL{T}) where {T <: AbstractFloat} = obs_pos_ECEF(x.rcvr) +# TODO: avoid sv_ecef_to_ecef overload by defining proper product between DCMs and Taylor1/TaylorN +# method below has been adapted from SatelliteToolboxTransformations.jl, MIT-licensed +# https://github.com/JuliaSpace/SatelliteToolboxTransformations.jl function sv_ecef_to_ecef( sv::OrbitStateVector, T_ECEF1::Val{:ITRF}, @@ -56,6 +59,9 @@ function sv_ecef_to_ecef( return OrbitStateVector(sv.t, r_ecef, v_ecef, a_ecef) end +# TODO: avoid sv_ecef_to_eci overload by defining proper product between DCMs and Taylor1/TaylorN +# method below has been adapted from SatelliteToolboxTransformations.jl, MIT-licensed +# https://github.com/JuliaSpace/SatelliteToolboxTransformations.jl function sv_ecef_to_eci( sv::OrbitStateVector, T_ECEF::Union{Val{:PEF}, Val{:TIRS}}, diff --git a/src/orbit_determination/osculating.jl b/src/orbit_determination/osculating.jl index 3e86c7c3..ac3c370b 100644 --- a/src/orbit_determination/osculating.jl +++ b/src/orbit_determination/osculating.jl @@ -120,7 +120,7 @@ end Return cartesian state vector of orbit `osc` at time `t` (Julian day). """ -function (osc::OsculatingElements{T})(t::T) where {T <: Number} +function (osc::OsculatingElements{T})(t::U) where {T,U <: Number} # Mean motion n = PE.meanmotion(μ_S, osc.a) @@ -176,6 +176,6 @@ See https://doi.org/10.1016/j.icarus.2013.02.004. - `e`: eccentricity. - `μ_S`: mass parameter of the Sun. """ -function yarkp2adot(A2, a, e, μ_S) - return 2A2/(sqrt(a)*(1-e^2)*sqrt(μ_S)) +function yarkp2adot(A2, a, e; μ = μ_S) + return 2A2/(sqrt(a)*(1-e^2)*sqrt(μ)) end diff --git a/src/propagation/propagation.jl b/src/propagation/propagation.jl index fecda942..c44c52cb 100644 --- a/src/propagation/propagation.jl +++ b/src/propagation/propagation.jl @@ -29,98 +29,6 @@ function rvelea(dx, x, params, t) return true, (x[1]-xe[1])*(x[4]-xe[4]) + (x[2]-xe[2])*(x[5]-xe[5]) + (x[3]-xe[3])*(x[6]-xe[6]) end -@doc raw""" - save2jldandcheck(objname, sol) - -Save `sol` in a file `objname_jt.jld2`. - -See also [`__save2jldandcheck`](@ref). -""" -function save2jldandcheck(objname, sol) - # Name of the file - outfilename = string(objname, "_jt.jld2") - # Save sol in outfilename - return __save2jldandcheck(outfilename, sol) -end - -@doc raw""" - __save2jldandcheck(outfilename, sol) - -Savs `sol` in `outfilename`. -""" -function __save2jldandcheck(outfilename, sol) - println("Saving solution to file: $outfilename") - # Open file - JLD2.jldopen(outfilename, "w") do file - # Loop over solution variables - for ind in eachindex(sol) - # Name of the variable - varname = string(ind) - println("Saving variable: ", varname) - # Write the varaible - write(file, varname, sol[ind]) - end - end - # Check that saved solution is equal to the original - println("Checking that all variables were saved correctly...") - # Loop over solution variables - for ind in eachindex(sol) - # Name of the variable - varname = string(ind) - # Read varname from files and assign recovered variable to recovered_sol_i - recovered_sol_i = JLD2.load(outfilename, varname) - # Check that varname was recovered succesfully - @show recovered_sol_i == sol[ind] - end - println("Saved solution") - return outfilename -end - -@doc raw""" - taylor_minimum(pol::Taylor1{T}, x0::T; niters::Int=10) where {T<:Real} - -Return the minimum of the Taylor polynomial `pol` computed via Newton's method. `x0` is the -initial guess and `niters` is the number of iterations. -""" -function taylor_minimum(pol::Taylor1{T}, x0::T; niters::Int=10) where {T<:Real} - # First derivative - dpol = PE.ordpres_differentiate(pol) - # Second derivative - dpol2 = PE.ordpres_differentiate(dpol) - # Initial guess - xnewton::T = x0 - #@show xnewton - # Newton iteration - for i in 1:niters - # Newton update rule - xnewton -= dpol(xnewton)/dpol2(xnewton) - #@show xnewton, dpol(xnewton) - end - - return xnewton -end - -@doc raw""" - taylor_roots(pol::Taylor1{T}, x0::T; niters::Int=10) where {T<:Real} - -Return the root of the Taylor polynomial `pol` computed via Newton's method. `x0` is the -initial guess and `niters` is the number of iterations. -""" -function taylor_roots(pol::Taylor1{T}, x0::T; niters::Int=10) where {T<:Real} - # First derivative - dpol = PE.ordpres_differentiate(pol) - # Initial guess - xnewton::T = x0 - #@show xnewton - # Newton iteration - for i in 1:niters - # Newton update rule - xnewton -= pol(xnewton)/dpol(xnewton) - #@show xnewton, pol(xnewton) - end - return xnewton -end - @doc raw""" scaled_variables(c::Vector{T} = fill(1e-6, 6), names::String = "δx"; order::Int = 5) where {T <: Real} diff --git a/test/osculating.jl b/test/osculating.jl index db0a3d25..6cf5e078 100644 --- a/test/osculating.jl +++ b/test/osculating.jl @@ -12,5 +12,5 @@ using LinearAlgebra: norm qa = [-1.0506628055913627, -0.06064314196134998, -0.04997102228887035, 0.0029591421121582077, -0.01423233538611057, -0.005218412537773594] jda = datetime2julian(DateTime(2004,6,1)) kepa = pv2kep(qa) - @test norm(kep(jd2000) - q, Inf) < 1e-10 + @test norm(kepa(jd2000) - qa, Inf) < 1e-10 end \ No newline at end of file