Apogee

Preface, about the calculation method

This problem is of "didactic" nature - to implement certain improved algorithm for calculating the path of the moving body - satellite of a planet. So regretfully we are to start with the dull algorithm explanation :)

We have already calculated motion equations in the problems like Safe Landing and Billiard Ball. Let's recap the "general form" of such task:

1. We have some moving object(s) with position described by one or more coordinates (it could be just height or X, Y etc)
2. Object(s) has velocity(es), which tends to be constant, unless the object is affected by some forces (like gravity or friction).
3. Forces sum up to create acceleration, which can increase or decrease velocity, or change its direction by and by.
4. The difficulty is that forces themselves often are not constant, and may depend on position or velocity.

To "track" how the object(s) moves under these circumstances, we use very clever method named after Euler:

• work in iterations, making small steps in time, choosing each step duration dt (e.g. 1 second)
• on every iteration we calculate next position according to current velocity
• also we calculate next velocity according to current acceleration (e.g. accelerations arising from forces based on current position and current velocity)
• now reassign next position to become new current position and next velocity to become new current velocity.

Let's regard an example. The stone is thrown upwards with speed 20 m/s. It is affected by the gravity force creating constant acceleration g = -10m/s (minus sign means it is directed downward, counter to speed). What height the stone shall reach, when its speed becomes zero?

Iterating with dt = 1 sec we do two iterations, starting with height Hcur = 0 and Vcur = 20:

1st iteration:

    Hnext = Hcur + Vcur * dt = 0 + 20 * 1 = 20
    Vnext = Vcur + g * dt = 20 + (-10 * 1) = 10

    so at iteration end Hcur = 20, Vcur = 10

2nd iteration:

    Hnext = Hcur + Vcur * dt = 20 + 10 * 1 = 30
    Vnext = Vcur + g * dt = 10 + (-10 * 1) = 0

    so we reach peak point with Vcur = 0 and height of 30 meters

But this is quite not precise! Sure, because Vcur really decreases during iteration itself, it should be an overshot. Really, if we calculate this with dt = 0.5 we'll only reach height of 25, and decreasing iteration time steps further we'll soon figure out the real target height is only 20 meters!

How to deal with it? Obvious idea is to reduce calculation steps, but in some situations this increases the overall calculation time, which could be bad for "real-time" situations (e.g. control systems of real spacecrafts etc).

One cunning method, called Adams-Bashforth (feel free to look it up in google) is using the values not only from current step, but also from step one iteration ago! It adds "effect" of values of the current step with coefficient 1.5 and subtracts "effect" of values of the preceding step with coefficient 0.5:

Hnext = Hcur + Vcur * dt * 3/2 - Vprev * dt * 1/2
Vnext = Vcur + Acur * dt * 3/2 - Aprev * dt * 1/2

Here Acur and Aprev are values for total acceleration on the current and previous step. Applying it to our example we know acceleration is constant (it is -g) and so formula for Vnext is actually unchanged. However formula for height shall give different values.

We only need to decide what is Vprev on start. If we are going to calculate many steps, we can leave it equal to Vstart since error will be negligible. In our case the process is very short, so let's cleverly assume it should be Vstart - g * dt = 30.

1st iteration:

    Hnext = Hcur + Vcur * dt * 3/2 - Vprev * dt * 1/2 = 0 + 20 * 1.5 - 30 * 0.5 = 15
    Vnext = Vcur + g * dt = 20 + (-10 * 1) = 10
    Vprev = Vcur = 20

2nd iteration:

    Hnext = Hcur + Vcur * dt * 3/2 - Vprev * dt * 1/2 = 15 + 10 * 1.5 - 20 * 0.5 = 20
    Vnext = Vcur + g * dt = 10 + (-10 * 1) = 0

So magically we get correct value even with very large step. You can decrease step to dt = 0.5 or even 0.1 but the result is the same. This formula gives correct results for "second degree equations" (and ours is one of them).

Satellite trajectory calculation

