Skip to content

Commit

Permalink
Improve naming. Formating. Change default value. Estimate coordinates…
Browse files Browse the repository at this point in the history
… at the center of the time interval instead of the start.

Signed-off-by: Gabriel Emery <gabriel.emery@cta-consortium.org>
  • Loading branch information
gabemery committed Oct 4, 2024
1 parent 9e45af0 commit b554ca0
Show file tree
Hide file tree
Showing 3 changed files with 33 additions and 31 deletions.
10 changes: 6 additions & 4 deletions gammapy/makers/map.py
Original file line number Diff line number Diff line change
Expand Up @@ -39,8 +39,10 @@ class MapDatasetMaker(Maker):
Pad one bin in offset for 2d background map.
This avoids extrapolation at edges and use the nearest value.
Default is True.
fov_rotation_error_limit : `~astropy.units.Quantity`
fov_rotation_step : `~astropy.units.Quantity`, optional
Maximum error on the rotation angle between AltAz and RaDec frames during background evaluation.
Used only when the Background IRF has an AltAz alignement.
Default is 0.5 deg.
Examples
--------
Expand Down Expand Up @@ -109,12 +111,12 @@ def __init__(
background_oversampling=None,
background_interp_missing_data=True,
background_pad_offset=True,
fov_rotation_error_limit=0.5 * u.deg,
fov_rotation_step=1.0 * u.deg,
):
self.background_oversampling = background_oversampling
self.background_interp_missing_data = background_interp_missing_data
self.background_pad_offset = background_pad_offset
self.fov_rotation_error_limit = fov_rotation_error_limit
self.fov_rotation_step = fov_rotation_step
if selection is None:
selection = self.available_selection

Expand Down Expand Up @@ -242,7 +244,7 @@ def make_background(self, geom, observation):
ontime=observation.observation_time_duration,
bkg=bkg,
geom=geom,
fov_rotation_error_limit=self.fov_rotation_error_limit,
fov_rotation_step=self.fov_rotation_step,
oversampling=self.background_oversampling,
use_region_center=use_region_center,
obstime=observation.tmid,
Expand Down
12 changes: 6 additions & 6 deletions gammapy/makers/tests/test_utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -190,7 +190,7 @@ def test_map_background_2d(bkg_2d, fixed_pointing_info):
ontime="42 s",
bkg=bkg_2d,
geom=geom,
fov_rotation_error_limit=0.5 * u.deg,
fov_rotation_step=1.0 * u.deg,
)

assert_allclose(bkg.data[:, 1, 1], [1.869025, 0.186903], rtol=1e-5)
Expand All @@ -201,7 +201,7 @@ def test_map_background_2d(bkg_2d, fixed_pointing_info):
ontime="42 s",
bkg=bkg_2d,
geom=geom,
fov_rotation_error_limit=0.5 * u.deg,
fov_rotation_step=1.0 * u.deg,
obstime=obstime,
)
assert_allclose(bkg.data, bkg_fpi.data, rtol=1e-5)
Expand All @@ -215,7 +215,7 @@ def make_map_background_irf_with_symmetry(fpi, symmetry="constant"):
ontime="42 s",
bkg=bkg_3d_custom(symmetry),
geom=WcsGeom.create(npix=(3, 3), binsz=4, axes=[axis], skydir=fpi.fixed_icrs),
fov_rotation_error_limit=0.5 * u.deg,
fov_rotation_step=1.0 * u.deg,
obstime=obstime,
)

Expand Down Expand Up @@ -265,7 +265,7 @@ def test_make_map_background_irf(bkg_3d, pars, fixed_pointing_info):
ebounds=pars["ebounds"],
skydir=fixed_pointing_info.fixed_icrs,
),
fov_rotation_error_limit=0.5 * u.deg,
fov_rotation_step=1.0 * u.deg,
oversampling=10,
obstime=Time("2020-01-01T20:00"),
)
Expand Down Expand Up @@ -326,7 +326,7 @@ def test_make_map_background_irf_skycoord(fixed_pointing_info_aligned):
ontime="42 s",
bkg=bkg_3d_custom("asymmetric", "ALTAZ"),
geom=WcsGeom.create(npix=(3, 3), binsz=4, axes=[axis], skydir=position),
fov_rotation_error_limit=0.5 * u.deg,
fov_rotation_step=1.0 * u.deg,
)


