Comment 7 for bug 812906

Revision history for this message
Lance Simms (lance0381) wrote :

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 += dt

    #Recompute the solution with the SGP4
    CXBN.compute(EpDate)

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

    print RA_rad-RA_rad_prev
    #print DEC_rad-DEC_rad_prev

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

    #Append the differences
    RA_arr.append(RA_rad-RA_rad_prev)
    DEC_arr.append(DEC_rad-DEC_rad_prev)

    #Set the ra to the previous ra and dec to the previous dec
    RA_rad_prev = RA_rad
    DEC_rad_prev = DEC_rad

#Plot the RA and DEC differences to show the discontinuities
mplot.figure(0)
mplot.clf()
mplot.plot(RA_arr)
mplot.xlabel('Time Step')
mplot.ylabel('Difference in RA over Time Step')
mplot.savefig('./DiscontinuityInRA_wPyEphem.png')

mplot.figure(1)
mplot.clf()
mplot.plot(DEC_arr)
mplot.xlabel('Time Step')
mplot.ylabel('Difference in DEC over Time Step')
mplot.savefig('./DiscontinuityInDEC_wPyEphem.png')