Neil's Qbasic Programs

Qbasic programs
by Neil H. Jacoby, Jr.

Laplace.bas Program

This program computes a heliocentric orbit of an object say a planet or an asteroid,or even a comet; given three observations of its right ascension,alpha, and its declination,delta. Each of these observations must be nearly equally spaced in time as much as possible and also the solar coordinants must be at the same time the observations are made. These solar coordinants are obtained from the Astronomical Almanac which is published by the U.S. Naval Observatory. In all these orbit determination programs the observed object's right ascension,alpha,and its declination,delta are all based on the earth's equator system of coordinates as well as the solar coordinates,X,Y and Z; whereas the orientation angles i,Node and Omega of the orbit of the observed object are all based on the earth's ecliptic system of coordinates. The user of all these orbit determination programs should always remember that due to the earth's precession, the equator and vernal equinox changes through time. Therefore the object's right ascension and declination referred to a particular equator and vernal equinox of date must agree with the solar coordinates X,Y,Z referred to the same equator and vernal equinox of that date. If the orbit of the object is geocentric, then alpha,delta,X,Y,Z, and the orientation angles i,Node,and omega are all based on the earth's equator system of coordinates. Of course, the x-axis always points toward the vernal equinox in both the equator and ecliptic coordinate systems. The observation times are then measured in minutes. In addition, the times of each of these observations must be specified and the units are in days. The output of this program are the orbital elements which are a,the semi major axis, e , the eccentricity of the orbit and the orientation angles, which are ,i, the inclination of the orbit to the ecliptic plane,the ascending node OMEGA, and the argument of perihelion, omega. The final element, T, which is the time of the last epoch passage is not computed in this program but will be added later. The time interval between each of the three observations must be small enough to permit use of the f and g series so as to minimize truncation error but large enough for accurate orbit determination. Next the input, which this program asks the user, is an smart initial guess for r, which is the distance in astronomical units (a.u.'s) from the sun to the observed object. The final value for r is determined from Newton's method of iteration where delta(r) is a correction to r in the Laplace method. This iteration continues until the ABS(delta(r)) is less than or equal to eps, where eps is the tolerance or limit of convergence. Usually eps is less than or equal to 0.000001 a.u., which the user must specify. On this same line of input the user must also specify eps2 which is the tolerance (in degrees) for the representation of the observations of alpha and delta. Usually eps2 is specified to be equal to .00001 of a degree of the actual alpha and delta. eps2 is used in the Leuschner differential correction portion of this program using three observations and the one for five observations. Finally, the program asks the user to specify the maximum number of iterations,i.e. imax, and mu, where mu is defined as the ratio of the mass of the sun plus the observed object to the mass of the sun. Since the mass of the observed object, say an asteroid or comet, is negligible compared to mass of the sun, then mu can assumed to be equal to 1.0. In practice it is best to set imax = 10. An smart guess for r can be found as follows. First, most of the objects of concern to us in our Solar System are the asteroids which lie between the orbits of Mars and Jupiter, where r ranges between 1.52 a.u.'s and 5.23 a.u.'s from the sun approximately. Second, as long as the eccentricity of the observed asteroid's orbit is not too large, then r can be approximated by knowing the asteroid's synodic period, which can be observed or approximated from the observations. The synodic period is defined as the time interval between two successive passages of the observed asteroid relative to the same stellar background as observed from earth. From knowing this synodic period, the observed asteroid's sidereal period can be computed from the following formula. First,let E = 1 year,sidereal period of the earth. Then let Psy = observed synodic period of asteroid,in years. Therefore the observed asteroid's sidereal period,P,is found from P = Psy/(Psy-1),where P is in years. Then knowing the observed asteroid's sidereal period, r can then be approximated by applying Kepler's Third Law which is P^2 = r^3, where P is in years and r is in a.u.'s. This program uses the Laplace method of preliminary orbit determination and the Leuschner differential correction method. These methods use the f and g series which is carried up to and including the sixth time derivative. The final output of all these orbit determination programs are the elements of the preliminary orbit and the final elements after the Leuschner differential correction is completed. If the three or the five observations are good, then convergence in the Leuschner differential corrections, should achieve the desired tolerance within say 0.0001 of a degree of the actual observed alpha and delta and within ten iterations. Five Observation Laplace5.bas Program This program works the same way as the three observation-Laplace, except five observations instead of three must be inputted as well as the solar coordinates which are obtained at the same time the observations are made. As before, these solar coordinants are obtained from the Astronomical Almanac. It was found that five observations instead of the usual three, that overall greater accuracy was obtained in determining the orbital elements, i.e., the semi major axis, the eccentricity, the inclination, the ascending node, and omega, the argument of perihelion, for low to moderate eccentricity orbits. In addition, this five-observation program outputs T, the time of last perihelion passage. As before, the orientation angles, the inclination i, angle of the ascending node, called node, and omega, the argument of perihelion, are all based on the earth's ecliptic coordinate system. A very important finding so far, in this research on the use of the Laplace and Leuschner methods of orbit determination, is that each of the three or five observations should be chosen such that the object's right ascension and declination be nearly equally spaced as much as possible in alpha and delta as well as their observation times for the best accuracy. It was found that if any two of the three or five observations of alpha and delta are spaced too close together, poor determination of the orbital elements resulted and no convergence in the differential corrections. This unfortunate situation can occur if any two or more of the observations is at or very near a point when the observed object is starting or ending retrograde motion as observed from earth. Therefore it is best to choose the observations during the time when the observed object is totally in either direct or retrograde motion and not at the beginning or end of retrograde motion. Also is was found that the most accurate heliocentric orbit determinaton resulted when the time interval between each succesive observation was 15 to 25 days for the asteroids between the orbits of Mars and Jupiter. For heliocentric objects whose orbits are beyond Jupiter's orbit this time interval should range between 25 to 50 days. Another very important finding in this research was that the most accurate determination of the orbital elements occurred when the middle observation of the three and five observations was at or very near the object's opposition point,i.e.,its elongation angle was 180 degrees. (The elongation angle is defined as the angle measured from the direction to the sun to the direction to the object as observed from earth.) This was the finding for low to low-moderate eccentricity heliocentric orbits. Also at or near the opposition point the object's right ascension and declination change most rapidly, where the object's motion is nearly perpendicular to the line of sight. In the case of high eccentricity orbits, such a comet's orbit, the three or five observations should be chosen such that the comet's change in right ascension and declination are maximum. This maximum change may not necessarily be at or near the opposition point however. I conducted further research in determinating the elements of the orbits of comets using three and five observations of right ascension and declination and the important result of this research was that using three observations instead of five gave a far more accurate determination of the elements of the comet's orbit. The reason for this is that there was not only rapid convergence within ten iterations in the differential corrections to the desired tolerance, but also a very accurate prediction of the comet's right ascension and declination at some future time. This prediction accuracy was found to be within 0.01 of degree of the actual comet's right ascension and declination in almost all cases. On the other hand, when five observations of right ascension and declination were used, poor results were obtained because there was no convergence in the differential corrections no matter how many iterations were allowed by the user. Apparently, the additional two observations corrupted the final results because of the very high eccentricity of the comet's orbit. These unfortunate results occurred in almost all cases. Therefore it is strongly advised that the user use three observations and not five in determining the elements of very high eccentricity orbits of comets. The orbits of Comet Hale-Bopp and Comet S4 Linear were used in this research because their orbits have a very high eccentricity,i.e., 0.9 < e < 1.0. Finally in cases of orbits ranging from low to high eccentricity heliocentric orbits, all of the three or five observations of the object should have an elongation angle greater than 120 degrees east of the sun and greater that 120 degrees west of the sun for the best visibility of the object being observed. Another Laplace3.bas Position This program computes the heliocentric object's right ascension ,alpha, and declination, delta, given the object's orbital elements, the semi-major access a, the eccentricty e, the inclination i,the angle of the ascending node, called node, the argument of perihelion, and the mean anomaly, M, at the time of the reference observation. The angles i, Node, and omega all expressed in degrees. If the 3 observation program is used, the reference observation is the second observation, and if the 5 observation program is used, the reference observation is the 3rd observation. The desired time of the future observation must be specified and the solar coordinants at this desired time, must also be specified. As before, the solar coordinants are obtained from the Astronomical Almanac. The desired time of the future observation is measured from the time of the second observation if 3 observations are used and if 5 observations are used, this desired time is measured from the time of the 3rd observation. As before, the time is measured in days. The user of this program, should remember, that the orientation angles, i, node, and omega, are based on the ecliptic system of coordinants and the right ascension and declination are based on the equator system of coordinates. Finally, always remember that the orientation angles for heliocentric orbits are always based on the ecliptic system of coordinants and the orientation angles for geocentric orbits are always based on the earth's equator coordinate system. Kepler.bas Program This program solves Kepler's equation when the Mean Anomaly,M, and eccentricity, e , are given and one needs to solve for the eccentric anomaly,E. E, is solved by applying Newton's method of iteration. In most cases, M and e are given and not the eccentric anomaly,E. Kepler's Equation, M=E-e*sinE, is also used in the closed form of the f and g series. Kepler's Equation is useful when one is given the orbital elements and one needs to predict the object's, right ascension and declination at some given future or past time. When running this program, the program will ask the user to input the eccentricity,e,and the mean anomaly M in radians, and the specified limit of convergence or tolerance,eps. Usually eps is 0.000001 of a radian. If e is between 0 and 1 Kepler's equation will be directed to the elliptical part of the program, and if e is greater than 1 Kepler's equation will be directed to the hyperbolic part of the program. The case of parabolic orbits, where e = 1 and the semi-major axis a is infinity, is not presented here because these types of orbits exist only in theory and not in reality. Here's a Demonstration of the Solution to Kepler's Equation. In all these heliocentric orbit determination programs, which I developed, it should be noted that these programs involve only Keplerian (two-body) orbits and that no perturbations are included. The orbital elements and the ephemeris of a heliocentric object, which are given in the Astronomical Almanac, are determined by including the perturbations of the nearby planets in our solar system. Therefore the object's orbit is not Keplerian. Even with this Keplerian approximation, however, it was found that as long as the observations are good, the results of all these orbit determination and alpha,delta position-prediction programs ranged from good to excellent. Transfer.bas program This program computes a space vehicle's spin axis right ascension(SARA), and its spin axis declination (SADEC), when one specifies the drift rate in degrees per day and the space vehicle's minimum injection altitude above the earth's surface (sea level). This minimum injection altitude places the space vehicle at perigee in the transfer orbit. SARA and SADEC are essential because they give the direction the space vehicle must be fired for injection into the drift orbit at the apogee point of the transfer orbit. In addition, the transfer orbit's inclination,i, its angle of the ascending node and the drift orbit's inclination,i, and the angle of the ascending node for the drift orbit must be specified in degrees. Finally, delta v must be specified so as to place the space vehicle into a near geosynchronous orbit. This delta v is in feet per second and it must be less than or equal to the minimum delta v required. When this delta v is specified, this program automatically computes the minimum delta v required and the resulting outputs of this program are the spin axis right ascension(SARA), and declination(SADEC), and the elements for the near geosynchronous orbit. Often, the geosynchronous(24 hour) orbit cannot be exactly achieved, so therefore the term drift orbit is often used. The final resulting output of this program will be the semi-major axis, the eccentricity, and the ascending node. This drift orbit is very nearly circular, i.e., its eccentricity is almost zero which means omega, the argument of perigee, has really no meaning. This final output will verify that the space vehicle has achieved the drift orbit. The near geosynchronous or drift orbit is vitaly important for our weather and communication satellites. For more info about this, mail Neil at PATJACOBY@PRODIGY.NET.