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
#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.)
#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')
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.
============== bin/python2. 5
#!/usr/
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 11:13]) 14:16]) 17:19])
Year = int(DateTime[0:4])
Month = int(DateTime[5:7])
Day = int(DateTime[8:10])
Hours = int(DateTime[
Mins = int(DateTime[
Secs = int(DateTime[
DateTime = datetime(Year, Month, Day, Hours, Mins, Secs)
#Convert this into a PyEphem Date (# days since noon on 1899 December 31) DateTime) 0.01*98000
EpDate_0 = ephem.date(
EpDate = EpDate_0
EpDate = EpDate + ephem.second*
#Form the satellite object tle1,tle2, tle3)
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(
#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 compute( EpDate)
CXBN.
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 compute( EpDate)
CXBN.
#First get the RA and DEC from satellite object elevation+ 6738100. )
RA_rad = float(CXBN.ra)
DEC_rad = float(CXBN.dec)
Radius = float(CXBN.
print RA_rad-RA_rad_prev DEC_rad_ prev
#print DEC_rad-
#Now convert to ECI math.cos( DEC_rad) *math.cos( RA_rad) math.cos( DEC_rad) *math.sin( RA_rad) math.sin( DEC_rad)
x_eci = Radius*
y_eci = Radius*
z_eci = Radius*
#Append the differences arr.append( RA_rad- RA_rad_ prev) arr.append( DEC_rad- DEC_rad_ prev)
RA_
DEC_
#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 'Difference in RA over Time Step') './Discontinuit yInRA_wPyEphem. png')
mplot.figure(0)
mplot.clf()
mplot.plot(RA_arr)
mplot.xlabel('Time Step')
mplot.ylabel(
mplot.savefig(
mplot.figure(1) 'Difference in DEC over Time Step') './Discontinuit yInDEC_ wPyEphem. png')
mplot.clf()
mplot.plot(DEC_arr)
mplot.xlabel('Time Step')
mplot.ylabel(
mplot.savefig(