Planets and Minor planets ephemeris, and matching to sources

This page provides examples on how to generate minor planets' ephemeris and how to match Astrocatalog objects to known asteroids. Before you start, it is recommended to read the celestial.INPOP and celestial.OrbitalEl pages.

Another relevant class is the MovingSource class.

Generate ephemeris of asteroids

There are several ways to get asteroids ephemeris. You can either query the JPL horizons web page, or use the celestial.OrbitalEl object to calculate their position. The Earth position (and planetary perturbations) are calculated using the celestial.INPOP class.

Getting data from JPL horizons:

% To get a table of ephemeris for asteroid 299 - today:
[T] = celestial.SolarSys.getJPL_ephem('299;','EPHEM_TYPE','OBSERVER','TimeScale','TT'); 

% ephmeris of asteroid 9804 in some range of dates:
% see help for many more options
[T] = celestial.SolarSys.getJPL_ephem('299;','EPHEM_TYPE','OBSERVER','TimeScale','TT','StartTime',[1 1 2023],'StopTime',[1 2 2023]); 

% To get the J2000.0 Ecliptic Cartesian barycentric positions (every 1 hr)
T = celestial.SolarSys.getJPL_ephem('299;','EPHEM_TYPE','VECTORS','TimeScale','TT','CENTER','500@0','StartTime',[1 1 2023],'StopTime',[3 1 2023],'StepUnits','h','StepSize',1); 
 
% to get heliocentric orbital elements, for today epoch
[T] = celestial.SolarSys.getJPL_ephem('299;','EPHEM_TYPE','ELEMENTS','TimeScale','TDB');

If you like to interpolate the JPL ephemeris to other dates:

% get ephemeris from JPL:
[T] = celestial.SolarSys.getJPL_ephem('299;','EPHEM_TYPE','OBSERVER','TimeScale','TT','StartTime',[1 1 2000],'StopTime',[2 1 2000],'StepSize',10,'StepUnits','m','OutType','astrocatalog'); 
% interpolate to two new data:
IT=T.interpCoo(2451545.6+[0;0.001],'ColJD','JDTT');

Other related tools exist. For example, to retrieve the list of current close approaches of NEOs:

[Data,Ephem]=celestial.SolarSys.getJPL_CloseApproachNEO

To make some predictions of TNO occultations using the JPL ephemeris:

% The ElementsIndex 1 is for using all the numbered asteroids (By default only asteroids with semi-major axis >9 au are used):
[Result, Table] = celestial.conjunctions.conjunctionsSearchStarMP('ElementsIndex',1,'StartDate',[1 12 2023],'EndDate',[1 1 2024]);

Using the celestial.OrbiralEl class

Search for known Asteroid in a specific position and time

The first step is to prepare a celestial.OrbitalEl object containing the orbital elements of all known asteroids.

% Download the most updated orbital elements files from JPL for numbered asteroids, unnumbered asteroids, and comets:
I = Installer;
% This will download the 3 files to ~/matlab/data/SolarSystem/MinorPlanets/
% numbered asteroids, unnumbered, comets:
I.install('MinorPlanets'); 
% and ~/matlab/data/SolarSystem/MinorPlanetsCT/
% which include a single list of numbered and selected unnumbered asteroids with a common epoch:
I.install('MinorPlanetsCT');


% Install INPOP
I.install('INPOP');

% load and populate INPOP
IN = celestial.INPOP.init;

% Read the files into an OrbitalEl object:
% The following will read all the common epoch asteroids:
OrbEl= celestial.OrbitalEl.loadSolarSystem('merge');

% alternatively, read all the 3 files (numbered, unnumbered, comets):
OrbEl= celestial.OrbitalEl.loadSolarSystem

% Merge the three elements into a single element OrbitalEl object.
% Note it is best to remove asteroids with old epochs, as their orbits are highly uncertain.
El = merge(OrbitalEl, 'MinEpoch',2451545, 'MaxEccen',0.9999);

% For the sake of the example, you can also read only the numbered asteroid list:
El = celestial.OrbitalEl.loadSolarSystem('num');

