astropy / astropy/astroplan

Rise/Meridian/Set calculations - interesting effect of low precision

Open
#453 7 comments 0 reactions 0 assignees View on GitHub
Dominant language
Python
Stars
225
Forks
121
Avg merge
28m
Merged PRs (30d)
1

Description

Hi,

We've been testing the new precision setting in the AstroPlan observer object’s target_rise_time, target_set_time and target_meridian_transit_time functions.

Generally, setting a lower precision drastically reduces the time it takes to run the calculations with only a little effect on the rise/meridian/set times. Our setup used to take over 11 hours and it ran in 3 with a precision of 10 n_grid_points. We only found this difference in precision caused a difference in the rise/set times of a couple of minutes art most. The difference in the meridian could be a little higher but this is understandable with a 'gentle curve' of the target as goes from rising to setting.

What we did find though is quite a large difference in calculations when the declination approaches 90 degrees:
```
--------------------------------------------------------------------------------------
PRECISION RISE SET MERIDIAN
--------------------------------------------------------------------------------------
(RA: 1 DEC: 5 DATETIME: 2020-01-01 00:00:00+00:00)
10 2020-01-01 14:55:49.947 2020-01-01 18:33:05.796 2020-01-01 22:10:22.254
20 2020-01-01 14:54:08.678 2020-01-01 18:33:49.447 2020-01-01 22:11:48.077
50 2020-01-01 14:53:58.219 2020-01-01 18:32:56.463 2020-01-01 22:12:00.698
150 2020-01-01 14:53:54.445 2020-01-01 18:32:59.678 2020-01-01 22:12:04.434
1000 2020-01-01 14:53:53.993 2020-01-01 18:32:59.425 2020-01-01 22:12:04.845
--------------------------------------------------------------------------------------
(RA: 1 DEC: 85 DATETIME: 2020-01-01 00:00:00+00:00)
10 2020-01-01 17:06:32.347 2020-01-01 19:13:02.999 2020-01-01 20:04:58.071
20 2020-01-01 16:40:57.439 2020-01-01 18:33:48.087 2020-01-01 20:23:25.087
50 2020-01-01 16:37:16.543 2020-01-01 18:33:04.621 2020-01-01 20:28:52.744
150 2020-01-01 16:37:07.810 2020-01-01 18:33:04.453 2020-01-01 20:29:04.200
1000 2020-01-01 16:37:04.750 2020-01-01 18:33:04.452 2020-01-01 20:29:04.130
--------------------------------------------------------------------------------------
```
See the larger differnce when Dec is 85. This seems to only start occurring when the Dec is 80deg or above and happens regardless of RA and date.

We think this is caused by the gentle rise/set arc that the target makes at higher declinations. It's not too much of an issue, we can raise the precision to compensate for this if needed and we only use a lower precision for simulations to see how the scheduler runs.

