Defining an equally spaced range of times using only years argument gives a sawtooth

[davidmikolas]

The following script causes the Moon to jump, likely because the the years argument doesn't like tiny spaces.

1. Is this the best behavior for float years?
2. Is there a better way to get an evenly distributed series of times over decades with sub-day step size?

from skyfield.api import Loader, Topos
import numpy as np
import matplotlib.pyplot as plt

load = Loader('~/Documents/fishing/SkyData') # avoid multiple copies of large files
ts = load.timescale()
eph = load('de421.bsp')

earth, sun, moon = [eph[x] for x in ('earth', 'sun', 'moon')]
years = 2020 + np.arange(0, 20, 0.0001)
times = ts.utc(years)

antarctica = earth + Topos('90.0 S', '0.0 E')

moon_apparent = antarctica.at(times).observe(moon).apparent()
sun_apparent = antarctica.at(times).observe(sun).apparent()

malt, maz = [x.degrees for x in moon_apparent.altaz()[:2]]
salt, saz = [x.degrees for x in sun_apparent.altaz()[:2]]

fig = plt.figure()
limits = ((2019.9, 2040.1), (2019.995, 2021.005), (2019.9995, 2020.1005))
axes = [fig.add_subplot(3, 1, n) for n in (1, 2, 3)]
for ax, lims in zip(axes, limits):
    ax.plot(years, malt)
    ax.set_xlim(*lims)
    ax.get_xaxis().get_major_formatter().set_useOffset(False)
plt.suptitle('Elevation of Moon from South Pole')
plt.show()

[brandon-rhodes]

I will have to think about this one. The sensible move would be to outlaw non-integral UTC years, as UTC does not define any notion (so far as I'm aware) of what a fractional UTC year would mean. The exception could point users to a new method that could be added, ts.J(floating_point_year), that would be the inverse of the existing routine t.J:

python-skyfield/skyfield/timelib.py

Line 673 in f2fe099

@reify

But that could cause existing Skyfield programs to die with an error. It might be worth the cost but I'll have to think about it.

The only other reasonable alternative I can think of would be to interpret the fraction as fractions of whatever length that particular year is — so 2011.5 would mean "2011 plus half of 365 days" while 2012.5 would mean "2012 plus half of 366 days". All programs that currently work would hopefully continue to work (but would give somewhat different answers than today, as the stair-stepping you see would disappear), but the computation would be a bit more costly — and would introduce a new Skyfield-specific non-uniform time scale "UTC years" that might generate confusion since the fraction would mean different amounts in different years.

Hmm.

[davidmikolas]

What where the uses you did expect and plan for?

For example:

• were arrays of hours from say 0 to 1000 okay, or for days, or for months?
• was the use of float arrays in any field besides seconds not anticipated?

What did you envision exactly for each field?

I ask because you have excellent intuition, so it's likely you had a plan for this but since you've been working on so many aspects in parallel you might need to re-remember it.

[brandon-rhodes]

My basic vision was:

• Integers make sense in every field.
• Fractions make sense for seconds.

I'm not sure fractions make sense in any of the other fields — or at least they definitely won't create motion whose derivative is stable, because (a) the length of the minute, and thus of the hour and day, jumps upward temporarily whenever there is a leap second, and (b) the lengths of months and years are of course not equal.

[glangford]

@davidmikolas For your question 2

Is there a better way to get an evenly distributed series of times over decades with sub-day step size?

You could break the time generation into subparts, something along the lines of

for i in range(1, 365*20): # 20 year sample daily
    times = ts.utc(2020, 1, i, 0, 0, range(0, 60*60*24, 300)) # single day sample, every 5 minutes
    ...

and append the computations inside the loop as appropriate (note the index i in the creation of times). I use this method in a different plot and it works well. In addition,

times._nutation_angles = iau2000b_radians(times)

will be quicker.

ps. this looks a bit odd but I think it will work because days can count beyond a month, eg.

>>> times = ts.utc(2020, 1, 954, 0, 0, range(1, 100)) # 954 is random
>>> times
<Time tt=[2459802.500812315 ... 2459802.501946574] len=99>
>>> times[0].utc_iso()
'2022-08-11T00:00:01Z'

[davidmikolas]

@glangford Thank you for the suggestions. My understanding of time (in Skyfield, in Astronomy, and even in the real world) is limited! My goal here was to get a time array that spanned 20 years, spaced at roughly 1 hour intervals.

Moving the fractions from year to hour seems to work, checking several derivatives of malt shows no more discontinuities or staircasing effects. This is good enough for me, for now, since I'm just making a plot to look at.

dyears = np.arange(0, 20, 0.0001)
years = dyears + 2020.
hours = dyears * 365.2564 * 24

times = ts.utc(2020, 1, 1, hours)

@brandon-rhodes But it would be great to have a time generator in Skyfield that simply accepts a start time object and an end time object and either a number of points or a spacing that worked just like np.linspace(start, stop, number) or np.arrange(start, stop, step) except that the start and stop parameters were Skyfield time objects, and perhaps step should be hard-wired as seconds to keep things simple.

I suppose they would need one more argument that specified in which time scale they should be equally spaced or perhaps there is one obvious choice (perhaps TDB) and it would simply use that.

[brandon-rhodes]

Having thought about it for a few days, I am not quite ready yet to duplicate or wrap linspace() or arange() as I would still prefer for Skyfield's docs to teach good NumPy techniques instead of duplicating them over again in its own API. But I have indeed gone ahead and added Julian years to the Timescale object so that a range of years can be created through ts.J(range(2000, 2020)) and so forth. The feature should come out in the next release, hopefully by this weekend. I'll also see about updating docs to use it in examples when I need a range of years!

[davidmikolas]

Thank you very much for the speedy disposition, conclusion, and solution!