Skip to content

Commit

Permalink
reverted the pointing API to the old behaviour, documentation updated
Browse files Browse the repository at this point in the history
  • Loading branch information
paganol committed Nov 5, 2024
1 parent 83a1364 commit c6c6021
Show file tree
Hide file tree
Showing 9 changed files with 55 additions and 34 deletions.
35 changes: 25 additions & 10 deletions docs/source/scanning.rst
Original file line number Diff line number Diff line change
Expand Up @@ -168,15 +168,15 @@ similar to what is going to be used for LiteBIRD:

.. testoutput::

Shape: (1, 600, 3)
Shape: (600, 3)
Pointings:
[[[ 2.182 -0. -1.571]
[[ 2.182 -0. -1.571]
[ 2.182 -0.006 -1.576]
[ 2.182 -0.012 -1.582]
...
[ 0.089 -2.967 -1.738]
[ 0.088 -3.021 -1.687]
[ 0.087 -3.075 -1.635]]]
[ 0.087 -3.075 -1.635]]

All the details in this code are explained in the next sections, so
for now just keep in mind the overall shape of the code:
Expand All @@ -201,12 +201,27 @@ for now just keep in mind the overall shape of the code:
the example above, this is done by the function
:func:`.get_pointings`.

3. The method :meth:`.Observation.get_pointings` returns a ``(D, N, 3)``
matrix, where D represents the detector index, N the index of the sample
and the three final columns contain the colatitude :math:`\theta`,
the longitude :math:`\phi`, and the orientation angle :math:`\psi`,
all expressed in radians. These angles are expressed in the Ecliptic
Coordinate System, where the Equator is aligned with the Ecliptic Plane of
3. The method :meth:`.Observation.get_pointings` returns an array with
either 2 or 3 fields depending on the argument passed:

- if an integer is passed, this is interpreted as the index of the
detector in the observation, and a ``(N, 3)`` matrix is returned
where the first column contains the colatitude :math:`\theta`,
the second column the longitude :math:`\phi`, and the third column
the orientation angle :math:`\psi`, all expressed in radians.

- if a list containing indices is passed, this is interpreted as
a list of detectors in the observation for which we want to compute
the pointing. It returns a ``(D, N, 3)`` matrix where D represents
the detector index, N the index of the sample and the three final
columns are the same described in the first case.

- if the string "all" is passed then a ``(D, N, 3)`` matrix is returned
containig the pointing information for all the detectors in the
observation.

These angles are expressed in the Ecliptic Coordinate
System, where the Equator is aligned with the Ecliptic Plane of
the Solar System.


Expand Down Expand Up @@ -361,7 +376,7 @@ split in several blocks inside the :class:`.Observation` class.
unlikely to be relevant.

Once all the quaternions have been computed at the proper sampling
rate, the direction of the detector on the sky and its o]rientation
rate, the direction of the detector on the sky and its orientation
angle can be computed via a call to :meth:`.Observation.get_pointings`.
The calculation works as follows:

