Comment 13 for bug 812906

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. #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
>
> To manage notifications about this bug go to:
> https://bugs.launchpad.net/pyephem/+bug/812906/+subscriptions
>