skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.38k stars 208 forks source link

calculate earth subpoint (long/lat), given earth velocity vector at time (barycentric, from skyfield) #890

Closed beluej123 closed 6 months ago

beluej123 commented 10 months ago

Please help me figure out how to calculate the subpoint (in long/lat) of the earth velocity vector (or other skyfield supported object). Using the following I can get position and velocity vectors, but I do not know how to get the subpoint of the velocity vector:

from skyfield.api import load t = load.timescale( ).now( ) planets = load( "de421.bsp" ) planet_position=planets[ "earth" ].at( t ).position.au planet_velocity=planets[ "earth" ].at( t ).velocity.km_per_s

Thank you.

brandon-rhodes commented 10 months ago

That's an interesting question, because it requires treating a velocity as though it's really a position. Maybe:

from skyfield.api import load, wgs84
from skyfield.positionlib import ICRF

t = load.timescale().now()
planets = load('de421.bsp')
earth = planets['Earth']
position = earth.at(t)
planet_velocity = planets["earth"].at(t).velocity.km_per_s

# Make a really big vector in the direction of the velocity.
v = position.velocity.au_per_d * 1e10

# Pretend that it's really a position in space, and ask for its subpoint.
p = ICRF(v, t=t, center=399)
g = wgs84.geographic_position_of(p)

print(g)
print(g.latitude)
print(g.longitude)

Does that give results you might expect? For me at this moment it gives:

WGS84 latitude +18.5121 N longitude -75.1915 E elevation 2.5e+19 m
18deg 30' 43.5"
-75deg 11' 29.5"

See if that gives you reasonable numbers or not! If not, then see if you can find a specific time t for which you know what the lat/lon should be, and we'll see if we can get Skyfield to return those exact numbers.

beluej123 commented 10 months ago

My motivation for this calculation, i.e. earth subpoint of earth velocity vector (at a given time); The best viewing of meteor shower is between 1:am->dawn (Perseid shower); ignoring light sources such as moon, local light pollution? As reported in the general astro-press (sky&telescope and others).

I have not seen a calculation to substantiate that claim (yes, I realize there is a "statistical" cloud of debris the earth passes thru; and viewing in the dark is a must). Never-the-less, I think this is a worth while and practical calculation "exorcise" since the Perseid shower peak just finished for this year.

Two logical questions seem to follow: (1) given a time, what earth location (long/lat) faces the earth velocity vector (flight path) and, (2) when is my location closest to the velocity vector for a given day/night?

*** Background for calculation support ** (1.1) Treating the velocity vector as a position vector seems reasonable since this calculation only interested direction (flight path). The velocity/position magnitude does not matter. Since the position of the velocity vector defines the flight path, this is the leading edge of the earth facing the meteor debris.

(1.2) use skyfield to translate ICRF coordinates to earth coordinates (long/lat) CONSIDERING earth rotation also causes velocity subpoint to change long/lat.

(2.1) I am not sure how to setup a skyfield procedure to find when my location is closest to the velocity subpoint vector. For calculating closest subpoint to my location, I believe finding minimum (subpointLong - myLong) answers the question of when my location is closest to the velocity vector.


BTW: I do not know how to add labels or include code like you do above... Can you tell me where I can go to figure out how?


I made a couple minor changes to your code, above:

// calculate earth subpoint (long/lat) at earth velocity vector, given time from skyfield.api import load, wgs84 from skyfield.positionlib import ICRF

// https://rhodesmill.org/skyfield/time.html //t = load.timescale().now() ts = load.timescale() t = ts.tt(2023, 8, 13, 5, 50)

planets = load('de421.bsp') earth = planets['Earth'] position = earth.at(t) planet_velocity = planets["earth"].at(t).velocity.km_per_s

// Make a really big vector in the direction of the velocity. // v = position.velocity.au_per_d * 1e10 # not sure why the multiplier, 1e10 v = position.velocity.au_per_d # treat velocity vector as a position vector

// Pretend velocity vector is a position in space, and ask for its earth subpoint. p = ICRF(v, t=t, center=399) g = wgs84.geographic_position_of(p)

print("given time:", t.utc_strftime()) print("calculated subpoint:", g) print("v-vector subpoint latitude:", g.latitude) print("v-vector subpoint longitude:", g.longitude)

brandon-rhodes commented 10 months ago

BTW: I do not know how to add labels or include code like you do above... Can you tell me where I can go to figure out how?

Here's how GitHub lets you mark up text as code:

https://docs.github.com/en/get-started/writing-on-github/getting-started-with-writing-and-formatting-on-github/basic-writing-and-formatting-syntax#quoting-code

You should be able to press Edit on your post and add the markup!

(And I'll answer the more technical astronomy comments when I get some time, hopefully later this week.)

brandon-rhodes commented 6 months ago

@beluej123 — Have you had the chance to follow the above link, which should let you get your source code in your previous comment into shape? If you don't have the time right now for further work on this issue, then I'll be happy to temporarily close it until you have more free time; I'll wait a few days, though, in case you want to start this thread back up. Thanks!

brandon-rhodes commented 6 months ago

I'll temporarily close this issue, but please feel free to comment further, even with the issue closed, if you have more questions.