Discontinuity in Satellite Orbit

[Lance Simms]

Hello,

I am trying to simulate the orbit of a CubeSat satellite with the following TLE File (stored as TLE_File):

CXBN
1 99999U 11158.66666667 .00001767 00000-0 83150-4 0 00005
2 99999 060.0472 001.9175 0286684 357.7762 001.3936 14.89594857000013

I use the following code to start and propagate the orbit:

  DateTime = '2012/01/01 20:20:00'
  #Format the time appropriately
  Year = int(DateTime[0:4])
  Month = int(DateTime[5:7])
  Day = int(DateTime[8:10])
  Hours = int(DateTime[11:13])
  Mins = int(DateTime[14:16])
  Secs = int(DateTime[17:19])
  DateTime = datetime(Year, Month, Day, Hours, Mins, Secs)

  #Convert this into a PyEphem Date (# days since noon on 1899 December 31)
  EpDate_0 = ephem.date(DateTime)
  EpDate = float(EpDate)

  #Call getsat to return an ephemeris satellite object based on the TLE
  CXBN = getsat(TLE_File, EpDate)

  #Set the time increment for the propogation
  dt = double(0.01) #Time difference in seconds
  t_final = 1*1. #Final time in seconds
  NumTimes = t_final/dt
  Times = arange(NumTimes)

  #Do the progagation
  for i in Times:

    #Inrement the time
    EpDate = double(EpDate) + ephem.second*dt
    Time += dt
    #Recompute the solution with the SGP4
    CXBN.compute(EpDate)

    #Get the position of the satellite in Earth-Centered Inertial Coordinates
    r_eci = getSatECI(CXBN, c)

------
Here are the functions used above

def getSatECI(sat, c):

  #First get the RA and DEC from satellite object
  RA_rad = float(sat.ra)
  DEC_rad = float(sat.dec)
  Radius = float(sat.elevation+c.EARTH_R)

  #Now convert to ECI
  x_eci = Radius*cos(DEC_rad)*cos(RA_rad)
  y_eci = Radius*cos(DEC_rad)*sin(RA_rad)
  z_eci = Radius*sin(DEC_rad)

  #Form a three vector and return it
  r_eci = vec3(x_eci, y_eci, z_eci)

def getsat(tlefile,date):
    "Get the a pyephem object representing a earth satellite at a particular date, given a TLE file describing the orbit"
    f = open(tlefile,'r')
    tle = f.readlines()
    f.close()
    sat = ephem.readtle(tle[0],tle[1],tle[2])
    sat.compute(date)
    return sat

------
The problem is that there are large jumps in the CXBN.elevation, and this causes jumps in the components of r_eci. I am starting to look into the SGP4 code you are using, but this I don't yet understand why this is happening. Is this a problem with the code's treatment of the orbit around the point of perigee?

Lance

[Lance Simms (lance0381)]

As an additional note of information, it seems that the discontinuity is actually in the RA and DEC of the satellite object. Here you can see a plot of the difference in cosine of the declination between the current and precious time step

[Lance Simms (lance0381)]

I don't mean to be impatient, but is anyone looking at this website? I can't see any activity whatsoever in the last 10 days...

[Brandon Rhodes (brandon-rhodes)]

Lance — I apologize that you have been waiting! Yes — I, the author of PyEphem, am indeed watching this website. The email alerting me about your bug is, in fact, marked "important" in my incoming email folder. But not until I have given my talks this weekend at PyOhio will I be able to dig into the math of the satellite function, and try to figure out what is causing your discontinuity. I will let you know what I find — and thanks for your interest in PyEphem!

[Lance Simms (lance0381)]

Brandon, I apologize. I did not realize it was just you maintaining this package. I have been trying to figure out where in the code the error is coming from, but I still don't have a firm conclusion. I appreciate what you're doing with PyEphem. It's a really wonderful tool.

[Brandon Rhodes (brandon-rhodes)]

Lance —

I am now happily ensconced at PyOhio and have spare cycles to check out your issue! I have not yet been able to reproduce your problem after a few tries, however. First, your code does not run when pasted into a Python file — for one thing, the statement "EpDate = float(EpDate)" is trying to access the "EpDate" before assigning it a value for the first time. Another problem is that your code has many symbols that are never defined — but two friendly fluid dynamics engineers at the table suggested that they looked like NumPy names ("double", "arange"), so maybe I am missing some import statements from the top of your program? Finally, the numbers I got out looked nothing like yours; could it be that the transitions you are making to and from NumPy floats are losing some precision somewhere?

