diff --git a/pygdsm/base_observer.py b/pygdsm/base_observer.py index 8013aa4..fa422b1 100644 --- a/pygdsm/base_observer.py +++ b/pygdsm/base_observer.py @@ -40,10 +40,11 @@ def _setup(self): self._pix0 = None self._mask = None + self._horizon_elevation = 0.0 self._observed_ra = None self._observed_dec = None - def generate(self, freq=None, obstime=None): + def generate(self, freq=None, obstime=None, horizon_elevation=None): """ Generate the observed sky for the observer, based on the GSM. Parameters @@ -52,6 +53,8 @@ def generate(self, freq=None, obstime=None): Frequency of map to generate, in units of MHz (default). obstime: astropy.time.Time Time of observation to generate + horizon_elevation: float + Elevation of the artificial horizon (default 0.0) Returns ------- @@ -75,8 +78,20 @@ def generate(self, freq=None, obstime=None): self._time = Time(obstime) # This will catch datetimes, but Time() object should be passed self.date = obstime.to_datetime() + # Match pyephem convertion -- string is degrees, int/float is rad + horizon_elevation = ephem.degrees(horizon_elevation or 0.0) + if self._horizon_elevation == horizon_elevation: + horizon_has_changed = False + else: + self._horizon_elevation = horizon_elevation + horizon_has_changed = True + + # Checking this separately encase someone tries to be smart and sets both the attribute and kwarg + if self._horizon_elevation < 0: + raise ValueError(f"Horizon elevation must be greater or equal to 0 degrees ({np.rad2deg(horizon_elevation)=}).") + # Rotation is quite slow, only recompute if time or frequency has changed, or it has never been run - if time_has_changed or self.observed_sky is None: + if time_has_changed or self.observed_sky is None or horizon_has_changed: # Get RA and DEC of zenith ra_zen, dec_zen = self.radec_of(0, np.pi/2) sc_zen = SkyCoord(ra_zen, dec_zen, unit=('rad', 'rad')) @@ -89,7 +104,7 @@ def generate(self, freq=None, obstime=None): # Generate below-horizon mask using query_disc mask = np.ones(shape=self._n_pix, dtype='bool') - pix_visible = hp.query_disc(self._n_side, vec=vec_zen, radius=np.pi/2) + pix_visible = hp.query_disc(self._n_side, vec=vec_zen, radius=np.pi/2 - self._horizon_elevation) mask[pix_visible] = 0 self._mask = mask diff --git a/tests/test_gsm2016.py b/tests/test_gsm2016.py index ca9cf70..ec314e6 100644 --- a/tests/test_gsm2016.py +++ b/tests/test_gsm2016.py @@ -45,6 +45,29 @@ def test_observer_test(): d = ov.view(logged=True) plt.show() + ov = GSMObserver16() + ov.lon = longitude + ov.lat = latitude + ov.elev = elevation + ov.date = datetime(2000, 1, 1, 23, 0) + + horizon_elevation = 85.0 + ov.generate(1400, horizon_elevation=str(horizon_elevation)) + d_85deg_horizon = ov.view(logged=True) + + ov = GSMObserver16() + ov.lon = longitude + ov.lat = latitude + ov.elev = elevation + ov.date = datetime(2000, 1, 1, 23, 0) + + ov.generate(1400, horizon_elevation=np.deg2rad(horizon_elevation)) + d_85deg2rad_horizon = ov.view(logged=True) + + assert np.all(d_85deg_horizon == d_85deg2rad_horizon), "The two methods for calculating the artificial horizon do not match." + assert np.ma.count_masked(d_85deg_horizon).sum() == 12558953 + plt.show() + def test_interp(): f = np.arange(40, 80, 5) for interp in ('pchip', 'cubic'): diff --git a/tests/test_gsm_observer.py b/tests/test_gsm_observer.py index 4382d42..482d9a6 100644 --- a/tests/test_gsm_observer.py +++ b/tests/test_gsm_observer.py @@ -5,6 +5,7 @@ import numpy as np import os from astropy.time import Time +import pytest def test_gsm_observer(show=False): """ Test GSMObserver() is working @@ -42,6 +43,10 @@ def test_gsm_observer(show=False): ov.generate(50) ov.view(logged=True) ov.view_observed_gsm(logged=True) + + with pytest.raises(ValueError) as e: + ov.generate(horizon_elevation=-1e-3) + if show: plt.show() @@ -57,9 +62,11 @@ def test_observed_mollview(): obs = [] if not os.path.exists('generated_sky'): os.mkdir('generated_sky') + freq = 50 + horizon_elevation = '85.0' for ii in range(0, 24, 4): ov.date = datetime(2000, 1, 1, ii, 0) - ov.generate(50) + ov.generate(freq) sky = ov.view_observed_gsm(logged=True, show=False, min=9, max=20) plt.savefig('generated_sky/galactic-%02d.png' % ii) plt.close() @@ -76,9 +83,15 @@ def test_observed_mollview(): plt.savefig('generated_sky/ortho-%02d.png' % ii) plt.close() + ov.generate(freq=freq, horizon_elevation=horizon_elevation) + ov.view(logged=True, show=False, min=9, max=20) + plt.savefig('generated_sky/ortho_85deg_horizon-%02d.png' % ii) + plt.close() + print(ii) os.system('convert -delay 20 generated_sky/ortho-*.png ortho.gif') + os.system('convert -delay 20 generated_sky/ortho_85deg_horizon-*.png ortho_85deg_horizon.gif') os.system('convert -delay 20 generated_sky/galactic-*.png galactic.gif') os.system('convert -delay 20 generated_sky/ecliptic-*.png ecliptic.gif') os.system('convert -delay 20 generated_sky/equatorial-*.png equatorial.gif') @@ -106,6 +119,10 @@ def test_generate_with_and_without_args(): ov.generate(obstime=now) ov.generate(obstime=now, freq=53) ov.generate(obstime=now, freq=52) + ov.generate(obstime=now, freq=53, horizon_elevation=0.0) + ov.generate(obstime=now, freq=52, horizon_elevation='0.0') + ov.generate(obstime=now, freq=53, horizon_elevation=np.deg2rad(85.0)) + ov.generate(obstime=now, freq=52, horizon_elevation='85.0') if __name__ == "__main__": test_gsm_observer(show=True) diff --git a/tests/test_haslam.py b/tests/test_haslam.py index a8aaf26..23c00c6 100644 --- a/tests/test_haslam.py +++ b/tests/test_haslam.py @@ -42,6 +42,29 @@ def test_observer_test(): d = ov.view(logged=True) plt.show() + ov = GSMObserver() + ov.lon = longitude + ov.lat = latitude + ov.elev = elevation + ov.date = datetime(2000, 1, 1, 23, 0) + + horizon_elevation = 85.0 + ov.generate(1400, horizon_elevation=str(horizon_elevation)) + d_85deg_horizon = ov.view(logged=True) + + ov = GSMObserver() + ov.lon = longitude + ov.lat = latitude + ov.elev = elevation + ov.date = datetime(2000, 1, 1, 23, 0) + + ov.generate(1400, horizon_elevation=np.deg2rad(horizon_elevation)) + d_85deg2rad_horizon = ov.view(logged=True) + + assert np.all(d_85deg_horizon == d_85deg2rad_horizon), "The two methods for calculating the artificial horizon do not match." + assert np.ma.count_masked(d_85deg_horizon).sum() == 3139749 + plt.show() + def test_cmb_removal(): g = HaslamSkyModel(freq_unit='MHz', include_cmb=False) sky_no_cmb = g.generate(400) diff --git a/tests/test_lfsm.py b/tests/test_lfsm.py index 1735c27..cd9e0a1 100644 --- a/tests/test_lfsm.py +++ b/tests/test_lfsm.py @@ -44,6 +44,29 @@ def test_observer_test(): d = ov.view(logged=True) plt.show() + ov = LFSMObserver() + ov.lon = longitude + ov.lat = latitude + ov.elev = elevation + ov.date = datetime(2000, 1, 1, 23, 0) + + horizon_elevation = 85.0 + ov.generate(200, horizon_elevation=str(horizon_elevation)) + d_85deg_horizon = ov.view(logged=True) + + ov = LFSMObserver() + ov.lon = longitude + ov.lat = latitude + ov.elev = elevation + ov.date = datetime(2000, 1, 1, 23, 0) + + ov.generate(200, horizon_elevation=np.deg2rad(horizon_elevation)) + d_85deg2rad_horizon = ov.view(logged=True) + + assert np.all(d_85deg_horizon == d_85deg2rad_horizon), "The two methods for calculating the artificial horizon do not match." + assert np.ma.count_masked(d_85deg_horizon).sum() == 784951 + plt.show() + def test_cmb_removal(): g = LowFrequencySkyModel(freq_unit='MHz', include_cmb=False) sky_no_cmb = g.generate(400)