Satellite altitude as a function of time?

In order to obtain a more accurate radiation estimate in this question, I am trying to approximate and/or determine the altitude of a satellite as a function of time. The satellite has an apogee of 32190 km above earth, and perigee of 320 km above earth.

I assume the following properties for simplicity:

• No orbit inclination, so in the equatorial plane,
• No J2 effects.
• No atmospheric drag/other drag forces.

If you see any mistakes please let me know :)

Approaches:

1. Find closed form analytical solution
2. Find iterative analytical solution
3. Find orbit simulation software to retrieve orbit altitude
4. Find a dataset of orbit altitudes as a function of time of equivalent satellites.

Summary of the attempt outcomes:

1. According to this attempt there is no (simple) closed form. I was trying to be effective, not to re-invent the wheel so I will switch to a numerical method to get a function of the orbit altitude.

2. Worked out in answers: The Orbital Dynamics Part 18 -- Formula for True Anomaly lecture explains how the the eccentric anomaly can be found by an iterative attempt. I was curious how this relates to the differential nature of the equation.

This solution is worked out in the answers.

1. I tried:

https://www.thanassis.space/gravity.html

http://en.homasim.com/orbitsimulation.php

https://nl.mathworks.com/matlabcentral/fileexchange/57132-satellite-orbit-transfer-simulation

The first 2 pages of these models:

http://www.winsite.com/satellite/satellite+simulator+orbit/

https://github.com/pytroll/pyorbital

https://nl.mathworks.com/matlabcentral/fileexchange/57132-satellite-orbit-transfer-simulation

4. Worked out in answers: Worked, though midnight resulted in a bug, so I restarted it and awaited a half orbit (from apogee to perigee or the other way round), since symmetry yields the entire orbit. The solution is worked out in the answers.

I am not sure if this is appropriate, but thanks everyone for the critical and supportive feedback! :)

• Orbit inclination ? Do you require J2 perturbation? – Prakhar May 2 '18 at 12:38
• Could you provide a link to the other question? Right now it just says [1] with no link. – barrycarter May 2 '18 at 16:51
• I'm voting to close this question as off-topic because it is extremely long-winded way to discover that there isn't a closed-form solution for the Mean Anomaly M (...is that mean anomaly? I forget) which also fails to understand that you can't just substitute in whatever you want for E(t). – Erin Anne May 3 '18 at 22:45
• @a.t. I can see you're putting in quite a substantial amount of work and prior research into this question! It's growing quite large now and is getting harder and harder to read. I propose you simplify it. To get the total dose per orbit with a given 1D dose vs altitude profile, all you need is the distribution of time spent at each altitude, in other words $dt/dr$. That might be a little simpler, and I think the math is interesting. – uhoh May 4 '18 at 9:19

Approach 4 (Succesfull approach)

At this point in time, I am not able to pursue the analytical approach any further, so I found a satellite (debris) via http://stuffin.space/?intldes=2012-008N that had an orbit similar to the one I am trying to determine, wrote a small excel script that crawls the data from http://www.satview.org/forec.php?sat_id=43273U once every minute, so that tomorrow, I can sort the data and get the altitude approximation as a function of time.

The script is:

Sub absorb_data_from_vdb()
Application.Wait (Now + TimeValue("0:00:10")) 'H:MM:SS
Dim ie As Object
Dim objelement As Object
Dim c As Integer
Dim i
Dim strImageSource As String
Dim img As HTMLhtmlElement 'tools =>reference Microsoft Html Object Library
Set ie = CreateObject("InternetExplorer.Application")

Dim electric_bullet_count As Integer
Dim gas_bullet_count As Integer

With ie
.Visible = True
.navigate "http://www.satview.org/forec.php?sat_id=43273U"
DoEvents
Loop
Application.Wait (Now + TimeValue("0:00:02")) 'H:MM:SS
whole_day = 2
Do While whole_day < 10
.Refresh
DoEvents
Loop
Application.Wait (Now + TimeValue("0:00:02")) 'H:MM:SS

Set elements0e = ie.document.getElementsByClassName("link_mudar") 'Get innertext to get value
If elements0e.Length <> 0 Then
'MsgBox (elements0e.item(2).innerText)
elements0e.item(2).Focus
elements0e.item(2).Click
End If