However, if you've any thoughts then we'd be interested to hear them. I've added my test code below. Feel free to have a play:
```python
import datetime
import pytz
from astral import *
import astropy.units as u
from astropy.time import Time
from astropy.coordinates import SkyCoord
from astropy.coordinates import EarthLocation
from astroplan import FixedTarget
from astroplan import Observer

#global variables
INT_Lon = -17.8748333
INT_Lat = 28.76027778
INT_height = 2336
INT_solar_depression = 'astronomical'
INT_sidereal_time_type = 'apparent'
INT_lowest_altitude = 33
utc = pytz.UTC
INT_location = EarthLocation(lat=INT_Lat*u.deg, lon=INT_Lon*u.deg, height=INT_height*u.m)
INT = Observer(location=INT_location, name='INT', timezone='utc')

def get_time(rise_set_meridian, precision, targ, dt_now):
""" get rise/set/meridian time """

time = Time(dt_now)

try:
rise_time = INT.target_rise_time(time, targ, which="next", horizon=INT_lowest_altitude * u.deg, n_grid_points=precision)
meridian_transit_time = INT.target_meridian_transit_time(rise_time, targ, which="next", n_grid_points=precision)
set_time = INT.target_set_time(meridian_transit_time, targ, which="next", horizon=INT_lowest_altitude * u.deg, n_grid_points=precision)

if rise_set_meridian == 'rise':
t = rise_time
elif rise_set_meridian == 'meridian':
t = meridian_transit_time
elif rise_set_meridian == 'set':
t = set_time

except ValueError:
print('Target never rises so no rise, set or meridian ' + str(dt_now))

except TypeError:
print('Gridpoint error ' + str(dt_now))

return t.iso

def print_times(ra_decimal, dec_decimal, dt_now):
""" print times """

# create skycoord object for our target
targ = FixedTarget(coord=SkyCoord(ra=ra_decimal * u.deg, dec=dec_decimal * u.deg))

rise_10 = get_time('rise', 10, targ, dt_now)
rise_20 = get_time('rise', 20, targ, dt_now)
rise_50 = get_time('rise', 50, targ, dt_now)
rise_150 = get_time('rise', 150, targ, dt_now)
rise_1000 = get_time('rise', 1000, targ, dt_now)

meridian_10 = get_time('meridian', 10, targ, dt_now)
meridian_20 = get_time('meridian', 20, targ, dt_now)
meridian_50 = get_time('meridian', 50, targ, dt_now)
meridian_150 = get_time('meridian', 150, targ, dt_now)
meridian_1000 = get_time('meridian', 1000, targ, dt_now)

set_10 = get_time('set', 10, targ, dt_now)
set_20 = get_time('set', 20, targ, dt_now)
set_50 = get_time('set', 50, targ, dt_now)
set_150 = get_time('set', 150, targ, dt_now)
set_1000 = get_time('set', 1000, targ, dt_now)

print('10'.ljust(11) + str(rise_10).ljust(26) + str(meridian_10).ljust(26) + str(set_10).ljust(26))
print('20'.ljust(11) + str(rise_20).ljust(26) + str(meridian_20).ljust(26) + str(set_20).ljust(26))
print('50'.ljust(11) + str(rise_50).ljust(26) + str(meridian_50).ljust(26) + str(set_50).ljust(26))
print('150'.ljust(11) + str(rise_150).ljust(26) + str(meridian_150).ljust(26) + str(set_150).ljust(26))
print('1000'.ljust(11) + str(rise_1000).ljust(26) + str(meridian_1000).ljust(26) + str(set_1000).ljust(26))

# get datetime
dt = datetime.datetime.strptime("2020-01-01 00:00:00", "%Y-%m-%d %H:%M:%S").astimezone(utc)
# dt = datetime.datetime.strptime("2020-07-01 00:00:00", "%Y-%m-%d %H:%M:%S").astimezone(utc)
ra = 1

print('--------------------------------------------------------------------------------------')
print('PRECISION'.ljust(11) + 'RISE'.ljust(26) + 'SET'.ljust(26) + 'MERIDIAN'.ljust(26))
print('--------------------------------------------------------------------------------------')

# check the rise/meridian/set times as a function of DEC
for n in range(5,90,10): # increment DEC from 5-85 using steps of 10
print('(RA: ' + str(ra) + ' DEC: ' + str(n) + ' DATETIME: ' + str(dt) + ')')
print_times(ra, n, dt)
print('--------------------------------------------------------------------------------------')
```

Contributor guide

No contributing guide indexed for this repository

Research direction

Start with the supplied Python reproduction and the Observer.target_rise_time, target_set_time, and target_meridian_transit_time entry points. Run the precision comparison across declinations, especially 80–85 degrees, then determine whether the large differences are expected or indicate a calculation issue. Done means reaching and recording a clear conclusion, with a focused regression test if a correction is warranted.

Written by the indexing model from the issue text.

Assessment

Tech stack
python
Domain
backend
Issue type
Bug
Difficulty
4/5
Estimated time
3-5 days
Activity status
Stale
Clarity
Needs clarification
Newbie friendliness
25/100

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.