From 360f073891225ffe913d86cffa253db22eb49ba8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sigurd=20Kirkevold=20N=C3=A6ss?= <sigurdkn@astro.uio.no>
Date: Tue, 18 Jun 2024 05:21:20 -0400
Subject: [PATCH 1/4] Plans for arbitrary detector response stuff

---
 src/Projection.cxx | 52 ++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 52 insertions(+)

diff --git a/src/Projection.cxx b/src/Projection.cxx
index 0a9ae34..ae7951d 100644
--- a/src/Projection.cxx
+++ b/src/Projection.cxx
@@ -121,6 +121,58 @@ inline bool isNone(const bp::object &pyo)
 // shape [1][ndet][nrange]. More complicated schemes could be used, but this
 // seems to work well enough.
 
+// 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.
+//
+// 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.
+//
+// Probably best to leave GetCoords alone, and instead augment spin_proj_factors with
+// T and P support. For example:
+//
+// 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);
+// }
+//
+// 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
+
 // State the template<CoordSys> options
 class ProjQuat;
 class ProjFlat;

From b790f7c924e281ba15f9fded5654965292680c47 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sigurd=20Kirkevold=20N=C3=A6ss?= <sigurdkn@astro.uio.no>
Date: Sat, 22 Jun 2024 17:10:24 -0400
Subject: [PATCH 2/4] Detector response now supported in c++, and represented
 in a redone FocalPlane class in python

---
 demos/demo_proj1.py   |  21 +++----
 demos/demo_proj2.py   |  14 +++--
 include/Projection.h  |   8 +--
 python/proj/coords.py | 126 +++++++++++++++---------------------------
 python/proj/wcs.py    |  14 ++---
 src/Projection.cxx    |  72 ++++++++++++++++--------
 6 files changed, 123 insertions(+), 132 deletions(-)

diff --git a/demos/demo_proj1.py b/demos/demo_proj1.py
index 81337ce..00b6c18 100644
--- a/demos/demo_proj1.py
+++ b/demos/demo_proj1.py
@@ -47,6 +47,7 @@
 pe = test_utils.get_proj(system, 'TQU', pxz, tiled=args.tiled)
 ptg = test_utils.get_boresight_quat(system, x, y)
 ofs = test_utils.get_offsets_quat(system, dx, dy, polphi)
+resp= np.ones((n_det,2), np.float32)
 
 sig = np.ones((1,n_det,n_t), 'float32') * .5
 
@@ -119,7 +120,7 @@ def map_delta(a, b):
     print('Get spin projection factors, too.',
           end='\n ... ')
     with Timer() as T:
-        pix2, spin_proj = pe.pointing_matrix(ptg, ofs, None, None)
+        pix2, spin_proj = pe.pointing_matrix(ptg, ofs, resp, None, None)
 
     del pix, pix_list, pix3, pix2, spin_proj
 
@@ -139,30 +140,30 @@ def map_delta(a, b):
     else:
         map1 += np.array([1,0,0])[:,None,None]
     with Timer() as T:
-        sig1 = pe.from_map(map1, ptg, ofs, None)
+        sig1 = pe.from_map(map1, ptg, ofs, resp, None)
 
 if 1:
     print('Project TOD-to-map (TQU)', end='\n ... ')
     map0 = pe.zeros(3)
     sig_list = [x for x in sig[0]]
     with Timer() as T:
-        map1 = pe.to_map(map0,ptg,ofs,sig_list,None,None)
+        map1 = pe.to_map(map0,ptg,ofs,resp,sig_list,None,None)
 
 if 1:
     print('TOD-to-map again but with None for input map', end='\n ... ')
     with Timer() as T:
-        map1 = pe.to_map(None,ptg,ofs,sig_list,None,None)
+        map1 = pe.to_map(None,ptg,ofs,resp,sig_list,None,None)
 
 if 1:
     print('Project TOD-to-weights (TQU)', end='\n ... ')
     map0 = pe.zeros((3, 3))
     with Timer() as T:
-        map2 = pe.to_weight_map(map0,ptg,ofs,None,None)
+        map2 = pe.to_weight_map(map0,ptg,ofs,resp,None,None)
 
 if 1:
     print('TOD-to-weights again but with None for input map', end='\n ... ')
     with Timer() as T:
-        map2 = pe.to_weight_map(None,ptg,ofs,None,None)
+        map2 = pe.to_weight_map(None,ptg,ofs,resp,None,None)
 
 print('Compute thread assignments (OMP prep)... ', end='\n ... ')
 with Timer():
@@ -171,12 +172,12 @@ def map_delta(a, b):
 if 1:
     print('TOD-to-map with OMP (%s): ' % n_omp, end='\n ... ')
     with Timer() as T:
-        map1o = pe.to_map(None,ptg,ofs,sig_list,None,threads)
+        map1o = pe.to_map(None,ptg,ofs,resp,sig_list,None,threads)
 
 if 1:
     print('TOD-to-weights with OMP (%s): ' % n_omp, end='\n ... ')
     with Timer() as T:
-        map2o = pe.to_weight_map(None,ptg,ofs,None,threads)
+        map2o = pe.to_weight_map(None,ptg,ofs,resp,None,threads)
 
     print('Checking that OMP and non-OMP forward calcs agree: ', end='\n ... ')
     assert map_delta(map1, map1o) == 0
@@ -187,7 +188,7 @@ def map_delta(a, b):
 if 1:
     print('Cache pointing matrix.', end='\n ...')
     with Timer() as T:
