Skip to content

Commit

Permalink
Better docstrings and comments
Browse files Browse the repository at this point in the history
  • Loading branch information
@amaurea
amaurea committed Jan 14, 2025
1 parent c396f10 commit 4624992
Show file tree
Hide file tree
Showing 2 changed files with 105 additions and 51 deletions.
93 changes: 87 additions & 6 deletions python/proj/coords.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -214,7 +214,7 @@ def coords(self, fplane=None, output=None):
Arguments:
fplane: A FocalPlane object representing the detector
offsets and responses, or None
output: An optional structure for holding the results. For
output: An optional structure for holding the results. For
that to work, each element of output must support the
buffer protocol.
Expand All @@ -240,11 +240,33 @@ class FocalPlane:
"""This class represents the detector positions and intensity and
polarization responses in the focal plane.
This used to be a subclass of OrderedDict, but this was hard to
generalize to per-detector polarization efficiency.
Members:
quats: Array of float64 with shape [ndet,4] representing the
pointing quaternions for each detector
resps: Array of float64 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
Expand All @@ -263,10 +285,57 @@ def __init__(self, quats=None, resps=None, dets=None):
self._detmap = {name:i for i,name in enumerate(self._dets)}
@classmethod
def boresight(cls):
"""Construct a FocalPlane representing a single detector with a
responsivity of one pointing along the telescope boresight"""
quats = np.array([[1,0,0,0]],np.float64)
return cls(quats=quats)
@classmethod
def 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).
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. T < P is valid, even though it
wouldn't make sense for most detector types.
* 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
Expand All @@ -290,10 +359,24 @@ def __repr__(self):
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 1d slice 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 quat.G3VectorQuat(self.quats)[self._dets.index(sel)]
return FocalPlane(quats=self.quats[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(quat.G3VectorQuat(self.quats), self.resps):
yield q, resp
# FIXME: old sotodlib compat - remove later
Expand Down Expand Up @@ -356,9 +439,7 @@ def attach(cls, sight_line, fplane):
@classmethod
def for_boresight(cls, sight_line):
"""Returns an Assembly where a single detector serves as a dummy for
the boresight.
"""
the boresight."""
self = cls(collapse=True)
self.Q = sight_line.Q
self.fplane = FocalPlane.boresight()
Expand Down
63 changes: 18 additions & 45 deletions src/Projection.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -125,55 +125,28 @@ inline bool isNone(const bp::object &pyo)

// Arbitrary detector response support
// -----------------------------------
// Useful to support detectors with polarization response < 1. Also occasionally
// useful with detectors with 0 total intensity response. Currently the detector
// response is only encoded indirectly through the detector angle on the sky, part
// of the per-detector quaternion. Hard to cram the rest of the response into this,
// so must be handled with a new argument.
// It's useful to support detectors with polarization response < 1. It's also
// sometimes useful with detectors with 0 total intensity response. This is
// implemented as follows.
//
// spin_proj_factors: cos(2psi),sin(2psi) → detector on-sky TQU response
// psi is here the projected detector angle, which is computed by Pointer::GetCoords
// It makes sense for this to deal with angles since it's the result of a quaternion
// pointing calculation.
// At the top level, the functions where the responsivity is relevant, like
// from_map, to_map etc. take a bp::object reponse argument, representing
// the [ndet,2] detector responsivities.
//
// Probably best to leave GetCoords alone, and instead augment spin_proj_factors with
// T and P support. For example:
// This is then converted into a BufferWrapper<FSIGNAL>, before each detector's
// responsivity is extracted using get_response(BufferWrapper<FSIGNAL> & buf, int i_det),
// which returns a Response { FSIGNAL T, P; } struct.
//
// void spin_proj_factors<SpinTQU>(const double * coords, const double * response, FSIGNAL * projfacs) {
// const double c = coords[2];
// const double s = coords[3];
// projfacs[0] = response[0];
// projfacs[1] = response[1]*(c*c - s*s);
// projfacs[2] = response[1]*(2*c*s);
// }
// Finally, this Response struct is passed to spin_proj_factors, e.g.
//
// None of the arguments to to_map and its relatives are a natural place to put response.
// It could be shoehorned into pofs by making it length 6, but the first four components
// would hold quaternion coefficients making it unnatural to have the last 2 be responses.
// It's probably better to modify these functions to take an additional response argument,
// e.g.
//
// template<typename C, typename P, typename S>
// bp::object ProjectionEngine<C,P,S>::to_map(
// bp::object map, bp::object pbore, bp::object pofs, bp::object det_response,
// bp::object signal, bp::object det_weights, bp::object thread_intervals)
//
// where these would then be extracted using
//
// auto _det_response = BufferWrapper<FSIGNAL>("det_response", det_response, false, vector<int>{n_det,2});
//
// Could make it optional, but easier to handle this on the python side.
//
// List of functions needing modification:
// * ProjectionEngine<C,P,S>::pointing_matrix
// * ProjectionEngine<C,P,S>::from_map
// * ProjectionEngine<C,P,S>::to_map
// * ProjectionEngine<C,P,S>::to_weight_map
// * to_map_single_thread
// * to_weight_map_single_thread
//
// * python/wcs.py/Projectionist
// * python/mapthreads.py/get_threads_domdir
// void spin_proj_factors<SpinTQU>(const double* coords, const Response & response, FSIGNAL *projfacs)
// {
// const double c = coords[2];
// const double s = coords[3];
// projfacs[0] = response.T;
// projfacs[1] = response.P*(c*c - s*s);
// projfacs[2] = response.P*(2*c*s);
// }

// State the template<CoordSys> options
class ProjQuat;
Expand Down

0 comments on commit 4624992

Please sign in to comment.