astroplan
astroplan copied to clipboard
astropy
Rise/Meridian/Set calculations - interesting effect of low precision
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:
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('--------------------------------------------------------------------------------------')
PS. Apologies for the poor formatting above, I'll have a look at this later.
Thanks for the note! That's really interesting, and I think your interpretation makes sense. Do you think it'd be a reasonable solution for this use-case to use a higher n_grid_points value for the stars at higher declinations, and lower n_grid_points values for, say, targets with dec < 80 deg?
Hi Brett,
Yes that might work. I currently use a global variable to set which precision level should be used in the code as we have calculated several levels of precision per target to test the difference it makes. However, it would be trivial to amend this slightly to use different levels of precision depending on the declination of the target.
Thanks for your help.
Best regards
Ian
From: "Brett M. Morris" [email protected] Reply to: astropy/astroplan [email protected] Date: Friday, 17 January 2020 at 10:43 To: astropy/astroplan [email protected] Cc: "Wellaway, Ian" [email protected], Author [email protected] Subject: Re: [astropy/astroplan] Rise/Meridian/Set calculations - interesting effect of low precision (#453)
Thanks for the note! That's really interesting, and I think your interpretation makes sense. Do you think it'd be a reasonable solution for this use-case to use a higher n_grid_points value for the stars at higher declinations, and lower n_grid_points values for, say, targets with dec < 80 deg?
— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fastropy%2Fastroplan%2Fissues%2F453%3Femail_source%3Dnotifications%26email_token%3DAAMX4SBALUBOAUEQ4O2BYI3Q6GDFTA5CNFSM4KIGENDKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEJHIUAQ%23issuecomment-575572482&data=02%7C01%7CI.J.Wellaway%40exeter.ac.uk%7C528c136d28844d64753208d79b3a1354%7C912a5d77fb984eeeaf321334d8f04a53%7C0%7C0%7C637148546039433927&sdata=mosn85bcFbmAcSIXM8wXkJIoca9RKmLL4aNzlb%2BvZ94%3D&reserved=0, or unsubscribehttps://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FAAMX4SCWEZXNL2UNBWXWHELQ6GDFTANCNFSM4KIGENDA&data=02%7C01%7CI.J.Wellaway%40exeter.ac.uk%7C528c136d28844d64753208d79b3a1354%7C912a5d77fb984eeeaf321334d8f04a53%7C0%7C0%7C637148546039438922&sdata=UYEkZ4PxBiQu5yKrTAcKaU7J%2BqaxCVNPHNHyFdavy2w%3D&reserved=0.
I'm not an astrophysicist, so I hope I am not making a dumb mistake or assumption. but I noticed different times being computed for solar noon, and wrote this little test program:
today = Time('2020-02-05')
sun = get_sun(today)
for lat in [-89, -80, -70, -60, -50, -40, -30, -20, -10, 0, 10,20,30,40,50,60,70,80,89, 89.9]:
obs1 = Observer(longitude=0, latitude=lat)
noon = obs1.target_meridian_transit_time(today, sun, n_grid_points=1000)
noon.format = 'datetime'
print(f'{lat: 4.1f} {noon} {noon.scale}')
Which yields
-89.0 2020-02-04 12:03:03.850783 utc
-80.0 2020-02-04 12:12:50.643175 utc
-70.0 2020-02-04 12:13:23.874014 utc
-60.0 2020-02-04 12:13:35.416659 utc
-50.0 2020-02-04 12:13:41.559173 utc
-40.0 2020-02-04 12:13:45.568566 utc
-30.0 2020-02-04 12:13:48.543645 utc
-20.0 2020-02-04 12:13:50.951276 utc
-10.0 2020-02-04 12:13:53.094335 utc
0.0 2020-02-04 12:13:55.104302 utc
10.0 2020-02-04 12:13:57.107670 utc
20.0 2020-02-04 12:13:59.238538 utc
30.0 2020-02-04 12:14:01.661216 utc
40.0 2020-02-04 12:14:04.633117 utc
50.0 2020-02-04 12:14:08.637722 utc
60.0 2020-02-04 12:14:14.774603 utc
70.0 2020-02-04 12:14:26.311093 utc
80.0 2020-02-04 12:14:59.536581 utc
89.0 2020-02-04 12:24:46.493044 utc
89.9 2020-02-04 14:07:09.982392 utc
Should these not all be the same value?
@rlabbe in fact, they shouldn’t all be the same time. The UTC time of solar noon moves around a bit because the Earth’s orbit is not circular. If you look for photos of the analemma you can see a striking visualisation of this.
@iwellaway using a two pass approach mostly fixes this issue. I’ve got a Code branch where I’m working on it, but there are a few subtleties I want to iron out first.
@StuartLittlefair – just a ping to see how progress is coming on the two-pass solution. Hope all is well 😄
Hi @bmorris3 ! All well here thanks.
The status of this is that I've a working solution that is really nice, except that it doesn't handle cases where objects don't rise/set at the moment. The issue is that the first pass returns nan for the rise-set time, so the approach of making a fine grid around the time suggested by the first pass doesn't work...
I'll hopefully get some time to work on this over the next few weeks...