12832805671

Committed 17 Jan 2025 04:16PM UTC coverage: 50.45% (+0.4%) from 50.055%

Build # 12832805671

Build Type

push

github

Committed by

web-flow

Commit Message

Stokes response (#193)

* Plans for arbitrary detector response stuff

* Detector response now supported in c++, and represented in a redone FocalPlane class in python

* Arbitrary response support

* Passes tests after adding .items(), updating xieta arguments and replacing asm.dets with asm.fplane.quats

* Temporary compatibility with old sotodlib

* Default FocalPlane is not empty, so make a separate constructor for making one with just the boresight

* Better docstrings and comments

* Working on new documentation

* More documentation. Updated numbers in examples which seem to have rotted at some point

* Matthews comments, plus go to G3VectorQuat instead of numpy array

* Ready for new review

* Fixed last-minute typo

* Addressed last round of commments

Run Details

69 of 87 new or added lines in 3 files covered. (79.31%)

1 existing line in 1 file now uncovered.

1401 of 2777 relevant lines covered (50.45%)

0.5 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

76.47

/python/proj/coords.py

import so3g
from . import quat
from .weather import weather_factory

from collections import OrderedDict

import numpy as np

DEG = np.pi / 180.


class EarthlySite:
    def __init__(self, lon, lat, elev, typical_weather=None):
        """Arguments in degrees E, degrees N, meters above sea level.  The
        typical_weather argument should be an so3g Weather object, or
        None.

        """
        self.lon, self.lat, self.elev = lon, lat, elev
        self.typical_weather = typical_weather

    def ephem_observer(self):
        """Return an ephem.Observer corresponding to this location."""
        import ephem
        site = ephem.Observer()
        site.lat =  str(self.lat)
        site.long = str(self.lon)
        site.elev = self.elev
        return site

    def skyfield_site(self, spice_kernel):
        """Return a skyfield VectorSum for this location on 'earth' (starting
        from wgs84).  The spice_kernel is the thing one might get from
        jpllib.SpiceKernel(emphemeris_filename).

        """
        from skyfield.api import wgs84
        earth = spice_kernel['earth']
        return earth + wgs84.latlon(self.lat, self.lon, self.elev)


def _debabyl(deg, arcmin, arcsec):
    return deg + arcmin/60 + arcsec/3600

SITES = {
    'act': EarthlySite(-67.7876, -22.9585, 5188.,
                       typical_weather=weather_factory('toco')),
    # SO coords taken from SO-SITE-HEF-003A on 2020-06-02; altitude
    # set to same as ACT.
    'so_lat': EarthlySite(-_debabyl(67,47,15.68), -_debabyl(22,57,39.47), 5188.,
                          typical_weather=weather_factory('toco')),
    'so_sat1': EarthlySite(-_debabyl(67,47,18.11), -_debabyl(22,57,36.38), 5188.,
                           typical_weather=weather_factory('toco')),
    'so_sat2': EarthlySite(-_debabyl(67,47,17.28), -_debabyl(22,57,36.35), 5188.,
                           typical_weather=weather_factory('toco')),
    'so_sat3': EarthlySite(-_debabyl(67,47,16.53), -_debabyl(22,57,35.97), 5188.,
                           typical_weather=weather_factory('toco')),
}
# Take LAT as default.  It makes most sense to use a single position
# for all telescopes and absorb the discrepancy into the pointing
# model as a few seconds of base tilt.
SITES['so'] = SITES['so_lat']
SITES['_default'] = SITES['so']
DEFAULT_SITE = 'so'  # deprecated, use SITES['_default']

# These definitions are used for the naive horizon -> celestial
# conversion.
ERA_EPOCH = 946684800 + 3600 * 12  # Noon on Jan 1 2000.
ERA_POLY = np.array([6.300387486754831,
                     4.894961212823756])  # Operates on "days since ERA_EPOCH".


class CelestialSightLine:
    """Carries a vector of celestial pointing data.

    Instantiation should result in the pointing quaternion being
    stored in self.Q.

    """
    @staticmethod
    def decode_site(site=None):
        """Convert site argument to a Site object.  The argument must be one
        of:

        - an EarthlySite object
        - a string, corresponding to the name of an internally known
          site
        - None, in which casethe default site info will be loaded.

        """
        if site is None:
            site = DEFAULT_SITE
        if isinstance(site, EarthlySite):
            return site
        if site in SITES:
            return SITES[site]
        raise ValueError("Could not decode %s as a Site." % site)

    @classmethod
    def naive_az_el(cls, t, az, el, roll=0., site=None, weather=None):
        """Construct a SightLine from horizon coordinates, az and el (in
        radians) and time t (unix timestamp).

        This will be off by several arcminutes... but less than a
        degree.  The weather is ignored.

        """
        site = cls.decode_site(site)
        assert isinstance(site, EarthlySite)

        self = cls()

        J = (t - ERA_EPOCH) / 86400
        era = np.polyval(ERA_POLY, J)
        lst = era + site.lon * DEG

        self.Q = (
            quat.euler(2, lst) *
            quat.euler(1, np.pi/2 - site.lat * DEG) *
            quat.euler(2, np.pi) *
            quat.euler(2, -az) *
            quat.euler(1, np.pi/2 - el) *
            quat.euler(2, roll)
            )
        return self

    @classmethod
    def az_el(cls, t, az, el, roll=None, site=None, weather=None, **kwargs):
        """Construct a SightLine from horizon coordinates.  This uses
        high-precision pointing. kwargs are passed to qpoint.Qpoint.

        """
        import qpoint  # https://github.com/arahlin/qpoint

        site = cls.decode_site(site)
        assert isinstance(site, EarthlySite)

        if isinstance(weather, str):
            if weather == 'typical':
                weather = site.typical_weather
            else:
                weather = weather_factory(weather)

        if weather is None:
            raise ValueError('High-precision pointing requires a weather '
                             'object.  Try passing \'toco\', \'typical\', or '
                             '\'vacuum\'.')

        self = cls()
        # Note that passing num_threads explicitly to qpoint will
        # cause openmp_set_thread_count to be called!
        qp = qpoint.QPoint(accuracy='high', fast_math=True, mean_aber=True,
                           rate_ref='always', **weather.to_qpoint(), **kwargs)

        az, el, t = map(np.asarray, [az, el, t])
        Q_arr = qp.azel2bore(az / DEG, el / DEG, None, None,
                             lon=site.lon, lat=site.lat, ctime=t)
        self.Q = quat.G3VectorQuat(Q_arr)

        # Apply boresight roll manually (note qpoint pitch/roll refer
        # to the platform rather than to something about the boresight
        # axis).  Regardless we need a pi rotation because of our
        # definition of boresight coordinates.
        if roll is None:
            self.Q *= quat.euler(2, np.pi)
        if roll is not None:
            self.Q *= quat.euler(2, np.pi + roll)

        return self

    @classmethod
    def for_lonlat(cls, lon, lat, psi=0):
        """Construct a SightLine directly from lonlat angles representing
        celestial coordinates.

        This is appropriate if you want a SightLine that takes focal
        plane coordinates directly to some celestial pointing.  In
        that case, lon and lat are RA and dec, and psi is the
        parallactic rotation of the focal plane on the sky at the
        target position.  psi=0 will correspond to focal plane "up"
        (+eta axis) mapping parallel to lines of longitude; psi > 0 is
        a clockwise rotation of the focal plane on the sky
        (i.e. opposite IAU Position Angle direction).

        """
        self = cls()
        self.Q = quat.euler(2, lon) * quat.euler(1, np.pi/2 - lat) * quat.euler(2, psi)
        return self

    @classmethod
    def for_horizon(cls, t, az, el, roll=None, site=None, weather=None):
        """Construct the trivial SightLine, where "celestial" coordinates are
        taken to simply be local horizon coordinates (without any
        effects of atmosphere).  Although this accepts site and
        weather arguments, they are ignored.

        """
        self = cls()
        az, el, t = map(np.asarray, [az, el, t])
        if roll is None:
            roll = 0.
        else:
            roll = np.asarray(roll)
        self.Q = (
            quat.euler(2, -az) *
            quat.euler(1, np.pi/2 - el) *
            quat.euler(2, roll)
            )
        return self

    def coords(self, fplane=None, output=None):
        """Get the celestial coordinates of each detector at each time.

        Arguments:
          fplane: A FocalPlane object representing the detector
            offsets and responses, or None
        output: An optional structure for holding the results.  For
            that to work, each element of output must support the
            buffer protocol.

        Returns:
          If fplane is None, then the result will be
          [n_samp,{lon,lat,cos2psi,sin2psi}]. Otherwise it will
          be [n_det,n_samp,{lon,lat,cos2psi,sin2psi}]
        """
        # Get a projector, in CAR.
        p = so3g.ProjEng_CAR_TQU_NonTiled((1, 1, 1., 1., 1., 1.))
        # Pre-process the offsets
        collapse = (fplane is None)
        if collapse:
            fplane = FocalPlane.boresight()
            if output is not None:
                output = output[None]
        output = p.coords(self.Q, fplane.quats, output)
        if collapse:
            output = output[0]
        return output

class FocalPlane:
    """This class represents the detector positions and intensity and
    polarization responses in the focal plane.

    Attributes:
     quats: G3VectorQuat representing the pointing quaternions for
      each detector. Can be turned into a numpy array of coefficients
      with np.array(). Or call .coeffs() to get them directly.
     resps: Array of float32 with shape [ndet,2] representing the
      total intensity and polarization responsivity of each detector
     ndet: The number of detectors (read-only)

    (This used to be a subclass of OrderedDict, but this was hard to
    generalize to per-detector polarization efficiency.)
    """
    # FIXME: old sotodlib compat - remove dets argument later
    def __init__(self, quats=None, resps=None, dets=None):
        """Construct a FocalPlane from detector quaternions and responsivities.

        Arguments:
         quats: Detector quaternions. Either:
           * An array-like of floats with shape [ndet,4]
           * An array-like of so3g.proj.quat.quat with shape [ndet]
           * An so3g.proj.quat.G3VectorQuat
           * None, which results in an empty focalplane with no detectors
         resps: Detector responsivities. Either:
           * An array-like of floats with shape [ndet,2], where the first and
             second number in the last axis are the total intensity and
             polarization response respectively
           * None, which results in a T and P response of 1 for all detectors.
         dets: Deprecated argument temporarily present for backwards
           compatibility."""
        # Building them this way ensures that
        # quats will be an quat coeff array-2 and resps will be a numpy
        # array with the right shape, so we don't need to check
        # for this when we use FocalPlane later
        if quats is None: quats = []
        # Asarray needed because G3VectorQuat doesn't handle list of lists,
        # which we want to be able to accept
        self.quats = quat.G3VectorQuat(np.asarray(quats))
        self.resps = np.ones((len(self.quats),2),np.float32)
        if resps is not None:
            self.resps[:] = resps
        if np.any(~np.isfinite(self.quats)):
            raise ValueError("nan/inf values in detector quaternions")
        if np.any(~np.isfinite(self.resps)):
            raise ValueError("nan/inf values in detector responses")
        # FIXME: old sotodlib compat - remove later
        self._dets   = list(dets) if dets is not None else []
        self._detmap = {name:i for i,name in enumerate(self._dets)}
    def coeffs(self):
        """Return an [ndet,4] numpy array representing the quaternion
        coefficients of ``.quats``. Useful for passing the detector
        pointing to C++ code"""
        return np.array(self.quats, dtype=np.float64)
    @classmethod
    def boresight(cls):
        """Construct a FocalPlane representing a single detector with a
        responsivity of one pointing along the telescope boresight"""
        return cls(quats=[[1,0,0,0]])
    # FIXME: old sotodlib compat - expand to actual argument list later
    @classmethod
    def from_xieta(cls, *args, **kwargs):
        """
        ``from_xieta(cls, xi, eta, gamma=0, T=1, P=1, Q=1, U=0, hwp=False)``
        Construct a FocalPlane from focal plane coordinates (xi,eta).

        For backwards compatibility, the signature
        ``from_xieta(names, xi, eta, gamma=0)`` is also supported; but
        this will be removed in the future.

        These are Cartesian projection plane coordinates. When looking
        at the sky along the telescope boresight, xi is parallel to
        increasing azimuth and eta is parallel to increasing elevation.
        The angle gamma, which specifies the angle of polarization
        sensitivity, is measured from the eta axis, increasing towards
        the xi axis.

        Arguments:
         xi:  An array-like of floats with shape [ndet]
         eta: An array-like of floats with shape [ndet]
         gamma: The detector polarization angles. [ndet], or a scalar
         T, P: The total intensity and polarization responsivity.
            Array-like with shape [ndet], or scalar.
         Q, U: The Q- and U-polarization responsivity.
            Array-like with shape [ndet], or scalar.
            Q,U are alternatives to P,gamma for specifying the
            polarization responsivity and angle.

        So there are two ways to specify the polarization angle and
        responsivity:
          1. gamma and P
          2. Q and U

        Examples, assuming ndet = 2
         * ``from_xieta(xi, eta, gamma=[0,pi/4])``
           Constructs a FocalPlane with T and P responsivity of 1
           and polarization angles of 0 and 45 degrees, representing
           a Q-sensitive and U-sensitive detector.
         * ``from_xieta(xi, eta, gamma=[0,pi/4], P=0.5)``
           Like the above, but with a polarization responsivity of
           just 0.5.
         * ``from_xieta(xi, eta, gamma=[0,pi/4], T=[1,0.9], P=[0.5,0.6])``
           Like above, but with a detector-dependent intensity and
           polarization responsivity. There is no restriction that
           T > P. For the pseudo-detector timestreams one gets after
           HWP demodulation, one would have T=0 for the cos-modulated
           and sin-modulated timestreams, for example.
         * ``from_xieta(xi, eta, Q=[1,0], U=[0,1])``
           Construct the FocalPlane with explicit Q and U responsivity.
           This example is equivalent to example 1.

        Usually one would either use gamma,P or Q,U. If they are
        combined, then ``gamma_total = gamma + arctan2(U,Q)/2`` and
        ``P_tot = P * (Q**2+U**2)**0.5``.
        """
        # The underlying code wants polangle gamma and the T and P
        # response, but we support speifying these as the T, Q and U
        # response too. Writing it like this handles both cases, and
        # as a bonus(?) any combination of them
        # FIXME: old sotodlib compat - remove later
        xi, eta, gamma, T, P, Q, U, hwp, dets = cls._xieta_compat(*args, **kwargs)
        gamma = gamma + np.arctan2(U,Q)/2
        P     = P * (Q**2+U**2)**0.5
        if hwp: gamma = -gamma
        # Broadcast everything to 1d
        xi, eta, gamma, T, P, _ = np.broadcast_arrays(xi, eta, gamma, T, P, [0])
        quats = quat.rotation_xieta(xi, eta, gamma)
        resps = np.ones((len(quats),2))
        resps[:,0] = T
        resps[:,1] = P
        # FIXME: old sotodlib compat - remove dets argument later
        return cls(quats, resps=resps, dets=dets)
    def __repr__(self):
        return "FocalPlane(quats=%s,resps=%s)" % (repr(self.coeffs()), repr(self.resps))
    @property
    def ndet(self): return len(self.quats)
    def __len__(self): return self.ndet
    def __getitem__(self, sel):
        """Slice the FocalPlane with slice sel, resulting in a new
        FocalPlane with a subset of the detectors. Only a 1d slice
        or boolean mask is supported, not integers or multidimensional
        slices. Example: ``focal_plane[10:20]`` would make a
        sub-FocalPlane with detector indices 10,11,...,19.

        Deprecated: Temporarily also supports that sel is a detector name,
        in which case an  ``spt3g.core.quat`` is returned for that detector.
        This is provided for backwards compatibility."""
        # FIXME: old sotodlib compat - remove later
        if isinstance(sel, str): return self.quats[self._dets.index(sel)]
        # We go via .coeffs() here to get around G3VectorQuat's lack
        # of boolean mask support
        return FocalPlane(quats=self.coeffs()[sel], resps=self.resps[sel])
    def items(self):
        """Iterate over detector quaternions and responsivities. Yields
        ``(spt3g.core.quat, array([Tresp,Presp]))`` pairs. Unlke the raw
        entries in the .quats member, which are just numpy arrays with
        length 4,  ``spt3g.core.quat`` are proper quaternon objects that
        support quaternion math.
        """
        for q, resp in zip(self.quats, self.resps):
            yield q, resp
    # FIXME: old sotodlib compat - remove later
    @staticmethod
    def _xieta_compat(*args, **kwargs):
        # Accept the alternative format (names,xi,eta,gamma)
        def helper(xi, eta, gamma=0, T=1, P=1, Q=1, U=0, hwp=False, dets=None):
            return xi, eta, gamma, T, P, Q, U, hwp, dets
        if isinstance(args[0][0], str):
            return helper(*args[1:], dets=args[0], **kwargs)
        else:
            return helper(*args, **kwargs)
    # FIXME: old sotodlib compat - remove later
    def __setitem__(self, name, q):
        # Make coords/pmat.py _get_asm work in old sotodlib. It
        # expects to be able to build up a focalplane by assigning
        # quats one at a time
        if name in self._detmap:
            i = self._detmap[i]
            self.quats[i] = q
        else:
            self._dets.append(name)
            self._detmap[name] = len(self._dets)-1
            self.quats.append(q)
            # This is inefficient, but it's temporary
            self.resps = np.concatenate([self.resps,np.ones((1,2),np.float32)])

class Assembly:
    """This class groups together boresight and detector offset
    information for use by the Projectionist.

    The basic implementation simply attachs some number of detectors
    rigidly to a boresight.  But this abstraction layer could
    facilitate more complex arrangements, eventually.

    """
    def __init__(self, collapse=False):
        self.collapse = collapse

    @classmethod
    def attach(cls, sight_line, fplane):
        """Create an Assembly based on a CelestialSightLine and a FocalPlane.

        Args:
            sight_line (CelestialSightLine): The position and
                orientation of the boresight.  This can just be a
                G3VectorQuat if you don't have a whole
                CelestialSightLine handy.
            fplane (FocalPlane): The "offsets" of each detector
                relative to the boresight, and their response to
                intensity and polarization
        """
        self = cls()
        if isinstance(sight_line, quat.G3VectorQuat):
            self.Q = sight_line
        else:
            self.Q = sight_line.Q
        self.fplane = fplane
        return self

    @classmethod
    def for_boresight(cls, sight_line):
        """Returns an Assembly where a single detector serves as a dummy for
        the boresight."""
        self = cls(collapse=True)
        self.Q = sight_line.Q
        self.fplane = FocalPlane.boresight()
        return self
    # FIXME: old sotodlib compat - remove later
    @property
    def dets(self): return self.fplane

1	import so3g	1✔
2	from . import quat	1✔
3	from .weather import weather_factory	1✔
4
5	from collections import OrderedDict	1✔
6
7	import numpy as np	1✔
8
9	DEG = np.pi / 180.	1✔
10
11
12	class EarthlySite:	1✔
13	def __init__(self, lon, lat, elev, typical_weather=None):	1✔
14	"""Arguments in degrees E, degrees N, meters above sea level. The
15	typical_weather argument should be an so3g Weather object, or
16	None.
17
18	"""
19	self.lon, self.lat, self.elev = lon, lat, elev	1✔
20	self.typical_weather = typical_weather	1✔
21
22	def ephem_observer(self):	1✔
23	"""Return an ephem.Observer corresponding to this location."""
24	import ephem	×
25	site = ephem.Observer()	×
26	site.lat = str(self.lat)	×
27	site.long = str(self.lon)	×
28	site.elev = self.elev	×
29	return site	×
30
31	def skyfield_site(self, spice_kernel):	1✔
32	"""Return a skyfield VectorSum for this location on 'earth' (starting
33	from wgs84). The spice_kernel is the thing one might get from
34	jpllib.SpiceKernel(emphemeris_filename).
35
36	"""
37	from skyfield.api import wgs84	×
38	earth = spice_kernel['earth']	×
39	return earth + wgs84.latlon(self.lat, self.lon, self.elev)	×
40
41
42	def _debabyl(deg, arcmin, arcsec):	1✔
43	return deg + arcmin/60 + arcsec/3600	1✔
44
45	SITES = {	1✔
46	'act': EarthlySite(-67.7876, -22.9585, 5188.,
47	typical_weather=weather_factory('toco')),
48	# SO coords taken from SO-SITE-HEF-003A on 2020-06-02; altitude
49	# set to same as ACT.
50	'so_lat': EarthlySite(-_debabyl(67,47,15.68), -_debabyl(22,57,39.47), 5188.,
51	typical_weather=weather_factory('toco')),
52	'so_sat1': EarthlySite(-_debabyl(67,47,18.11), -_debabyl(22,57,36.38), 5188.,
53	typical_weather=weather_factory('toco')),
54	'so_sat2': EarthlySite(-_debabyl(67,47,17.28), -_debabyl(22,57,36.35), 5188.,
55	typical_weather=weather_factory('toco')),
56	'so_sat3': EarthlySite(-_debabyl(67,47,16.53), -_debabyl(22,57,35.97), 5188.,
57	typical_weather=weather_factory('toco')),
58	}
59	# Take LAT as default. It makes most sense to use a single position
60	# for all telescopes and absorb the discrepancy into the pointing
61	# model as a few seconds of base tilt.
62	SITES['so'] = SITES['so_lat']	1✔
63	SITES['_default'] = SITES['so']	1✔
64	DEFAULT_SITE = 'so' # deprecated, use SITES['_default']	1✔
65
66	# These definitions are used for the naive horizon -> celestial
67	# conversion.
68	ERA_EPOCH = 946684800 + 3600 * 12 # Noon on Jan 1 2000.	1✔
69	ERA_POLY = np.array([6.300387486754831,	1✔
70	4.894961212823756]) # Operates on "days since ERA_EPOCH".
71
72
73	class CelestialSightLine:	1✔
74	"""Carries a vector of celestial pointing data.
75
76	Instantiation should result in the pointing quaternion being
77	stored in self.Q.
78
79	"""
80	@staticmethod	1✔
81	def decode_site(site=None):	1✔
82	"""Convert site argument to a Site object. The argument must be one
83	of:
84
85	- an EarthlySite object
86	- a string, corresponding to the name of an internally known
87	site
88	- None, in which casethe default site info will be loaded.
89
90	"""
91	if site is None:	1✔
92	site = DEFAULT_SITE	×
93	if isinstance(site, EarthlySite):	1✔
94	return site	1✔
95	if site in SITES:	1✔
96	return SITES[site]	1✔
97	raise ValueError("Could not decode %s as a Site." % site)	×
98
99	@classmethod	1✔
100	def naive_az_el(cls, t, az, el, roll=0., site=None, weather=None):	1✔
101	"""Construct a SightLine from horizon coordinates, az and el (in
102	radians) and time t (unix timestamp).
103
104	This will be off by several arcminutes... but less than a
105	degree. The weather is ignored.
106
107	"""
108	site = cls.decode_site(site)	1✔
109	assert isinstance(site, EarthlySite)	1✔
110
111	self = cls()	1✔
112
113	J = (t - ERA_EPOCH) / 86400	1✔
114	era = np.polyval(ERA_POLY, J)	1✔
115	lst = era + site.lon * DEG	1✔
116
117	self.Q = (	1✔
118	quat.euler(2, lst) *
119	quat.euler(1, np.pi/2 - site.lat * DEG) *
120	quat.euler(2, np.pi) *
121	quat.euler(2, -az) *
122	quat.euler(1, np.pi/2 - el) *
123	quat.euler(2, roll)
124	)
125	return self	1✔
126
127	@classmethod	1✔
128	def az_el(cls, t, az, el, roll=None, site=None, weather=None, **kwargs):	1✔
129	"""Construct a SightLine from horizon coordinates. This uses
130	high-precision pointing. kwargs are passed to qpoint.Qpoint.
131
132	"""
133	import qpoint # https://github.com/arahlin/qpoint	1✔
134
135	site = cls.decode_site(site)	1✔
136	assert isinstance(site, EarthlySite)	1✔
137
138	if isinstance(weather, str):	1✔
139	if weather == 'typical':	1✔
140	weather = site.typical_weather	×
141	else:
142	weather = weather_factory(weather)	1✔
143
144	if weather is None:	1✔
145	raise ValueError('High-precision pointing requires a weather '	×
146	'object. Try passing \'toco\', \'typical\', or '
147	'\'vacuum\'.')
148
149	self = cls()	1✔
150	# Note that passing num_threads explicitly to qpoint will
151	# cause openmp_set_thread_count to be called!
152	qp = qpoint.QPoint(accuracy='high', fast_math=True, mean_aber=True,	1✔
153	rate_ref='always', weather.to_qpoint(), kwargs)
154
155	az, el, t = map(np.asarray, [az, el, t])	1✔
156	Q_arr = qp.azel2bore(az / DEG, el / DEG, None, None,	1✔
157	lon=site.lon, lat=site.lat, ctime=t)
158	self.Q = quat.G3VectorQuat(Q_arr)	1✔
159
160	# Apply boresight roll manually (note qpoint pitch/roll refer
161	# to the platform rather than to something about the boresight
162	# axis). Regardless we need a pi rotation because of our
163	# definition of boresight coordinates.
164	if roll is None:	1✔
165	self.Q *= quat.euler(2, np.pi)	1✔
166	if roll is not None:	1✔
167	self.Q *= quat.euler(2, np.pi + roll)	×
168
169	return self	1✔
170
171	@classmethod	1✔
172	def for_lonlat(cls, lon, lat, psi=0):	1✔
173	"""Construct a SightLine directly from lonlat angles representing
174	celestial coordinates.
175
176	This is appropriate if you want a SightLine that takes focal
177	plane coordinates directly to some celestial pointing. In
178	that case, lon and lat are RA and dec, and psi is the
179	parallactic rotation of the focal plane on the sky at the
180	target position. psi=0 will correspond to focal plane "up"
181	(+eta axis) mapping parallel to lines of longitude; psi > 0 is
182	a clockwise rotation of the focal plane on the sky
183	(i.e. opposite IAU Position Angle direction).
184
185	"""
186	self = cls()	×
187	self.Q = quat.euler(2, lon) * quat.euler(1, np.pi/2 - lat) * quat.euler(2, psi)	×
188	return self	×
189
190	@classmethod	1✔
191	def for_horizon(cls, t, az, el, roll=None, site=None, weather=None):	1✔
192	"""Construct the trivial SightLine, where "celestial" coordinates are
193	taken to simply be local horizon coordinates (without any
194	effects of atmosphere). Although this accepts site and
195	weather arguments, they are ignored.
196
197	"""
198	self = cls()	×
199	az, el, t = map(np.asarray, [az, el, t])	×
200	if roll is None:	×
201	roll = 0.	×
202	else:
203	roll = np.asarray(roll)	×
204	self.Q = (	×
205	quat.euler(2, -az) *
206	quat.euler(1, np.pi/2 - el) *
207	quat.euler(2, roll)
208	)
209	return self	×
210
211	def coords(self, fplane=None, output=None):	1✔
212	"""Get the celestial coordinates of each detector at each time.
213
214	Arguments:
215	fplane: A FocalPlane object representing the detector
216	offsets and responses, or None
217	output: An optional structure for holding the results. For
218	that to work, each element of output must support the
219	buffer protocol.
220
221	Returns:
222	If fplane is None, then the result will be
223	[n_samp,{lon,lat,cos2psi,sin2psi}]. Otherwise it will
224	be [n_det,n_samp,{lon,lat,cos2psi,sin2psi}]
225	"""
226	# Get a projector, in CAR.
227	p = so3g.ProjEng_CAR_TQU_NonTiled((1, 1, 1., 1., 1., 1.))	1✔
228	# Pre-process the offsets
229	collapse = (fplane is None)	1✔
230	if collapse:	1✔
231	fplane = FocalPlane.boresight()	1✔
232	if output is not None:	1✔
233	output = output[None]	×
234	output = p.coords(self.Q, fplane.quats, output)	1✔
235	if collapse:	1✔
236	output = output[0]	1✔
237	return output	1✔
238
239	class FocalPlane:	1✔
240	"""This class represents the detector positions and intensity and
241	polarization responses in the focal plane.
242
243	Attributes:
244	quats: G3VectorQuat representing the pointing quaternions for
245	each detector. Can be turned into a numpy array of coefficients
246	with np.array(). Or call .coeffs() to get them directly.
247	resps: Array of float32 with shape [ndet,2] representing the
248	total intensity and polarization responsivity of each detector
249	ndet: The number of detectors (read-only)
250
251	(This used to be a subclass of OrderedDict, but this was hard to
252	generalize to per-detector polarization efficiency.)
253	"""
254	# FIXME: old sotodlib compat - remove dets argument later
255	def __init__(self, quats=None, resps=None, dets=None):	1✔
256	"""Construct a FocalPlane from detector quaternions and responsivities.
257
258	Arguments:
259	quats: Detector quaternions. Either:
260	* An array-like of floats with shape [ndet,4]
261	* An array-like of so3g.proj.quat.quat with shape [ndet]
262	* An so3g.proj.quat.G3VectorQuat
263	* None, which results in an empty focalplane with no detectors
264	resps: Detector responsivities. Either:
265	* An array-like of floats with shape [ndet,2], where the first and
266	second number in the last axis are the total intensity and
267	polarization response respectively
268	* None, which results in a T and P response of 1 for all detectors.
269	dets: Deprecated argument temporarily present for backwards
270	compatibility."""
271	# Building them this way ensures that
272	# quats will be an quat coeff array-2 and resps will be a numpy
273	# array with the right shape, so we don't need to check
274	# for this when we use FocalPlane later
275	if quats is None: quats = []	1✔
276	# Asarray needed because G3VectorQuat doesn't handle list of lists,
277	# which we want to be able to accept
278	self.quats = quat.G3VectorQuat(np.asarray(quats))	1✔
279	self.resps = np.ones((len(self.quats),2),np.float32)	1✔
280	if resps is not None:	1✔
281	self.resps[:] = resps	1✔
282	if np.any(~np.isfinite(self.quats)):	1✔
NEW 283	raise ValueError("nan/inf values in detector quaternions")	×
284	if np.any(~np.isfinite(self.resps)):	1✔
NEW 285	raise ValueError("nan/inf values in detector responses")	×
286	# FIXME: old sotodlib compat - remove later
287	self._dets = list(dets) if dets is not None else []	1✔
288	self._detmap = {name:i for i,name in enumerate(self._dets)}	1✔
289	def coeffs(self):	1✔
290	"""Return an [ndet,4] numpy array representing the quaternion
291	coefficients of ``.quats``. Useful for passing the detector
292	pointing to C++ code"""
293	return np.array(self.quats, dtype=np.float64)	1✔
294	@classmethod	1✔
295	def boresight(cls):	1✔
296	"""Construct a FocalPlane representing a single detector with a
297	responsivity of one pointing along the telescope boresight"""
298	return cls(quats=[[1,0,0,0]])	1✔
299	# FIXME: old sotodlib compat - expand to actual argument list later
300	@classmethod	1✔
301	def from_xieta(cls, args, *kwargs):	1✔
302	"""
303	``from_xieta(cls, xi, eta, gamma=0, T=1, P=1, Q=1, U=0, hwp=False)``
304	Construct a FocalPlane from focal plane coordinates (xi,eta).
305
306	For backwards compatibility, the signature
307	``from_xieta(names, xi, eta, gamma=0)`` is also supported; but
308	this will be removed in the future.
309
310	These are Cartesian projection plane coordinates. When looking
311	at the sky along the telescope boresight, xi is parallel to
312	increasing azimuth and eta is parallel to increasing elevation.
313	The angle gamma, which specifies the angle of polarization
314	sensitivity, is measured from the eta axis, increasing towards
315	the xi axis.
316
317	Arguments:
318	xi: An array-like of floats with shape [ndet]
319	eta: An array-like of floats with shape [ndet]
320	gamma: The detector polarization angles. [ndet], or a scalar
321	T, P: The total intensity and polarization responsivity.
322	Array-like with shape [ndet], or scalar.
323	Q, U: The Q- and U-polarization responsivity.
324	Array-like with shape [ndet], or scalar.
325	Q,U are alternatives to P,gamma for specifying the
326	polarization responsivity and angle.
327
328	So there are two ways to specify the polarization angle and
329	responsivity:
330	1. gamma and P
331	2. Q and U
332
333	Examples, assuming ndet = 2
334	* ``from_xieta(xi, eta, gamma=[0,pi/4])``
335	Constructs a FocalPlane with T and P responsivity of 1
336	and polarization angles of 0 and 45 degrees, representing
337	a Q-sensitive and U-sensitive detector.
338	* ``from_xieta(xi, eta, gamma=[0,pi/4], P=0.5)``
339	Like the above, but with a polarization responsivity of
340	just 0.5.
341	* ``from_xieta(xi, eta, gamma=[0,pi/4], T=[1,0.9], P=[0.5,0.6])``
342	Like above, but with a detector-dependent intensity and
343	polarization responsivity. There is no restriction that
344	T > P. For the pseudo-detector timestreams one gets after
345	HWP demodulation, one would have T=0 for the cos-modulated
346	and sin-modulated timestreams, for example.
347	* ``from_xieta(xi, eta, Q=[1,0], U=[0,1])``
348	Construct the FocalPlane with explicit Q and U responsivity.
349	This example is equivalent to example 1.
350
351	Usually one would either use gamma,P or Q,U. If they are
352	combined, then ``gamma_total = gamma + arctan2(U,Q)/2`` and
353	``P_tot = P * (Q2+U2)**0.5``.
354	"""
355	# The underlying code wants polangle gamma and the T and P
356	# response, but we support speifying these as the T, Q and U
357	# response too. Writing it like this handles both cases, and
358	# as a bonus(?) any combination of them
359	# FIXME: old sotodlib compat - remove later
360	xi, eta, gamma, T, P, Q, U, hwp, dets = cls._xieta_compat(args, *kwargs)	1✔
361	gamma = gamma + np.arctan2(U,Q)/2	1✔
362	P = P * (Q2+U2)**0.5	1✔
363	if hwp: gamma = -gamma	1✔
364	# Broadcast everything to 1d
365	xi, eta, gamma, T, P, _ = np.broadcast_arrays(xi, eta, gamma, T, P, [0])	1✔
366	quats = quat.rotation_xieta(xi, eta, gamma)	1✔
367	resps = np.ones((len(quats),2))	1✔
368	resps[:,0] = T	1✔
369	resps[:,1] = P	1✔
370	# FIXME: old sotodlib compat - remove dets argument later
371	return cls(quats, resps=resps, dets=dets)	1✔
372	def __repr__(self):	1✔
NEW 373	return "FocalPlane(quats=%s,resps=%s)" % (repr(self.coeffs()), repr(self.resps))	×
374	@property	1✔
375	def ndet(self): return len(self.quats)	1✔
376	def __len__(self): return self.ndet	1✔
377	def __getitem__(self, sel):	1✔
378	"""Slice the FocalPlane with slice sel, resulting in a new
379	FocalPlane with a subset of the detectors. Only a 1d slice
380	or boolean mask is supported, not integers or multidimensional
381	slices. Example: ``focal_plane[10:20]`` would make a
382	sub-FocalPlane with detector indices 10,11,...,19.
383
384	Deprecated: Temporarily also supports that sel is a detector name,
385	in which case an ``spt3g.core.quat`` is returned for that detector.
386	This is provided for backwards compatibility."""
387	# FIXME: old sotodlib compat - remove later
388	if isinstance(sel, str): return self.quats[self._dets.index(sel)]	1✔
389	# We go via .coeffs() here to get around G3VectorQuat's lack
390	# of boolean mask support
391	return FocalPlane(quats=self.coeffs()[sel], resps=self.resps[sel])	1✔
392	def items(self):	1✔
393	"""Iterate over detector quaternions and responsivities. Yields
394	``(spt3g.core.quat, array([Tresp,Presp]))`` pairs. Unlke the raw
395	entries in the .quats member, which are just numpy arrays with
396	length 4, ``spt3g.core.quat`` are proper quaternon objects that
397	support quaternion math.
398	"""
399	for q, resp in zip(self.quats, self.resps):	1✔
400	yield q, resp	1✔
401	# FIXME: old sotodlib compat - remove later
402	@staticmethod	1✔
403	def _xieta_compat(args, *kwargs):	1✔
404	# Accept the alternative format (names,xi,eta,gamma)
405	def helper(xi, eta, gamma=0, T=1, P=1, Q=1, U=0, hwp=False, dets=None):	1✔
406	return xi, eta, gamma, T, P, Q, U, hwp, dets	1✔
407	if isinstance(args[0][0], str):	1✔
NEW 408	return helper(args[1:], dets=args[0], *kwargs)	×
409	else:
410	return helper(args, *kwargs)	1✔
411	# FIXME: old sotodlib compat - remove later
412	def __setitem__(self, name, q):	1✔
413	# Make coords/pmat.py _get_asm work in old sotodlib. It
414	# expects to be able to build up a focalplane by assigning
415	# quats one at a time
NEW 416	if name in self._detmap:	×
NEW 417	i = self._detmap[i]	×
NEW 418	self.quats[i] = q	×
419	else:
NEW 420	self._dets.append(name)	×
NEW 421	self._detmap[name] = len(self._dets)-1	×
NEW 422	self.quats.append(q)	×
423	# This is inefficient, but it's temporary
NEW 424	self.resps = np.concatenate([self.resps,np.ones((1,2),np.float32)])	×
425
426	class Assembly:	1✔
427	"""This class groups together boresight and detector offset
428	information for use by the Projectionist.
429
430	The basic implementation simply attachs some number of detectors
431	rigidly to a boresight. But this abstraction layer could
432	facilitate more complex arrangements, eventually.
433
434	"""
435	def __init__(self, collapse=False):	1✔
436	self.collapse = collapse	1✔
437
438	@classmethod	1✔
439	def attach(cls, sight_line, fplane):	1✔
440	"""Create an Assembly based on a CelestialSightLine and a FocalPlane.
441
442	Args:
443	sight_line (CelestialSightLine): The position and
444	orientation of the boresight. This can just be a
445	G3VectorQuat if you don't have a whole
446	CelestialSightLine handy.
447	fplane (FocalPlane): The "offsets" of each detector
448	relative to the boresight, and their response to
449	intensity and polarization
450	"""
451	self = cls()	1✔
452	if isinstance(sight_line, quat.G3VectorQuat):	1✔
453	self.Q = sight_line	1✔
454	else:
455	self.Q = sight_line.Q	1✔
456	self.fplane = fplane	1✔
457	return self	1✔
458
459	@classmethod	1✔
460	def for_boresight(cls, sight_line):	1✔
461	"""Returns an Assembly where a single detector serves as a dummy for
462	the boresight."""
UNCOV 463	self = cls(collapse=True)	×
464	self.Q = sight_line.Q	×
NEW 465	self.fplane = FocalPlane.boresight()	×
466	return self	×
467	# FIXME: old sotodlib compat - remove later
468	@property	1✔
469	def dets(self): return self.fplane	1✔

simonsobs / so3g / 12832805671

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous