Subsolar/sublunar points calculation

[ARK-ARadi]

Hello,

It probably is my ignorance, but it's been 3 days that I'm trying to calculate subsolar and sublunar points, at no avail.

I do not seem to be able to produce a geocentric object from the lib. Would you please help?

And please forgive me if it were a simple line of code, but the following (which I thought would be the most intuitive) did not work:

earth.at(ts.now()).observe(sun).subpoint() earth.at(ts.now()).observe(moon).subpoint()

Thanks a lot in anticipation, and kind regards.

[ghost]

I vaguely remember this being a huge issue at one point, but can't seem to find that much on it now. This not-very-helpful link refers to Brandon's other project, pyephem, which is now sadly deprecated:

https://astronomy.stackexchange.com/questions/7936/given-a-date-obtain-latitude-and-longitude-where-is-the-sun-zenith

[ARK-ARadi]

Thank you barrycarter for your response.

When I stumbled upon subpoint() I thought it was god-sent, but then I thoroughly failed to generate geocentric objects with skyfield. Indeed, observing a satellite is the only way I found to generate a geocentric object.

Still, the following "ugly" code works:

`from skyfield.api import load, pi

planets = load('de421.bsp') earth = planets['earth'] sun = planets['sun'] ts = load.timescale()

def lat(dec): lat = dec.radians / pi * 180 return lat

def lon(ra, t): lon = (ra.radians / pi 12 - t.gmst) 15 if lon > 180: lon -= 360 elif lon < -180: lon += 360 return lon

t = ts.tt(1964, 1, 1) ras, decs, dist = earth.at(t).observe(sun).radec() print(lat(decs), lon(ras, t)) `

[brandon-rhodes]

Good afternoon! Could you provide a working script that tries calling "subpoint()" as described in your question, along with the exception or error you get as a result? Those two items should make it easy for one of us to jump in and figure out how to get your code working or get Skyfield fixed. Thanks!

[ARK-ARadi]

Hello Brandon,

Calling subpoint() as follows fails, letting me know that an Astrometric object has no attribute subpoint().

from skyfield.api import load

ts = load.timescale()

planets = load('de421.bsp') earth = planets['earth'] sun = planets['sun']

t = ts.now()

a = earth.at(t).observe(sun)

print(a.subpoint())

Also, a.center and a.target are 'None'.

The code in my second post works fine, it just doesn't look good. Would you please help me fix it if you see room for improvement?

[brandon-rhodes]

One issue, I think, is that an astrometric position is a bit of a fiction — it's not a literal position where the object is really located, but is where the object appears to have been based on light-travel time. To get a subpoint, I think you'd have to use vector math. Maybe:

p = (sun - earth).at(t)
print(p.subpoint())

Try that and see if it works!

[ghost]

Not sure how helpful this is, but https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/spkezp_c.html#Particulars (the particulars section) explains how CSPICE deals with the issue of light travel time. I vaguely recall reading somewhere that light travel time has a surprisingly small effect on position, because most of the angular movement you see from Earth is due to the Earth's own rotation and the correction for light time is extremely small compared to this. However, I can't find a source.

[ARK-ARadi]

Thanks a lot Brandon. The following code not only works, it is also by far more elegant:

from skyfield.api import load

ts = load.timescale()

planets = load('de421.bsp') earth = planets['earth'] sun = planets['sun']

t = ts.now() a = (sun - earth).at(t)

print(t, a.subpoint(), sep = '\n')

I wish I could help tune a.subpoint().dstr a bit, as the output from the code above reads as follows: Topos 01deg 11' 11.8" N -148deg 26' 38.9" E But that would be preposterous, given my Python skills :-(

Now, if I understand right, 'a' above describes the actual zenith location of the sun, NOT taking into account the time light needs to travel, right?

[brandon-rhodes]

If you check the documentation, I think you'll find routines that will help you print the numbers in whatever format you want if you dislike the default. :)

