Russian version of the page

Description of the site "Numerical simulation of the motion of a small body of the Solar system"

This site consists of several computational programs with web-interfaces which are intended for simulating the motion of asteroids and comets. The site contains now the folowing three programs: Simulation of the motion of a small body, Transformation of the small body's coordinates and Transformation of the date. It is planned to add soon a program of determination of the orbit of a small body using its astromertric, radar and satellite observations.

The specificity of the site is a possibility of computing perturbations using different planet's ephemerides: DE405, DE408, DE414, DE421, DE422, DE423, DE424, DE425, DE430, DE431, DE432, DE433, DE434, DE435, DE436, DE437, DE438, DE440, DE441, [1] EPM2011, EPM2015 [2]. Besides, the site allows one to transform the motion parameters of a small body from one kind to another (for example, from rectangular coordinates and velocities to osculating orbital elements), to compute these parameters for another epoch by means of the Keplerian formulae of the non-perturbed motion, to estimate the accuracy of numerical integration by means of comparison of the results of the forward and backward integration, to transform a date from one time scale to another. For a convenience of usability of the site it is available to select both a way of output and a format of output for all computed parameters.

Below there is a description of the web-interface of the basic program of the site

Simulation of the motion of a small body

This program is mainly indended to compute different motion parameters of a small body at specified time moments on the base of pre-known orbital elements. To use the program one should input initial motion papameters into the field with the same name. The format of parameters' location is an arbitrary: you can input all parameters into a single line or into different lines or in another way. You just have to type at least one space between the parameters, however, a new line symbol is not a space, so if your parameters are located in several lines then at least one space must be typed either before or after each parameter. Moreover, a set of the orbital elements and their order must strictly match the description in the column "The way of input". If your elements have another structure then you can use the auxiliary program "Transformation of the small body's coordinates". This program can transform an arbitrary set of Keplerian elements into any other set, else it can transfrom orbital elements into rectangular coordinates and velocities and vice versa, and also it can recompute the motion parameters to another epoch by means of the formulae of non-perturbed motion. The rectangular coordinates and velocities of a small body may also be used as its initial motion parameters. You should specify with the switch above which set of initial data you want to use.

The epoch of the initial parameters should be input in the field "Epoch, JD". It must be a Julian date. If your epoch is of another kind then use the program "Transformation of the date" which can transform any anno domini date from one kind to another.