Expand All @@ -344,7 +344,7 @@ def test_make_map_background_irf_AltAz_align(fixed_pointing_info):
axes=[axis],
skydir=fixed_pointing_info.get_icrs(obstime),
),
fov_rotation_error_limit=0.5 * u.deg,
fov_rotation_step=1.0 * u.deg,
obstime=obstime,
)

Expand Down
42 changes: 21 additions & 21 deletions gammapy/makers/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,8 @@

log = logging.getLogger(__name__)

MINIMUM_TIME_STEP = 1 * u.s # Minimum time step used to handle FoV rotations


def _get_fov_coords(pointing, irf, geom, use_region_center=True, obstime=None):
# TODO: create dedicated coordinate handling see #5041
Expand Down Expand Up @@ -98,17 +100,17 @@ def _get_fov_coords(pointing, irf, geom, use_region_center=True, obstime=None):
return coords


def time_limit_rotation(maximum_rotation, pointing_altaz):
def compute_rotation_time_step(rotation, pointing_altaz):
"""
Compute the maximum time such that the FoV associated to a fixed RaDec position rotates by less than
'maximum_rotation' in AltAz frame. It assumes that the rotation rate at the provided pointing is a good estimate
of the rotation rate over the full output duration (first order approximation).
Compute the time such that the FoV associated to a fixed RaDec position rotates by 'rotation' in AltAz frame.
It assumes that the rotation rate at the provided pointing is a good estimate of the rotation rate over the full
output duration (first order approximation).
Parameters
----------
maximum_rotation : `~astropy.units.Quantity`
Maximum rotation angle.
rotation : `~astropy.units.Quantity`
Rotation angle.
pointing_altaz : `~astropy.coordinates.SkyCoord`
Pointing direction.
Expand All @@ -117,17 +119,15 @@ def time_limit_rotation(maximum_rotation, pointing_altaz):
duration : `~astropy.units.Quantity`
Time associated with the requested rotation.
"""
earth_rotation_deg_per_second = 0.004167
earth_rotation = 360 * u.deg / u.day
denom = (
earth_rotation_deg_per_second
earth_rotation
* np.cos(pointing_altaz.location.lat.rad)
* np.cos(pointing_altaz.az.rad)
* np.abs(np.cos(pointing_altaz.az.rad))
)
if denom == 0:
if denom.value == 0:
return 3600 * u.s
return (
maximum_rotation.to_value(u.deg) * np.cos(pointing_altaz.alt.rad) / denom
) * u.s
return rotation * np.cos(pointing_altaz.alt.rad) / denom


def make_map_exposure_true_energy(
Expand Down Expand Up @@ -223,7 +223,7 @@ def evaluate_bkg(pointing, ontime, bkg, geom, use_region_center, obstime, d_omeg
irf=bkg,
geom=geom,
use_region_center=use_region_center,
obstime=obstime,
obstime=obstime + ontime / 2,
)
coords["energy"] = broadcast_axis_values_to_geom(geom, "energy", False)

Expand All @@ -236,7 +236,7 @@ def make_map_background_irf(
ontime,
bkg,
geom,
fov_rotation_error_limit,
fov_rotation_step,
oversampling=None,
use_region_center=True,
obstime=None,
Expand All @@ -260,7 +260,7 @@ def make_map_background_irf(
Background rate model.
geom : `~gammapy.maps.WcsGeom`
Reference geometry.
fov_rotation_error_limit : `~astropy.units.Quantity`
fov_rotation_step : `~astropy.units.Quantity`
Maximum rotation error to create sub-time bins if the irf is 3D and in AltAz.
oversampling : int
Oversampling factor in energy, used for the background model evaluation.
Expand Down Expand Up @@ -302,17 +302,17 @@ def make_map_background_irf(
endtime = obstime + ontime
data = np.zeros(geom.data_shape)
while obstime < endtime:
# Evaluate the step time needed to have FoV rotation of less than fov_rotation_error_limit
# Evaluate the step time needed to have a FoV rotation of fov_rotation_step
# Minimum step of 1 second to avoid infinite computation very close to zenith
steptime = max(
time_limit_rotation(
fov_rotation_error_limit, pointing.get_altaz(obstime)
compute_rotation_time_step(
fov_rotation_step, pointing.get_altaz(obstime)
),
1 * u.s,
MINIMUM_TIME_STEP,
)
steptime = (
steptime
if steptime < endtime - obstime - 1 * u.s
if steptime < endtime - obstime - MINIMUM_TIME_STEP
else endtime - obstime
)
data += evaluate_bkg(
Expand Down

0 comments on commit b554ca0

Please sign in to comment.