Skip to content

Commit

Permalink
Ready for new review
Browse files Browse the repository at this point in the history
  • Loading branch information
amaurea committed Jan 16, 2025
1 parent 63f92f1 commit 5efc1f7
Show file tree
Hide file tree
Showing 4 changed files with 22 additions and 35 deletions.
11 changes: 2 additions & 9 deletions demos/demo_proj_hp.py
Original file line number Diff line number Diff line change
Expand Up @@ -191,13 +191,6 @@ def simulate_detectors(ndets, rmax, seed=0):
eta = rr * np.sin(phi)
return xi, eta, gamma

def make_fp(xi, eta, gamma):
fp = so3g.proj.quat.rotation_xieta(xi, eta, gamma)
so3g_fp = so3g.proj.FocalPlane()
for i, q in enumerate(fp):
so3g_fp[f'a{i}'] = q
return so3g_fp

def extract_coord(sight, fp, isignal=0, groupsize=100):
coord_out = []
for ii in range(0, len(fp), groupsize):
Expand All @@ -208,14 +201,14 @@ def extract_coord(sight, fp, isignal=0, groupsize=100):

def make_signal(ndets, rmax, sight, isignal, seed=0):
xi, eta, gamma = simulate_detectors(ndets, rmax, seed)
fp = so3g.proj.quat.rotation_xieta(xi, eta, gamma)
fp = so3g.proj.FocalPlane.from_xieta(xi, eta, gamma)
signal = extract_coord(sight, fp, isignal)
return signal

def make_tod(time, az, el, ndets, rmax, isignal, seed=0):
sight = so3g.proj.CelestialSightLine.naive_az_el(time, az, el)
signal = make_signal(ndets, rmax, sight, isignal, seed)
so3g_fp = make_fp(*simulate_detectors(ndets, rmax, seed))
so3g_fp = so3g.proj.FocalPlane.from_xieta(*simulate_detectors(ndets, rmax, seed))
asm = so3g.proj.Assembly.attach(sight, so3g_fp)
return signal, asm

Expand Down
20 changes: 7 additions & 13 deletions docs/proj.rst
Original file line number Diff line number Diff line change
Expand Up @@ -125,22 +125,16 @@ the quaternion for the first detector, while ``fp.resps[0]``
gives its total intensity and polarization responsivity.::

>>> fp.quats[2]
array([ 0.86601716, 0.00377878, -0.00218168, 0.49999524])
spt3g.core.quat(0.866017,0.00377878,-0.00218168,0.499995)
>>> fp.resps[2]
array([1., 1.], dtype=float32)

