Get intermediate points of otbit propagation

[tmamedzadeh]

I'm propagating the orbit for 2 days time interval. Is there a way to get the intermediate points of propagation with 1 second periodicity?

[brandon-rhodes]

Yes. The most straightforward approach would be a for loop that called propagate() repeatedly, increasing the value for second from 0 to 1, 2, 3, … up to the number of seconds in two days (60 60 24 * 2 = 172800). There's no need to stop and turn the seconds into minutes and hours and dates, as the fields of propagate are simply added together without regard for whether the minutes and seconds are in range or not.

[tmamedzadeh]

That's what I'm doing now:

  while start_time <= end_time:
      date= datetime.datetime(2000, 1, 1, 12, 0) + datetime.timedelta(days=start_time - 2451545)
      p,v= tle.propagate(date.year,date.month, date.day, date.hour, date.minute, date.second) 
      start_time += 1

You mean, I could do:

  date= datetime.datetime(2000, 1, 1, 12, 0) + datetime.timedelta(days=start_time - 2451545)
  seconds=date.second
  while start_time <= end_time:
      p,v= tle.propagate(date.year,date.month, date.day, date.hour, date.minute, seconds) 
      seconds+=1
      start_time += 1

That's OK. However, it doesn't make calculations quicker. It would be better to propagate once in 2 days interval, storing the intermediate points of calculations. Therefore, calling the function only once, not 60* 60* 24* 2 times.

[tmamedzadeh]

BTW, my code works longer due to conversion to lat/long:

  date= datetime.datetime(2000, 1, 1, 12, 0) + datetime.timedelta(days=start_time - 2451545)
  seconds=date.second
  while start_time <= end_time:
      p,v= tle.propagate(date.year,date.month, date.day, date.hour, date.minute, seconds) 

      # Conversion to ITRS    
      p,v= sgp4lib.TEME_to_ITRF(start_time,np.asarray(p),np.asarray(v)*86400)
      v=v/86400
      p=p*1000

      # Conversion to spherical    
      ecef = pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84')
      lla = pyproj.Proj(proj='latlong', ellps='WGS84', datum='WGS84')
      lon, lat, alt = pyproj.transform(ecef, lla, p[0], p[1], p[2], radians=False)

      json_arr.append({'id':id,'name':name, 'time':str(datetime.datetime(2000, 1, 1, 12, 0) + datetime.timedelta(days=start_time - 2451545)),'lat':lat,'lon':lon,'alt':alt})
      start_time += 1/86400 #It shows days count
      seconds += 1

Would you suggest a better way?

[brandon-rhodes]

I'll try to find time this weekend to run your sample code!

[tmamedzadeh]

I'm willing to use Skyfield for this purpose, if I could find an example of calculation of lat/lon/alt from the propagated TLE data.

Now, I use the following to get the result in GCRS:

geocentric = satellite.at(t)
print(geocentric.position.km)

Also, I use alt, az, distance = topocentric.altaz() to get az/el.

However, in the documentation I couldn't find any information about getting the latitude and longitude of ITRS frame. I would appreciate if you could share some examples.

Replaced to https://github.com/skyfielders/python-skyfield/issues/168