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

[iwellaway]

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('--------------------------------------------------------------------------------------')

[iwellaway]

PS. Apologies for the poor formatting above, I'll have a look at this later.

[bmorris3]

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?

[iwellaway]

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" <notifications@github.com> Reply to: astropy/astroplan <reply@reply.github.com> Date: Friday, 17 January 2020 at 10:43 To: astropy/astroplan <astroplan@noreply.github.com> Cc: "Wellaway, Ian" <I.J.Wellaway@exeter.ac.uk>, Author <author@noreply.github.com> 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 GitHub<https://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 unsubscribe<https://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>.

[rlabbe]

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?

[StuartLittlefair]

@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.

[bmorris3]

@StuartLittlefair – just a ping to see how progress is coming on the two-pass solution. Hope all is well 😄

[StuartLittlefair]

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...