Orbital Simulation

Program: RMorbits.z80

I decided to write a simulation of a spacecraft in motion around a planet, because I thought it would be interesting to model numerically what can seem counter-intuitive at first: i.e. an orbiting object is always falling, but does so without hitting the ground! I had already played around with this on a programmable calculator I'd borrowed, a little while before I got my Spectrum+, but this had just produced numbers, and I wanted to see a picture! The above image shows the result of one simulation where the orbit is supposed to be circular. Because of accumulated errors during the simulation, the start point and end point do not quite meet here, but this is configurable. The green blob at the centre is supposed to represent a planet! 

The user is prompted for the starting x, y co-ordinates of the space craft, and the x, y components of the orbital velocity (note - use 'crtl-M' on a PC keyboard to input a decimal point).

I've also provided thrusters so that the pilot can adjust the path taken. A prompt is given for the (constant) acceleration from the thrusters (labelled somewhat confusingly 'FR', 'FT', suggesting I was thinking about force and not acceleration when I wrote the program), although I have chosen to prompt for these in radial co-ordinates, one component being outwards along a line from the planet to the spacecraft and the other being along a line at right angles to the radial thrust in an anticlockwise direction, the idea being that the spacecraft is always suitably oriented for this.

Finally, the user is prompted for an interval of time which elapses before acceleration, a new position and a new velocity are calculated*. The smaller this interval, the more accurate the simulation. In the above image, this was set to 0.01 - setting it to 0.001 would have given just about a perfect circle.

Once the simulation starts, the position, velocity, and time are updated at the right of the display.

Note that the closer the craft gets to the planet, the faster it goes and the greater the distance it travels in each time interval. As a result, the simulation becomes less accurate. For this reason, if the craft gets too close, then it is deemed to have crashed and the program stops.

Thrusters are activated during the simulation by pressing the 'f' key (for 'fire') and you can interrupt the program to change the thruster settings by pressing the 'i' key. Also, you can stop the program by pressing the 's' key.

One point of note about this program is that all units are normalized - i.e. the planet's mass is 1 unit, the gravitational constant G is 1.

On reflection, the program could be improved by simulating the motion of the spacecraft just when the thrusters are fired, and at other times using the position and velocity to determine the shape of the orbit, using Newtonian mechanics.

Below is an image after about 1.5 orbits. The starting point was (1,0) with velocity (0,1) as before, although this time I chose a smaller time interval, 0.001, for a more accurate depiction (although that also contributed to the time display overflowing onto the next line!) After half an orbit, the (tangential) thruster was fired for a short while. This has changed the path of the craft, so that after one full orbit, the height of the craft is perhaps 25% higher than before. After travelling on, at 1.5 orbits, the craft has returned to its original height, so we have an elliptical orbit with apogee ~1.25 and perigee of 1. The elliptical orbit between smaller and larger circular orbits is known as a Hohmann transfer orbit. To raise the perigee of the orbit, the thrusters will need to be fired again, but this time at apogee.

Note:

• * ignoring the thrusters, acceleration, position and velocity are determined for each time interval dt. Acceleration is calculated first as e.g. (ax) = -x/(r cubed) where r is the distance to the planet (remember we are using normalized units, so there is no gravitational constant (G), nor mass of the planet to consider). Then the change in position is given by e.g. delta x = (vx)*dt + 0.5* (ax)*dt*dt, where vx and ax are the components of velocity and acceleration respectively along the x axis. Then the change in velocity is given by e.g. delta (vx) = (ax)*dt.

• unfortunately, a lot of lines towards the end were missing from my tape backup, but luckily, I still had a print-out of the program. Rather than loading the program and laboriously keying in the missing lines using the Spectrum editor, I entered the text in an ASCII file and then converted it to a TAP image file using a utility called bas2tap. Then I used the MERGE statement on a Spectrum emulator to append the new text to the old program. In doing so, I even attempted to improve somewhat upon the old program - my first bit of Spectrum programming in the new millennium!

• in the picture below, produced by an earlier version of the program, the time display should show '6.44' units (using an interval of 0.01 units), but instead shows '6.4399999' and has flowed onto the following line at some point (actually at 3.27/8 units). (I should point out that exactly the same behaviour occurred on the actual machine when I first wrote the program). All the program does here is add the time interval to an accumulator multiple times. I put this strange behaviour down to a coding quirk in the Spectrum's ROM (the z80 processor having no native support for numbers with decimal places), and have since changed the code slightly to work-around the issue. The following program demonstrates the issue:

100 LET tt=0

110 LET t=0.01

120 PRINT tt

130 LET tt=tt+t

140 GOTO 110