Anyway, to move forward I think I'll need your full .py file; making random changes to what you pasted above in an effort to produce working code is not producing anything like your graph. I'll keep watching this issue this evening and tomorrow, and start up the investigation again if I see you attach a file.

[Lance Simms (lance0381)]

Brandon,

I am sorry for posting an incomplete script. Not only did I forget to include all my import statements; I forgot to do things as critical as initializing EpDate! In the code that I was running, EpDate was indeed initialized, so that's not the problem.

I have created a bare-bones script that includes a minimal amount of packages to show the error. Assuming you have matplotlib, you should be able to plot the RA and DEC difference arrays and see the large jumps like the ones shown in the plots I have attached. If not, you will have to comment these statements out.

I've included the script as an attachment. You can see the only thing I've done is convert the ra and dec using the float() statement. But this is something that is done in the PyEphem examples, so I'm guessing this is not the problem...unless it's the way the Mac processor treats floats, but I can't imagine that would be the case.

I looked through the C source code, and I see how there could be a problem around 200km where the orbit is evaluated using different terms, but my orbit is nowhere near that altitude. So I'm still a bit puzzled.

Thanks for looking into this,

Lance

Note that I also run this out of the pyraf environment, so pyraf works interactively with matplotlib. I'm not sure if the plots will behave nicely from the command line, but I imagine you have your own way of plotting arrays.

==============
#!/usr/bin/python2.5

import matplotlib.pylab as mplot
import os
import ephem
import math
from datetime import *
from array import *

# Main routine if we are called as a script
if __name__ == "__main__":

  #DEFAULT VALUES FOR KEYWORDS
  DateTime = '2012/01/01 20:20:00'

  #Format the time appropriately
  Year = int(DateTime[0:4])
  Month = int(DateTime[5:7])
  Day = int(DateTime[8:10])
  Hours = int(DateTime[11:13])
  Mins = int(DateTime[14:16])
  Secs = int(DateTime[17:19])
  DateTime = datetime(Year, Month, Day, Hours, Mins, Secs)

  #Convert this into a PyEphem Date (# days since noon on 1899 December 31)
  EpDate_0 = ephem.date(DateTime)
  EpDate = EpDate_0
  EpDate = EpDate + ephem.second*0.01*98000

  #Form the satellite object
  tle1 = 'CXBN'
  tle2 = '1 99999U 11158.66666667 .00001767 00000-0 83150-4 0 00005'
  tle3 = '2 99999 060.0472 001.9175 0286684 357.7762 001.3936 14.89594857000013'
  CXBN = ephem.readtle(tle1,tle2,tle3)

  #Set the time increment for the propogation
  dt = 0.01 #Time difference in seconds
  t_final = 1*100. #Final time in seconds
  NumTimes = t_final/dt
  Times = range(NumTimes)

  #Set the initial time to t=0
  Time=0

  #Get the initial RA and DEC
  CXBN.compute(EpDate)
  RA_rad_prev = float(CXBN.ra)
  DEC_rad_prev = float(CXBN.dec)

  #Set up arrays
  RA_arr = array('f')
  DEC_arr = array('f')

  #############################################################################################################
  #PROPAGATE THE SATELLITE IN ITS ORBIT
  #############################################################################################################
  for i in Times:

    #Inrement the time
    EpDate = EpDate + ephem.second*dt
    Time ...

[Brandon Rhodes (brandon-rhodes)]

The script raised the exception "TypeError: range() integer end argument expected, got float." but then I threw an int() call into the expression "Times = range(int(NumTimes))" and the script ran! (I commented out the "print" statement since it was otherwise taking forever to run and actually draw the graphs.) So I can now indeed see the continuity problem on my own screen, and will drill in and see if I can find its cause. Thanks for packing the script in an executable form for me!

[Brandon Rhodes (brandon-rhodes)]

Lance, I think I have found the problem — the routine in "sgp4.c", beneath the comment "SOLVE KEPLERS EQUATION", has a loop that narrows down its estimate of the satellite's position. It is willing to run 10 iterations (!) to get a close-enough answer, but is willing to stop early if the answer seems to have already gotten close enough. And in this case "close enough" is defined by the constant E6A:

 if(fabs(EPW - TEMP2) <= E6A)
             break;

Which, higher up in the file, is defined as:

    #define E6A (1.E-6)

And this is the source of your problem: the region which you are exploring is one in which the loop usually does 2 iterations, but on occasion is happy enough with its first answer, and ends after only 1 iteration. And this causes your discontinuities: you are seeing the difference between where 1 iteration puts the answer and the more refined value that 2 iterations produces.

