This is a program to calculate the torque on a Keplerian elliptic orbit (test orbit) by another Keplerian elliptic orbit (field orbit). I wrote it for a configuration where there are many field orbits but only one test orbit. I modified the code a little to make it more general, but still the design reflects this origin. The code is written in C and has some subroutines which may be useful on their own. They are described below.
Download the tarball, untar it into a directory and run make. This should create an executable kep_torque. You can run it and see what it does. Please see the terms for usage and distribution below. Note that in its current form, this program is not very useful. I expect other people to use this code as part of their programs rather than by itself.
The method is described in Gürkan and Hopman (2007). I gave a talk about it in the Near Keplerian Systems workshop in Leiden. Here are links to OO Impress (with OOoLaTeX) and PDF files of that talk.
Briefly the method works as follows. The test orbit is in the x-y plane with its epicenter on positive x-axis. The field orbit has a random orientation. I represent the orbits with solid wires. The density of a given portion of these wires is inversely proportional to the time spent by the star at that portion. I discreticize these wires by equal mass points, distributed equidistantly in the mean anomaly of the orbit. Initially 64 points on each orbit are used to calculate the torque. This is done twice, with two sets of points that are in the middle (in mean anomaly) of each other (see the figure below). If the relative difference of the torques calculated by these sets is more than the required tolerance, the number of points is doubled and the calculation is repeated. At most 65536 points are used on each orbit. The computation is terminated at this point even if the desired accuracy is not reached.
The main data structure of the code, orbits, are stored as structs which contain doubles for mass, shape and orientation, an unsigned int for number of points in the orbit and a pointer to pointer to double to hold the position of points in the orbit.
typedef struct {
double m; /* mass */
double a, e; /* shape */
double theta, phi, gamma; /* orientation */
unsigned Np; /* number of points */
double **p; /* the points (each a triplet) */
} Orbit;
The points do not traverse the orbit linearly. The first point (index 0 in C) has mean anomaly M=0; the second point has M=π; the third and fourth have M=π/2 and M=3π/2, respectively and so and so forth. The mean anomaly values (divided by 2π) corresponding to the points are:
point index | mean anomaly / 2pi -----------------+----------------------- 0 | 0 1 | 1/2 2, 3 | 1/4, 3/4 4, 5, 6, 7 | 1/8, 3/8, 5/8, 7/8 8, 9, 10, ... | 1/16, 3/16, 5/16, ... 16, 17, ... | 1/32, 3/32, ... . . | . . . . | . . . . | . .
In this arrangement, the points with indices k to 2k have angles in between the points with indices 0 to k (Note that I use C convention for referring to indices: "0 to k" means integers from 0 up to but not including k). The mean anomaly values for this indexing trick are calculated and stored in an array Ma by the routine precompute_Ma(), which is called near the beginning of the program. This is only nontrivial part of the program. The rest is pretty straightforward (maybe too much so, see notes on efficiency below). It consists of initialization routines and a main loop.
Initialization: Here we assign values to various variables used in the code, initialize the random number generator, prepare an array containing mean anomaly values to be used (see above), allocate memory for orbits, and assign values to coordinates of the points that form the test orbit. In this program, we calculate the torque from a single orbit and it is possible that not all the points of the test orbit will be necessary. However, when there are many field orbits it is simpler to determine the orbit once at the beginning of the main loop.
Main loop: We first set parameters for the field orbit. We can either assign values to variables and make a call to set_orbit_params(), similar to test orbit; or make a call to set_rnd_orb_params(). The latter method makes further calls to other routines to pick random semimajor axis, eccentricity and orientation angle values.
Then, we make a call to calc_torque_on_test_orbit(). Here is a breakdown of this routine:
The call to set_orbit_points() for the field orbit above will involve two rotations. First, we choose a direction for angular momentum, determined by two angles: φ (azimuthal angle) and θ (polar angle). The first rotation is by an angle θ around the vector n0= (cos(φ), sin(φ), 0). This sets the inclination and right ascension of the ascending node. Then we rotate the orbit around its angular momentum vector (i.e., in the orbital plane) by ω. This sets the argument of pericentre. For a random orbit, φ and ω are sampled from a uniform distribution in the interval [0, 2π) and cosine of θ is sampled from a uniform distribution in [-1, 1]. Note that the rotations given here are different from the rotations given in Gürkan and Hopman (2007), I think these are easier to understand.
The calculation of torque on a bunch of points from a bunch of points is done by the routine calc_arc_torque(). It calls calc_f_from_arc() for each test orbit point to calculate the force from a bunch of field orbit points and then takes the cross product to calculate the torque and adds things together. calc_f_from_arc() is a little optimized while sacrifiying clarity. If you try to make this more clear, make sure you do not make it much slower.
E_from_M_e(M, e): This routine solves Kepler's equation for an elliptical orbit, M=E-e sin(E). That is, finds eccentric anomaly given the mean anomaly and eccentricity. It is based on equation (6.6.7) in Danby (1992). If this routine is your bottleneck, see Fukushima (1999) and references therein.
rotate3_xy(x, y, theta, phi, omega, *newpos): This routine takes a point in x-y plane with coordinates x and y; carries out the rotations described above; and stores the new x, y, z values in newpos[0], newpos[1] and newpos[2], respectively. I used the formula given in MathWorld and MAXIMA to obtain the formulae for rotated coordinates. You can find the expressions in the source file.
There is a lot of room for improvement here. The code can be made both more efficient and more robust. Here are some ideas:
I did not try any of the above and would be curious to hear if you did. Also, if you have any ideas for further improvement, I would like to hear those as well.
This software is provided to help you get the job done. Hopefully to solve a problem, but possibly as a nice pastime. There are no restrictions whatsoever on the usage of this software. You can modify and adapt it to any use you see fit. However, I cannot be held responsible for any damage arising from using this software, including but not limited to wasting your time. Also, distributing the software, even parts of it as part of another software, is not usage. Such actions are governed by "Terms for copying and distribution", as indicated below.
I appreciate, but do not require, attribution. If you are unsure about what reference to use for attribution, feel free to ask. For torque calculation work, the appropiate reference would be Gürkan and Hopman (2007). Note that, if you simply use this software to obtain some results to write a paper, you are not required to offer me co-authorship or such. Only if I provide extensive help to you for using the software, and hence I contribute significant amount of my time specifically for your project, I would like to be offered co-authorship. Do not let this to discourage you from asking questions, but please read the documentation (at least this page and possibly our paper) before doing so.
This software is copyright (c) 2007 by Atakan Gürkan, except where noted in the source files.
This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License version 2 as published by the Free Software Foundation.
Note that I leave out the phrase "... or any later version", that is, you cannot merge this software with software under different versions of GPL. If you need to do this, please contact me and if the version you propose is reasonable for me, I will make an addition or exception to these terms.
This software is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
Parts of this software may be provided with licenses different from but compatible with GNU General Public License version 2. All such parts are clearly marked in the source code files. For those parts, you are free to use GNU General Public License version 2 or the particular license indicated in the source code.
If you make an improvement to this software, I would appreciate if you let me know about it. I also hope that you would allow me to further distribute your enhancements as part of this software, under the terms of this license. Note that, you do not need to transfer your copyright to me to make your improvements part of this software.