Highly precise python code to compute sat positions

From: Giuseppe via Seesat-l <seesat-l_at_satobs.org>
Date: Tue, 21 May 2019 09:51:48 +0000
Dear all,
please, let me show a very concise software, which uses the powerful
python libraries available here https://rhodesmill.org/skyfield/earth-satellites.html
and easy uploadable through (pip install skyfield) command. Their accuracy
is there described, as follows:
-------------
Skyfield computes positions for the stars, planets, and satellites in orbit
around the Earth. Its results should agree with the positions generated by
the United States Naval Observatory and their Astronomical Almanac to within
0.0005 arcseconds (which equals half a “mas” or milliarcsecond).
-------------

Results of the example below can be compared with the screenshot
ISS-visible.png catchable from my GitHub repo
https://github.com/sunrise125/Satellites-Orbits/find/master

In details, (elevation, azimuth, distance) values.
* Figure: (80.69 degs, 348.07 degs, 415.5 km)
* Software python (outputs below):
  Altitude=  81deg 06' 39.2"
  Azimuth =  342deg 02' 03.9"
  Distance=  415.500 km

In full harmony

Regards,
Giuseppe

------------------------------------------
Python code to fully determine sat coords
------------------------------------------
# --------- fullsat.py --------- Apr.27-May.07, 2019 --------------------
from skyfield.api import EarthSatellite, Topos, load
from numpy import around

ts = load.timescale()
timestring= '2019, 5, 21, 18, 22, 50'
t = ts.utc(2019, 5, 21, 18, 22, 50)   # datetime selection
# TLE twoline dbase
line1= '1 25544U 98067A   19141.19851322  .00000663  00000-0  18020-4 0  9997'
line2= '2 25544  51.6411 135.2638 0001720  14.4836 102.9344 15.52691834171066'
#
loc = Topos('37.0328 N', '15.0650 E')    # location coords
satellite = EarthSatellite(line1, line2)

# Geocentric
geometry = satellite.at(t)

# Geographic point beneath satellite
subpoint = geometry.subpoint()
latitude = subpoint.latitude
longitude = subpoint.longitude
elevation = subpoint.elevation

# Topocentric
difference = satellite - loc
geometry = difference.at(t)
topoc= loc.at(t)
#
topocentric = difference.at(t)
geocentric = satellite.at(t)
# ------ Start outputs -----------
print ('\n Ephemeris time:', timestring)
print (' JD time: ',t)
print ('',loc)
print ('\n Subpoint Longitude= ', longitude )
print (' Subpoint Latitude = ', latitude )
print (' Subpoint Elevation=  {0:.3f}'.format(elevation.km),'km')
# ------ Step 1: compute sat horizontal coords ------
alt, az, distance = topocentric.altaz()
if alt.degrees > 0:
    print('\n',satellite, '\n is above the horizon')
print ('\n Altitude= ', alt )
print (' Azimuth = ', az )
print (' Distance=  {0:.3f}'.format(distance.km), 'km')
#
# ------ Step 2: compute sat RA,Dec [equinox of date] ------
ra, dec, distance = topocentric.radec(epoch='date')
print ('\n Right Ascension RA= ', ra )
print (' Declination     de= ', dec )
#
# ------ Step 3: compute sat equatorial coordinates  --------
print ('\n Vectors to define sat position: r = R + rho')
print('    Obs. posit.(R): ',topoc.position.km,'km')
print(' Topocentric (rho): ',topocentric.position.km,'km')
print(' -----------------')
print('    Geocentric (r): ',geocentric.position.km,'km')
#
# ------ Step 4: sat equatorial coordinates roundoff 3 decimals  --------
sho1= around(topoc.position.km, decimals=3)
sho2= around(topocentric.position.km, decimals=3)
sho3= around(geocentric.position.km, decimals=3)
print ('\n Rounded Vectors to define sat position: r = R + rho')
print('    Obs. posit.(R): ',sho1,'km')
print(' Topocentric (rho): ',sho2,'km')
print(' -----------------')
print('    Geocentric (r): ',sho3,'km')
# EOF: ----- fullsat.py -----------

Running command
---------------
python fullsat.py


# ------- OUTPUT -------------------------------------
C:\Training>python fullsat.py
[#################################] 100% deltat.data
[#################################] 100% deltat.preds
[#################################] 100% Leap_Second.dat

 Ephemeris time: 2019, 5, 21, 18, 22, 50
 JD time:  <Time tt=2458625.2666572225>
 Topos 37deg 01' 58.1" N 15deg 03' 54.0" E

 Subpoint Longitude=  14deg 51' 16.0"
 Subpoint Latitude =  37deg 32' 58.3"
 Subpoint Elevation=  410.814 km

 EarthSatellite number=25544 epoch=2019-05-21T04:45:52Z
 is above the horizon

 Altitude=  81deg 06' 39.2"
 Azimuth =  342deg 02' 03.9"
 Distance=  415.500 km

 Right Ascension RA=  11h 03m 55.09s
 Declination     de=  +45deg 25' 43.4"

 Vectors to define sat position: r = R + rho
    Obs. posit.(R):  [-5007.45284283   917.56772044  3829.58886589] km
 Topocentric (rho):  [-282.0581147    71.84247777  296.51702639] km
 -----------------
    Geocentric (r):  [-5289.51095753   989.41019821  4126.10589229] km

 Rounded Vectors to define sat position: r = R + rho
    Obs. posit.(R):  [-5007.453   917.568  3829.589] km
 Topocentric (rho):  [-282.058   71.842  296.517] km
 -----------------
    Geocentric (r):  [-5289.511   989.41   4126.106] km
------------------End OUTPUT -------------------------------------------


_______________________________________________
Seesat-l mailing list
http://mailman.satobs.org/mailman/listinfo/seesat-l
Received on Tue May 21 2019 - 04:52:42 UTC

This archive was generated by hypermail 2.3.0 : Tue May 21 2019 - 09:52:42 UTC