skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.43k stars 213 forks source link

Satellites: rounding-off Geoc/Topocentric matrices to 3 decimal places #252

Closed sunrise125 closed 5 years ago

sunrise125 commented 5 years ago

The code below shows the proper way to cut some variables up to 3 decimals, namely subpoint elevation, distance, but the same command ({0:.3f}) fails when applied to vectors (rho), (R), (r). May I have the appropriate drift to perform them? Thanks. Edit: sorry, unable to put tags for the code.

`from skyfield.api import EarthSatellite, Topos, load ts = load.timescale() timestring= '2019, 4, 27, 4, 18, 0' t = ts.utc(2019, 4, 27, 4, 18, 0) # line1 = '1 25544U 98067A 19116.21512024 .00001545 00000-0 32153-4 0 9998' line2 = '2 25544 51.6416 259.5314 0000745 223.5003 196.2377 15.52590521167188' # loc = Topos('37.0328 N', '15.0650 E') satellite = EarthSatellite(line1, line2, name='ISS (ZARYA)')

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')

EOF: ----- fullsat.py -----------`

sunrise125 commented 5 years ago

`` from skyfield.api import EarthSatellite, Topos, load

ts = load.timescale() timestring= '2019, 4, 27, 4, 18, 0' t = ts.utc(2019, 4, 27, 4, 18, 0) # line1 = '1 25544U 98067A 19116.21512024 .00001545 00000-0 32153-4 0 9998' line2 = '2 25544 51.6416 259.5314 0000745 223.5003 196.2377 15.52590521167188' # loc = Topos('37.0328 N', '15.0650 E') satellite = EarthSatellite(line1, line2, name='ISS (ZARYA)')

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')

EOF: ----- fullsat.py -----------

``

brandon-rhodes commented 5 years ago

Can you find a guide to writing GitHub markdown? I think if you'll edit the comments you'll find that three back-ticks in a row let you create code sections that I'll be able to cut and paste to try your examples.

Note that this might not be a question that Skyfield can answer, if you're asking in general about how to get NumPy arrays to display to a certain number of decimal places?

sunrise125 commented 5 years ago

Good! Your code provides proper cutoffs as fòllows Subpoint Elevation= 412.682 km Distance= 454.442 km but not on computing next irrealistic precision [1e-5 meters] Obs. posit.(R): [ 2095.83987655 -4650.28825894 3816.35162127] km being devoted the matter to ... Numpy? Very strange

ghost commented 5 years ago

I think this is a question of preserving significant digits, and the question of accuracy versus precision. It's a tricky issue in general (even Mathematica struggles with it), and I don't think it's something that can or should be solved in SkyField. Perhaps you could create or find another project that deals with numerical accuracy, precision, and the preservation of significant digits.

brandon-rhodes commented 5 years ago

I am not aware of any astronomy libraries that keep track of significant digits, so I have not thought through all of the complex issues that would arise when every computation had to process both the quantities and also their error bars. If you have seen a library before that keeps track of significant digits, I would be interested in your sharing a link so I could take a look at their interface. I have not seen it done before in a Python library.

sunrise125 commented 5 years ago

Just solved as follows: Code from numpy import around sho2= around(topocentric.position.km, decimals=3) print(' Topocentric (rho): ',sho2,'km') Output Topocentric (rho): [ -33.244 -290.066 348.244] km

brandon-rhodes commented 5 years ago

Good, I'm glad you found a solution in NumPy's documentation!