Closed bradfordrobbins closed 1 month ago
brandon-rhodes
commented
1 month ago Could you provide your TLE data, along with the coordinates that you expected instead of these? Your code looks fine, so I don't see any way to proceed without a bit more information. See if you can provide a small Python script that runs, including its import statements and the data it needs. Thanks!
bradfordrobbins
commented
1 month ago My apologies, I had a bug that caused the discrepancy. With the code below, it works. You can close this.
However, I would value your advice. I am trying to understand how celetrak socrates detects satellites that may collide without adjustment. The TLE info in the code below is when they predicted a collision (a conjunction) 5/20 GMT with a distance < 1 km, however, I get a much bigger separate. Any suggestions where I could look to learn more about why I'm getting a different answer? Thanks.
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
To demonstrate puzzle Sat distance vs. Celetrak
https://celestrak.org/SOCRATES/table-socrates.php?NAME=,&ORDER=MAXPROB&MAX=10
it predicted a conjuction at tclosest for the two satellites below.
Note the orbits were adjusted and the satellites will not meet any longer
However at the time listed I find the satellites are not that close so I
trying to determine what I am doing wrong.
@author: brad
"""
from skyfield.api import load, EarthSatellite, wgs84
from datetime import datetime
from pytz import timezone
import numpy as np
# get orbital period of satellite in seconds
def get_period_sec(sat):
mean_motion = sat.model.no_kozai # in radians per minute
p = 2*3.141592*60/mean_motion # orbital period in second
return p
# print out information about a satellite
def print_sat_info(sat, t):
p = get_period_sec(sat)/60 # orbital period in minutes
geocentric = sat.at(t)
pos = geocentric.position.km
# basically length = (pos[0]**2+pos[1]**2+pos[2]**2)**.5
length = np.linalg.norm(pos)
lat, lon = wgs84.latlon_of(geocentric)
h = wgs84.height_of(geocentric).km
print(f'{sat.name:15}', f" Period (min): {p:,.0f}")
print(f" Epoch: {sat.epoch.astimezone(gmt).strftime('%m-%d'):}",
f" Time: {t.astimezone(gmt).strftime('%m-%d %H:%M:%S'):}")
print(f" Position: [{pos[0]:.0f}, {pos[1]:.0f}, {pos[2]:.0f}]",
f" Vector length: {length:,.0f}")
print(f" Latitude: {lat}",f" Longitude: {lon}",f" Height km: {h:,.0f}")
return
# print details about the closest pass
def print_details_closest_pass(sat, t):
geocentric = sat.at(t)
lat, lon = wgs84.latlon_of(geocentric)
height = wgs84.height_of(geocentric)
geo_pos = wgs84.geographic_position_of(geocentric)
g_lon = geo_pos.longitude.degrees
g_lat = geo_pos.latitude.degrees
print(f"{sat.name:15} at time: {t.astimezone(gmt).strftime('%m-%d %h:%m')}")
print(f" Position km: {geocentric.position.km}")
print(f"Velocity km/s: {geocentric.velocity.km_per_s}")
print(f" Time: {geocentric.t.astimezone(gmt).strftime('%m-%d %H:%M:%S')}",
f" Center: {geocentric.center}",
f" Target: {geocentric.target}")
print('Space track',f" Lat: {lat.degrees:.4f}",
f" Lon: {lon.degrees:.4f}",
f" Height km: {height.km:,.0f}")
print('Ground track',f" Lat: {g_lat:,.4f}",
f" Lon: {g_lon:,.4f}")
return
#--------------------- Main line --------------------------------------
# Create a timescale
gmt = timezone('GMT')
us_pacific = timezone('US/Pacific')
now = datetime.now(gmt)
ts = load.timescale()
tclosest = ts.utc(2024,5,20,3,54,20.206)
t = tclosest
# hard code the values to avoid calling celetrak too many times!
# https://celestrak.org/SOCRATES/data.php?CATNR=55042,58381
# these are predicted to collide unless orbits are adjusted
# the time of collision or conjuction is tclosest
print("---------------------------------------------------------------------")
line1 = '1 58381U 23178H 24137.89624224 .00009589 00000+0 35847-3 0 9990'
line2 = '2 58381 53.1277 197.3288 0001546 98.6590 261.4585 15.28012846 28461'
satellite1 = EarthSatellite(line1, line2, 'STARLINK-30923', ts)
print_sat_info(satellite1, t)
print("---------------------------------------------------------------------")
line1 = '1 55042U 23001AK 24137.43384556 .00023799 00000+0 88384-3 0 9990'
line2 = '2 55042 97.4269 199.3699 0008343 272.7584 87.2700 15.27528193 75935'
satellite2 = EarthSatellite(line1, line2, 'FLOCK 4Y-33', ts)
print_sat_info(satellite2, t)
#tstart, min = find_closest(satellite1, satellite2, t0)
#print(tstart.astimezone(gmt).strftime('%m-%d %H:%M'), f' dist: {min:,.0f}')
print("---------------------------------------------------------------------")
d = (satellite1.at(t) - satellite2.at(t)).distance().km
print(f'At predicted closest time {tclosest.astimezone(gmt):%m-%d %h:%m} distance between: {d:,.0f}')
# print out details of closest pass
print("---------------------------------------------------------------------")
print_details_closest_pass(satellite1,tclosest)
print("---------------------------------------------------------------------")
print_details_closest_pass(satellite2,tclosest)
print("---------------------------------------------------------------------")Thanks for the update, I'll go ahead and close! My first guess about your lack-of-collision problem would be that you are using different TLEs than Celestrak was using when it made the prediction; are you sure that you selected the TLEs that would have been current as of the moment that they made that determination?
bradfordrobbins
commented
1 month ago The Celetrak website lists the TLE data and time of conjunction, but I’ll double check that. Thanks.
For example 
On May 18, 2024, at 8:00 PM, Brandon Rhodes @.***> wrote:
Thanks for the update, I'll go ahead and close! My first guess about your lack-of-collision problem would be that you are using different TLEs than Celestrak was using when it made the prediction; are you sure that you selected the TLEs that would have been current as of the moment that they made that determination?
— Reply to this email directly, view it on GitHub https://github.com/skyfielders/python-skyfield/issues/967#issuecomment-2119079656, or unsubscribe https://github.com/notifications/unsubscribe-auth/AML45TY6AAX73FJLCPBB2JLZDAIVTAVCNFSM6AAAAABH5XKYB6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCMJZGA3TSNRVGY. You are receiving this because you authored the thread.
I ran the calculation for several satellites and found the Starlink satellites ground reference location appears to be wrong, 10-30 degrees. Both the geographic_position_of and subpoint methods gave the same result. The code snippet I ran each satellite through is shown below as well as the output.
p.s. you have written a wonderful set of code!
--------------------------------------------------------------------- STARLINK-30923 at time: 05-20 May:05 Position km: [-6453.87903664 1270.69855276 -1958.48059866] Velocity km/s: [-2.58059505 -4.3538169 5.69866407] Time: 05-20 03:54:20 Center: 399 Target: -158381 Space track Lat: -26.3513 Lon: -90.1225 Height km: 491 Ground track Lat: -16.8114 Lon: -127.7222 --------------------------------------------------------------------- FLOCK 4Y-33 at time: 05-20 May:05 Position km: [-5521.47038643 -2743.78729126 -3018.84064098] Velocity km/s: [-3.45530927 -0.50134564 6.77193602] Time: 05-20 03:54:20 Center: 399 Target: -155042 Space track Lat: -26.3513 Lon: -90.1225 Height km: 491 Ground track Lat: -26.3513 Lon: -90.1225 ---------------------------------------------------------------------