a is the point on Earth that the sun is directly above. It takes several minutes for the light to reach Earth so the sun would appear overhead somewhere slightly different. I'll think about whether there should be a way to at least have Skyfield estimate the location from which an object appears overhead. (You could force Skyfield to do it today by directly creating a Geocentric position with the x,y,z position of the Astrometric one you calculated earlier; but maybe there should be a way to do that more directly.)

[ARK-ARadi]

Hello barrycarter. Given the distance between the sun and earth = 150*10^6 km, and the speed of light in vacuum = 300*10^3 km/s, the time needed for light to travel from the sun to earth would be 500 seconds, around 8 minutes and 20 seconds. Now, whether 8 minutes are negligible or considerable totally depends on the particular application. For example, the sun's altitude wouldn't differ much in 8 minutes around noon, but would differ a lot in 8 minutes around sunset.

[ARK-ARadi]

If you check the documentation, I think you'll find routines that will help you print the numbers in whatever format you want if you dislike the default. :)

a is the point on Earth that the sun is directly above. It takes several minutes for the light to reach Earth so the sun would appear overhead somewhere slightly different. I'll think about whether there should be a way to at least have Skyfield estimate the location from which an object appears overhead. (You could force Skyfield to do it today by directly creating a Geocentric position with the x,y,z position of the Astrometric one you calculated earlier; but maybe there should be a way to do that more directly.)

Thanks again Brandon. The default works just fine :-)

Now to create a geocentric position from the astrometric one, would it be a straight forward vector subtraction like before? Because if it is that straight forward already then you should not waste much time to improve on it. I'll try and be back for more help if my ignorance gets in my way again.

[brandon-rhodes]

Once you'd created the astrometric position a, I think you'd just do: g = Geocentric(a.position) or somesuch to directly copy the x, y, z numbers into a Geocentric object.

[ghost]

Hello barrycarter. Given the distance between the sun and earth = 15010^6 km, and the speed of light in vacuum = 30010^3 km/s, the time needed for light to travel from the sun to earth would be 500 seconds, around 8 minutes and 20 seconds. Now, whether 8 minutes are negligible or considerable totally depends on the particular application. For example, the sun's altitude wouldn't differ much in 8 minutes around noon, but would differ a lot in 8 minutes around sunset.

Here's how I worked it out. Over a year, the Earth completes a full 360 revolution around the Sun. Therefore, in 7 minutes, the Sun's angular position for a geocentric observer would change by about 17 seconds of arc. I would argue that's negligible.

In comparison, the Sun's movement on the equinox for a topographical observer at the equator would change by 1.75 degree, about 365.2425 times as much (the fact that this is the length of the year in days is not a coincidence)

[ARK-ARadi]

Apologies for the rather long delay.

The only way I found to generate a 'complete' geocentric object (for example \) is to follow your earlier example: (sun - earth).at(ts.now()).

Geocentric(earth.at(ts.now()).observe(sun)) complains that an 'Astrometric' object has no attribute 'shape'.

Geocentric(earth.at(ts.now()).observe(sun).apparent()) complains that an 'Apparent' object has no attribute 'shape'.

Geocentric(earth.at(ts.now()).observe(sun).to_skycoord()).subpoint() breaks with the error 'you can only ask for the geographic subpoint of a position measured from Earth's center'.

However, the closest was the object Geocentric(earth.at(ts.now()).observe(sun).to_skycoord()). It is indeed a geocentric object missing time, center, and target information: \<Geocentric position and velocity>

[ghost]

I don't know what I'm talking about, but I'm bored.

Is the subsolar/lunar point actually well defined? The most obvious definition would be the point where a line from the geocenter to the Sun's center intersects the Earth's surface.

However, because of the Earth's ellipsoidity, the Sun isn't necessarily overhead there. In other words, the surface normal vector doesn't necessarily point away from the geocenter.

[brandon-rhodes]

@barrycarter Skyfield uses geodetic latitude and longitude, which indeed does not go through the center of the Earth (usually). Check out the code to see the routine that converges those values.

[ARK-ARadi]

Just wanted to thank you Brandon for the great work. I wish I were a better python programmer to lend you a helping hand, but who knows perhaps I'll be back in a few months ;-)

Kind regards