brandon-rhodes / python-sgp4

Python version of the SGP4 satellite position library
MIT License
371 stars 87 forks source link

Difference between SGP4 1.4 and Pyhem #11

Closed powersoft closed 5 years ago

powersoft commented 8 years ago

Hello,

I started to play with tracking the ISS space station. Have approached it with Pyhem and SGP4 1.4. When I do the calculation I get a different answer for the longitude. I have checked the results with a tracker found on the internet, and the value of the calculation with Pyephem are the same!

I use the ECEF routine to convert the (x,y,z) to lat,lon and alt.

Thanks for any support,or help.

Jan Kromhout Hellevoetsluis-NL

import urllib.request from sgp4.earth_gravity import wgs84 from sgp4.io import twoline2rv from datetime import datetime from math import sqrt, atan2, sin, cos, pi import ephem

def ecef(x, y, z): #Convert (x,y,z) to (lat,lon,alt) a = 6378137 e = 8.1819190842622e-2 b = sqrt(a * 2 * (1 - e * 2)) ep = sqrt((a * 2 - b * 2) / (b * 2)) p = sqrt(x * 2 + y * 2) th = atan2(a * z, b * p) lon = atan2(y,x) lat = atan2((z + ep * 2 * b * sin(th) * 3), (p - e * 2 * a * cos(th) * 3)) n = a / sqrt(1 - e * 2 * sin(lat) * 2) alt = p / cos(lat) - n lon = lon % (2 \ pi) return lat_180/pi, lon_180/pi, alt

def update_tle(): data=[]

Schrijf data tijdelijk in file

urllib.request.urlretrieve('http://www.celestrak.com/NORAD/elements/stations.txt','test.dat') f=open('test.dat','r') data=f.read() data=data.split('\n') return data[0],data[1],data[2] #Gegevens International Space Station

print('Read stations data') name,line1,line2=update_tle()

print(name) print(line1) print(line2)

satellite = twoline2rv(line1, line2, wgs84) now=datetime.utcnow() print ('UTC time ',now) print('Calculate with sgp4 1.4') position, velocity = satellite.propagate(now.year,now.month,now.day,now.hour,now.minute,now.second) lat,lon,alt=ecef(position[0]_1000,position[1]_1000,position[2]_1000) speed=sqrt(velocity[0]2+velocity[1]2+velocity[2]__2)_3600 print('epochyr ',satellite.epochyr) print('epochdays ',satellite.epochdays) print('jdsatepoch ',satellite.jdsatepoch) print('epoch ',satellite.epoch) print('position (km) ',position) print('velocity (km/s) ',velocity) print('orbital speed (km/h) ',speed) print() print('Position of ISS at time ',now) print('Latitude (deg) ',lat) print('Longitude(deg) ',lon) print('Height ',alt)

home = ephem.Observer() home.lon = '52' # +E home.lat = '05' # +N

iss= ephem.readtle(name, line1, line2) iss.compute(now) lat=ephem.degrees(iss.sublat) lon=ephem.degrees(iss.sublong) print('Calculate with Ephem') print('Latitude (deg) ',float(lat)_180/pi) print('Longitude (deg) ',float(lon)_180/pi)

dcajacob commented 8 years ago

SGP4 returns position and velocity in an inertial frame, not the rotating ECEF frame, so calculating lat/lon directly from position (without properly doing a coordinate transformation to ECEF first) is meaningless.

powersoft commented 8 years ago

Thanks for the reply. Are these transformations available in SGP4?

dcajacob commented 8 years ago

They probably are at a low level, but they are not exposed.

powersoft commented 8 years ago

Maybe you can help me with this problem.

I want whit a given RA and DEC calculate the position of the moon perpendicular to the earth (lat- and longitude) so I can plot this on my world map. Have done some search on the web to look for a solution but without result.

brandon-rhodes commented 5 years ago

There wasn’t an easy answer back when you asked, but Skyfield now has a routine that will give you the latitude and longitude beneath any object in the sky:

https://rhodesmill.org/skyfield/api-position.html#skyfield.positionlib.Geocentric.subpoint

(I realize that it’s probably far too late for this to help the project you were working on it 2016, but I mention it here in case it helps you or other readers in future projects!)