As you can see, the default responsivity is 1 for both total
intensty and polarization. Note that ``.quats`` contains
quaternion *coefficients*, not quaternion *objects*. To do
quaternion math, you need to convert them to actual quaternion
objects, e.g. ``q = fp.quats[0]``, or convert
them all at once with ``qs = spt3g.core.G3VectorQuat(fp.quats)```.

To represent detectors with responsivity different from 1,
use the ``T`` and ``P`` arguments to :func:`from_xieta()`
to set the total intensity and polarization responsivity
respectively. These can be either single numbers or
array-likes with lengths ``ndet``.::
intensty and polarization. To represent detectors with
responsivity different from 1, use the ``T`` and ``P``
arguments to :func:`from_xieta()` to set the total intensity
and polarization responsivity respectively. These can be
either single numbers or array-likes with lengths ``ndet``.::

xi = np.array([-0.5, 0.0, 0.5]) * DEG
eta = np.zeros(3)
Expand All @@ -157,7 +151,7 @@ specify these directly. The example above is equivalent to::
xi = np.array([-0.5, 0.0, 0.5]) * DEG
eta = np.zeros(3)
gamma = np.array([0,30,60]) * DEG
T = np.array([1.0, 0.9])
T = np.array([1.0, 0.9, 1.1])
Q = np.array([0.5, 0.3, -0.025])
U = np.array([0.0, 0.51961524, 0.04330127])
fp2 = so3g.proj.FocalPlane.from_xieta(xi, eta, T=T, Q=Q, U=U)
Expand Down
2 changes: 1 addition & 1 deletion python/proj/coords.py
Original file line number Diff line number Diff line change
Expand Up @@ -231,7 +231,7 @@ def coords(self, fplane=None, output=None):
fplane = FocalPlane.boresight()
if output is not None:
output = output[None]
output = p.coords(self.Q, fplane.coeffs(), output)
output = p.coords(self.Q, fplane.quats), output)
if collapse:
output = output[0]
return output
Expand Down
24 changes: 12 additions & 12 deletions python/proj/wcs.py
Original file line number Diff line number Diff line change
Expand Up @@ -183,7 +183,7 @@ def get_pixels(self, assembly):
"""
projeng = self.get_ProjEng('TQU')
q_native = self._cache_q_fp_to_native(assembly.Q)
return projeng.pixels(q_native, assembly.fplane.coeffs(), None)
return projeng.pixels(q_native, assembly.fplane.quats, None)

def get_pointing_matrix(self, assembly):
"""Get the pointing matrix information, which is to say both the pixel
Expand All @@ -199,7 +199,7 @@ def get_pointing_matrix(self, assembly):
"""
projeng = self.get_ProjEng('TQU')
q_native = self._cache_q_fp_to_native(assembly.Q)
return projeng.pointing_matrix(q_native, assembly.fplane.coeffs(), assembly.fplane.resps, None, None)
return projeng.pointing_matrix(q_native, assembly.fplane.quats, assembly.fplane.resps, None, None)

def get_coords(self, assembly, use_native=False, output=None):
"""Get the spherical coordinates for the provided pointing Assembly.
Expand All @@ -225,7 +225,7 @@ def get_coords(self, assembly, use_native=False, output=None):
q_native = self._cache_q_fp_to_native(assembly.Q)
else:
q_native = assembly.Q
return projeng.coords(q_native, assembly.fplane.coeffs(), output)
return projeng.coords(q_native, assembly.fplane.quats, output)

def get_planar(self, assembly, output=None):
"""Get projection plane coordinates for all detectors at all times.
Expand All @@ -240,7 +240,7 @@ def get_planar(self, assembly, output=None):
"""
q_native = self._cache_q_fp_to_native(assembly.Q)
projeng = self.get_ProjEng('TQU')
return projeng.coords(q_native, assembly.fplane.coeffs(), output)
return projeng.coords(q_native, assembly.fplane.quats, output)

def zeros(self, super_shape):
"""Return a map, filled with zeros, with shape (super_shape,) +
Expand Down Expand Up @@ -272,7 +272,7 @@ def to_map(self, signal, assembly, output=None, det_weights=None,
comps = self._guess_comps(self._get_map_shape(output))
projeng = self.get_ProjEng(comps)
q_native = self._cache_q_fp_to_native(assembly.Q)
map_out = projeng.to_map(output, q_native, assembly.fplane.coeffs(),
map_out = projeng.to_map(output, q_native, assembly.fplane.quats,
assembly.fplane.resps, signal, det_weights, threads)
return map_out

Expand All @@ -298,7 +298,7 @@ def to_weights(self, assembly, output=None, det_weights=None,
comps = self._guess_comps(shape[1:])
projeng = self.get_ProjEng(comps)
q_native = self._cache_q_fp_to_native(assembly.Q)
map_out = projeng.to_weight_map(output, q_native, assembly.fplane.coeffs(),
map_out = projeng.to_weight_map(output, q_native, assembly.fplane.quats,
assembly.fplane.resps, det_weights, threads)
return map_out

Expand All @@ -318,7 +318,7 @@ def from_map(self, src_map, assembly, signal=None, comps=None):
comps = self._guess_comps(self._get_map_shape(src_map))
projeng = self.get_ProjEng(comps)
q_native = self._cache_q_fp_to_native(assembly.Q)
signal_out = projeng.from_map(src_map, q_native, assembly.fplane.coeffs(),
signal_out = projeng.from_map(src_map, q_native, assembly.fplane.quats,
assembly.fplane.resps, signal)
return signal_out

Expand Down Expand Up @@ -362,7 +362,7 @@ def assign_threads(self, assembly, method='domdir', n_threads=None):
if method == 'simple':
projeng = self.get_ProjEng('T')
q_native = self._cache_q_fp_to_native(assembly.Q)
omp_ivals = projeng.pixel_ranges(q_native, assembly.fplane.coeffs(), None, n_threads)
omp_ivals = projeng.pixel_ranges(q_native, assembly.fplane.quats, None, n_threads)
return wrap_ivals(omp_ivals)

elif method == 'domdir':
Expand Down Expand Up @@ -402,7 +402,7 @@ def assign_threads_from_map(self, assembly, tmap, n_threads=None):
projeng = self.get_ProjEng('T')
q_native = self._cache_q_fp_to_native(assembly.Q)
n_threads = mapthreads.get_num_threads(n_threads)
omp_ivals = projeng.pixel_ranges(q_native, assembly.fplane.coeffs(), tmap, n_threads)
omp_ivals = projeng.pixel_ranges(q_native, assembly.fplane.quats, tmap, n_threads)
return wrap_ivals(omp_ivals)

def get_active_tiles(self, assembly, assign=False):
Expand Down Expand Up @@ -441,7 +441,7 @@ def get_active_tiles(self, assembly, assign=False):
projeng = self.get_ProjEng('T')
q_native = self._cache_q_fp_to_native(assembly.Q)
# This returns a G3VectorInt of length n_tiles giving count of hits per tile.
hits = np.array(projeng.tile_hits(q_native, assembly.fplane.coeffs()))
hits = np.array(projeng.tile_hits(q_native, assembly.fplane.quats))
tiles = np.nonzero(hits)[0]
hits = hits[tiles]
if assign is True:
Expand All @@ -460,7 +460,7 @@ def get_active_tiles(self, assembly, assign=False):
# print(f"Warning: Threads poorly balanced. Max/mean hits = {max_ratio}")

# Now paint them into Ranges.
R = projeng.tile_ranges(q_native, assembly.fplane.coeffs(), group_tiles)
R = projeng.tile_ranges(q_native, assembly.fplane.quats, group_tiles)
R = wrap_ivals(R)
return {
'active_tiles': list(tiles),
Expand Down Expand Up @@ -795,7 +795,7 @@ def get_coords(self, assembly, use_native=False, output=None):
q_native = self._cache_q_fp_to_native(assembly.Q)
else:
q_native = assembly.Q
return projeng.coords(q_native, assembly.fplane.coeffs(), output)
return projeng.coords(q_native, assembly.fplane.quats, output)

def from_map(self, src_map, assembly, signal=None, comps=None):
"""De-project from a map, returning a Signal-like object.
Expand Down

0 comments on commit 5efc1f7

Please sign in to comment.