Try editing your own copy of PyEphem — are you able to download it and install it from source code? — and make this E6A symbol be a smaller number, by editing the line where it is defined and putting something like this:

#define E6A (1.E-12)

I think this will make the discontinuity small enough that you will no longer be bothered by it. If so, let me know, and I'll make that change to my copy to. But if not, then please experiment with even smaller values for E6A until you are happy with the result, and let me know what you wind up with. Thanks!

[Lance Simms (lance0381)]

Wow, great de-bugging-tective work! This makes sense as the root of the
problem. I installed PyEphem with the easy install tool, but I should be
able to compile it on my own. I will let you know once I have done so.
Thanks again.

On Mon, Aug 1, 2011 at 2:04 PM, Brandon Craig Rhodes <
<email address hidden>> wrote:

> Lance, I think I have found the problem — the routine in "sgp4.c",
> beneath the comment "SOLVE KEPLERS EQUATION", has a loop that narrows
> down its estimate of the satellite's position. It is willing to run 10
> iterations (!) to get a close-enough answer, but is willing to stop
> early if the answer seems to have already gotten close enough. And in
> this case "close enough" is defined by the constant E6A:
>
> if(fabs(EPW - TEMP2) <= E6A)
> break;
>
> Which, higher up in the file, is defined as:
>
> #define E6A (1.E-6)
>
> And this is the source of your problem: the region which you are
> exploring is one in which the loop usually does 2 iterations, but on
> occasion is happy enough with its first answer, and ends after only 1
> iteration. And this causes your discontinuities: you are seeing the
> difference between where 1 iteration puts the answer and the more
> refined value that 2 iterations produces.
>
> Try editing your own copy of PyEphem — are you able to download it and
> install it from source code? — and make this E6A symbol be a smaller
> number, by editing the line where it is defined and putting something
> like this:
>
> #define E6A (1.E-12)
>
> I think this will make the discontinuity small enough that you will no
> longer be bothered by it. If so, let me know, and I'll make that change
> to my copy to. But if not, then please experiment with even smaller
> values for E6A until you are happy with the result, and let me know what
> you wind up with. Thanks!
>
> --
> You received this bug notification because you are subscribed to the bug
> report.
> https://bugs.launchpad.net/bugs/812906
>
> Title:
> Discontinuity in Satellite Orbit
>
> Status in PyEphem astronomy library for Python:
> New
>
> Bug description:
> Hello,
>
> I am trying to simulate the orbit of a CubeSat satellite with the
> following TLE File (stored as TLE_File):
>
> CXBN
> 1 99999U 11158.66666667 .00001767 00000-0 83150-4 0 00005
> 2 99999 060.0472 001.9175 0286684 357.7762 001.3936 14.89594857000013
>
> I use the following code to start and propagate the orbit:
>
> DateTime = '2012/01/01 20:20:00'
> #Format the time appropriately
> Year = int(DateTime[0:4])
> Month = int(DateTime[5:7])
> Day = int(DateTime[8:10])
> Hours = int(DateTime[11:13])
> Mins = int(DateTime[14:16])
> Secs = int(DateTime[17:19])
> DateTime = datetime(Year, Month, Day, Hours, Mins, Secs)
>
> #Convert this into a PyEphem Date (# days since noon on 1899 December
> 31)
> EpDate_0 = ephem.date(DateTime)
> EpDate = float(EpDate)
>
> #Call getsat to return an ephemeris satellite object based on the TLE
> CXBN = getsat(TLE_File, EpDate)
>
> #Set the time increment for the propogation
> dt = double(0.01) #Time difference in seconds
> t_final = 1*1. ...

[Lance Simms (lance0381)]

Brandon,

The suggestion you made helps out quite a bit. It's removed the major kinks that I showed in the previous plot. There is still one discontinuity near apogee, even with E6A=1.E-14. I have to make sure this is due to a jump in RA or DEC, though.

I would say that setting E6A to 1.E-12 is a good start for now. I will let you know if lowering the value further helps get rid of this last kink.

Lance

[Dave Gordon (python-bugz)]

Apropos #11, that's probably down to a change in Python 2.7. Up until 2.6, if you passed a float to a function that was defined as taking an int, the float would automatically be converted, though possibly with a DeprecationWarnng. In 2.7 (and Python 3000) it's now a TypeError. So, you have to put in an explicit conversion to int wherever this occurs.

There are quite a lot of packages affected by this!

Dave

[Brandon Rhodes (brandon-rhodes)]

The increased e-12 precision will be part of PyEphem 3.7.5.1, which I am releasing this evening. Enjoy!