Stokes response

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')

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...
 

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;

python/proj/coords.py

@@ -208,93 +208,80 @@ 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, 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.
     """
+    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
+        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, 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 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
+        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])
+    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
@@ -305,38 +292,28 @@ 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:
             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.
-            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)
+        self = cls()
         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 +324,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