If ThisWorkbook.Worksheets("Sheet1").Range("A1").Value - 2 < 2 Then
ThisWorkbook.Worksheets("Sheet1").Range("A1").Value = 2
End If
For Iteration = ThisWorkbook.Worksheets("Sheet1").Range("A1").Value - 2 To ThisWorkbook.Worksheets("Sheet1").Range("A1").Value - 2 + 100
Set elements0e = ie.document.getElementsByClassName("texto_track2") 'Get innertext to get value
If elements0e.Length <> 0 Then
ThisWorkbook.Worksheets("Sheet1").Range("A" & 2 + Iteration).Value = elements0e.item(0).innerHTML
ThisWorkbook.Worksheets("Sheet1").Range("F" & 2 + Iteration).Value = Now()
'MsgBox ("found from")
End If
Next Iteration
For Iteration = ThisWorkbook.Worksheets("Sheet1").Range("A1").Value - 2 To ThisWorkbook.Worksheets("Sheet1").Range("A1").Value - 2 + 100
Set elements0e = ie.document.getElementsByClassName("texto_track2") 'Get innertext to get value
If elements0e.Length <> 0 Then
ThisWorkbook.Worksheets("Sheet1").Range("G" & 2 + Iteration).Value = elements0e.item(1).innerHTML
'MsgBox ("found from")
End If
If ThisWorkbook.Worksheets("Sheet1").Range("A1").Value < Iteration + 2 Then
ThisWorkbook.Worksheets("Sheet1").Range("A1").Value = Iteration + 2
End If

Next Iteration

proceed_boolean = False
Do While proceed_boolean = False
If Left(Right(Now(), 5), 2) <> Left(Right(ThisWorkbook.Worksheets("Sheet1").Range("F" & 2 + Iteration - 1).Value, 5), 2) Then
'MsgBox ("A minute has passed")
'MsgBox (Left(Right(Now(), 5), 2))
'MsgBox (Left(Right(ThisWorkbook.Worksheets("Sheet1").Range("F" & 2 + Iteration - 1).Value, 5), 2))
proceed_boolean = True
Application.Wait Now + #12:00:02 AM# 'This waits for 5 seconds
End If
Loop

Loop
End With
MsgBox ("Done")

End Sub


After two runs, the following graph is plotted for the equivalent orbit:

Visual inspection yields a clear orbit function, after removing the outliers. The following sanity check are performed:

1. the velocity at the perigee (lowest altitude) which should be highest, meaning the gradient of the orbit altitude function should be highest.
2. The velocity at apogee (highest altitude) should be lowest, meaning the gradient of the orbit altitude should be lowest near the perigee.
3. The gradient of the orbit altitude should equal 0 (a horizontal line at the apogee).
4. The Velocity profile is plotted separately and checked on continuity.
5. The velocity profile is plotted separately and visually checked on matching with the orbit altitude gradient.

All 5 verification are completed and successful, meaning no errors have been found in this approach.

Since the orbit has an axis of symmetry only half the velocity profile is needed. As one can see, two recordings were made, the right recording contains a complete half symmetric orbit, so it yields the complete orbit altitude function after a short data manipulation. This is the orbit data that will represent the orbit altitude as a function of time, for the radiation estimation purposes.

Next the data was manually filtered (by colourmarking the magnitude of the differences between the altitude datapoints and removing the largest differences untill I was satisfied with the differences.) yielding:

After filtering I used the following software/excel sheet to make a least squares fit: http://www.jkp-ads.com/articles/leastsquares.asp

The following polynomial function was used to approach the function**:

$y=ax^4+bx^3+c(x-f)^2+dx+e$

where x represented the time in minutes, and y was computed and compared with $y_{actual}$. The difference between the two y's was squared and summed for all datapoints (times and altitudes)scraped by the first excel). The boundary conditions (b.c.) were determined initially for the Evolutionary algorithm, by selecting the maximum value for which the highest altitude could be obtained by just that 1 coefficient multiplied by it's x, and then multiplying that value with 10. This was done to ensure the solution was obtainable, but limited the option space.

It yielded the following function**:

$$\begin{split} const_a = 4.52765E-06\\ const_b = -0.002389915\\ const_c = -0.074298946\\ const_d = 252.5709003\\ const_e = 260.8620904\\ const_f = -108.1781644 \end{split}$$

$h(t)$ is computed in $km$ $$h=Const_a*xValues^4+Const_b*xValues^3+Const_c*(xValues-Const_f)^2+Const_d*xValues+Const_e \label{altitude_orbit_as_func_time}$$

The Excel script can be further improved by adding an errorhandler closing the ie browser and restarting the script to ensure a higher level of automation, at the risk of greater data errors.