% The next step is to ensure that all the orbital elements in the list have a common epoch.
% For the full list this may take a few hours, but for the numbered list it should take seconds: 
El.propagate2commonEpoch('INPOP',IN);

Next, we will perform the search. The accuracy of asteroids ephemeris depends on the method of calculation. By solving the Kepler equation, the accuracy may degrade very fast if the time at which the coordinates are evaluated is far away from the epoch of the orbital elements. Typically, the accuracy may degrade to the level of a few arcsecond after 100 days. Therefore, in the function celestial.OrbitalEl/searchMinorPlanetsNearPosition, the search is done in two steps. In the first step, we calculate the position of all the asteroids specified in the OrbitalEl object using the Kepler equation. Next, for all the asteroids found within the search radius plus a search buffer from the search position, we evaluate their position using direct integration of the orbit. This provides more accurate coordinates.

% Define the position and coordinates you would like to search:
JD = celestial.time.julday([1 1 2023 17 57 12]);
RA = 180.0; % deg
Dec = 0; % deg
% search within 1000 arcsec from specified location:
[Result] = searchMinorPlanetsNearPosition(El, JD, RA, Dec, 1000, 'CooUnits','deg', 'SearchRadiusUnits','arcsec', 'INPOP',IN, 'ConeSearch',true);
% The output is an AstroCatalog object with the list of asteroids found within the search radius.

Search for known asteroids among the entries of an AstroCatalog object

Given an (AC) AstroImage populated with an Astrocatalog, or an AstroCatalog, you can search the AstroCatalog for known asteroids:

% Loading - see previous example for details:
OrbEl= celestial.OrbitalEl.loadSolarSystem('num');
OrbEl.propagate2commonEpoch;
IN = celestial.INPOP;
IN.populateAll;

% Search, with 1" search radius
[OnlyMP, AstCat, AC] = imProc.match.match2solarSystem(AC, 'JD',JD, 'GeoPos',[], 'OrbEl',OrbEl, 'SearchRadius',1);

% OnlyMP is an AstroCatalog containing the sources selected from AC which were matched with a known asteroid.
% a few columns are added (e.g., angular dist to asteroid and its name).

% AstCat, contains the full AstroCatalog of asteroid in this search.
% to save time you can use it in the next search (assuming the same JD), using:
[OnlyMP, AstCat, AC1] = imProc.match.match2solarSystem(AC, 'JD',JD, 'GeoPos',[], 'AstCat',AstCat, 'SearchRadius',1);

% AC is the updated AC input, with additional columns (e.g., DistMP - the distance to the nearest known asteroid; NaN if no asteroid within the search radius

Plotting the position of known asteroids in ds9

You can use the DS9analysis tool to plot the position of known asteroids on top of a ds9 image:

% Put your image file name here:
FileName = 'LAST.01.02.01_20230425.184452.001_clear_185-02_000_001_002_sci_coadd_Image_1.fits';
% create a DS9analysis object:
D9=DS9analysis;               
% load an image and display in frame=1 (the image will be stored in the 'Images' property.                                                            
D9.load(FileName);

% you can also load the image into a different frame:
D9.load(FileName, 'Frames',2)

% To plot the position of all known asteroids on top of the current image:
OrbEl = celestial.OrbitalEl.loadSolarSystem('merge');  % loading this time takes time, but you have to do it once
INPOP = celestial.INPOP;
INPOP.populateAll;
% Mark the positions of all known asteroids (note that the position maybe off by a few arcsec):
[KA] = D9.plotKnownAst('OrbEl',OrbEl, 'INPOP',INPOP);

% to delete the ds9 markers (regions)
ds9.delete_region

% Mark the positions of all known asteroids up to some mag limit
[KA] = D9.plotKnownAst('OrbEl',OrbEl, 'INPOP',INPOP, 'MagLimit',21);

Forced photometry of moving target

The forced photometry utilities are applicable for moving targets and examples are given here.

Select sources near the trajectory of a moving target