Overview

Newtonian (plain ol' inverse square) part

Asteroid perturbers

Natural satellite (moon, Galileans, etc.) perturbers

Non-spherical planet effects

Atmospheric drag

Non-gravitational effects

General relativity

Forces not modelled (at least yet)

Discontinuities in the force model

When numerically integrating an object's trajectory, there are a huge range of forces acting on it. In reality, any object in space is being acted upon by all the forces listed above. In practice, if Find_Orb evaluated, at each step of an integration, the forces exerted by the 300 asteroids it knows about, all the planetary satellites, general relativistic effects, and the complicated effects caused by planets not being perfect spheres, the program would drag almost to a halt.

As a practical matter, it has to be bright enough to apply the needed
forces at needed moments. If you have tracked a main-belt asteroid for a few
days, it's silly to include all these minor effects. The orbit is so poorly
determined (usually) that the errors in it are far greater than even planetary
effects. (The
Minor Planet Center, for example, routinely produces
unperturbed, "two-body" orbits for which the *only* force comes from the
sun. It can almost always do that for objects short- and medium-length arcs of
observation, and get results that wouldn't really be improved much by adding in
everything else. The errors caused by neglecting small effects are much, much
less than the errors caused by the data not being absolutely perfect.)

The following describes which forces Find_Orb can model, and the circumstances under which they are included.

The sun is always included as a perturber. Depending on which check boxes you've ticked, planets and the moon may be included as perturbers. Find_Orb will turn on some planet perturbers automatically when an object is loaded, making some choices based on how close the object happens to be to a planet. Its judgment of which will be needed is usually pretty good, but you can sometimes find that turning on another planet or two and doing a 'full step' helps significantly.

It's possible to override Find_Orb's automated judgment by
editing the `PERTURBERS` parameter in `environ.dat`. That
file contains a fair number of obscure settings, including this one.
You basically can set a value saying that objects X, Y, and Z are to
be included in the force model at all times. (Remarks in the file
provide details on how to do this.)

The mass of planets not selected as perturbers are "thrown in" to the sun when
outside the orbit of that planet. (This is a common trick in numerical
integration of planetary systems. If you're out in the asteroid belt and
don't want to include the perturbations of Mercury and Venus, it's a decent
approximation -- at least for many purposes -- to simply add their mass in
to that of the sun's.) The process is a gradual one, because
we want to avoid discontinuities. If you are more
than 120% of the planet's mean semimajor axis from the sun, the entire
mass is thrown in. If you're less than 100%, no mass is thrown in. In
between, we ramp up the mass linearly. For details, see the
`include_thrown_in_planets()` function in `runge.cpp` in the
Find_Orb source code.

If you tick the earth but not the moon, Find_Orb computes the position and mass of the earth-moon barycenter. As a result, including the moon as a separate perturber usually doesn't matter all that much, unless you're within a million or so kilometers of the earth/moon system.

Up to 300 asteroids can be included. (These are essentially the largest
asteroids, or at least those whose perturbing effects are most noticeable,
as described in
the BC-405 paper. Essentially, Baer and Chesley -- the 'BC' in BC-405 --
determined masses and orbits for these 300 objects.) Including all 300 asteroids,
all the time, can result in very slow integrations. Most simply won't matter
most of the time, since your object may be several AU from many tiny rocks.
Editing `environ.dat` and resetting the `BC405_ASTEROIDS` and/or
`ASTEROID_THRESH` parameters can help a lot, in the (usual) situations
where you are willing to ignore distant asteroids. You can
click here for more information
on asteroid perturbers.

For points within 0.03 AU of Saturn and Jupiter, the effects of the satellites of each of those planets are included separately. (This includes all four Galileans for Jupiter, and Tethys, Dione, Rhea, Titan, and Iapetus for Saturn.) Beyond that distance, the masses of these satellites are "thrown in" to their primaries.

The equatorial "bulge" (J2) terms, and the smaller J3 and J4 zonal terms, are included for the Earth, Mars, and all four gas giants, when within 0.015 AU of the planets in question. (J3 accounts for a planet being a little bit "pear-shaped". I don't know of a simple, non-mathematical way to explain J4.) This is, however, not the "whole story"; there are fairly elaborate models for the gravitational fields of all four inner planets and the Earth's moon, running (in some cases) to tens of thousands of terms. Such models also exist for Ceres and Vesta.

The earth's J2 (the oblateness term) is, on rare occasion, important for asteroids. (Quite rare. We had to wait until 2011 for an object close enough and moving slowly enough for J2 to have an observable effect.) J2 is almost always important for earth-orbiting objects. I've never seen higher-terms matter for an asteroid, but they often matter for artificial satellites. Find_Orb (as of January 2017) handles these terms, adding in more terms as an object gets closer to the earth, using some empirically derived methods to figure out just how many terms are needed at a given distance.

For the moon and non-Earth planets, Find_Orb only handles some lower-order zonal terms (J2, J3, J4). It actually wouldn't be all that hard to add in higher-order spherical harmonics for other objects (now that it's been done for the earth), but I'd want to have an object for which those effects would be measurable. Otherwise, I could add the force in with reversed sign or something like that and have no way of testing it. Quite good models exist for the gravitational fields of Mercury, Venus, Mars, and the moon, as well as for Vesta and Ceres, and the sun's J2 value has been determined.

This was added in October 2015, in anticipation of the re-entry into Earth's atmosphere of WT1190F. Only the earth's atmosphere is currently modelled, though it would not be too difficult to generalize the code to include the atmospheres of Mars, Saturn, Venus, and Titan (because models of the average atmospheric density as a function of altitude are available for those objects).

The key parameter for determining deceleration due to drag is the area/mass ratio (AMR). Fortunately, this is exactly the parameter that is determined when measuring the effects of solar radiation pressure. If you're turned SRP on in Find_Orb's force model, and have determined (or constrained) a particular AMR, Find_Orb will turn on atmospheric drag using that value.

For meteors, one can turn on SRP and do a seven-parameter fit. The AMR will then be determined from the object's observed deceleration. This has only been tested for 1964-054A = OGO-1, an artsat that dipped into the earth's atmosphere during several successive perigees in August 2020, eventually re-entering over Tahiti. We got quite a bit of data on it, showing the extent to which atmospheric drag slowed it, and enabling at atmospheric-drag based determination of its area/mass ratio.

The same trick really ought to be useful for meteors, but I haven't had a chance to give that a try yet.

If you have a newly observed artsat and don't know the AMR from observations (usually the case), you can enter a constraint such as "A=0.006" (m^2/kg), do a full step, and then compute ephemerides based on that assumed AMR. I've done that on occasions where we had some idea what the object was, and I could look up data on its mass and dimensions.

You *can* have an object that is repeatedly passing through
the upper atmosphere at perigee, and also being kicked around by solar
radiation pressure throughout its orbit (at least, the sunlit part).
In such a case, the "effective" area/mass ratio due to drag may be
quite different from the AMR determined from solar radiation pressure,
and we really ought to be doing an eight-parameter fit rather than a
seven-parameter one. I haven't gotten around to trying to fit both
parameters, though I do have data for a couple of objects where the
difference in AMRs might actually be detectable.

The simplest non-gravitational effect is solar radiation pressure (SRP). SRP can be expressed as a simple radial force proportional to the sun's gravity. In other words, the object bounces some sunlight back toward the sun, cancelling out some of the sun's gravity. If it were "fluffy" enough, i.e., had a really high area/mass ratio, it would bounce back enough sunlight to cancel out the sun's gravity and travel in a straight line relative to the sun (except for planetary perturbations, of course.) However, for a one-kilogram object, you'd need about 1300 square meters of "sail area" to manage that cancellation. There's a reason we aren't making lots of solar sails.

For comets, two parameters (sometimes three or four) are used. The comet parameters A1 and A2 (or A1, A2, and A3, and sometimes DT) follow the standard Marsden-Sekanina model. In this model, there is almost no cometary force at 2.6 AU or more from the sun (because almost no outgassing occurs), and then immense forces close to the sun as the comet starts boiling away.

For rocks and artificial satellites, the force will be a simple
inverse square one. That matches the physics for rocks such as 2009 BD
and (101955) Bennu, where there is an A1 outward force
similar to SRP and an A2 tangential force that is accelerating or
decreasing the object's speed (i.e., increasing or decreasing the
object's orbital period/semimajor axis). You could consider this to be
a "poor man's Yarkovsky" model. (In general, a *very* poor man's
Yarkovsky. It should work quite well if the object's pole is nearly
perpendicular to its orbit; for example, it results in an excellent
improvement in the fit for (101955) Bennu, which presumably has a spin
axis perpendicular to the plane of its orbit. If its spin axis was
closer to being in the plane of the orbit, like that of Uranus, it
wouldn't model the actual physics well.)

For 2009 BD, an object only a few meters across for which we have a good bit of data, both A1 and A2 are significant. For larger objects, you can set the constraint A1=0 and be on safe ground; any A1 you determined for Bennu, for example, would be spurious.

Solar radiation pressure used to be treated as 'always on', even
when an object is going through the earth's shadow. Since early 2020,
Find_Orb has decreased such forces if the sun is partially eclipsed by the earth,
and they go to zero if the sun is totally eclipsed. For certain
ETBOs ("empty trash bag
objects", probably bits of spacecraft insulation), this turned out
to be quite significant; the residuals in the observations decreased
a *lot* once shadowing effects were considered.

Note that, at present, no other shadowing effects (due to the
moon or other planets) are included. And the earth is treated as
round and airless. I don't have any observations yet sensitive
enough to show such effects. Limb darkening (the fact that you get
more photons per square arcsecond near the middle of the sun than
you do from the edges) *is* taken into account, though.

In some cases, the comet activity lags the comet's actual position by days or months. It just takes the comet a while to warm up enough, such that maximum non-grav forces only happen a bit after perihelion. This is fitted in Find_Orb with the DT parameter.

GR is included only for the sun. Click here to read about the mathematical method used to compute GR in Find_Orb.

At present, this first-order approximation is more than ample.
We'd need *much* better observations to notice any shortcomings.
But if Gaia truly produces results as wondrous as claimed, in the milli-
or even micro-arcsecond realm, it's possible that some smaller effects
will be noticeable. GR from planets and second-order terms from the sun
may have to be included. Most likely, I would use the "parameterized
post-Newtonian" (PPN) method described in the
* Explanatory Supplement to the
Astronomical Almanac * on page 281. This would be quite slow, but
I'd keep it as an option, and would determine which of the many terms involved
can be dropped in normal use without meaningfully affecting the results.

As of Gaia DR2, the effects of relativistic light-bending have become (just barely) noticeable, and are now included.

The forces not yet modelled are mostly effects I can't even detect yet: that is to say, they are so small that if I included them, the fit to the orbit wouldn't change noticeably. That makes them difficult to add them in; I can't tell if I've done it correctly or not.

One exception, which I *can* detect, is the Yarkovsky effect.
(And possibly YORP.) This looks to be quite tricky to model. You need
to know a good deal about the object in question, such as its spin
axis and rotation rate.

For some objects where we don't have anywhere near that level of
knowledge, the Yarkovsky effect is modelled simply in terms of a
gentle acceleration or deceleration in the direction of motion.
Basically, solar radiation pressure with an added along-orbit
component. You can do that by turning on the two-parameter SRP
model in Find_Orb, and then setting the constraint A1=0. So
this isn't a *completely* unmodelled, or at least unmodellable,
force in Find_Orb.

I don't expect to add a "real" Yarkovsky or YORP model soon. (Though there is a sort of "poor man's Yarkovsky" model in four parameters that I'd like to try. In this model, three of the parameters are Euler angles; the idea is that the incoming solar radiation is absorbed, the object rotates by that amount, and then re-radiates, with the magnitude of that re-emission being given by the fourth parameter. This could apply to both rocks and artsats, though I think it may not work well with artsats; the polar axis of an artsat, it seems to me, is apt to precess much more rapidly than that of a rock.)

There are three situations I know of (and probably more I don't
know of) where adding an extra perturber might help. Jupiter's
satellite Himalia has a small,
but noticeable, effect on Elara, sufficient to allow the mass of the
former to be determined. There are satellites of Saturn that may
react to the mass of Phoebe (though this is uncertain, and it's true
that we have a much better mass for Phoebe from the
*Cassini* flyby on 2004 June 11; trying to determine the mass
from satellite astrometry isn't going to help, though it would serve as
a sanity check on Find_Orb's ability to determine such masses.) And
I've heard from 'Planet X' researchers who would like to be able to
model the effects of hypothetical masses on observed TNOs. To
accommodate these situations, I would allow one to first generate an
ephemeris for Himalia or Phoebe or the hypothetical perturber, and then
tell Find_Orb to use that object to perturb Elara or a different
Saturnian satellite or on a TNO.

As described in the section on natural satellite perturbations above, only ten planetary satellites are modelled: the moon, the four Galileans, and five satellites of Saturn. (In reality, the remaining satellites of all the planets could be approximated, quite well at large distances, as an increase in their J2 (oblateness) factor. I should look into that...)

As described in the above section on Earth and other irregular objects, Find_Orb handles all effects that have been observable thus far. There is some possibility of solar J2 and maybe even J4 being detectable (and even measurable) with radar astrometry, so I may have to add that in. I do treat the irregular parameters as static; there's a seasonal change in the Earth's J2 as ice melts and freezes, for example, and longer-term variations as some of the ice fails to re-freeze.

General relativistic effects are currently handled only for the Sun, and only to first order. Someday, we may have astrometry good enough to measure such effects for other planets, or higher-order effects for the Sun. (At present, we aren't even close.)

Some care was taken in the force model to ensure that the accelerations and their derivatives would be both finite and continuous, even when the objects pass through planets. It's true that real objects never do that. But hypothetical ones do, and orbits in the initial stages of determination may do so (before, we hope, converging on an above-the-surface solution). Numerical integration becomes simpler if you can be confident that there are no discontinuities.

The only real problem in doing this was handling situations where objects pass inside of planets. In the real world, this causes a very abrupt discontinuity, i.e., an impact. But again, I really wanted to avoid discontinuities, and therefore didn't want the acceleration to increase to infinity as the object approached the center of the planet.

The solution was this. Accelerations are computed in their 'normal' manner and left unchanged for distances greater than 90% of the planet radius. Within about 60% of the remainder, accelerations are zero; i.e., there is a "bubble" inside each object within which the object has no acceleration at all relative to the sun.

Between those limits, accelerations are multiplied by a factor computed
in the `compute_accel_multiplier()` function in `runge.cpp` (see
the source code for Find_Orb).
This factor is exactly 1.0 at the 90% limit, and zero at the inner "bubble"
limit. Its derivative at both points is zero, i.e., both this function and
its derivative are continuous across the outer and inner limits. It is
computed as a cubic spline. Thus, instead of having accelerations climb to
infinity as you approach a planet center, they instead reach a peak, then
drop smoothly to zero.

All this weirdness has no end effect on what orbit is determined; the only thing it does is to make sure the program doesn't hang, or go off on weird tangents, when determining the orbit. It can admittedly lead to odd effects if you ask Find_Orb to determine ephemerides for times after an object hits a planet. The object may continue through the planet until it reaches the "bubble", then continue on a straight line until it reaches the other side of the "bubble", then climb up through the planet and come out the other side. For the earth, with atmospheric drag turned on, you may see the object slow down due to drag and then bounce around inside the "bubble".