Error discussion: As was pointed out in the comments, the 4th degree polynomial fit will limit the accuracy of the orbit fit. (This is visible if you look at the start of the actual orbit measurement, it curves towards/around -t when going down, whereas the fitting polynomial only appears to have 1 curving radius (towards/around +t). This made the fit a bit challenging, since several iterations were required before a decent fit was found. However, this degree of accuracy was considered sufficient, by means of visual inspection. The procedure can be improved by a computation of the actual error squared.

**I think a better approach would indeed be to approximate the elliptical orbit altitude as a function of time with a sinusoidal.

• "As a next solution attempt I will compute the surface of the ellipse and the surface of the circular orbit, and then apply a constant angular area sweep rate to each of them separately." ...but that won't work. The angular sweep rate absolutely isn't constant for the ellipse. The constant, per Kepler's equal area law, is the swept area per time. – Erin Anne May 3 '18 at 6:14
• This isn't an answer to the question. While it seems helpful to use answer posts as scratch pads to show failed attempts, it isn't really allowed in Stack Exchange. – uhoh May 4 '18 at 9:20
• Thank you @uhoh, the result has been posted, so that it answers the question. I will simplify/clean up the results after my deadline. – a.t. May 5 '18 at 14:40
• Wow! Okay that's great news! Ya you should clean this up at some point; while it would be fun if Stack Exchange had scratch-pad posts, it really only has questions and answers. I'll give this a read tomorrow. In the mean time I can switch my vote to up. Thanks for your message! – uhoh May 5 '18 at 15:07
• A polynomial really isn't the right solution for this problem. The altitude of a satellite is a sinusoid. I'm worried about this being misleading for others who come across it; you rejected the standard iterative solution in favor of one that doesn't appear to give the apogee you selected. – Erin Anne May 10 '18 at 22:03

Aproach 2 (Succesfull approach)

The orbit is described in the following figure:

After a long, tedious, yet wonderful derivation the following 3 equations are applicable to the elliptical orbit:

$M= \frac{t}{T}2 \pi$

$M = E-\epsilon sin{E}$

$r = a(1-\epsilon cos {E})$

$\upsilon = \theta = 2 arctan (\sqrt{\frac{1+\epsilon}{1-\epsilon}} tan{\frac{E}{2}})$

$T=2\pi\sqrt{\frac{a^3}{\mu}}$ with $\mu_{earth} = 398600.4415 \frac{{km}^3}{s^2}$

$r = \frac{a(1-e^2)}{1+e cos{\theta}} = a(1-\epsilon cos{E})$

$M(t)$ = Mean anomaly (=Angle between green line, center of ellipse, dashed line between center of ellipse and perigee)

$E(t)$ = Eccentric anomaly (=Angle between red line, center of ellipse, line between center of ellipse and perigee)

$\upsilon(t) = \theta(t)$ = True anomaly (=Angle between black line, center of ellipse, line between center of ellipse and perigee)

$a$ = semi-major axis

$r(t)$ = altitude of satellite

$\epsilon$ = eccentricity

The anomalies are visualized in the following orbit animation. Note that, as was pointed out in the comments, M(t) is constant and unequal to E(t):

Furthermore, it is necessary to compute the orbital period $T$, semi major axis $a$ and orbit eccentricity $\epsilon$ once.

Given that the perigee height = 320 km, and the apogee height 32190 km, the semi major axis is found by using $r_e = 6371 km$:

$a = \frac{320+r_e+32190+r_e}{2}= 25971.5 km$

Now the eccentricity is found using $\theta = 0 rad$:

$r = \frac{a(1-e^2)}{1+e cos{\theta}}$

$r = \frac{a(1-e^2)}{1+e} = a(1-e)$

$e =1-\frac{r}{a} = 1-\frac{r_e+320}{a} =1-\frac{6371+320}{25971.5}=0.74237144562308$

The orbital period is found with $\mu_{earth} = 398600.4415 \frac{{km}^3}{s^2}$:

$T=2\pi\sqrt{\frac{a^3}{\mu}}=41685.4 s$

So to determine the orbit at a point in time t (e.g. 0.5 $s$), the following computation steps are made for every point t at which we want to know the satellite altitude:

The mean anomaly M(t) is computed. $M= \frac{t}{T}2 \pi$

$M= \frac{0.5}{41685.4}2 \pi=0.0000753643.. rad$

Then M is substituted in:

$M = E-\epsilon sin{E}$

This leaves 1 unknown: the eccentric anomaly $E$. This is find with an iterative process (Just assume a first E, then see what you compute for $M$, check if it equals the actual computed $M$ (probably not), so then lower or increase your estimate for $E$, and see if you stray further away from $M$, or get closer to it, untill you are satisfied with the degree of accuracy of $M$, (Which implicates the degree of accuracy of $E$).

Once E(t) is determined with a sufficient degree of accuracy, it is substituded in:

$r = a(1-\epsilon cos{E})$ to find $r{t}$.

That is the orbit altitude at time t = 0.5 s. The procedure is repeated for another time t by recomputing $M(t) again. On top of the altitude, one might also want to know the location of the satellite, relative to earth, by knowing the true anomaly. To determine the true anomaly at time t, one substitutes the iteratively found$E(t)$in the following formula:$\upsilon = \theta = 2 arctan (\sqrt{\frac{1+\epsilon}{1-\epsilon}} tan{\frac{E}{2}})$This is the documentation of the procedure explained in 28:55 to 30:16 of: . (I did not write a script for this {iterative loop of determining/guessing$E\$} (yet) so it is not semi-automated, unlike the answer of approach 4.)

• Congratulations on your progress! It looks like you've lost your images. Also at some point you can consider accepting this answer instead other, if you feel this better answers your question. – uhoh May 8 '18 at 17:13
• Yes, thank you, I must say it was pretty cool to find/understand the analytical iterative solution after all. If I manage to automate the guessing of the eccentric anomaly E(t) I will accept this answer, but right now, the least effort is using solution 4. Especially with the conversion of measured data points to function by means of least squares included. – a.t. May 8 '18 at 18:44