Expand Down
2 changes: 1 addition & 1 deletion litebird_sim/dipole.py
Original file line number Diff line number Diff line change
Expand Up @@ -225,7 +225,7 @@ def add_dipole(
if type(pointings) is np.ndarray:
theta_phi_det = pointings[detector_idx, :, :]
else:
theta_phi_det = pointings(detector_idx)[0][0, :, 0:2]
theta_phi_det = pointings(detector_idx)[0][:, 0:2]

add_dipole_for_one_detector(
tod_det=tod[detector_idx],
Expand Down
1 change: 0 additions & 1 deletion litebird_sim/hwp_sys/hwp_sys.py
Original file line number Diff line number Diff line change
Expand Up @@ -1256,7 +1256,6 @@ def fill_tod(
cur_point, cur_hwp_angle = cur_obs.get_pointings(
detector_idx=idet, pointings_dtype=dtype_pointings
)
cur_point = cur_point.reshape(-1, 3)
else:
cur_point = ptg_list[idx_obs][idet, :, :]
cur_hwp_angle = hwp_angle_list[idx_obs]
Expand Down
1 change: 0 additions & 1 deletion litebird_sim/mapmaking/common.py
Original file line number Diff line number Diff line change
Expand Up @@ -205,7 +205,6 @@ def _compute_pixel_indices(
curr_pointings_det = pointings[idet, :, :]
else:
curr_pointings_det, hwp_angle = pointings(idet)
curr_pointings_det = curr_pointings_det.reshape(-1, 3)

if hwp_angle is None:
hwp_angle = 0
Expand Down
15 changes: 12 additions & 3 deletions litebird_sim/observations.py
Original file line number Diff line number Diff line change
Expand Up @@ -680,20 +680,30 @@ def get_pointings(
# All the detectors are included
pointings, hwp_angle = cur_obs.get_pointings("all")
# This returns ``(N_det, N_samples, 3)`` array
# Only the first two detectors are included
pointings, hwp_angle = cur_obs.get_pointings([0, 1])
# This returns ``(2, N_samples, 3)`` array
# Only the first detector is used
pointings, hwp_angle = cur_obs.get_pointings(0)
# This returns ``(N_samples, 3)`` array
# NB if a list of indices is passed an array of dimension ``(N_det, N_samples, 3)``
# is always returned
# For example this
pointings, hwp_angle = cur_obs.get_pointings([0])
# returns ``(1, N_samples, 3)`` array
The return value is a pair containing (1) the pointing matrix and (2) the
HWP angle. The pointing matrix is a NumPy array with shape ``(N_det, N_samples, 3)``,
where ``N_det` is the number of detectors and ``N_samples`` is the number of
samples in the TOD (the field ``Observation.n_samples``). The last dimension
spans the three angles θ (colatitude, in radians), φ (longitude, in radians),
and ψ (orientation angle, in radians). *Important*: if you ask for just *one*
detector, the shape of the pointing matrix will always be ``(N_samples, 3)``.
detector passing the index of the detector, the shape of the pointing matrix
will always be ``(N_samples, 3)``.
The HWP angle is always a vector with shape ``(N_samples,)``, as it does
not depend on the list of detectors.
Expand All @@ -712,7 +722,7 @@ def get_pointings(
(detector_idx >= 0) and (detector_idx < self.n_detectors)
), f"Invalid detector index {detector_idx}, it must be a number between 0 and {self.n_detectors - 1}"

pointing, hwp = self.pointing_provider.get_pointings(
return self.pointing_provider.get_pointings(
detector_quat=self.quat[detector_idx],
start_time=self.start_time,
start_time_global=self.start_time_global,
Expand All @@ -722,7 +732,6 @@ def get_pointings(
hwp_buffer=hwp_buffer,
pointings_dtype=pointings_dtype,
)
return pointing.reshape(1, -1, 3), hwp

# More complex case: we need all the detectors
if isinstance(detector_idx, str):
Expand Down
1 change: 0 additions & 1 deletion litebird_sim/scan_map.py
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,6 @@ def scan_map(
curr_pointings_det = pointings[detector_idx, :, :]
else:
curr_pointings_det, hwp_angle = pointings(detector_idx)
curr_pointings_det = curr_pointings_det.reshape(-1, 3)

if hwp_angle is None:
hwp_angle = 0
Expand Down
10 changes: 5 additions & 5 deletions test/test_coordinates.py
Original file line number Diff line number Diff line change
Expand Up @@ -49,14 +49,14 @@ def test_coordinates():

r = hp.Rotator(coord=["E", "G"])

pointings_gal_hp = np.empty_like(pointings[0])
pointings_gal_hp = np.empty_like(pointings)

pointings_gal_hp[:, 0:2] = r(pointings[0, :, 0], pointings[0, :, 1]).T
pointings_gal_hp[:, 2] = pointings[0, :, 2] + r.angle_ref(
pointings[0, :, 0], pointings[0, :, 1]
pointings_gal_hp[:, 0:2] = r(pointings[:, 0], pointings[:, 1]).T
pointings_gal_hp[:, 2] = pointings[:, 2] + r.angle_ref(
pointings[:, 0], pointings[:, 1]
)

pointings_gal_lbs = lbs.coordinates.rotate_coordinates_e2g(pointings[0])
pointings_gal_lbs = lbs.coordinates.rotate_coordinates_e2g(pointings)

np.testing.assert_allclose(
pointings_gal_hp, pointings_gal_lbs, rtol=1e-6, atol=1e-6
Expand Down
12 changes: 6 additions & 6 deletions test/test_pointing_sys.py
Original file line number Diff line number Diff line change
Expand Up @@ -112,7 +112,7 @@ def test_PointingSys_add_single_offset_to_FP(
pointings, hwp_angle = cur_obs.get_pointings(
det_idx, pointings_dtype=np.float32
)
pointings_list.append(pointings[0].tolist())
pointings_list.append(pointings.tolist())

if func.__name__ not in results_dict:
results_dict[func.__name__] = {}
Expand Down Expand Up @@ -150,7 +150,7 @@ def test_PointingSys_add_multiple_offsets_to_FP(
pointings, hwp_angle = cur_obs.get_pointings(
det_idx, pointings_dtype=np.float32
)
pointings_list.append(pointings[0].tolist())
pointings_list.append(pointings.tolist())

if func.__name__ not in results_dict:
results_dict[func.__name__] = {}
Expand Down Expand Up @@ -196,7 +196,7 @@ def test_PointingSys_add_uncommon_disturb_to_FP(
pointings, hwp_angle = cur_obs.get_pointings(
det_idx, pointings_dtype=np.float32
)
pointings_list.append(pointings[0].tolist())
pointings_list.append(pointings.tolist())

if func.__name__ not in results_dict:
results_dict[func.__name__] = {}
Expand Down Expand Up @@ -241,7 +241,7 @@ def test_PointingSys_add_common_disturb_to_FP(
pointings, hwp_angle = cur_obs.get_pointings(
det_idx, pointings_dtype=np.float32
)
pointings_list.append(pointings[0].tolist())
pointings_list.append(pointings.tolist())

if func.__name__ not in results_dict:
results_dict[func.__name__] = {}
Expand Down Expand Up @@ -279,7 +279,7 @@ def test_PointingSys_add_single_offset_to_spacecraft(
pointings, hwp_angle = cur_obs.get_pointings(
det_idx, pointings_dtype=np.float32
)
pointings_list.append(pointings[0].tolist())
pointings_list.append(pointings.tolist())

if func.__name__ not in results_dict:
results_dict[func.__name__] = {}
Expand Down Expand Up @@ -323,7 +323,7 @@ def test_PointingSys_add_common_disturb_to_spacecraft(
pointings, hwp_angle = cur_obs.get_pointings(
det_idx, pointings_dtype=np.float32
)
pointings_list.append(pointings[0].tolist())
pointings_list.append(pointings.tolist())

if func.__name__ not in results_dict:
results_dict[func.__name__] = {}
Expand Down
12 changes: 6 additions & 6 deletions test/test_simulations.py
Original file line number Diff line number Diff line change
Expand Up @@ -465,7 +465,7 @@ def test_smart_pointings_consistency_with_hwp(tmp_path):

for det_idx in range(obs.n_detectors):
(pointings, hwp_angle) = obs.get_pointings(det_idx)
assert pointings.shape == (1, obs.n_samples, 3)
assert pointings.shape == (obs.n_samples, 3)
assert hwp_angle.shape == (obs.n_samples,)


Expand All @@ -479,7 +479,7 @@ def test_smart_pointings_consistency_without_hwp(tmp_path):

for det_idx in range(obs.n_detectors):
(pointings, hwp_angle) = obs.get_pointings(det_idx)
assert pointings.shape == (1, obs.n_samples, 3)
assert pointings.shape == (obs.n_samples, 3)
assert hwp_angle is None


Expand All @@ -497,7 +497,7 @@ def test_smart_pointings_angles(tmp_path):
# function present in LBS 0.12.0

np.testing.assert_allclose(
actual=pointings[0, 0:10, :],
actual=pointings[0:10, :],
desired=np.array(
[
[1.6580627894, 0.0000000000, 1.3977241805],
Expand Down Expand Up @@ -617,7 +617,7 @@ def test_compute_pointings_for_one_detector(dtype, tmp_path):
pointings, hwp_angle = cur_obs.get_pointings(0, pointings_dtype=dtype)

assert pointings.dtype == dtype
assert pointings.shape == (1, cur_obs.n_samples, 3)
assert pointings.shape == (cur_obs.n_samples, 3)

assert hwp_angle.dtype == dtype
assert hwp_angle.shape == (cur_obs.n_samples,)
Expand Down Expand Up @@ -648,7 +648,7 @@ def test_store_pointings_for_two_detectors(dtype, tmp_path):
cur_pointings, _ = cur_obs.get_pointings(abs_det_idx)
np.testing.assert_allclose(
actual=pointings[rel_det_idx, :, :],
desired=cur_pointings[0, :, :],
desired=cur_pointings[:, :],
rtol=1e-6,
)

Expand Down Expand Up @@ -677,6 +677,6 @@ def test_smart_pointings_for_all_detectors(dtype, tmp_path):
assert cur_pointings.dtype == dtype
np.testing.assert_allclose(
actual=pointings[det_idx, :, :],
desired=cur_pointings[0, :, :],
desired=cur_pointings[:, :],
rtol=1e-6,
)

0 comments on commit c6c6021

Please sign in to comment.