diff --git a/sotodlib/coords/helpers.py b/sotodlib/coords/helpers.py index 231e390db..022eac3f3 100644 --- a/sotodlib/coords/helpers.py +++ b/sotodlib/coords/helpers.py @@ -739,6 +739,34 @@ def to_g3(self): return so3g.proj.quat.G3VectorQuat(temp) raise ValueError("Can only convert 1- or 2-d arrays to G3.") +def get_deflection_quat(aman, wobble_meta=None): + ''' + Derives the deflection quaternion. + + Parameters + ---------- + aman : AxisManager + + wobble_meta : AxisManager + + ''' + if wobble_meta is None: + wobble_meta = aman.get('wobble_params') + elif isinstance(wobble_meta, str): + wobble_meta = aman.get(wobble_meta) + # the amp and phase are the same for a given wafer, so we can take any of them, in this case for detector index 0 + # !!!!! this won't work for mixing more than one wafer. + # the metadata has amplitudes in arcmin, and phases in radians + amp = wobble_meta.amp[0]/60.*np.pi/180.0 + phase = wobble_meta.phase[0] + + dxi = amp * np.cos(aman.hwp_angle - phase) + deta = -amp * np.sin(aman.hwp_angle - phase) + deflq = so3g.proj.quat.rotation_xieta(xi=dxi, eta=deta) + + return deflq + + def get_deflected_sightline(aman, wobble_meta=None, sight=None, site='so', weather='typical'): """ Constructs a deflected CelestialSightLine using HWP-synchronous @@ -781,19 +809,8 @@ def get_deflected_sightline(aman, wobble_meta=None, sight=None, site='so', weath if len(wafer_slots) != 1 or len(bands) != 1: raise ValueError("Detectors span multiple wafer_slots or bands.") - if wobble_meta is None: - wobble_meta = aman.get('wobble_params') - elif isinstance(wobble_meta, str): - wobble_meta = aman.get(wobble_meta) - # the amp and phase are the same for a given wafer, so we can take any of them, in this case for detector index 0 - # !!!!! this won't work for mixing more than one wafer. - # the metadata has amplitudes in arcmin, and phases in radians - amp = wobble_meta.amp[0]/60.*np.pi/180.0 - phase = wobble_meta.phase[0] - - dxi = amp * np.cos(aman.hwp_angle - phase) - deta = -amp * np.sin(aman.hwp_angle - phase) - deflq = so3g.proj.quat.rotation_xieta(xi=dxi, eta=deta) + + deflq = get_deflection_quat(aman, wobble_meta) if sight is None: sight = so3g.proj.CelestialSightLine.az_el( diff --git a/sotodlib/coords/planets.py b/sotodlib/coords/planets.py index 69fa1c92c..7541f7029 100644 --- a/sotodlib/coords/planets.py +++ b/sotodlib/coords/planets.py @@ -773,7 +773,8 @@ def make_map(tod, center_on=None, scan_coords=True, thread_algo=False, filename=None, source_flags=None, cuts=None, data_splits=None, low_pass=None, n_modes=10, - eigentol=1e-3, info={}): + eigentol=1e-3, info={}, + rot_bs=None, rot_fp=None): """Make a compact source map from the TOD. Specify filename to write things to disk; this should be a format string, for example '{obs_id}_{map}.fits', where 'map' will be given values of @@ -815,6 +816,14 @@ class MmTimer(coords.Timer): cuts=cuts, threads=thread_algo, wcs_kernel=wcsk) + + # apply an optional boresight rotation or focal plane rotation + # e.g. HWP deflection is rot_fp + if rot_bs is not None: + P.sight.Q = rot_bs * P.sight.Q + if rot_fp is not None: + P.sight.Q = P.sight.Q * rot_fp + with MmTimer('get_proj_threads'): P._get_proj_threads()