This month’s post is long and quite complex. So I will split it in two parts: Input & Output. I also decided to write an Extension post to the series, where I will discuss binary and multiple star systems, and things related.
Designing an extremely realistic planet (along with the whole system) is not a simple task, so if you have no interest in it, I suggest you don’t read math & tech parts of this article (I will mark them as such). We’ll be doing some math here and there. Also, advanced level of computer knowledge is required. If you don’t know what a command line is, skip those parts and enjoy the rest. And last but not least, this task will require a few days of your time (mostly computer time though). So, if you have something to do, like mow grass on the football field, now is a perfect opportunity. With all this said I’m moving onto the main thing.
Countless planetary systems exist throughout the vastness of space. Some of them are small, consisting only of one planet per star; some of them are big, like our Solar System with 8 planets and lots of other stuff orbiting around the Sun. In fact, all planetary systems discovered up until now consist of smaller number of planets than the Solar System, with 2 to 4 planets being most common.
Before designing your own system it might be a good idea to look through descriptions of those systems. There are numerous catalogues with data. For example:
For more planet stuff you can also read the Systemic blog, which is written by Greg Laughlin, Professor of Astronomy and Astrophysics at the University of California, Santa Cruz. Systemic reports recent developments in the field of extrasolar planets, with a particular focus on observational and theoretical astronomical research work supported by the NASA Astrobiology Institute, and by the Foundational Questions Institute (FQXi).
Now, there are two ways of getting your system up: top-to-bottom or bottom-to-top methods.
The top-to-bottom method suggests using some of the accretion algorithms to generate your system – you’ll get all stuff automatically, nothing to invent there. And I’m not talking about “planetary system generators” found on the web (e.g. StarGen), they don’t produce realistic systems. They tend to produce systems “similar to our own”; the truth is planetary systems don’t follow this “similarity” plan. A truly alien world will be something else, as you can see from worlds discovered up to date. Personally, at some point I’ve been using AstroSynthesis (v.3.0 will be released somewhere after September 30th). It does cool things, try it. I recommend it to those who don’t want to go deep into calculations.
State-of-the-art planetary accretion simulations typically use N-body integrations, but these are too computationally expensive and time consuming. So we’ll skip that part and get to the time when all planets are already in their respectful places. Trust your intuition; it’s better than some generator.
The bottom-to-top method implies you first imagine your main planet and do all the necessary calculations for its parameters (I will talk about those in detail below). The rest of the system is rather sketchy: all you decide is where to put other planets (semimajor axes) and how many there will be in total (you can also add some of their details as well, e.g. mass, density and orbital elements). This decision should be based on your star of choice and its parameters.
Planetary parameters
I assume you have already decided on mass, density and calculated radius. If not, check the previous post here.
Another two parameters required for our computations are axial tilt and rotation rate. The tilt can be anything you want from 0 to 180°. Note that object with tilt below 90° has a retrograde rotation, like Venus does. Also, planets usually move on their orbits in counterclockwise direction, with clockwise direction being retrograde orbital motion.
There are three cases of planetary rotation around its axis: normal (counterclockwise), retrograde and tidally locked (synchronous rotation or orbital resonance, e.g. Mercury’s 3:2 spin-orbit resonance).
When planets are born, they have their primordial spins. During their lifetimes planetary spins are affected by solar tides and tides inflicted by another bodies of the system. Planets and moons slow down, spin up (like Mars and Saturn did), or become tidally locked if too close to a star or any other massive body (e.g. Mercury, the Moon). Earth’s initial spin rate was about 6.5 hours comparing to modern 24 hour day.
Math Box 1 – Initial rotation of a planet
Approximate initial planetary spin can be estimated from its mass.
First, you calculate initial angular velocity using Dole’s equation (output is in rad/s). M is the mass of the planet; k is the moment of inertia ratio (I/MR^2), a value of 0.33 is taken for a terrestrial planet (you can look at planet moments of inertia in these fact sheets); r is planet’s radius, and j is 1.46E-20 m²/s²kg.
ω = [(2*j*M) / (k*r^2)] ^ 0.5
Then, initial rotation rate in seconds is calculated as t = (2*PI)/ω
Now, from this point onward things start to slow down, because a tidal deceleration torque from the primary is applied. (On some occasions things might actually spin up.)
Here (dωe/dt) = (-1.3E-6 /10^9) rad s-1/yrs; k is as mentioned above for terrestrials; m is a planetary mass; index letter e stands for Earth and p for planet.
You can estimate planetary angular velocity at any age of the system. In the following formula system age is in seconds. The new angular velocity then becomes
ωnew = Initial angular velocity ω + (system age (s) * tidal deceleration torque dω/dt)
Rotation rate in seconds is calculated >> t = (2*PI)/ωnew
In hours >> th = t/3600
Thus you can determine how your planet day becomes longer. However, the primary isn’t the only one influencing rotation rate of a planet. There’s also atmospheric drag, and a small influence of the system’s bodies. (Some of the formulas above are taken from Fogg, 1992.)
Note: If you have a moon or two, things get more complicated. You’ll have to take them into account as well.
Note 2: I suggest you calculate rotation rate this way instead of integrating the whole system for considerable amount of time; that could take days on an average machine, especially if precision and long timespan (more than 100×10^8 to X × 10^9 years) are chosen.
Angular velocity also gives us a clue to planetary flattening and difference between polar and equatorial radii.
Math Box 3 – Planet flattening (dynamical ellipticity) and radii
The radius you determined from mass and density is the mean radius, R. If not, the time is now:
R=0.5*[(mass/density)*(6/PI)]^(1/3)
or, if given mass and log g, through relation
Mstar/Msun=((Rstar/Rsun)^2)*10^(loggstar-loggsun).
Ellipticity e = SQRT(f*(2-f)) is related to flattening, which can then be expressed as f = 1-SQRT(1-e^2) or f = (1/SQRT(1-e^2))-1; or, if you have only radii, flattening becomes f=(Re-Rp)/R.
Flattening f is also calculated as f* = (5/4)*(((ω^2)*(R^3))/(M*G)); ω – from calculations above (math box 1 & 2); R – mean radius; M – planetary mass; G = 6.67259E-11. Equatorial radius thus becomes Re = R*(1+(f/3)), and polar Rp = R*(1-((2*f)/3)). The difference between radii is expressed as Re-Rp, or simply as mean R*f.
*Note: We’ll get somewhat overestimated degree of rotational flattening of the planet because the formula takes planet’s interior as incompressible fluid of uniform density. In reality terrestrial planet’s core is much denser than its crust. Or, in case of gaseous planets, an overestimate will rise from the fact that they have a mass distribution which is strongly concentrated at their cores.
Now, when we have our rotation rates, we can move onto the next step.
Orbits and orbital elements
The traditional orbital elements are the Keplerian elements. They describe the orbit of the body in space.
Keplerian elements can be obtained from and converted to orbital state vectors (x-y-z coordinates for position and velocity) by manual transformations or with computer software. For our integration Keplerian elements should be enough. If you want something else, the body-centered, inertial rectangular components of the radius vector can be determined from the classical orbital elements (see math box 4).
Math Box 4 – Orbital State Vectors
The mean longitude λ = M + ω + Ω where M is the mean anomaly, ω is the argument of perigee, Ωis the longitude of the ascending node.
The true anomaly ν is found from the expression ν=λ-Ω-ω
The body-centered, inertial rectangular components of the radius vector can be determined from the classical orbital elements as follows:
Position formulasVelocity formulas
Here a is the semimajor axis, e is the eccentricity, i is the inclination; in these equations p is called the semiparameter of the orbit and is calculated from p=a*(1-e^2). η is the gravitational constant of the primary or central body.
Math Box 5 – The spin angular momentum*
*this will be necessary for the main integration part to calculate obliquities and spin rates. If you won’t do integration, skip this part.
SAM is a vector, it has 3 components. It is calculated as H = H * Ĥ
H = 0.4 * M * (R^2) * v
Here, R is planet’s equatorial radius, M is planetary mass, and v is a spin rate (angular velocity).
The initial coordinates of unit vector Ĥ in the reference frame are:
X = sin ψ * sin ε
Y = cos ψ * sin ε
Z = cos ε
Where ψ is precession in longitude and ε is axial tilt of a planet.
The precession in longitude ψ is defined by ψ = L – Ω, where Ω is the longitude of the ascending node, and L is the inclination of ecliptic. The angle between the planet’s equatorial plane and inclination of ecliptic is the obliquity ε.
Your both precession angle ψ and obliquity ε should be in radians.
**You can calculate your precession rate as well. Note that precession formulas use semimajor axis (a) and eccentricity (e). You can start integration at the beginning of the system age (if you have time and patience) to get all things neatly; or you can skip to the necessary age and use the “invented” data as initial point. Whatever works for you.
The final step is to multiply each component of Ĥ by the quantity H.
Hx = H * X
Hy = H * Y
Hz = H * Z
Integration will require these quantities in AU^2/day, so take mass as Mass/Mass of the star, R of a planet as R/149597870.7, and angular velocity in radians/day.
TECH TOOLBOX: Environment and Integrator
Note: Numerous programs exist to perform integrations of n-body systems. If you want to learn more, google something like “n-body integrator”.
Mercury is written in Fortran77 in a slightly non-standard code and before you can use it, you have to compile it. So, first of all, you need a compiler. I used GNU gfortran for this purpose. Depending on your operating system there are some options:
Assuming you installed all necessary things, we now can move onto our main goal – Mercury package.
The package has the manual, the file called mercury6.man and can be opened with any text editor of your choice. Read it CAREFULLY and follow all the instructions.
If using gfortran, compile the drivers typing lines
gfortran -o mercury6 mercury6_2.for
gfortran -o element6 element6.for
gfortran -o close6 close6.for
To run what you have compiled, just type
./mercury6
and hit enter.
Pretty simple.
If you have questions regarding installation and usage, feel free to ask. But please note that I’m running all this under Windows XP and might not know about any issues with other operating systems.
So far I had no issues with any software mentioned here.
Integration of orbits
When integrating the equations of motion in classical or relativistic celestial mechanics, one has to know the values of the parameters, in particular the masses and the initial conditions (orbital parameters or state vectors). I won’t be going into the n-body problem details here. If you are interested I suggest reading some good literature on the subject, e.g. Gravitational N-body Simulations: Tools and Algorithms by Sverre J. Aarseth (a cheap copy can be found on abebooks.com or on any of their international sites).
So far all necessary info on required parameters is in the Mercury manual.
Are planetary systems stable? Probably. You can’t tell just by looking at them; the invisible force out there guides all the interactions between the bodies. That’s why scientists perform numerical integrations of planetary orbits.
#For the next part I’ll do a small sample system integration.
JENO MARZ is a science fiction writer from Latvia, Northern Europe, with background in electronics engineering and computer science. She is the author of two serial novels, Falaha’s Journey: A Spacegirl’s Account in Three Movements and Falaha’s Journey into Pleasure. Marz is current at work on a new SF trilogy. All her fiction is aimed at an adult audience. View all posts by Jeno Marz →
