Skip to content

Commit

Permalink
More documentation. Updated numbers in examples which seem to have ro…
Browse files Browse the repository at this point in the history
…tted at some point
  • Loading branch information
amaurea committed Jan 14, 2025
1 parent bac61f9 commit c381023
Show file tree
Hide file tree
Showing 2 changed files with 26 additions and 29 deletions.
29 changes: 12 additions & 17 deletions docs/proj.rst
Original file line number Diff line number Diff line change
Expand Up @@ -174,7 +174,7 @@ As expected, these coordinates are close to the ones computed before,
for the boresight::

>>> print(ra / DEG, dec / DEG)
[13.06731917] [-53.12633433]
[14.73387143] [-53.1250149]

But the more expedient way to get pointing for multiple detectors is
to call :func:`coords() <CelestialSightLine.coords>` with the
Expand All @@ -185,11 +185,6 @@ FocalPlane object as first argument::
array([[ 0.24261138, -0.92726871, 0.86536489, 0.5011423 ]]),
array([[ 0.22806774, -0.92722945, 0.50890634, 0.86082189]])]

OrderedDict([('a', array([[ 0.22806774, -0.92722945, -0.9999468 ,
0.01031487]])), ('b', array([[ 0.24261138, -0.92726871, -0.86536489,
-0.5011423 ]])), ('c', array([[ 0.25715457, -0.92720643, -0.48874018,
-0.87242939]]))])

To be clear, ``coords()`` now returns a dictionary whose keys are the
detector names. Each value is an array with shape (n_time,4), and at
each time step the 4 elements of the array are: ``(lon, lat,
Expand Down Expand Up @@ -256,17 +251,17 @@ from the pixels in pmap, which are the coordinates of the centers of
the pixels::

>>> [x/DEG for x in pix_ra]
[array([13.04], dtype=float32), array([13.88], dtype=float32), array([14.72], dtype=float32)]
[array([14.740001], dtype=float32), array([13.900001], dtype=float32), array([13.059999], dtype=float32)]
>>> [x/DEG for x in pix_dec]
[array([-53.100002], dtype=float32), array([-53.100002], dtype=float32), array([-53.100002], dtype=float32)]
[array([-53.12], dtype=float32), array([-53.12], dtype=float32), array([-53.12], dtype=float32)]

If you are not getting what you expect, you can grab the pixel indices
inferred by the projector -- perhaps your pointing is taking you off
the map (in which case the pixel indices would return value -1)::

>>> p.get_pixels(asm)
[array([[ 45, 148]], dtype=int32), array([[ 45, 106]], dtype=int32),
array([[45, 64]], dtype=int32)]
[array([[44, 63]], dtype=int32), array([[ 44, 105]], dtype=int32),
array([[ 44, 147]], dtype=int32)]

Let's project signal into an intensity map using
:func:`Projectionist.to_map`::
Expand All @@ -281,9 +276,9 @@ Inspecting the map, we see our signal values occupy the three non-zero
pixels:

>>> map_out.nonzero()
(array([0, 0, 0]), array([45, 45, 45]), array([ 64, 106, 148]))
(array([0, 0, 0]), array([44, 44, 44]), array([ 63, 105, 147]))
>>> map_out[map_out!=0]
array([100., 10., 1.])
array([ 1., 10., 100.])


If we run this projection again, but pass in this map as a starting
Expand All @@ -292,7 +287,7 @@ point, the signal will be added to the map:
>>> p.to_map(signal, asm, output=map_out, comps='T')
array([[[0., 0., 0., ..., 0.]]])
>>> map_out[map_out!=0]
array([200., 20., 2.])
array([ 2., 20., 200.])

If we instead want to treat the signal as coming from
polarization-sensitive detectors, we can request components
Expand All @@ -305,8 +300,8 @@ according to the projected detector angle on the sky::

>>> map_pol.shape
(3, 100, 200)
>>> map_pol[:,45,106]
array([10. , 4.97712803, 8.673419 ])
>>> map_pol[:,44,105]
array([10., 4.97712803, 8.673419])

For the most basic map-making, the other useful operation is the
:func:`Projectionist.to_weights` method. This is used to compute
Expand All @@ -330,7 +325,7 @@ the upper diagonal has been filled in, for efficiency reasons...::

>>> weight_out.shape
(3, 3, 100, 200)
>>> weight_out[...,45,106]
>>> weight_out[...,44,105]
array([[1. , 0.49771279, 0.86734194],
[0. , 0.24771802, 0.43168718],
[0. , 0. , 0.75228202]])
Expand Down Expand Up @@ -393,7 +388,7 @@ Inspecting::

>>> threads
RangesMatrix(4,3,1)
>>> map_pol2[:,45,106]
>>> map_pol2[:,44,105]
array([10. , 4.97712898, 8.67341805])

The same ``threads`` result can be passed to ``p.to_weights``.
Expand Down
26 changes: 14 additions & 12 deletions python/proj/coords.py
Original file line number Diff line number Diff line change
Expand Up @@ -320,24 +320,26 @@ def from_xieta(cls, xi, eta, gamma=0, T=1, P=1, Q=1, U=0, hwp=False):
2. Q and U
Examples, assuming ndet = 2
* from_xieta(xi, eta, gamma=[0,pi/4])
* ``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)
* ``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])
* ``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])
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.
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
Expand All @@ -364,20 +366,20 @@ 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]
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.
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
``(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
length 4, ``spt3g.core.quat`` are proper quaternon objects that
support quaternion math.
"""
for q, resp in zip(quat.G3VectorQuat(self.quats), self.resps):
Expand Down

0 comments on commit c381023

Please sign in to comment.