From fdc1ec3bba4603e5b4c08105be722d08ac4332b4 Mon Sep 17 00:00:00 2001 From: Michael Coughlin Date: Wed, 21 Jun 2023 13:02:45 +0200 Subject: [PATCH 1/7] Add hour angle constraint --- astroplan/constraints.py | 52 ++++++++++++++++++++++++++++- astroplan/tests/test_constraints.py | 17 +++++++++- 2 files changed, 67 insertions(+), 2 deletions(-) diff --git a/astroplan/constraints.py b/astroplan/constraints.py index 61a4cc0d..cca44e7f 100644 --- a/astroplan/constraints.py +++ b/astroplan/constraints.py @@ -35,7 +35,8 @@ "LocalTimeConstraint", "PrimaryEclipseConstraint", "SecondaryEclipseConstraint", "Constraint", "TimeConstraint", "observability_table", "months_observable", "max_best_rescale", - "min_best_rescale", "PhaseConstraint", "is_event_observable"] + "min_best_rescale", "PhaseConstraint", "is_event_observable", + "HourAngleConstraint"] _current_year = time.localtime().tm_year # needed for backward compatibility _current_year_time_range = Time( # needed for backward compatibility @@ -936,6 +937,55 @@ def compute_constraint(self, times, observer=None, targets=None): return mask +class HourAngleConstraint(Constraint): + """ + Constrain the hour angle of a target. + + Parameters + ---------- + min : float or `None` + Minimum hour angle of the target. `None` indicates no limit. + max : float or `None` + Maximum hour angle of the target. `None` indicates no limit. + """ + + def __init__(self, min=-5.5, max=5.5): + self.min = min + self.max = max + + def compute_constraint(self, times, observer, targets): + if times.isscalar: + jds = np.array([times.jd]) + else: + jds = np.array([t.jd for t in times]) + GMST = 18.697374558 + 24.06570982441908 * (jds - 2451545) + GMST = np.mod(GMST, 24) + + lon = observer.location.lon.value / 15 + if targets.size == 1: + lst = np.mod(GMST + lon, 24) + ras = np.tile([targets.ra.hour], len(jds)) + else: + lst = np.tile(np.mod(GMST + lon, 24), (len(targets), 1)) + ras = np.tile([target.ra.hour for target in targets], (len(jds), 1)).T + has = np.mod(lst - ras, 24) + + # Use hours from -12 to 12 + idx = np.where(has > 12)[0] + has[idx] = has[idx] - 24 + + if self.min is None and self.max is not None: + mask = has <= self.max + elif self.max is None and self.min is not None: + mask = self.min <= has + elif self.min is not None and self.max is not None: + mask = (self.min <= has) & (has <= self.max) + else: + raise ValueError("No max and/or min specified in " "HourAngleConstraint.") + + return np.squeeze(mask) + + def is_always_observable(constraints, observer, targets, times=None, time_range=None, time_grid_resolution=0.5*u.hour): """ diff --git a/astroplan/tests/test_constraints.py b/astroplan/tests/test_constraints.py index 0bdbde4c..b7b63802 100644 --- a/astroplan/tests/test_constraints.py +++ b/astroplan/tests/test_constraints.py @@ -20,7 +20,7 @@ TimeConstraint, LocalTimeConstraint, months_observable, max_best_rescale, min_best_rescale, PhaseConstraint, PrimaryEclipseConstraint, SecondaryEclipseConstraint, - is_event_observable) + HourAngleConstraint, is_event_observable) from ..periodic import EclipsingSystem from ..exceptions import MissingConstraintWarning @@ -163,6 +163,21 @@ def test_galactic_plane_separation(): assert np.all(is_constraint_met == [False, True, True]) +def test_hour_angle_constraint(): + time = Time('2003-04-05 06:07:08') + apo = Observer.at_site("APO") + targets = [vega, rigel, polaris] + + constraint = HourAngleConstraint() + + is_constraint_met = constraint(times=time, observer=apo, targets=targets) + assert np.all(is_constraint_met == [False, False, False]) + + time = Time('2003-10-05 06:07:08') + is_constraint_met = constraint(times=time, observer=apo, targets=targets) + assert np.all(is_constraint_met == [True, True, True]) + + # in astropy before v1.0.4, a recursion error is triggered by this test @pytest.mark.skipif('APY_LT104') def test_sun_separation(): From d764cacbfcf8dab1b82bd9d53f3705da4450900f Mon Sep 17 00:00:00 2001 From: Michael Coughlin Date: Wed, 21 Jun 2023 13:16:48 +0200 Subject: [PATCH 2/7] Treat times of length 1 different from multiple times --- astroplan/constraints.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/astroplan/constraints.py b/astroplan/constraints.py index cca44e7f..3a629d1d 100644 --- a/astroplan/constraints.py +++ b/astroplan/constraints.py @@ -966,8 +966,12 @@ def compute_constraint(self, times, observer, targets): lst = np.mod(GMST + lon, 24) ras = np.tile([targets.ra.hour], len(jds)) else: - lst = np.tile(np.mod(GMST + lon, 24), (len(targets), 1)) - ras = np.tile([target.ra.hour for target in targets], (len(jds), 1)).T + if len(jds) == 1: + lst = np.array([np.mod(GMST + lon, 24)] * len(targets)).flatten() + ras = np.array([target.ra.hour for target in targets]) + else: + lst = np.tile(np.mod(GMST + lon, 24), (len(targets), 1)) + ras = np.tile([target.ra.hour for target in targets], len(jds)) has = np.mod(lst - ras, 24) # Use hours from -12 to 12 From 120be20ed537b43be5d6cc42eceb5c9bdab08974 Mon Sep 17 00:00:00 2001 From: Michael Coughlin Date: Wed, 21 Jun 2023 16:24:42 +0200 Subject: [PATCH 3/7] linting --- astroplan/constraints.py | 4 ++-- astroplan/tests/test_constraints.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/astroplan/constraints.py b/astroplan/constraints.py index 3a629d1d..b2de29ef 100644 --- a/astroplan/constraints.py +++ b/astroplan/constraints.py @@ -944,8 +944,8 @@ class HourAngleConstraint(Constraint): Parameters ---------- min : float or `None` - Minimum hour angle of the target. `None` indicates no limit. - max : float or `None` + Minimum hour angle of the target. `None` indicates no limit. + max : float or `None` Maximum hour angle of the target. `None` indicates no limit. """ diff --git a/astroplan/tests/test_constraints.py b/astroplan/tests/test_constraints.py index b7b63802..99055124 100644 --- a/astroplan/tests/test_constraints.py +++ b/astroplan/tests/test_constraints.py @@ -165,7 +165,7 @@ def test_galactic_plane_separation(): def test_hour_angle_constraint(): time = Time('2003-04-05 06:07:08') - apo = Observer.at_site("APO") + apo = Observer.at_site("APO") targets = [vega, rigel, polaris] constraint = HourAngleConstraint() From 20c6263aaae1a1aa9e3a137e1a12478ea5f0f478 Mon Sep 17 00:00:00 2001 From: Michael Coughlin Date: Fri, 23 Jun 2023 10:22:54 +0200 Subject: [PATCH 4/7] Corner cases --- astroplan/constraints.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/astroplan/constraints.py b/astroplan/constraints.py index b2de29ef..ea495817 100644 --- a/astroplan/constraints.py +++ b/astroplan/constraints.py @@ -968,10 +968,13 @@ def compute_constraint(self, times, observer, targets): else: if len(jds) == 1: lst = np.array([np.mod(GMST + lon, 24)] * len(targets)).flatten() - ras = np.array([target.ra.hour for target in targets]) + ras = np.array([target.ra.hour for target in targets]).flatten() else: lst = np.tile(np.mod(GMST + lon, 24), (len(targets), 1)) - ras = np.tile([target.ra.hour for target in targets], len(jds)) + ras = np.tile( + np.array([target.ra.hour for target in targets]).flatten(), + (len(jds), 1), + ).T has = np.mod(lst - ras, 24) # Use hours from -12 to 12 From da48f8d28ba01001a208880b157ffe0d4e296543 Mon Sep 17 00:00:00 2001 From: Michael Coughlin Date: Fri, 23 Jun 2023 18:06:50 +0200 Subject: [PATCH 5/7] Use Astropy by default --- astroplan/constraints.py | 49 +++++++++++++++++++++++----------------- 1 file changed, 28 insertions(+), 21 deletions(-) diff --git a/astroplan/constraints.py b/astroplan/constraints.py index ea495817..1d6fed13 100644 --- a/astroplan/constraints.py +++ b/astroplan/constraints.py @@ -947,35 +947,42 @@ class HourAngleConstraint(Constraint): Minimum hour angle of the target. `None` indicates no limit. max : float or `None` Maximum hour angle of the target. `None` indicates no limit. + use_astropy : boolean + Use Astropy for hour angle calculation """ - def __init__(self, min=-5.5, max=5.5): + def __init__(self, min=-5.5, max=5.5, use_astropy=True): self.min = min self.max = max + self.use_astropy = use_astropy def compute_constraint(self, times, observer, targets): - if times.isscalar: - jds = np.array([times.jd]) - else: - jds = np.array([t.jd for t in times]) - GMST = 18.697374558 + 24.06570982441908 * (jds - 2451545) - GMST = np.mod(GMST, 24) - - lon = observer.location.lon.value / 15 - if targets.size == 1: - lst = np.mod(GMST + lon, 24) - ras = np.tile([targets.ra.hour], len(jds)) + + if self.use_astropy: + has = np.array(observer.target_hour_angle(times, targets).hour) else: - if len(jds) == 1: - lst = np.array([np.mod(GMST + lon, 24)] * len(targets)).flatten() - ras = np.array([target.ra.hour for target in targets]).flatten() + if times.isscalar: + jds = np.array([times.jd]) + else: + jds = np.array([t.jd for t in times]) + GMST = 18.697374558 + 24.06570982441908 * (jds - 2451545) + GMST = np.mod(GMST, 24) + + lon = observer.location.lon.value / 15 + if targets.size == 1: + lst = np.mod(GMST + lon, 24) + ras = np.tile([targets.ra.hour], len(jds)) else: - lst = np.tile(np.mod(GMST + lon, 24), (len(targets), 1)) - ras = np.tile( - np.array([target.ra.hour for target in targets]).flatten(), - (len(jds), 1), - ).T - has = np.mod(lst - ras, 24) + if len(jds) == 1: + lst = np.array([np.mod(GMST + lon, 24)] * len(targets)).flatten() + ras = np.array([target.ra.hour for target in targets]).flatten() + else: + lst = np.tile(np.mod(GMST + lon, 24), (len(targets), 1)) + ras = np.tile( + np.array([target.ra.hour for target in targets]).flatten(), + (len(jds), 1), + ).T + has = np.mod(lst - ras, 24) # Use hours from -12 to 12 idx = np.where(has > 12)[0] From a146ffd80843fb583dd8f9709327095d5afa61fc Mon Sep 17 00:00:00 2001 From: Michael Coughlin Date: Sat, 24 Jun 2023 09:40:06 +0200 Subject: [PATCH 6/7] Add squeeze to prevent extra dimensions --- astroplan/constraints.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/astroplan/constraints.py b/astroplan/constraints.py index 1d6fed13..14bb3cdf 100644 --- a/astroplan/constraints.py +++ b/astroplan/constraints.py @@ -984,6 +984,8 @@ def compute_constraint(self, times, observer, targets): ).T has = np.mod(lst - ras, 24) + # ensures no extra dimensions + has = np.squeeze(has) # Use hours from -12 to 12 idx = np.where(has > 12)[0] has[idx] = has[idx] - 24 From dd81fd2ff313d2c5aea364c4bf5d6066e1a4b1cd Mon Sep 17 00:00:00 2001 From: Michael Coughlin Date: Mon, 26 Jun 2023 18:46:16 +0200 Subject: [PATCH 7/7] Direct HA computation suggestion --- astroplan/constraints.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/astroplan/constraints.py b/astroplan/constraints.py index 14bb3cdf..54a586c7 100644 --- a/astroplan/constraints.py +++ b/astroplan/constraints.py @@ -949,17 +949,26 @@ class HourAngleConstraint(Constraint): Maximum hour angle of the target. `None` indicates no limit. use_astropy : boolean Use Astropy for hour angle calculation + sidereal_model : str + Sidereal time model. + See astropy.time.Time.sidereal_time docs for options. """ - def __init__(self, min=-5.5, max=5.5, use_astropy=True): + def __init__(self, min=-5.5, max=5.5, use_astropy=True, + sidereal_model=None): self.min = min self.max = max self.use_astropy = use_astropy + self.sidereal_model = sidereal_model def compute_constraint(self, times, observer, targets): if self.use_astropy: - has = np.array(observer.target_hour_angle(times, targets).hour) + lst = times.sidereal_time('mean', + longitude='tio', + model=self.sidereal_model) + \ + observer.location.lon + has = np.mod([(lst - target.ra).hour for target in targets], 24) else: if times.isscalar: jds = np.array([times.jd])