Newtonian (plain ol' inverse square) part
Natural satellite (moon, Galileans, etc.) perturbers
Non-spherical planet effects
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, and so forth, it 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 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 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. 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.
The mass of planets not ticked 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, 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 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) situationse 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 the Earth and Mars and moon. People dealing with low-orbiting satellites of those objects have to include a lot of spherical harmonic terms to model satellite motion correctly. Find_Orb (at least as of October 2015) doesn't handle those terms yet. (I've been working on it, but there seem to be some strange things going on, probably due to bugs in the code.)
So far, I haven't dealt with objects where those smaller terms mattered. For asteroids and high-orbiting satellites, you can't even tell these effects exist. (For these objects, J2 is often important, and J3 and J4 are somewhat noticeable.) But if I start dealing with lower orbiting earth satellites, that situation will change: higher-order terms and spherical harmonics will be required. (I've received some data for a low-orbiting satellite to use as a test case.)
This has only recently (October 2015) been added, and is somewhat experimental. 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.
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. (I should note that I haven't had an opportunity to test this. I'd love to see some astrometry for meteors and give it a try.) For WT1190F, we have a satellite which has been observed long enough for the AMR to be determined from solar radiation effects. If you have no idea what the AMR might be 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.
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) are used. The comet parameters A1 and A2 (or A1, A2, and A3) 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, one can edit environ.dat and adjust the 2009BD parameter to have the force dependence 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. Otherwise, it doesn't model the actual physics well.)
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 nay 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.
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. I don't expect to add Yarkovsky or YORP soon.
For low-orbiting satellites, the handling of non-spherical planet effects will have to include higher-order terms. This hasn't been an issue yet, but it would be nice to be able to handle low-orbiting objects.
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...)
General relativistic effects are currently handled only for the Sun, to first order. Someday, we may have astrometry good enough to measure such effects for other planets, or higher-order effects for the Sun.
Some care was taken in the force model to ensure that the accelerations and their derivatives would be both finite and continuous. The real problem in that department was situations where objects pass inside of planets. In the real world, this causes a very abrupt discontinuity, i.e., an impact. But in the initial stages of orbit determination, it can make sense to compute orbits passing inside planets. Also, abrupt discontinuities anywhere would be something of a headache for numerical integration.
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 planet exerts no force at all.
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".