Below initial motion parameters the mass of the considered bodye should be input (in the units of the Sun's mass which is 1.989 · 1030 kg). If the body is quite small a zero mass may be input.

The integration interval is defined by the start and finish ephemeris times which must be input in the form of Julian date. The field "Output step" defines the time interval between nieghboring outputs in the computed ephemeris. This step may be specified in days, hours or minutes (this is defined by a switch). During the integration process computed data will be output by means of exits from the integrator with specified output step.

There is also a button "Forward and backward integration" which runs an integration at first in forward direction - from the start ephemeris time to the finish one (or in contrary direction if the initial epoch lies after integration interval), and then in backward direction. During the forward integration all coordinates and velocities calculated for the output time moments are stored in the computer's memory and during the backward integration at these moments the differences in position (dr, AU) and velocity (dv, AU/day) are calculated and output in a new window. If the initial epoch is located between the start and finish ephemeris time then the described process is performed in two stages: at first forward upto finish ephemeris time and back, and then backward to start ephemeris time and back.
The time moments for dr and dv are output in the same manner as when computing the ephemeris. It means that if you want to specify a way and format of output of the time moments you should scroll the page down and do this in the lowest line of the table "What to output in the ephemeris".

A following switch defines the basic plane (equator 2000.0 or ecliptic 2000.0) of the reference frame of the initial motion parameters.

The following two switches are intended to include relativistic and nongravitational effects into the model of the motion. Relativistic effect noticeably influence on the body's motion just if its perihelion distance is quite small (< 0.5 AU). Non-gravitational forces should be taken into account just for comets. To do this turn on a corresponing swith and type values of the nongravitational parameters A1 and A2.

Further there is a section of computing of perturbations from planets and the Moon. With the help of turning off or on of selected perturbations one may research their structure, i.e. a contribution of the gravitation of some panet or several planets into the change of small body's orbital elements. When computing these perturbations the coordinates of planets and the Moon are interpolated with the help of Chebyshev's polinomials. The coefficients for these polinomials are extracted from the ephemeris selected by a user (DE405,..., DE441, EPM2011, EPM2015). In the brackets near the name of each ephemeris a time span covered by this ephemeris is presented in the form of a pair of Julian dates.

Here is also a switch which permits one to use the selected ephemeris in the form smoothed at the junctions of neighboring interpalation intervals. This smoothing has been done so that the interpolated planets' coordinates and their first derivatives keep continuity at the junctions upto 33-34 decimal places, i.e. they are completely continuous at the junctions in quad-precision calculations. The algorithm of such smoothing and some its applications are described, for example, in the papers [3, 4]. The smoothed variants of ephemerides it has meaning to use just in quad-precision calculations, i.e. when the 128-bit floating-point numbers' capacity is selected (this capacity corresponds to the decimal precision of 34 digits).

Besides, here is a switch to turn on or off a so called correction of the integration step according to ephemeris interpolation intervals. This correction is realized so that when some junction of interpolation intervals falls within some integration step then this step is finished at this junction. Such a correction [3, 4] allows to avoid a loss of accuracy which occurs at the junctions due to the breaks of the second- and higher-order derivatives of interpaloted coordinates. Note that this step correction in the realized here computing program may be applied just with a variable integration step (i.e. with a positive value of the LL parameter) and without computing of approaches to the planets.

Further following parameters should be specified: an initial integration step (in days); the order of the Everhart's method [5] of numerical integration of odinary differential equations; LL parametr which define a behavior of the integration step:
if a negative value is specified for LL then the integration step will always be constant and equal to the initial step,
if a positive value is specified for LL then the integration step will be variable and its value will be defined from the condition of the local error not to exceed the value 10-LL.
Hence, as it is clear from the practice of application of the Everhart's method the real accuracy of the method is always higher than this value, so the optimal value of LL is approximately equal to the two thirds of the order of numerical integration method.

As for the order itself, it is empirically established that its optimal value is one less than used decimal precision, i.e., for example, when using 16-digit decimal precision (64-bit floating-point capacity) the optimal order is 15. However, as the realized here program permits to use only several pre-defined values of the order (from 7 upto 39 with a step of 4) then for the other two precisions (19-digit and 34-digit) respectively the 19th and the 27th order is used.

After the integrator's parameters a bit capacity of the used floating-point numbers should be specified. The choice of this parameter is small and entirely determined by the capabilities of the compiler used. Since all programs that perform calculations on the site are written in the Fortran language modern compilers of which support only four floating-point number sizes (32, 64, 80 and 128 bit) then the last three sizes has been used. The size of the numbers of 32 bit (this is so called single precision numbers that contain only 7 decimal digits) in the site is not used because such numbers do not provide the necessary accuracy of astronomical calculations.

The last input parameters in the table are swithes and fields intended for the computing of the small body's approaches to the planets. When the computing of approaches to some planets is enabled on the right you should type maximal distances to the centers of corresponding planets. Approaches would be computed just if the distance to the planet's center is less than this specified value. Approaches are calculated by means of squared distance interpolation by three consecutive integration steps (i.e. by four time moments) with the help of a Lagrange polynomial. The derivative of this polinomial is equated to zero and the obtained quadratic equation is solved. From two roots the one for which the distance is less is chosen. Again, note that when the computing of approaches is enabled it is impossible to enable the correction of integration step according to ephemeris interpolation intervals.

Further there is a section of output information. At the request of the user as the output parameters it's possible to choose: right ascension and declination, rectangular heliocentric coordinates and velocities, orbital elements, distances to the planets. The units and the format of numerical data are specified by switches. When the computation of right ascension and declination is enabled the delay of a light signal and daily parallax are taken into account. For the computing of daily parallax it is necessary to specify a place on the Earth surface for which the ephemeris must be calculated. As these places just observatories or the geocenter may be used. So you should type a corresponding Observatory code. Right ascension and declination are output in the reference frame of the standard epoch 2000.0.

After clicking the button "Compute" a process of calculations is run and the selected above parameters are output in separate columns. The results are output in this page below the button "Compute". This is done so (output into the same window but not into a new one) because under such organization of the processing script the all data (the positions of switches, number data, etc.) input by the user are saved. So the user can save the page in his web-browser and all input changes of the page also will be saved. This permits to open the page later and its state will be the same as at the moment of saving. The user only needs to run a computing before saving because all input data are being stored just at the moment of calling to the server, i.e. when the button "Compute" is being clicked. Besides, when saving the page it's necessary to select the option "The entire (whole) web-page" but not the option "HTML only". After opening of the page saved in a such mannner this page may be changed and re-used as a page of the site: when clicking the button "Compute" the server solar.tsu.ru will be called and a new computation will be done, the new results will replace the former ones. The obtained page may again be saved etc.

References

1. ftp://ssd.jpl.nasa.gov/pub/eph/planets
2. ftp://quasar.ipa.nw.ru/incoming/EPM
3. Baturin, A.P., Influence of the methods of constructing ephemerides of major planets and the Moon on the accuracy of predicting motion of asteroids, Russ. Phys. J., 2015, vol. 57, no. 9, pp. 1230–1238.
4. Baturin, A.P., Enhancing the accuracy of numerical integration of the equations of asteroid motion with perturbations from major planets and the Moon from the DE ephemerides, Sol. Syst. Res., 2018, vol. 52, no. 4, pp. 360–363.
5. Everhart E. An efficient integrator that uses Gauss-Radau spacings // Dynamics of comets: their origin and evolution // Proc. 83rd IAU Colloq. Rome, 11–15 June 1984 / Eds Carusi A., Valsecchi G.B. Dordrecht: D. Reidel Publ. Co., 1985. P. 185–202.

Author, web-master: A.P.Baturin (TSU, Tomsk)