# python/asteria/asteria_api.py
from dataclasses import dataclass
from importlib.resources import files
import requests
import math
import numpy as np
import matplotlib.pyplot as plt
[docs]
@dataclass
class Simulation:
"""
Monte Carlo simulation of Yarkovsky effect thermal parameters for a
near-Earth asteroid.
The class wraps the ASTERIA Fortran core (``_asteria``) and exposes a
Python-friendly interface for setting up, running, and reading the results
of a thermal-parameter inversion based on a measured semi-major axis
drift (da/dt).
The simulation proceeds in two stages driven by :meth:`run`:
1. **Distribution generation** (``gen_distrib``) — draws ``sample_size``
realisations of all physical and orbital parameters from the supplied
prior distributions.
2. **Monte Carlo inversion** (``gamma_est_mc``) — iterates up to
``max_iter`` times to find thermal conductivity *K* and thermal
inertia *Gamma* solutions consistent with the observed da/dt.
Physical parameters (albedo, density, diameter, surface density,
obliquity) default to ``-1.0``, which is a **sentinel value** meaning
*"draw from the built-in NEA population model"*. Set both the mean and
standard deviation to a positive value to use a Gaussian prior instead.
Parameters
----------
sample_size : int, optional
Number of samples drawn for the a-priori distributions. Default 10 000.
sma : float
Semi-major axis in au. Must be > 0.
ecc : float
Eccentricity. Must satisfy 0 <= ecc < 1.
inc : float
Inclination in degrees. Must satisfy 0 <= inc <= 180.
mean_H : float
Mean absolute magnitude H. Must be >= 0.
std_H : float, optional
1-sigma uncertainty on H. Default 0.5 mag.
mean_pv : float
Mean geometric albedo pv. Set to ``-1.0`` to use the population
model. Default ``-1.0``.
std_pv : float
1-sigma uncertainty on pv. Ignored when ``mean_pv = -1.0``. Default ``-1.0``.
mean_rho : float
Mean bulk density in kg/m^3. Set to ``-1.0`` for the population model.
Default ``-1.0``.
std_rho : float
1-sigma uncertainty on bulk density in kg/m^3. Default ``-1.0``.
mean_D : float
Mean diameter in metres. Set to ``-1.0`` to derive from H and pv via
the population model. Default ``-1.0``.
std_D : float
1-sigma uncertainty on diameter in metres. Default ``-1.0``.
mean_rhos : float
Mean surface (regolith) density in kg/m^3. Set to ``-1.0`` for the
population model. Default ``-1.0``.
std_rhos : float
1-sigma uncertainty on surface density in kg/m^3. Default ``-1.0``.
mean_gam : float
Mean spin-axis obliquity in degrees. Set to ``-1.0`` for the
population model. Default ``-1.0``.
std_gam : float
1-sigma uncertainty on obliquity in degrees. Default ``-1.0``.
mean_P : float
Mean rotation period in hours. Must be >= 0. Default ``None``.
std_P : float
1-sigma uncertainty on rotation period in hours. Must be >= 0.
Default ``None``.
a_to_b : float
Axis ratio a/b for the non-spherical correction (0 < a_to_b <= 1).
Required when calling ``run(ns_corr=1)``. Default ``None``.
delta_m : float
Lightcurve peak-to-peak amplitude in magnitudes (>= 0). Used for the
non-spherical correction when calling ``run(ns_corr=2)``.
Default ``None``.
mean_dadt : float
Measured semi-major axis drift da/dt in au/My. Default ``None``.
std_dadt : float
1-sigma uncertainty on da/dt in au/My. Must be >= 0. Default ``None``.
alpha : float
Solar radiation absorption coefficient. Default 1.0.
C : float
Specific heat capacity in J/kg/K. Default 800.0.
thermalCondMin : float
Lower bound of the thermal conductivity search range in W/m/K.
Must be > 0. Default 1e-5.
thermalCondMax : float
Upper bound of the thermal conductivity search range in W/m/K.
Must be > 0. Default 10.0.
emissiv : float
Bolometric thermal emissivity. Default 0.984.
method : str
Yarkovsky model variant. One of:
* ``'circular'`` — circular-orbit approximation (fastest).
* ``'elliptic'`` — full elliptic orbit, single-layer model.
* ``'elliptic_2l'`` — full elliptic orbit, two-layer thermal model.
Default ``'circular'``.
filename : str
Base name (without extension) for the output text file. Default
``'asteria_output'``.
max_iter : int
Maximum number of Monte Carlo iterations. Must be > 0. Default 100 000.
expo : float, optional
Exponent governing how thermal inertia scales with heliocentric
distance (Gamma proportional to r^expo). Default 1.0.
n_proc : int, optional
Number of parallel CPU cores used by the Fortran core. Default 1.
"""
# Default sample size for gen_distrib
sample_size: int = 10000
# Set null values for orbital parameters
sma: float = None
ecc: float = None
inc: float = None
# Set null values for absolute magnitude
mean_H: float = None
std_H: float = 0.5
# Set population model as a default for physical parameters
mean_pv: float = -1.0
std_pv: float = -1.0
mean_rho: float = -1.0
std_rho: float = -1.0
mean_D: float = -1.0
std_D: float = -1.0
mean_rhos: float = -1.0
std_rhos: float = -1.0
mean_gam: float = -1.0
std_gam: float = -1.0
# Set null values rotation period
mean_P: float = None
std_P: float = None
# Non-spherical effect effects
a_to_b: float = None
delta_m: float = None
# Set null values for Yarkovsky effect
mean_dadt: float = None
std_dadt: float = None
# Set alpha to 1 as default
alpha: float = 1.0
# Heat capacity, default to 800
C: float = 800.0
# Set default to mean value
emissiv: float = 0.984
# Thermal conductivity limits, set defaults
thermalCondMin: float = 1.e-5
thermalCondMax: float = 10.0
# Set default Yarkovsky method to 'circular'
method: str = 'circular'
# Output filename
filename: str = 'asteria_output'
# Set maximum iteration default
max_iter: int = 100000
# Exponent of Gamma as function of distance from Sun, default to 1
expo: float = 1.0
# Set number of processors, default to 1
n_proc: int = 1
def __post_init__(self):
if self.sma is not None:
if self.sma <= 0:
raise ValueError('Semi-major axis must be positive')
if self.ecc is not None:
if self.ecc < 0 or self.ecc >= 1:
raise ValueError('Eccentricity must be between 0 and 1')
if self.inc is not None:
if self.inc < 0 or self.ecc >= 180:
raise ValueError('Inclination must be between 0 and 180')
if self.sample_size <= 0:
raise ValueError("sample_size must be positive")
if self.max_iter <= 0:
raise ValueError("max_iter must be positive")
if self.thermalCondMin <= 0:
raise ValueError("thermalCondMin must be positive")
if self.thermalCondMax <= 0:
raise ValueError("thermalCondMax must be positive")
if self.mean_H is not None:
if self.mean_H < 0:
raise ValueError("mean_H must be non-negative")
if self.std_H is not None:
if self.std_H < 0:
raise ValueError("std_H must be non-negative")
if self.mean_P is not None:
if self.mean_P < 0:
raise ValueError("mean_P must be non-negative")
if self.std_P is not None:
if self.std_P is not None:
if self.std_P < 0:
raise ValueError("std_P must be non-negative")
if self.std_dadt is not None:
if self.std_dadt < 0:
raise ValueError("std_dadt must be non-negative")
if self.a_to_b is not None:
if self.a_to_b <= 0 or self.a_to_b > 1:
raise ValueError('a_to_b must be >0 and <1')
if self.delta_m is not None:
if self.delta_m < 0:
raise ValueError('Lightcurve amplitude delta_m must be positive')
if self.method != "circular" and self.method != "elliptic" and self.method != "elliptic_2l":
raise ValueError('method not valid')
[docs]
def set_neocc_asteroid(self, astname):
"""
Populate orbital and physical parameters by fetching data from the
ESA NEOCC database.
Downloads the ``.ke0`` Keplerian orbit file for *astname*, parses it
with :func:`read_neocc_orb`, and sets ``sma``, ``ecc``, ``inc``,
``mean_H``, ``mean_dadt``, and ``std_dadt`` on this instance.
Parameters
----------
astname : str or int
Asteroid designation or number as recognised by ESA NEOCC
(e.g. ``"2024YR4"``, ``101955``).
Returns
-------
a : float
Semi-major axis in au.
e : float
Eccentricity.
i : float
Inclination in degrees.
H : float
Absolute magnitude.
dadt : float
Semi-major axis drift da/dt in au/My.
dadt_std : float
1-sigma uncertainty on da/dt in au/My.
Raises
------
ValueError
If the asteroid is not found (HTTP status other than 200).
ValueError
If the orbit file does not contain a measured A2 parameter.
Notes
-----
Requires an active internet connection. The ``.ke0`` file is saved
to the current working directory.
"""
astname = str(astname)
orbfile = astname.strip() + '.ke0'
url = 'https://neo.ssa.esa.int/PSDB-portlet/download?file=' + astname.strip() + '.ke0'
response = requests.get(url)
if response.status_code == 200:
# Write to a local file
with open(orbfile, "wb") as f:
f.write(response.content)
a, e, i, H, dadt, dadt_std = read_neocc_orb(orbfile)
# If da/dt is nan, raise an error
if math.isnan(dadt):
raise ValueError('Semi-major axis drift not present for ' + astname)
# Set asteroid orbital parameters
self.sma = a
self.ecc = e
self.inc = i
self.mean_dadt = dadt
self.std_dadt = dadt_std
self.mean_H = H
return a, e, i, H, dadt, dadt_std
else:
raise ValueError('Asteroid ' + astname.strip() + ' does not exist.')
[docs]
def set_jpl_asteroid(self, astname):
"""
Populate orbital and physical parameters by querying the JPL Small
Body Database (SBDB) API.
Sends a request to ``https://ssd-api.jpl.nasa.gov/sbdb.api``,
extracts orbital elements, absolute magnitude H, and the
non-gravitational parameter A2, then sets ``sma``, ``ecc``, ``inc``,
``mean_H``, ``mean_dadt``, and ``std_dadt`` on this instance.
Parameters
----------
astname : str or int
Asteroid name, designation, or SPK-ID accepted by JPL SBDB
(e.g. ``"2024YR4"``, ``469219``).
Returns
-------
a : float
Semi-major axis in au.
e : float
Eccentricity.
inc : float
Inclination in degrees.
H : float or None
Absolute magnitude, or ``None`` if not in SBDB.
dadt : float
Semi-major axis drift da/dt in au/My.
dadt_std : float
1-sigma uncertainty on da/dt in au/My.
Raises
------
requests.HTTPError
If the SBDB API returns a non-2xx status code.
ValueError
If no orbital data are found for *astname*.
ValueError
If the entry does not contain a measured A2 parameter.
Notes
-----
Requires an active internet connection.
"""
astname = str(astname)
dadt = math.nan
dadt_std = math.nan
url = "https://ssd-api.jpl.nasa.gov/sbdb.api"
params = {
"sstr": astname,
"full-prec": "true",
"phys-par": "true"
}
response = requests.get(url, params=params)
response.raise_for_status()
data = response.json()
if "orbit" not in data:
raise ValueError("No orbital data found in JPL for object " + astname.strip())
orbit = data["orbit"]
# Extract orbital data
a = e = inc = None
for el in orbit.get("elements", []):
if el["name"] == "a":
a = float(el["value"])
elif el["name"] == "e":
e = float(el["value"])
elif el["name"] == "i":
inc = float(el["value"])
# Extract absolute magnitude value
H = None
for phys in data.get("phys_par", []):
if phys.get("name") == "H":
H = float(phys["value"])
break
# ---- A2 parameter ----
A2 = math.nan
A2std = math.nan
for par in orbit.get("model_pars", []):
if par["name"] == "A2":
A2 = float(par["value"])
A2std = float(par["sigma"]) if par.get("sigma") else None
if not math.isnan(A2):
dadt = A2_to_dadt(A2, a, e)
if not math.isnan(A2std):
dadt_std = A2_to_dadt(A2std, a, e)
if math.isnan(dadt):
raise ValueError('Semi-major axis drift not present for ' + astname)
# Set
self.sma = a
self.ecc = e
self.inc = inc
self.mean_dadt = dadt
self.std_dadt = dadt_std
self.mean_H = H
return a, e, inc, H, dadt, dadt_std
[docs]
def print_simulation_parameters(self):
"""
Print a formatted summary of the simulation configuration to stdout.
Displays sample size, thermal conductivity search range, Yarkovsky
model, thermal-inertia distance exponent, output filename, maximum
iteration count, and CPU count.
"""
print("========================")
print(" SIMULATION PARAMETERS ")
print("========================")
print(" Sample size of a priori distributions : ", self.sample_size)
print(" Min. K solutions (W/m/K) : ", self.thermalCondMin)
print(" Max. K solutions (W/m/K) : ", self.thermalCondMax)
print(" Yarkovsky effect model : ", self.method)
print(" TI exponent for eccentric Yarkovsky : ", self.expo)
print(" Output filename : ", self.filename)
print(" Maximum Monte Carlo iterations : ", self.max_iter)
print(" Number of CPUs : ", self.n_proc)
[docs]
def print_asteroid_parameters(self):
"""
Print a formatted summary of the asteroid's orbital and physical
parameters to stdout.
Physical parameters whose mean is still ``-1.0`` are reported as
``"Population"`` to indicate that the built-in NEA population model
will be used during the simulation.
"""
print("======================")
print(" ASTEROID PARAMETERS ")
print("======================")
print("Orbital parameters:")
print(" Semi-major axis (au) : ", self.sma)
print(" Eccentricity : ", self.ecc)
print(" Inclination (deg) : ", self.inc)
print(" da/dt (au/My) : ", self.mean_dadt, " +- ", self.std_dadt)
print("Physical parameters:")
print(" Absolute Magnitude : ", self.mean_H, " +- ", self.std_H)
if self.mean_pv < 0:
print(" Albedo : Population")
else:
print(" Albedo : ", self.mean_pv, " +- ", self.std_pv)
if self.mean_D < 0:
print(" Diameter : Population")
else:
print(" Diameter (m) : ", self.mean_D, " +- ", self.std_D)
if self.mean_rho < 0:
print(" Density : Population")
else:
print(" Density (kg/m^3) : ", self.mean_rho, " +- ", self.std_rho)
if self.mean_gam < 0:
print(" Obliquity : Population")
else:
print(" Obliquity (deg) : ", self.mean_gam, " +- ", self.std_gam)
print(" Rotation Period (h) : ", self.mean_P, " +- ", self.std_P)
if self.delta_m is not None:
print(" Lightcurve amplitude : ", self.delta_m)
if self.a_to_b is not None:
print(" Ratio a/b : ", self.a_to_b)
print(" Heat Capacity (J/kg/K) : ", self.C)
print(" Emissivity : ", self.emissiv)
print(" Absorption Coefficient : ", self.alpha)
[docs]
def run(self, ns_corr=0):
"""
Execute the ASTERIA Monte Carlo thermal-parameter inversion.
Optionally applies a non-spherical shape correction to the observed
da/dt, then calls the two-stage Fortran core: distribution generation
(``gen_distrib``) followed by the Monte Carlo solver
(``gamma_est_mc``). Results are written to ``<filename>.txt``.
Parameters
----------
ns_corr : {0, 1, 2}, optional
Non-spherical shape correction flag:
* ``0`` — No correction (spherical approximation).
* ``1`` — Correct da/dt using the axis ratio ``a_to_b``.
Requires ``a_to_b`` to be set.
* ``2`` — Correct da/dt using the lightcurve amplitude ``delta_m``.
Requires ``delta_m`` to be set.
The scale factor is xi = f^(-0.3), where f is either ``a_to_b``
or exp(delta_m / 2.5). Default is ``0``.
Raises
------
ValueError
If ``ns_corr=1`` and ``a_to_b`` is not set.
ValueError
If ``ns_corr=2`` and ``delta_m`` is not set.
ValueError
If ``ns_corr`` is not 0, 1, or 2.
ValueError
If any required field (``sma``, ``ecc``, ``inc``, ``mean_H``,
``mean_dadt``, ``std_dadt``, ``mean_P``, ``std_P``) is ``None``.
Notes
-----
The Fortran extension ``_asteria`` is imported lazily so the package
can be imported without a compiled binary (e.g. during documentation
builds).
"""
from . import _asteria as core
if ns_corr == 1:
try:
f = self.a_to_b
except:
raise ValueError('Axes ratio a/b not set.')
# Compute csi factor
csi = f**(-0.3)
# Correct semi-major axis drift
mean_dadt = self.mean_dadt/csi
std_dadt = self.std_dadt/csi
elif ns_corr == 2:
try:
f = np.exp(self.delta_m/2.5)
except:
raise ValueError('Lightcurve amplitude not set.')
# Compute csi factor
csi = f**(-0.3)
# Correct semi-major axis drift
mean_dadt = self.mean_dadt/csi
std_dadt = self.std_dadt/csi
elif ns_corr == 0:
mean_dadt = self.mean_dadt
std_dadt = self.std_dadt
else:
raise ValueError('Flag ns_corr not valid, ns_corr = ' + str(ns_corr))
DATA_PATH = files("asteria") / "dat"
gmb_file = DATA_PATH / "gmb_model.dat"
# Check if orbital elements are none
if self.sma is None:
raise ValueError("Semi-major axis value not valid.")
if self.ecc is None:
raise ValueError("Eccentricity value not valid.")
if self.inc is None:
raise ValueError("Inclination value not valid.")
if self.mean_H is None:
raise ValueError("Absolute magnitude value not valid.")
if self.mean_dadt is None:
raise ValueError("Semi-major axis drift value not valid.")
if self.std_dadt is None:
raise ValueError("Semi-major axis drift 1sigma uncertainty value not valid.")
if self.mean_P is None:
raise ValueError("Rotation period value not valid.")
if self.std_P is None:
raise ValueError("Rotation period 1sigma uncertainty value not valid.")
core.gen_distrib_mod.gen_distrib(self.sample_size, self.sma, self.ecc, self.inc,
self.mean_H, self.std_H, self.mean_pv, self.std_pv,
self.mean_rho, self.std_rho, self.mean_D, self.std_D,
self.mean_rhos, self.std_rhos, self.mean_gam, self.std_gam,
self.mean_P, self.std_P, mean_dadt, std_dadt, self.alpha, gmb_file)
if self.method == 'circular':
method = 1
elif self.method == 'elliptic':
method = 2
elif self.method == 'elliptic_2l':
method = 3
else:
raise ValueError("Method not valid")
core.gamma_est_mc_mod.gamma_est_mc(self.C, self.thermalCondMin, self.thermalCondMax, self.sma, self.ecc,
self.emissiv, method, self.filename, self.max_iter, self.expo, self.n_proc)
[docs]
def read_output_data(self):
"""
Read and parse the Monte Carlo output file produced by :meth:`run`.
The file ``<filename>.txt`` contains one accepted solution per line.
Column layout depends on ``method``:
* ``'circular'`` / ``'elliptic'``: K Gamma rho D gamma
* ``'elliptic_2l'``: K Gamma rho <unused> D gamma
Returns
-------
K : list of float
Thermal conductivity solutions in W/m/K.
Gamma : list of float
Thermal inertia solutions in J m^-2 K^-1 s^-0.5 (TIU).
rho : list of float
Bulk density solutions in kg/m^3.
diam : list of float
Diameter solutions in metres.
obli : list of float
Spin-axis obliquity solutions in degrees.
Raises
------
FileNotFoundError
If ``<filename>.txt`` does not exist in the working directory.
"""
filename = self.filename + '.txt'
# Set variables names
K = []
Gamma = []
rho = []
diam = []
obli = []
with open(filename, 'r') as f:
for line in f:
# The output is different in these two cases
if self.method == 'circular' or self.method == 'elliptic':
K.append(float(line.split()[0]))
Gamma.append(float(line.split()[1]))
rho.append(float(line.split()[2]))
diam.append(float(line.split()[3]))
obli.append(float(line.split()[4]))
else:
K.append(float(line.split()[0]))
Gamma.append(float(line.split()[1]))
rho.append(float(line.split()[2]))
diam.append(float(line.split()[4]))
obli.append(float(line.split()[5]))
return K, Gamma, rho, diam, obli
[docs]
def plot_distribution(data, bins=50, xlabel="", ylabel="Count", fontsize=13, figsize=(8, 5), log_bins=False, filename=None):
"""
Plot a normalized histogram of a 1D dataset.
Non-finite values (NaN, Inf) are silently removed before plotting.
When ``log_bins=True``, the histogram is weighted by bin center so that
the plotted quantity approximates a log-space probability density
(i.e. dN/d log x).
Parameters
----------
data : array-like
Input data. Must contain at least one finite value.
bins : int, optional
Number of histogram bins. Default is 50.
xlabel : str, optional
Label for the x-axis. Default is an empty string.
ylabel : str, optional
Label for the y-axis. Default is ``"Count"``.
fontsize : int or float, optional
Font size for axis labels and tick labels. Default is 13.
figsize : tuple of (float, float), optional
Width and height of the figure in inches. Default is ``(8, 5)``.
log_bins : bool, optional
If ``True``, bin edges are spaced logarithmically and the x-axis is
shown on a log scale. All data values must be strictly positive.
Default is ``False``.
filename : str or None, optional
If provided, the figure is saved to this path at 600 dpi before
being displayed. Default is ``None``.
Raises
------
ValueError
If *data* is empty after removing non-finite values.
ValueError
"""
data = np.asarray(data)
data = data[np.isfinite(data)]
if len(data) == 0:
raise ValueError("Input data is empty after removing invalid values.")
if log_bins:
if np.any(data <= 0):
raise ValueError("Log bins require all data to be > 0.")
bins_edges = np.logspace(np.log10(data.min()),np.log10(data.max()),bins)
else:
bins_edges = bins
plt.figure(figsize=figsize)
if log_bins:
hist, edges = np.histogram(data, bins=bins_edges, density=True)
centers = np.sqrt(edges[:-1] * edges[1:])
log_pdf = hist * centers
plt.step(centers, log_pdf, where="mid", linewidth=2)
plt.xscale("log")
else:
plt.hist(data, bins=bins_edges, histtype="step", linewidth=2, density=True)
plt.xlabel(xlabel, fontsize=fontsize)
plt.ylabel(ylabel, fontsize=fontsize)
plt.xticks(fontsize=fontsize)
plt.yticks(fontsize=fontsize)
plt.tight_layout()
if filename is not None:
plt.savefig(filename, dpi=600)
plt.show()
def A2_to_dadt(A2, a, e):
"""
Convert the Yarkovsky non-gravitational parameter A2 to a semi-major
axis drift rate da/dt in au/Myr.
The conversion uses the linearised secular Yarkovsky formula:
.. math::
\\frac{da}{dt} = \\frac{2 A_2 (1 - e^2)}{p^2 \\, n}
where *p = a(1 - e^2)* is the semi-latus rectum and *n* is the mean
motion. The result is expressed in au/My (astronomical units per
million years).
Parameters
----------
A2 : float
Non-gravitational transverse acceleration parameter in au/d^2.
a : float
Osculating semi-major axis in au.
e : float
Osculating eccentricity (0 <= e < 1).
Returns
-------
dadt : float
Semi-major axis drift rate in au/My.
Notes
-----
The Sun's gravitational constant used internally is
GM = 1.32712440018e20 m^3/s^2.
"""
# Sun gravitational constant in m an s
GM = 1.32712440018e20
# Conversion factors
au2m = 149597870700.0
s2d = 1.0/86400.0
# Compute the mean motion
nstd = math.sqrt(GM/(a*au2m)**3)/s2d
# Convert the value of A2
p = a*(1.0-e**2)
d = 2.0
d2y = 1.0/365.25
y2my = 1.0/10.0**6
d2my = d2y*y2my
# The conversion formula is
# dadt = 2 A2 (1-e^2)*(1/(a(1-e^2)))^2/n
dadt = 2.0*A2*(1.0-e**2)/p**d/nstd/d2my
return dadt
def read_neocc_orb(orbfile):
"""
Parse an ESA NEOCC Keplerian orbit file (.ke0) and return the parameters.
Parameters
----------
orbfile : str or path-like
Path to the local ``.ke0`` orbit file downloaded from the ESA NEOCC
database.
Returns
-------
a : float
Semi-major axis in au.
e : float
Eccentricity.
i : float
Inclination in degrees.
H : float
Absolute magnitude.
dadt : float
Semi-major axis drift rate in au/My, or ``math.nan`` if A2 is absent.
dadt_std : float
1-sigma uncertainty on da/dt in au/My, or ``math.nan`` if absent.
Notes
-----
A2 values are read in units of au/d^2 * 1e-10 and converted to da/dt
via :func:`A2_to_dadt`. If the ``NGR`` line is absent, both drift
quantities are returned as ``math.nan``.
"""
dadt = math.nan
A2 = math.nan
A2std = math.nan
dadt = math.nan
dadt_std = math.nan
AMR = 0.0
with open(orbfile, 'r') as f:
for line in f:
if 'KEP' in line:
a = float(line.split()[1])
e = float(line.split()[2])
i = float(line.split()[3])
if 'MAG' in line:
H = float(line.split()[1])
if 'NGR' in line:
# A2 in au/d^2
AMR = float(line.split()[1])
A2 = float(line.split()[2]) * 10**-10
# If AMR = 0, it was not determined
if 'RMS' in line and not math.isnan(A2) and AMR == 0.0:
A2std = float(line.split()[8]) * 10**-10
elif 'RMS' in line and not math.isnan(A2) and AMR != 0.0:
A2std = float(line.split()[9]) * 10**-10
if not math.isnan(A2):
dadt = A2_to_dadt(A2, a, e)
if not math.isnan(A2std):
dadt_std = A2_to_dadt(A2std, a, e)
return a, e, i, H, dadt, dadt_std