-        pix_idx, spin_proj = pe.pointing_matrix(ptg, ofs, None, None)
+        pix_idx, spin_proj = pe.pointing_matrix(ptg, ofs, resp, None, None)
     pp = test_utils.get_proj_precomp(args.tiled)
 
     print('Map-to-TOD using precomputed pointing matrix',
@@ -212,7 +213,7 @@ def map_delta(a, b):
 
     print('Checking that it agrees with on-the-fly',
           end='\n ...')
-    sig1f = pe.from_map(map1, ptg, ofs, None)
+    sig1f = pe.from_map(map1, ptg, ofs, resp, None)
     thresh = map1[0].std() * 1e-6
     assert max([np.abs(a - b).max() for a, b in zip(sig1f, sig1p)]) < thresh
     print('yes')
diff --git a/demos/demo_proj2.py b/demos/demo_proj2.py
index 839b9df..ae23207 100644
--- a/demos/demo_proj2.py
+++ b/demos/demo_proj2.py
@@ -82,7 +82,7 @@ def get_pixelizor(emap):
 
 ptg = test_utils.get_boresight_quat(system, x*np.pi/180, y*np.pi/180)
 ofs = test_utils.get_offsets_quat(system, dx, dy, polphi)
-
+resp= np.ones((n_det,2),np.float32)
 
 #
 # Projection tests.
@@ -91,7 +91,7 @@ def get_pixelizor(emap):
 pe = test_utils.get_proj(system, 'TQU')(pxz)
 
 # Project the map into time-domain.
-sig0 = pe.from_map(beam, ptg, ofs, None)
+sig0 = pe.from_map(beam, ptg, ofs, resp, None)
 sig0 = np.array(sig0)
 
 # Add some noise...
@@ -118,10 +118,10 @@ def get_pixelizor(emap):
 # Then back to map.
 print('Create timestream...', end='\n ... ')
 with Timer():
-    map1 = pe.to_map(None, ptg, ofs, sig1, det_weights, None)
+    map1 = pe.to_map(None, ptg, ofs, resp, sig1, det_weights, None)
 
 # Get the weight map (matrix).
-wmap1 = pe.to_weight_map(None, ptg, ofs, det_weights, None)
+wmap1 = pe.to_weight_map(None, ptg, ofs, resp, det_weights, None)
 wmap1[1,0] = wmap1[0,1]  # fill in unpopulated entries...
 wmap1[2,0] = wmap1[0,2]
 wmap1[2,1] = wmap1[1,2]
@@ -151,6 +151,8 @@ def get_pixelizor(emap):
 
 sigd = (sig1[0::2,:] - sig1[1::2,:]) / 2
 ofsd = (ofs[::2,...])
+respd= resp[::2]
+
 if not args.uniform_weights:
     det_weights = det_weights[::2]
 
@@ -160,12 +162,12 @@ def get_pixelizor(emap):
 # Bin the map again...
 print('to map...', end='\n ... ')
 with Timer():
-    map1d = pe.to_map(None, ptg, ofsd, sigd, det_weights, None)
+    map1d = pe.to_map(None, ptg, ofsd, respd, sigd, det_weights, None)
 
 # Bin the weights again...
 print('to weight map...', end='\n ... ')
 with Timer():
-    wmap1d = pe.to_weight_map(None, ptg, ofsd, det_weights, None)
+    wmap1d = pe.to_weight_map(None, ptg, ofsd, respd, det_weights, None)
 
 wmap1d[1,0] = wmap1d[0,1] # fill in unpopulated entries again...
 
diff --git a/include/Projection.h b/include/Projection.h
index 8890e09..d5c8c21 100644
--- a/include/Projection.h
+++ b/include/Projection.h
@@ -48,17 +48,17 @@ class ProjectionEngine {
     bp::object pixels(bp::object pbore, bp::object pofs, bp::object pixel);
     vector<int> tile_hits(bp::object pbore, bp::object pofs);
     bp::object tile_ranges(bp::object pbore, bp::object pofs, bp::object tile_lists);
-    bp::object pointing_matrix(bp::object pbore, bp::object pofs,
+    bp::object pointing_matrix(bp::object pbore, bp::object pofs, bp::object response,
                                bp::object pixel, bp::object proj);
     bp::object zeros(bp::object shape);
     bp::object pixel_ranges(bp::object pbore, bp::object pofs, bp::object map, int n_domain=-1);
     bp::object from_map(bp::object map, bp::object pbore, bp::object pofs,
-                        bp::object signal);
-    bp::object to_map(bp::object map, bp::object pbore, bp::object pofs,
+                        bp::object response, bp::object signal);
+    bp::object to_map(bp::object map, bp::object pbore, bp::object pofs, bp::object response,
                       bp::object signal, bp::object det_weights,
                       bp::object thread_intervals);
     bp::object to_weight_map(bp::object map, bp::object pbore, bp::object pofs,
-                             bp::object det_weights, bp::object thread_intervals);
+                  bp::object response, bp::object det_weights, bp::object thread_intervals);
 
     int comp_count() const;
     int index_count() const;
diff --git a/python/proj/coords.py b/python/proj/coords.py
index a4fea0e..06d816b 100644
--- a/python/proj/coords.py
+++ b/python/proj/coords.py
@@ -208,93 +208,65 @@ def for_horizon(cls, t, az, el, roll=None, site=None, weather=None):
             )
         return self
 
-    def coords(self, det_offsets=None, output=None):
+    def coords(self, fplane=None, output=None):
         """Get the celestial coordinates of each detector at each time.
 
         Arguments:
-          det_offset: A dict or list or array of detector offset
-            tuples.  If each tuple has 2 elements or 3 elements, the
-            typles are interpreted as X,Y[,phi] coordinates in the
-            conventional way.  If 4 elements, it's interpreted as a
-            quaternion.  If this argument is None, then the boresight
-            pointing is returned.
+          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 the det_offset was passed in as a dict, then a dict with the same
-          keys is returned.  Otherwise a list is returned.  For each
-          detector, the corresponding returned object is an array with
-          shape (n_samp, 4) where the four components correspond to
-          longitude, latitude, cos(psi), sin(psi).
-
+          If fplane is not 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 = (det_offsets is None)
+        collapse = (fplane is None)
         if collapse:
-            det_offsets = np.array([[1., 0., 0., 0.]])
+            fplane = FocalPlane()
             if output is not None:
                 output = output[None]
-        redict = isinstance(det_offsets, dict)
-        if redict:
-            keys, det_offsets = zip(*det_offsets.items())
-            if isinstance(det_offsets[0], quat.quat):
-                # Individual quat doesn't array() properly...
-                det_offsets = np.array(quat.G3VectorQuat(det_offsets))
-            else:
-                det_offsets = np.array(det_offsets)
-            if isinstance(output, dict):
-                output = [output[k] for k in keys]
-        if np.shape(det_offsets)[1] == 2:
-            # XY
-            x, y = det_offsets[:, 0], det_offsets[:, 1]
-            theta = np.arcsin((x**2 + y**2)**0.5)
-            v = np.array([-y, x]) / theta
-            v[np.isnan(v)] = 0.
-            det_offsets = np.transpose([np.cos(theta/2),
-                                        np.sin(theta/2) * v[0],
-                                        np.sin(theta/2) * v[1],
-                                        np.zeros(len(theta))])
-            det_offsets = quat.G3VectorQuat(det_offsets.copy())
-        output = p.coords(self.Q, det_offsets, output)
-        if redict:
-            output = OrderedDict(zip(keys, output))
+        output = p.coords(self.Q, fplane.quats, fplane.resps, 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.
 
-class FocalPlane(OrderedDict):
-    """This class collects the focal plane positions and orientations for
-    multiple named detectors.  The classmethods can be used to
-    construct the object from some common input formats.
-
+    This used to be a subclass of OrderedDict, but this was hard to
+    generalize to per-detector polarization efficiency.
     """
-    @classmethod
-    def from_xieta(cls, names, xi, eta, gamma=0):
-        """Creates a FocalPlane object for a set of detector positions
-        provided in xieta projection plane coordinates.
-
-        Args:
-            names: vector of detector names.
-            xi: vector of xi positions, in radians.
-            eta: vector of eta positions, in radians.
-            gamma: vector or scalar detector orientation, in radians.
-
-        The (xi,eta) coordinates are Cartesian projection plane
-        coordinates.  When looking at the sky along the telescope
-        boresight, xi 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.
-
-        """
-        qs = quat.rotation_xieta(np.asarray(xi), np.asarray(eta), np.asarray(gamma))
-        return cls([(n,q) for n, q in zip(names, qs)])
+    def __init__(self, quats=None, resps=None):
+        if quats is None: quats = np.array([[1.0,0.0,0.0,0.0]])
+        # 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
+        self.quats = np.array(quat.G3VectorQuat(quats))
+        self.resps = np.ones((len(self.quats),2),np.float32)
+        if resps is not None:
+            self.resps[:] = resps
 
+    @classmethod
+    def from_xieta(cls, xi, eta, gamma=0, T=1, P=1, Q=1, U=0):
+        # 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
+        gamma = gamma + np.arctan2(U,Q)/2
+        P     = P * (Q**2+U**2)**0.5
+        quats = quat.rotation_xieta(np.asarray(xi), np.asarray(eta), np.asarray(gamma))
+        resps = np.ones((len(quats),2))
+        resps[:,0] = T
+        resps[:,1] = P
+        return cls(quats, resps=resps)
 
 class Assembly:
     """This class groups together boresight and detector offset
@@ -305,12 +277,11 @@ class Assembly:
     facilitate more complex arrangements, eventually.
 
     """
-    def __init__(self, keyed=False, collapse=False):
-        self.keyed = keyed
+    def __init__(self, collapse=False):
         self.collapse = collapse
 
     @classmethod
-    def attach(cls, sight_line, det_offsets):
+    def attach(cls, sight_line, fplane):
         """Create an Assembly based on a CelestialSightLine and a FocalPlane.
 
         Args:
@@ -318,25 +289,16 @@ def attach(cls, sight_line, det_offsets):
                 orientation of the boresight.  This can just be a
                 G3VectorQuat if you don't have a whole
                 CelestialSightLine handy.
-            det_offsets (FocalPlane): The "offsets" of each detector
-                relative to the boresight.
-
+            fplane (FocalPlane): The "offsets" of each detector
+                relative to the boresight, and their response to
+                intensity and polarization
         """
-        keyed = isinstance(det_offsets, dict)
         self = cls(keyed=keyed)
         if isinstance(sight_line, quat.G3VectorQuat):
             self.Q = sight_line
         else:
             self.Q = sight_line.Q
-        if self.keyed:
-            self.keys = list(det_offsets.keys())
-            self.dets = [det_offsets[k] for k in self.keys]
-        else:
-            self.dets = det_offsets
-        # Make sure it's a numpy array.  This is dumb.
-        if isinstance(self.dets[0], quat.quat):
-            self.dets = quat.G3VectorQuat(self.dets)
-        self.dets = np.array(self.dets)
+        self.fplane = fplane
         return self
 
     @classmethod
@@ -347,5 +309,5 @@ def for_boresight(cls, sight_line):
         """
         self = cls(collapse=True)
         self.Q = sight_line.Q
-        self.dets = [np.array([1., 0, 0, 0])]
+        self.fplane = FocalPlane()
         return self
diff --git a/python/proj/wcs.py b/python/proj/wcs.py
index 49eac62..f8b3a4d 100644
--- a/python/proj/wcs.py
+++ b/python/proj/wcs.py
@@ -373,7 +373,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.dets, 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.
@@ -446,8 +446,8 @@ def to_map(self, signal, assembly, output=None, det_weights=None,
             comps = self._guess_comps(output.shape)
         projeng = self.get_ProjEng(comps)
         q_native = self._cache_q_fp_to_native(assembly.Q)
-        map_out = projeng.to_map(
-            output, q_native, assembly.dets, signal, det_weights, threads)
+        map_out = projeng.to_map(output, q_native, assembly.fplane.quats,
+            assembly.fplane.resps, signal, det_weights, threads)
         return map_out
 
     def to_weights(self, assembly, output=None, det_weights=None,
@@ -471,8 +471,8 @@ def to_weights(self, assembly, output=None, det_weights=None,
             comps = self._guess_comps(output.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.dets, det_weights, threads)
+        map_out = projeng.to_weight_map(output, q_native, assembly.fplane.quats,
+            assembly.fplane.resps, det_weights, threads)
         return map_out
 
     def from_map(self, src_map, assembly, signal=None, comps=None):
@@ -494,8 +494,8 @@ def from_map(self, src_map, assembly, signal=None, comps=None):
             comps = self._guess_comps(src_map.shape)
         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.dets, signal)
+        signal_out = projeng.from_map(src_map, q_native, assembly.fplane.quats,
+            assembly.fplane.resps, signal)
         return signal_out
 
     def assign_threads(self, assembly, method='domdir', n_threads=None):
diff --git a/src/Projection.cxx b/src/Projection.cxx
index ae7951d..d5869ae 100644
--- a/src/Projection.cxx
+++ b/src/Projection.cxx
@@ -22,6 +22,7 @@ using namespace std;
 #include "exceptions.h"
 #include "so_linterp.h"
 
+#include <stdio.h>
 
 // TRIG_TABLE_SIZE
 //
@@ -855,35 +856,47 @@ inline int Pixelizor2_Flat<Tiled, Bilinear>::GetPixels(int i_det, int i_time, co
     return iout;
 }
 
+// This struct is used for representing the total intensity and polarization
+// response of a single detector. Originally this was just a length-2 array,
+// but those can't easily be returned from functions. If I could ensure
+// the buffer had stride-1 in the last dimension, then one could simply pass
+// in the address to the first element for a detector, BufferWrapper doesn't make
+// that straightforward, so I think this approach is best.
+struct Response { FSIGNAL T, P; };
+Response get_response(BufferWrapper<FSIGNAL> & buf, int i_det) {
+    Response r = {*buf.ptr_2d(i_det,0), *buf.ptr_2d(i_det,1)};
+    return r;
+}
+
 template <typename SpinSys>
 static inline
-void spin_proj_factors(const double* coords, FSIGNAL *projfacs);
+void spin_proj_factors(const double* coords, const Response & response, FSIGNAL *projfacs);
 
 template <>
-inline void spin_proj_factors<SpinT>(const double* coords, FSIGNAL *projfacs)
+inline void spin_proj_factors<SpinT>(const double* coords, const Response & response, FSIGNAL *projfacs)
 {
-    projfacs[0] = 1.;
+    projfacs[0] = response.T;
 }
 
 template <>
 inline
-void spin_proj_factors<SpinQU>(const double* coords, FSIGNAL *projfacs)
+void spin_proj_factors<SpinQU>(const double* coords, const Response & response, FSIGNAL *projfacs)
 {
     const double c = coords[2];
     const double s = coords[3];
-    projfacs[0] = c*c - s*s;
-    projfacs[1] = 2*c*s;
+    projfacs[0] = response.P*(c*c - s*s);
+    projfacs[1] = response.P*(2*c*s);
 }
 
 template <>
 inline
-void spin_proj_factors<SpinTQU>(const double* coords, FSIGNAL *projfacs)
+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] = 1.;
-     projfacs[1] = c*c - s*s;
-     projfacs[2] = 2*c*s;
+     projfacs[0] = response.T;
+     projfacs[1] = response.P*(c*c - s*s);
+     projfacs[2] = response.P*(2*c*s);
 }
 
 
@@ -1068,7 +1081,7 @@ bp::object ProjectionEngine<C,P,S>::pixels(
 // sample
 template<typename C, typename P, typename S>
 bp::object ProjectionEngine<C,P,S>::pointing_matrix(
-    bp::object pbore, bp::object pofs, bp::object pixel, bp::object proj)
+    bp::object pbore, bp::object pofs, bp::object response, bp::object pixel, bp::object proj)
 {
     auto _none = bp::object();
 
@@ -1081,9 +1094,10 @@ bp::object ProjectionEngine<C,P,S>::pointing_matrix(
 
     auto pixel_buf_man = SignalSpace<int32_t>(
         pixel, "pixel", NPY_INT32, n_det, n_time, P::index_count);
-
     auto proj_buf_man = SignalSpace<FSIGNAL>(
         proj, "proj", FSIGNAL_NPY_TYPE, n_det, n_time, S::comp_count);
+    auto _response = BufferWrapper<FSIGNAL>(
+         "response", response, false, vector<int>{n_det,2});
 
 #pragma omp parallel for
     for (int i_det = 0; i_det < n_det; ++i_det) {
@@ -1093,12 +1107,13 @@ bp::object ProjectionEngine<C,P,S>::pointing_matrix(
         FSIGNAL* const proj_buf = proj_buf_man.data_ptr[i_det];
         const int step = pixel_buf_man.steps[0];
         int pixel_offset[P::index_count] = {-1};
+        Response resp = get_response(_response, i_det);
         for (int i_time = 0; i_time < n_time; ++i_time) {
             double coords[4];
             FSIGNAL pf[S::comp_count];
             pointer.GetCoords(i_det, i_time, (double*)dofs, (double*)coords);
             _pixelizor.GetPixel(i_det, i_time, (double*)coords, pixel_offset);
-            spin_proj_factors<S>(coords, pf);
+            spin_proj_factors<S>(coords, resp, pf);
 
             for (int i_dim = 0; i_dim < P::index_count; i_dim++)
                 pix_buf[i_time * pixel_buf_man.steps[0] +
@@ -1403,7 +1418,7 @@ bp::object ProjectionEngine<C,P,S>::zeros(bp::object shape)
 
 template<typename C, typename P, typename S>
 bp::object ProjectionEngine<C,P,S>::from_map(
-    bp::object map, bp::object pbore, bp::object pofs, bp::object signal)
+    bp::object map, bp::object pbore, bp::object pofs, bp::object response, bp::object signal)
 {
     auto _none = bp::object();
 
@@ -1419,6 +1434,8 @@ bp::object ProjectionEngine<C,P,S>::from_map(
     // Get pointers to the signal and (optional) per-det weights.
     auto _signalspace = SignalSpace<FSIGNAL>(
             signal, "signal", FSIGNAL_NPY_TYPE, n_det, n_time);
+    auto _response = BufferWrapper<FSIGNAL>(
+         "response", response, false, vector<int>{n_det,2});
 
     #pragma omp parallel for
     for (int i_det = 0; i_det < n_det; ++i_det) {
@@ -1427,10 +1444,11 @@ bp::object ProjectionEngine<C,P,S>::from_map(
         int pixinds[P::interp_count][P::index_count] = {-1};
         FSIGNAL pixweights[P::interp_count];
         FSIGNAL pf[S::comp_count];
+        Response resp = get_response(_response, i_det);
         for (int i_time = 0; i_time < n_time; ++i_time) {
             double coords[4];
             pointer.GetCoords(i_det, i_time, (double*)dofs, (double*)coords);
-            spin_proj_factors<S>(coords, pf);
+            spin_proj_factors<S>(coords, resp, pf);
             FSIGNAL *sig = (_signalspace.data_ptr[i_det] + _signalspace.steps[0]*i_time);
             int n_point = _pixelizor.GetPixels(i_det, i_time, coords, pixinds, pixweights);
             for(int i_point = 0; i_point < n_point; ++i_point)
@@ -1445,6 +1463,7 @@ bp::object ProjectionEngine<C,P,S>::from_map(
 template<typename C, typename P, typename S>
 static
 void to_map_single_thread(Pointer<C> &pointer,
+                          BufferWrapper<FSIGNAL> & _response,
                           P &_pixelizor,
                           const vector<RangesInt32> & ivals,
                           BufferWrapper<FSIGNAL> &_det_weights,
@@ -1462,13 +1481,14 @@ void to_map_single_thread(Pointer<C> &pointer,
         double coords[4];
         FSIGNAL pf[S::comp_count];
         pointer.InitPerDet(i_det, dofs);
+        Response resp = get_response(_response, i_det);
         // Pointing matrix interpolation stuff
         int pixinds[P::interp_count][P::index_count] = {-1};
         FSIGNAL pixweights[P::interp_count] = {0};
         for (auto const &rng: ivals[i_det].segments) {
             for (int i_time = rng.first; i_time < rng.second; ++i_time) {
                 pointer.GetCoords(i_det, i_time, (double*)dofs, (double*)coords);
-                spin_proj_factors<S>(coords, pf);
+                spin_proj_factors<S>(coords, resp, pf);
                 const FSIGNAL sig = _signalspace->data_ptr[i_det][_signalspace->steps[0]*i_time];
                 // In interpolated mapmaking like bilinear mampamking, each sample hits multipe
                 // pixels, each with its own weight.
@@ -1484,6 +1504,7 @@ void to_map_single_thread(Pointer<C> &pointer,
 template<typename C, typename P, typename S>
 static
 void to_weight_map_single_thread(Pointer<C> &pointer,
+                                 BufferWrapper<FSIGNAL> & _response,
                                  P &_pixelizor,
                                  vector<RangesInt32> ivals,
                                  BufferWrapper<FSIGNAL> &_det_weights)
@@ -1500,12 +1521,13 @@ void to_weight_map_single_thread(Pointer<C> &pointer,
         double coords[4];
         FSIGNAL pf[S::comp_count];
         pointer.InitPerDet(i_det, dofs);
+        Response resp = get_response(_response, i_det);
         int pixinds[P::interp_count][P::index_count] = {-1};
         FSIGNAL pixweights[P::interp_count] = {0};
         for (auto const &rng: ivals[i_det].segments) {
             for (int i_time = rng.first; i_time < rng.second; ++i_time) {
                 pointer.GetCoords(i_det, i_time, (double*)dofs, (double*)coords);
-                spin_proj_factors<S>(coords, pf);
+                spin_proj_factors<S>(coords, resp, pf);
 
                 int n_point = _pixelizor.GetPixels(i_det, i_time, coords, pixinds, pixweights);
                 // Is this enough? Do we need a double loop over i_point?
@@ -1592,8 +1614,8 @@ vector<vector<vector<RangesInt32>>> derive_ranges(
 
 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 signal, bp::object det_weights,
-    bp::object thread_intervals)
+    bp::object map, bp::object pbore, bp::object pofs, bp::object response,
+    bp::object signal, bp::object det_weights, bp::object thread_intervals)
 {
     //Initialize it / check inputs.
     auto pointer = Pointer<C>();
@@ -1613,6 +1635,8 @@ bp::object ProjectionEngine<C,P,S>::to_map(
             signal, "signal", FSIGNAL_NPY_TYPE, n_det, n_time);
     auto _det_weights = BufferWrapper<FSIGNAL>(
          "det_weights", det_weights, true, vector<int>{n_det});
+    auto _response = BufferWrapper<FSIGNAL>(
+         "response", response, false, vector<int>{n_det,2});
 
     // For multi-threading, the principle here is that we loop serially
     // over bunches, and then inside each block all threads loop over
@@ -1627,15 +1651,15 @@ bp::object ProjectionEngine<C,P,S>::to_map(
         // or if ivals.size() == 1
         #pragma omp parallel for
         for (int i_thread = 0; i_thread < ivals.size(); i_thread++)
-            to_map_single_thread<C,P,S>(pointer, _pixelizor, ivals[i_thread], _det_weights, &_signalspace);
+            to_map_single_thread<C,P,S>(pointer, _response, _pixelizor, ivals[i_thread], _det_weights, &_signalspace);
     }
     return map;
 }
 
 template<typename C, typename P, typename S>
 bp::object ProjectionEngine<C,P,S>::to_weight_map(
-    bp::object map, bp::object pbore, bp::object pofs, bp::object det_weights,
-    bp::object thread_intervals)
+    bp::object map, bp::object pbore, bp::object pofs, bp::object response,
+    bp::object det_weights, bp::object thread_intervals)
 {
     auto _none = bp::object();
 
@@ -1655,6 +1679,8 @@ bp::object ProjectionEngine<C,P,S>::to_weight_map(
     // Get pointer to (optional) per-det weights.
     auto _det_weights = BufferWrapper<FSIGNAL>(
          "det_weights", det_weights, true, vector<int>{n_det});
+    auto _response = BufferWrapper<FSIGNAL>(
+         "response", response, false, vector<int>{n_det,2});
 
     // For multi-threading, the principle here is that we loop serially
     // over bunches, and then inside each block all threads loop over
@@ -1669,7 +1695,7 @@ bp::object ProjectionEngine<C,P,S>::to_weight_map(
         // or if ivals.size() == 1
         #pragma omp parallel for
         for (int i_thread = 0; i_thread < ivals.size(); i_thread++)
-            to_weight_map_single_thread<C,P,S>(pointer, _pixelizor, ivals[i_thread], _det_weights);
+            to_weight_map_single_thread<C,P,S>(pointer, _response, _pixelizor, ivals[i_thread], _det_weights);
     }
     return map;
 }

From ac3e64a37b22063711a562d9225aa2cd477cb38c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sigurd=20Kirkevold=20N=C3=A6ss?= <sigurdkn@astro.uio.no>
Date: Sun, 23 Jun 2024 13:14:22 -0400
Subject: [PATCH 3/4] Arbitrary response support

---
 python/proj/coords.py     | 22 +++++++++++++++++-----
 python/proj/mapthreads.py | 24 ++++++++++++------------
 python/proj/wcs.py        | 20 ++++++++++----------
 3 files changed, 39 insertions(+), 27 deletions(-)

diff --git a/python/proj/coords.py b/python/proj/coords.py
index 06d816b..093504c 100644
--- a/python/proj/coords.py
+++ b/python/proj/coords.py
@@ -231,7 +231,7 @@ def coords(self, fplane=None, output=None):
             fplane = FocalPlane()
             if output is not None:
                 output = output[None]
-        output = p.coords(self.Q, fplane.quats, fplane.resps, output)
+        output = p.coords(self.Q, fplane.quats, output)
         if collapse:
             output = output[0]
         return output
@@ -253,20 +253,32 @@ def __init__(self, quats=None, resps=None):
         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")
     @classmethod
-    def from_xieta(cls, xi, eta, gamma=0, T=1, P=1, Q=1, U=0):
+    def from_xieta(cls, xi, eta, gamma=0, T=1, P=1, Q=1, U=0, hwp=False):
         # 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
         gamma = gamma + np.arctan2(U,Q)/2
         P     = P * (Q**2+U**2)**0.5
-        quats = quat.rotation_xieta(np.asarray(xi), np.asarray(eta), np.asarray(gamma))
+        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
         return cls(quats, resps=resps)
+    def __repr__(self):
+        return "FocalPlane(quats=%s,resps=%s)" % (repr(self.quats), repr(self.resps))
+    @property
+    def ndet(self): return len(self.quats)
+    def __getitem__(self, sel):
+        return FocalPlane(quats=self.quats[sel], resps=self.resps[sel])
 
 class Assembly:
     """This class groups together boresight and detector offset
@@ -293,7 +305,7 @@ def attach(cls, sight_line, fplane):
                 relative to the boresight, and their response to
                 intensity and polarization
         """
-        self = cls(keyed=keyed)
+        self = cls()
         if isinstance(sight_line, quat.G3VectorQuat):
             self.Q = sight_line
         else:
diff --git a/python/proj/mapthreads.py b/python/proj/mapthreads.py
index 0440ea6..92e71f2 100644
--- a/python/proj/mapthreads.py
+++ b/python/proj/mapthreads.py
@@ -19,8 +19,8 @@ def get_num_threads(n_threads=None):
         return so3g.useful_info()['omp_num_threads']
     return n_threads
 
-def get_threads_domdir(sight, offs, shape, wcs, tile_shape=None,
-                       active_tiles=None, n_threads=None, offs_rep=None,
+def get_threads_domdir(sight, fplane, shape, wcs, tile_shape=None,
+                       active_tiles=None, n_threads=None, fplane_rep=None,
                        plot_prefix=None):
     """Assign samples to threads according to the dominant scan
     direction.
@@ -37,7 +37,7 @@ def get_threads_domdir(sight, offs, shape, wcs, tile_shape=None,
 
     Arguments:
       sight (CelestialSightLine): The boresight pointing.
-      offs (array of quaternions): The detector pointing.
+      fplane (FocalPlane): The detector pointing
       shape (tuple): The map shape.
       wcs (wcs): The map WCS.
       tile_shape (tuple): The tile shape, if this should be done using
@@ -46,9 +46,9 @@ def get_threads_domdir(sight, offs, shape, wcs, tile_shape=None,
         will be computed along the way.
       n_threads (int): The number of threads to target (defaults to
         OMP_NUM_THREADS).
-      offs_rep (array of quaternions): A representative set of
+      fplane_rep (FocalPlane): A representative set of
         detectors, for determining scan direction (if not present then
-        offs is used for this).
+        fplane is used for this).
       plot_prefix (str): If present, pylab will be imported and some
         plots saved with this name as prefix.
 
@@ -67,8 +67,8 @@ def get_threads_domdir(sight, offs, shape, wcs, tile_shape=None,
                                       if _m is not None]
 
     n_threads = get_num_threads(n_threads)
-    if offs_rep is None:
-        offs_rep = offs
+    if fplane_rep is None:
+        fplane_rep = fplane
 
     if tile_shape is None:
         # Let's pretend it is, though; this simplifies logic below and
@@ -77,7 +77,7 @@ def get_threads_domdir(sight, offs, shape, wcs, tile_shape=None,
         active_tiles = [0]
 
     # The full assembly, for later.
-    asm_full = so3g.proj.Assembly.attach(sight, offs)
+    asm_full = so3g.proj.Assembly.attach(sight, fplane)
 
     # Get a Projectionist -- note it can be used with full or
     # representative assembly.
@@ -93,10 +93,10 @@ def get_threads_domdir(sight, offs, shape, wcs, tile_shape=None,
     # For the scan direction map, use the "representative" subset
     # detectors, with polarization direction aligned parallel to
     # elevation.
-    xi, eta, gamma = so3g.proj.quat.decompose_xieta(offs_rep)
-    offs_xl = np.array(so3g.proj.quat.rotation_xieta(xi, eta, gamma*0 + 90*so3g.proj.DEG))
-    asm_rep = so3g.proj.Assembly.attach(sight, offs_xl)
-    sig = np.ones((len(offs_xl), len(asm_rep.Q)), dtype='float32')
+    xi, eta, gamma = so3g.proj.quat.decompose_xieta(fplane_rep.quats)
+    fplane_xl = so3g.proj.FocalPlane.from_xieta(xi, eta, gamma*0+90*so3g.proj.DEG)
+    asm_rep = so3g.proj.Assembly.attach(sight, fplane_xl)
+    sig = np.ones((fplane_xl.ndet, len(asm_rep.Q)), dtype='float32')
     scan_maps = pmat.to_map(sig, asm_rep, comps='TQU')
 
     # Compute the scan angle, based on Q and U weights.  This assumes
diff --git a/python/proj/wcs.py b/python/proj/wcs.py
index f8b3a4d..4ed06bb 100644
--- a/python/proj/wcs.py
+++ b/python/proj/wcs.py
@@ -357,7 +357,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.dets, 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
@@ -399,7 +399,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.dets, 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.
@@ -414,7 +414,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.dets, 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,) +
@@ -538,11 +538,11 @@ 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.dets, None, n_threads)
+            omp_ivals = projeng.pixel_ranges(q_native, assembly.fplane.quats, None, n_threads)
             return wrap_ivals(omp_ivals)
 
         elif method == 'domdir':
-            offs_rep = assembly.dets[::100]
+            fplane_rep = assembly.fplane[::100]
             if (self.tiling is not None) and (self.active_tiles is None):
                 tile_info = self.get_active_tiles(assembly)
                 active_tiles = tile_info['active_tiles']
@@ -550,9 +550,9 @@ def assign_threads(self, assembly, method='domdir', n_threads=None):
             else:
                 active_tiles = None
             return mapthreads.get_threads_domdir(
-                assembly.Q, assembly.dets, shape=self.naxis[::-1], wcs=self.wcs,
+                assembly.Q, assembly.fplane, shape=self.naxis[::-1], wcs=self.wcs,
                 tile_shape=self.tile_shape, active_tiles=active_tiles,
-                n_threads=n_threads, offs_rep=offs_rep)
+                n_threads=n_threads, fplane_rep=fplane_rep)
 
         elif method == 'tiles':
             tile_info = self.get_active_tiles(assembly, assign=n_threads)
@@ -578,7 +578,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.dets, 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):
@@ -617,7 +617,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.dets))
+        hits = np.array(projeng.tile_hits(q_native, assembly.fplane.quats))
         tiles = np.nonzero(hits)[0]
         hits = hits[tiles]
         if assign is True:
@@ -631,7 +631,7 @@ def get_active_tiles(self, assembly, assign=False):
                 group_n[idx] += _n
                 group_tiles[idx].append(_t)
             # Now paint them into Ranges.
-            R = projeng.tile_ranges(q_native, assembly.dets, group_tiles)
+            R = projeng.tile_ranges(q_native, assembly.fplane.quats, group_tiles)
             R = wrap_ivals(R)
             return {
                 'active_tiles': list(tiles),

From a2eb4dd51d6bb782fe6abffd88b0121f73ac8ef5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sigurd=20Kirkevold=20N=C3=A6ss?= <sigurdkn@astro.uio.no>
Date: Sun, 1 Dec 2024 14:28:20 -0500
Subject: [PATCH 4/4] Passes tests after adding .items(), updating xieta
 arguments and replacing asm.dets with asm.fplane.quats

---
 python/proj/coords.py    |  3 +++
 test/test_proj_astro.py  | 11 +++++------
 test/test_proj_eng.py    |  5 ++---
 test/test_proj_eng_hp.py |  4 ++--
 4 files changed, 12 insertions(+), 11 deletions(-)

diff --git a/python/proj/coords.py b/python/proj/coords.py
index ddb2a6e..1e11395 100644
--- a/python/proj/coords.py
+++ b/python/proj/coords.py
@@ -279,6 +279,9 @@ def __repr__(self):
     def ndet(self): return len(self.quats)
     def __getitem__(self, sel):
         return FocalPlane(quats=self.quats[sel], resps=self.resps[sel])
+    def items(self):
+        for q, resp in zip(quat.G3VectorQuat(self.quats), self.resps):
+            yield q, resp
 
 class Assembly:
     """This class groups together boresight and detector offset
diff --git a/test/test_proj_astro.py b/test/test_proj_astro.py
index 35ebe94..3d488f3 100644
--- a/test/test_proj_astro.py
+++ b/test/test_proj_astro.py
@@ -195,12 +195,10 @@ def test_horizon(self):
 
     def test_focalplane(self):
         # Focal plane
-        names = ['a', 'b', 'c']
         xi  = np.array([1., 0., 0.]) * DEG
         eta = np.array([0., 0., 1.]) * DEG
         gamma = np.array([45, 0, -45]) * DEG
-        fp = so3g.proj.FocalPlane.from_xieta(names, xi, eta, gamma)
-        print(fp)
+        fp = so3g.proj.FocalPlane.from_xieta(xi, eta, gamma)
 
         # Make a CelestialSightLine with some known equatorial pointing.
         r = so3g.proj.quat.rotation_lonlat(35*DEG, 80*DEG, 15*DEG)
@@ -213,13 +211,14 @@ def test_focalplane(self):
         coords = csl.coords(fp)
 
         # Compare to a more direct computation, done here.
-        for k, roffset in fp.items():
-            print('Check detector %s:' % k)
+        for i, (roffset, resp) in enumerate(fp.items()):
+            print('Check detector %d:' % i)
             r1 = r * roffset
+            print("AAA", type(r1), r1)
             lon0, lat0, gamma0 = so3g.proj.quat.decompose_lonlat(r1)
             print('  - inline computation: ',
                   lon0 / DEG, lat0 / DEG, gamma0 / DEG)
-            lon1, lat1, cosg, sing = coords[k][0]
+            lon1, lat1, cosg, sing = coords[i][0]
             gamma1 = np.arctan2(sing, cosg)
             print('  - proj computation: ',
                   lon1 / DEG, lat1 / DEG, gamma1 / DEG)
diff --git a/test/test_proj_eng.py b/test/test_proj_eng.py
index b5c190f..b812d57 100644
--- a/test/test_proj_eng.py
+++ b/test/test_proj_eng.py
@@ -35,7 +35,7 @@ def get_scan():
 def get_basics(clipped=True):
     t, az, el = get_scan()
     csl = proj.CelestialSightLine.az_el(t, az, el, weather='vacuum', site='so')
-    fp = proj.FocalPlane.from_xieta(['a', 'b'], [0., .1*DEG], [0, .1*DEG])
+    fp = proj.FocalPlane.from_xieta([0., .1*DEG], [0, .1*DEG])
     asm = proj.Assembly.attach(csl, fp)
 
     # And a map ... of where?
@@ -73,8 +73,7 @@ def test_00_basic(self):
                      det_weights=np.array([0., 0.], dtype='float32'))[0]
         assert(np.all(m==0))
 
-        # Raise if pointing invalid.
-        asm.dets[1,2] = np.nan
+        asm.fplane.quats[1,2] = np.nan
         with self.assertRaises(ValueError):
            p.to_map(sig, asm, comps='T')
         with self.assertRaises(ValueError):
diff --git a/test/test_proj_eng_hp.py b/test/test_proj_eng_hp.py
index a7eb88c..6ef20e1 100644
--- a/test/test_proj_eng_hp.py
+++ b/test/test_proj_eng_hp.py
@@ -23,7 +23,7 @@ def get_scan():
 def get_basics():
     t, az, el = get_scan()
     csl = proj.CelestialSightLine.az_el(t, az, el, weather='vacuum', site='so')
-    fp = proj.FocalPlane.from_xieta(['a', 'b'], [0., .1*DEG], [0, .1*DEG])
+    fp = proj.FocalPlane.from_xieta([0., .1*DEG], [0, .1*DEG])
     asm = proj.Assembly.attach(csl, fp)
     return ((t, az, el), asm)
 
@@ -52,7 +52,7 @@ def test_00_basic(self):
         assert(np.all(m==0))
 
         # Raise if pointing invalid.
-        asm.dets[1,2] = np.nan
+        asm.fplane.quats[1,2] = np.nan
         with self.assertRaises(ValueError):
            p.to_map(sig, asm, comps='T')
         with self.assertRaises(ValueError):