This October (in 2022) we commemorate the beginning of the "Space Age" 75 years ago. Scientists and Engineers of the Soviet missile development team persuaded government to spare few rockets to launch some very naive device into the orbit around the Earth. It was back in 1957 and the device was using vacuum tubes, doing nothing besides giving out constant beep-beep-beep on frequencies available to radio-amateurs all over the world. Their example soon was followed by US scientists, launching some more clever satellite, packed with instrument payload, built upon semiconductor transistors (invented less then 10 years before). Thus people on both sides made great job in diverting military interests into more peaceful, constructive and scientific path of space research.

We are going to trace the satellite's path around the Earth! From Kepler's times we know it is going to be closed path in the shape of the ellipse. However, if we use Euler's method, we'll see the path won't "close"!

Let's go into details. Suppose the earth is perfectly round with radius of Re = 6371 km. The satellite has just been launched by the rocket to the height of 500 km and we start "tracing" from this point. It has the linear horizontal speed of 8000 km/s. No forces affect the satellite except gravity (engines exhausted the fuel and are shut down). The gravity has acceleration of g = -9.81 at the earth surface, but decreases with height according to the law of squares, e.g. at the height H it is g * (Re / (Re + H))^2.

For convenience we have (below) simple demo program in BASIC which implements Euler's method to calculate one turn over the globe. We shall see that in the end the calculated height is quite wrong, different from that at start.

Problem Statement

You'll get starting height and speed of the satellite - and the goal is to calculate its apogee height and when it is reached, using Adams method (recommended step is 1 sec). The "apogee" point of trajectory is the one with maximum height. Since we are supposed to start in "perigee", the calculation just should check height after every step and stop when the new height becomes less then previous.

Input data are just two values, starting height (in km) and velocity (in m/s).
Answer should be the apogee height after a hour of flight, rounded to the nearest integer (in km) and time in seconds to reach it (effectively the number of steps).

Example:

input:
200 7900

answer:
605 2775

This is not the only way to calculate answer - one may try Euler with reduced steps or calculate apogee analytically - but as the "checker" for the problem uses Adams itself, other methods may not exactly match the expected answer.

Sample code with Euler's method

Study this simple code, hopefully comments (starting with rem) may help:

re = 6371000 : h = 500000 : g = 9.81 : rem Earth radius, initial height and gravity
x = 0 : y = re + h                   : rem Start coordinates
vx = 8000 : vy = 0                   : rem Start speed, by its orthogonal components
dt = 10                              : rem Time-step size (10 sec)

loop:
r = sqr(x^2 + y^2)                   : rem Full radius to Earth center
a = -g*(re/r)^2                      : rem Respective acceleration due to gravity
ax = a*(x/r) : ay = a*(y/r)          : rem Acceleration split into orthogonal components

vxn = vx + ax * dt : vyn = vy + ay * dt : rem Euler step for Vnext, by components
xn = x + vx * dt : yn = y + vy * dt  : rem Euler step for Xnext and Ynext

xp = x : x = xn : y = yn             : rem Update current X and Y
vx = vxn : vy = vyn                  : rem Same to current Vx and Vy

rem print x, y, sqr(vx^2 + vy^2)     : rem Uncomment this to print position at every step

if xp < 0 and x >= 0 then goto done  : rem Stop when step started at negative X and ended in positive
goto loop                            : rem Otherwise continue iterating

done:
print "end height: ", sqr(y^2 + x^2) - re

Run this code (i.e. copy-paste to "solution" area and click "Basic" button below) - and you'll see the satellite ends at the height of 1200 km instead of 500 from which it started! Reduce time step to 3 and see it is still 727 km. Even with step of 1 sec we get 576 km.

The task is easy but subtle mistakes are possible. Here we provide a list of coordinates on few steps after start (for the example given above with h=200 v=7900) which we are assume should arise with correct calculations.

0 0 6571000
1 7900 6571000
2 15800 6570986.1671221
3 23699.97505414 6570963.1123707
4 31599.916847142 6570930.8357633
5 39499.814292059 6570889.3373525
6 47399.656302011 6570838.6172018
7 55299.431790219 6570778.6753887
8 63199.129670023 6570709.512005
9 71098.7388549 6570631.1271564
10 78998.248258485 6570543.5209628
11 86897.646794592 6570446.6935581