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

demos/demo_proj_hp.py

@@ -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):
@@ -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
 

docs/proj.rst

@@ -85,9 +85,9 @@ quaternions into rotation angles::
   >>> csl.Q
   spt3g.core.G3VectorQuat([(-0.0384748,0.941783,0.114177,0.313891)])
 
-  >>> csl.coords()   
+  >>> csl.coords()
   array([[ 0.24261138, -0.92726871, -0.99999913, -0.00131952]])
-  
+
 The :func:`coords() <CelestialSightLine.coords>` returns an array with
 shape (n_time, 4); each 4-tuple contains values ``(lon, lat,
 cos(gamma), sin(gamma))``.  The ``lon`` and ``lat`` are the celestial
@@ -110,48 +110,77 @@ Pointing for many detectors
 Create a :class:`FocalPlane` object, with some detector positions and
 orientations::
 
-  names = ['a', 'b', 'c']
-  x = np.array([-0.5, 0., 0.5]) * DEG
-  y = np.zeros(3)
+  xi  = np.array([-0.5, 0.0, 0.5]) * DEG
+  eta = np.zeros(3)
   gamma = np.array([0,30,60]) * DEG
-  fp = so3g.proj.FocalPlane.from_xieta(names, x, y, gamma)
+  fp  = so3g.proj.FocalPlane.from_xieta(xi, eta, gamma)
 
 This particular function, :func:`from_xieta()
 <FocalPlane.from_xieta>`, will apply the SO standard coordinate
-definitions and store a rotation quaternion for each detector.
-FocalPlane is just a thinly wrapped OrderedDict, where the detector
-name is the key and the value is the rotation quaternion::
-
-  >>> fp['c']
+definitions and stores quaternions (``.quats``) and
+responsivities (``.resps``) for each detector. These are a
+``G3VectorQuat``  wth length ``ndet`` and a numpy array
+with shape ``[ndet,2]`` respectively. ``fp.quats[0]`` gives
+the quaternion for the first detector, while ``fp.resps[0]``
+gives its total intensity and polarization responsivity.::
+
+  >>> fp.quats[2]
   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. 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)
+  gamma = np.array([0,30,60]) * DEG
+  T   = np.array([1.0, 0.9, 1.1])
+  P   = np.array([0.5, 0.6, 0.05])
+  fp2 = so3g.proj.FocalPlane.from_xieta(xi, eta, gamma, T=T, P=P)
+
+Together, gamma, T and P specify the full responsivity of a
+detector to the T, Q and U Stokes parameters in focal plane
+coordinates. But as an alternative, it's also possible to
+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, 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)
 
 At this point you could get the celestial coordinates for any one of
 those detectors::
 
-  # Get vector of quaternion pointings for detector 'a'
-  q_total = csl.Q * fp['a']
+  # Get vector of quaternion pointings for detector 0
+  q_total = csl.Q * fp.quats[0]
   # Decompose q_total into lon, lat, roll angles
   ra, dec, gamma = so3g.proj.quat.decompose_lonlat(q_total)
 
 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
 FocalPlane object as first argument::
 
   >>> csl.coords(fp)
-  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,
+  [array([[ 0.25715457, -0.92720643,  0.9999161 ,  0.01295328]]),
+   array([[ 0.24261138, -0.92726871,  0.86536489,  0.5011423 ]]),
+   array([[ 0.22806774, -0.92722945,  0.50890634,  0.86082189]])]
+
+So ``coords()`` returns a list of numpy arrays with shape (n_time,4),
+and at each time step the 4 elements of the array are: ``(lon, lat,
 cos(gamma), sin(gamma))``.
 
 
@@ -215,17 +244,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`::
@@ -240,9 +269,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
@@ -251,7 +280,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
@@ -264,8 +293,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
@@ -289,7 +318,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]])
@@ -352,7 +381,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``.

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;